C++ OsiSolverInterface::readMps方法代码示例

int
main(void)
{
   // Create a problem pointer.  We use the base class here.
   OsiSolverInterface *si;

   // When we instantiate the object, we need a specific derived class.
   si = new OsiClpSolverInterface;

   // Read in an mps file.  This one's from the MIPLIB library.
   si->readMps("../../Data/Sample/p0033");

   // Solve the (relaxation of the) problem
   si->initialSolve();

   // Check the solution
   if ( si->isProvenOptimal() ) { 
      std::cout << "Found optimal solution!" << std::endl; 
      std::cout << "Objective value is " << si->getObjValue() << std::endl;

      int n = si->getNumCols();
      const double *solution;
      solution = si->getColSolution();
      // We could then print the solution or examine it.
   } else {
      std::cout << "Didn't find optimal solution." << std::endl;
      // Could then check other status functions.
   }

   return 0;
}

//--------------------------------------------------------------------------
// test cut debugger methods.
void
OsiRowCutDebuggerUnitTest(const OsiSolverInterface * baseSiP, const std::string & mpsDir)
{
  
  CoinRelFltEq eq;
  
  // Test default constructor
  {
    OsiRowCutDebugger r;
    OSIUNITTEST_ASSERT_ERROR(r.integerVariable_ == NULL, {}, "osirowcutdebugger", "default constructor");
    OSIUNITTEST_ASSERT_ERROR(r.knownSolution_   == NULL, {}, "osirowcutdebugger", "default constructor");
    OSIUNITTEST_ASSERT_ERROR(r.numberColumns_   == 0,    {}, "osirowcutdebugger", "default constructor");
  }
  
  {
    // Get non trivial instance
    OsiSolverInterface * imP = baseSiP->clone();
    std::string fn = mpsDir+"exmip1";
    imP->readMps(fn.c_str(),"mps");
    OSIUNITTEST_ASSERT_ERROR(imP->getNumRows() == 5, {}, "osirowcutdebugger", "read exmip1");
    
    /*
      Activate the debugger. The garbled name here is deliberate; the
      debugger should isolate the portion of the string between '/' and
      '.' (in normal use, this would be the base file name, stripped of
      the prefix and extension).
    */
    imP->activateRowCutDebugger("ab cd /x/ /exmip1.asc");
    
    int i;
    
    // return debugger
    const OsiRowCutDebugger * debugger = imP->getRowCutDebugger();
    OSIUNITTEST_ASSERT_ERROR(debugger != NULL, {}, "osirowcutdebugger", "return debugger");
    OSIUNITTEST_ASSERT_ERROR(debugger->numberColumns_ == 8, {}, "osirowcutdebugger", "return debugger");
    
    const bool type[]={0,0,1,1,0,0,0,0};
    const double values[]= {2.5, 0, 1, 1, 0.5, 3, 0, 0.26315789473684253};
    CoinPackedVector objCoefs(8,imP->getObjCoefficients());
   
    bool type_ok = true;
#if 0
    for (i=0;i<8;i++)
      type_ok &= type[i] == debugger->integerVariable_[i];
    OSIUNITTEST_ASSERT_ERROR(type_ok, {}, "osirowcutdebugger", "???");
#endif
    
    double objValue = objCoefs.dotProduct(values);
    double debuggerObjValue = objCoefs.dotProduct(debugger->knownSolution_);
    OSIUNITTEST_ASSERT_ERROR(eq(objValue, debuggerObjValue), {}, "osirowcutdebugger", "objective value");
    
    OsiRowCutDebugger rhs;
    {
      OsiRowCutDebugger rC1(*debugger);

      OSIUNITTEST_ASSERT_ERROR(rC1.numberColumns_ == 8, {}, "osirowcutdebugger", "copy constructor");
      type_ok = true;
      for (i=0;i<8;i++)
      	type_ok &= type[i] == rC1.integerVariable_[i];
      OSIUNITTEST_ASSERT_ERROR(type_ok, {}, "osirowcutdebugger", "copy constructor");
      OSIUNITTEST_ASSERT_ERROR(eq(objValue,objCoefs.dotProduct(rC1.knownSolution_)), {}, "osirowcutdebugger", "copy constructor");
      
      rhs = rC1;
      OSIUNITTEST_ASSERT_ERROR(rhs.numberColumns_ == 8, {}, "osirowcutdebugger", "assignment operator");
      type_ok = true;
      for (i=0;i<8;i++)
      	type_ok &= type[i] == rhs.integerVariable_[i];
      OSIUNITTEST_ASSERT_ERROR(type_ok, {}, "osirowcutdebugger", "assignment operator");
      OSIUNITTEST_ASSERT_ERROR(eq(objValue,objCoefs.dotProduct(rhs.knownSolution_)), {}, "osirowcutdebugger", "assignment operator");
    }
    // Test that rhs has correct values even though lhs has gone out of scope
    OSIUNITTEST_ASSERT_ERROR(rhs.numberColumns_ == 8, {}, "osirowcutdebugger", "assignment operator");
    type_ok = true;
    for (i=0;i<8;i++)
    	type_ok &= type[i] == rhs.integerVariable_[i];
    OSIUNITTEST_ASSERT_ERROR(type_ok, {}, "osirowcutdebugger", "assignment operator");
    OSIUNITTEST_ASSERT_ERROR(eq(objValue,objCoefs.dotProduct(rhs.knownSolution_)), {}, "osirowcutdebugger", "assignment operator");

    OsiRowCut cut[2];
    
    const int ne = 3;
    int inx[ne] = { 0, 2, 3 };
    double el[ne] = { 1., 1., 1. };
    cut[0].setRow(ne,inx,el);
    cut[0].setUb(5.);
    
    el[1]=5;
    cut[1].setRow(ne,inx,el);
    cut[1].setUb(5);
    OsiCuts cs; 
    cs.insert(cut[0]);
    cs.insert(cut[1]);
    OSIUNITTEST_ASSERT_ERROR(!debugger->invalidCut(cut[0]), {}, "osirowcutdebugger", "recognize (in)valid cut");
    OSIUNITTEST_ASSERT_ERROR( debugger->invalidCut(cut[1]), {}, "osirowcutdebugger", "recognize (in)valid cut");
    OSIUNITTEST_ASSERT_ERROR(debugger->validateCuts(cs,0,2) == 1, {}, "osirowcutdebugger", "recognize (in)valid cut");
    OSIUNITTEST_ASSERT_ERROR(debugger->validateCuts(cs,0,1) == 0, {}, "osirowcutdebugger", "recognize (in)valid cut");
    delete imP;
  }
//.........这里部分代码省略.........

//--------------------------------------------------------------------------
// test the simple rounding cut generators methods.
void
CglSimpleRoundingUnitTest(
  const OsiSolverInterface * baseSiP,
  const std::string mpsDir )
{

  // Test default constructor
  {
    CglSimpleRounding cg;
  }

  // Test copy & assignment
  {
    CglSimpleRounding rhs;
    {
      CglSimpleRounding cg;
      CglSimpleRounding cgC(cg);
      rhs=cg;
    }
  }

  // Test gcd and gcdn
  {
    CglSimpleRounding cg;
    int v = cg.gcd(122,356);
    assert(v==2);
    v=cg.gcd(356,122);
    assert(v==2);
    v=cg.gcd(54,67);
    assert(v==1);
    v=cg.gcd(67,54);
    assert(v==1);
    v=cg.gcd(485,485);
    assert(v==485);
    v=cg.gcd(17*13,17*23);
    assert( v==17);
    v=cg.gcd(17*13*5,17*23);
    assert( v==17);
    v=cg.gcd(17*13*23,17*23);
    assert(v==17*23);

    int a[4] = {12, 20, 32, 400};
    v= cg.gcdv(4,a);
    assert(v== 4);
    int b[4] = {782, 4692, 51, 2754};
    v= cg.gcdv(4,b);
    assert(v== 17);
    int c[4] = {50, 40, 30, 10};
    v= cg.gcdv(4,c);
    assert(v== 10);
  }


  // Test generate cuts method on exmip1.5.mps
  {
    CglSimpleRounding cg;
    
    OsiSolverInterface * siP = baseSiP->clone();
    std::string fn = mpsDir+"exmip1.5.mps";
    siP->readMps(fn.c_str(),"");
    OsiCuts cuts;
    cg.generateCuts(*siP,cuts);

    // there should be 3 cuts
    int nRowCuts = cuts.sizeRowCuts();
    assert(nRowCuts==3);

    // get the last "sr"=simple rounding cut that was derived
    OsiRowCut srRowCut2 = cuts.rowCut(2); 
    CoinPackedVector srRowCutPV2 = srRowCut2.row();

    // this is what the last cut should look like: i.e. the "solution"
    const int solSize = 2;
    int solCols[solSize]={2,3};
    double solCoefs[solSize]={5.0, 4.0};
    OsiRowCut solRowCut;
    solRowCut.setRow(solSize,solCols,solCoefs);
    solRowCut.setLb(-COIN_DBL_MAX);
    solRowCut.setUb(2.0);

    // Test for equality between the derived cut and the solution cut

    // Note: testing two OsiRowCuts are equal invokes testing two
    // CoinPackedVectors are equal which invokes testing two doubles
    // are equal.  Usually not a good idea to test that two doubles are equal, 
    // but in this cut the "doubles" represent integer values. Also allow that
    // different solvers have different orderings in packed vectors, which may
    // not match the ordering defined for solRowCut.

    assert(srRowCut2.OsiCut::operator==(solRowCut)) ;
    assert(srRowCut2.row().isEquivalent(solRowCut.row())) ;
    assert(srRowCut2.lb() == solRowCut.lb()) ;
    assert(srRowCut2.ub() == solRowCut.ub()) ;

    delete siP;
  }

  // Test generate cuts method on p0033
//.........这里部分代码省略.........

//--------------------------------------------------------------------------
// test EKKsolution methods.
void
CglOddHoleUnitTest(
  const OsiSolverInterface * baseSiP,
  const std::string mpsDir )
{
  CoinRelFltEq eq(0.000001);

  // Test default constructor
  {
    CglOddHole aGenerator;
  }
  
  // Test copy & assignment
  {
    CglOddHole rhs;
    {
      CglOddHole bGenerator;
      CglOddHole cGenerator(bGenerator);
      rhs=bGenerator;
    }
  }


  // test on simple case
  {  
    const int nRows=3;
    const int nCols=3;
    const int nEls=6;
    const double elem[]={1.0,1.0,1.0,1.0,1.0,1.0};
    const int row[]={0,1,0,2,1,2};
    const CoinBigIndex start[]={0,2,4};
    const int len[]={2,2,2};
    CoinPackedMatrix matrix(true,nRows,nCols,nEls,elem,row,start,len);
    const double sol[]={0.5,0.5,0.5};
    const double dj[]={0,0,0};
    const int which[]={1,1,1};
    const int fixed[]={0,0,0};
    OsiCuts cs;
    CglOddHole test1;
    CglTreeInfo info;
    info.randomNumberGenerator=NULL;
    test1.generateCuts(NULL,matrix,sol,dj,cs,which,fixed,info,true);
    CoinPackedVector check;
    int index[] = {0,1,2};
    double el[] = {1,1,1};
    check.setVector(3,index,el);
    //assert (cs.sizeRowCuts()==2);
    assert (cs.sizeRowCuts()==1);
    // sort Elements in increasing order
    CoinPackedVector rpv=cs.rowCut(0).row();
    rpv.sortIncrIndex();
    assert (check==rpv);
  }
  
  // Testcase /u/rlh/osl2/mps/scOneInt.mps
  // Model has 3 continous, 2 binary, and 1 general
  // integer variable.
  {
    OsiSolverInterface  * siP = baseSiP->clone();
    
    std::string fn = mpsDir+"scOneInt";
    siP->readMps(fn.c_str(),"mps");
#if 0
    CglOddHole cg;
    int nCols=siP->getNumCols();
    
    // Test the siP methods for detecting
    // variable type
    int numCont=0, numBinary=0, numIntNonBinary=0, numInt=0;
    for (int thisCol=0; thisCol<nCols; thisCol++) {
      if ( siP->isContinuous(thisCol) ) numCont++;
      if ( siP->isBinary(thisCol) ) numBinary++;
      if ( siP->isIntegerNonBinary(thisCol) ) numIntNonBinary++;
      if ( siP->isInteger(thisCol) ) numInt++;
    }
    assert(numCont==3);
    assert(numBinary==2);
    assert(numIntNonBinary==1);
    assert(numInt==3);
    
    
    // Test initializeCutGenerator
    siP->initialSolve();
    assert(xstar !=NULL);
    for (i=0; i<nCols; i++){
      assert(complement[i]==0);
    }
    int nRows=siP->getNumRows();
    for (i=0; i<nRows; i++){
    int vectorsize = siP->getMatrixByRow()->getVectorSize(i);
    assert(vectorsize==2);
    }
    
    kccg.cleanUpCutGenerator(complement,xstar);
#endif  
    delete siP;
  }

//.........这里部分代码省略.........

void
CglTwomirUnitTest(const OsiSolverInterface *baseSiP,
		  const std::string mpsDir)
{
  // Test default constructor
  {
    CglTwomir aGenerator;
  }
  
  // Test copy & assignment
  {
    CglTwomir rhs;
    {
      CglTwomir bGenerator;
      CglTwomir cGenerator(bGenerator);
      rhs=bGenerator;
    }
  }

  // Test get/set methods
  {
    CglTwomir getset;
    
    int gtmin = getset.getTmin() + 1;
    int gtmax = getset.getTmax() + 1;
    getset.setMirScale(gtmin, gtmax);
    double gtmin2 = getset.getTmin();
    double gtmax2 = getset.getTmax();
    assert(gtmin == gtmin2);
    assert(gtmax == gtmax2);

    int gamax = 2 * getset.getAmax() + 1;
    getset.setAMax(gamax);
    int gamax2 = getset.getAmax();
    assert(gamax == gamax2);
  }

  // Test generateCuts
  {
    CglTwomir gct;
    OsiSolverInterface  *siP = baseSiP->clone();
    std::string fn = mpsDir+"capPlan1";
    std::string fn2 = mpsDir+"capPlan1.mps";
    FILE *in_f = fopen(fn2.c_str(), "r");
    if(in_f == NULL) {
      std::cout<<"Can not open file "<<fn2<<std::endl<<"Skip test of CglTwomir::generateCuts()"<<std::endl;
    }
    else {
      fclose(in_f);
      siP->readMps(fn.c_str(),"mps");
 
      siP->initialSolve();
      double lpRelax = siP->getObjValue();
      
      OsiCuts cs;
      gct.generateCuts(*siP, cs);
      int nRowCuts = cs.sizeRowCuts();
      std::cout<<"There are "<<nRowCuts<<" Twomir cuts"<<std::endl;
      assert(cs.sizeRowCuts() > 0);
      OsiSolverInterface::ApplyCutsReturnCode rc = siP->applyCuts(cs);
      
      siP->resolve();
      
      double lpRelaxAfter= siP->getObjValue(); 
      std::cout<<"Initial LP value: "<<lpRelax<<std::endl;
      std::cout<<"LP value with cuts: "<<lpRelaxAfter<<std::endl;
      assert( lpRelax < lpRelaxAfter );
      assert(lpRelaxAfter < 964);
    }
    delete siP;
  }

}

void
CglCliqueUnitTest(const OsiSolverInterface *baseSiP,
			    const std::string mpsDir)
{
  // Test default constructor
  {
    CglClique aGenerator;
  }

  // Test copy & assignment
  {
    CglClique rhs;
    {
      CglClique bGenerator;
      CglClique cGenerator(bGenerator);
      //rhs=bGenerator;
    }
  }

  // Test get/set methods
  {
    CglClique getset;
    // None to test
  }

  // Test generateCuts
  {
    CglClique gct;
    OsiSolverInterface  *siP = baseSiP->clone();
    std::string fn = mpsDir+"l152lav";
    std::string fn2 = mpsDir+"l152lav.mps";
    FILE *in_f = fopen(fn2.c_str(), "r");
    if(in_f == NULL) {
      std::cout<<"Can not open file "<<fn2<<std::endl<<"Skip test of CglClique::generateCuts()"<<std::endl;
    }
    else {
      fclose(in_f);
      siP->readMps(fn.c_str(),"mps");

      siP->initialSolve();
      double lpRelax = siP->getObjValue();

      OsiCuts cs;
      gct.generateCuts(*siP, cs);
      int nRowCuts = cs.sizeRowCuts();
      std::cout<<"There are "<<nRowCuts<<" Clique cuts"<<std::endl;
      assert(cs.sizeRowCuts() > 0);
      OsiSolverInterface::ApplyCutsReturnCode rc = siP->applyCuts(cs);

      siP->resolve();

      double lpRelaxAfter= siP->getObjValue();
      std::cout<<"Initial LP value: "<<lpRelax<<std::endl;
      std::cout<<"LP value with cuts: "<<lpRelaxAfter<<std::endl;
      assert( lpRelax < lpRelaxAfter );
      assert(lpRelaxAfter < 4722.1);
    }
    delete siP;
  }

}

void
CglResidualCapacityUnitTest(const OsiSolverInterface *baseSiP,
			    const std::string mpsDir)
{
  // Test default constructor
  {
    CglResidualCapacity aGenerator;
  }
  
  // Test copy & assignment
  {
    CglResidualCapacity rhs;
    {
      CglResidualCapacity bGenerator;
      CglResidualCapacity cGenerator(bGenerator);
      rhs=bGenerator;
    }
  }

  // Test get/set methods
  {
    CglResidualCapacity getset;
    
    double geps = 10 * getset.getEpsilon();
    getset.setEpsilon(geps);
    double geps2 = getset.getEpsilon();
    assert(geps == geps2);

    double gtol = 10 * getset.getTolerance();
    getset.setTolerance(gtol);
    double gtol2 = getset.getTolerance();
    assert(gtol == gtol2);

    int gpre = getset.getDoPreproc();
    gpre = (gpre + 1) % 3 - 1;
    getset.setDoPreproc(gpre);
    int gpre2 = getset.getDoPreproc();
    assert(gpre == gpre2);
  }

  // Test generateCuts
  {
    CglResidualCapacity gct;
    OsiSolverInterface  *siP = baseSiP->clone();
    std::string fn = mpsDir+"capPlan1";
    std::string fn2 = mpsDir+"capPlan1.mps";
    FILE *in_f = fopen(fn2.c_str(), "r");
    if(in_f == NULL) {
      std::cout<<"Can not open file "<<fn2<<std::endl<<"Skip test of CglResidualCapacity::generateCuts()"<<std::endl;
    }
    else {
      fclose(in_f);
      siP->readMps(fn.c_str(),"mps");
 
      siP->initialSolve();
      double lpRelax = siP->getObjValue();
      
      OsiCuts cs;
      gct.setDoPreproc(1); // Needed for DyLP
      gct.generateCuts(*siP, cs);
      int nRowCuts = cs.sizeRowCuts();
      std::cout<<"There are "<<nRowCuts<<" Residual Capacity cuts"<<std::endl;
      assert(cs.sizeRowCuts() > 0);
      OsiSolverInterface::ApplyCutsReturnCode rc = siP->applyCuts(cs);
      
      siP->resolve();
      
      double lpRelaxAfter= siP->getObjValue(); 
      std::cout<<"Initial LP value: "<<lpRelax<<std::endl;
      std::cout<<"LP value with cuts: "<<lpRelaxAfter<<std::endl;
      assert( lpRelax < lpRelaxAfter );
      assert(lpRelaxAfter < 964);
    }
    delete siP;
  }

}

void
CglMixedIntegerRoundingUnitTest(const OsiSolverInterface *baseSiP,
			    const std::string mpsDir)
{
  // Test default constructor
  {
    CglMixedIntegerRounding aGenerator;
  }

  // Test copy & assignment
  {
    CglMixedIntegerRounding rhs;
    {
      CglMixedIntegerRounding bGenerator;
      CglMixedIntegerRounding cGenerator(bGenerator);
      rhs=bGenerator;
    }
  }

  // Test get/set methods
  {
    CglMixedIntegerRounding getset;

    int gagg = 10 * getset.getMAXAGGR_();
    getset.setMAXAGGR_(gagg);
    int gagg2 = getset.getMAXAGGR_();
    assert(gagg == gagg2);

    bool gmult = !getset.getMULTIPLY_();
    getset.setMULTIPLY_(gmult);
    bool gmult2 = getset.getMULTIPLY_();
    assert(gmult == gmult2);

    int gcrit = getset.getCRITERION_();
    gcrit = (gcrit) % 3 + 1;
    getset.setCRITERION_(gcrit);
    int gcrit2 = getset.getCRITERION_();
    assert(gcrit == gcrit2);

    int gpre = getset.getDoPreproc();
    gpre = (gpre + 1) % 3 - 1;
    getset.setDoPreproc(gpre);
    int gpre2 = getset.getDoPreproc();
    assert(gpre == gpre2);
  }

  // Test generateCuts
  {
    CglMixedIntegerRounding gct;
    OsiSolverInterface  *siP = baseSiP->clone();
    std::string fn = mpsDir+"capPlan1";
    std::string fn2 = mpsDir+"capPlan1.mps";
    FILE *in_f = fopen(fn2.c_str(), "r");
    if(in_f == NULL) {
      std::cout<<"Can not open file "<<fn2<<std::endl<<"Skip test of CglMixedIntegerRounding::generateCuts()"<<std::endl;
    }
    else {
      fclose(in_f);
      siP->readMps(fn.c_str(),"mps");

      siP->initialSolve();
      double lpRelax = siP->getObjValue();

      OsiCuts cs;
      gct.generateCuts(*siP, cs);
      int nRowCuts = cs.sizeRowCuts();
      std::cout<<"There are "<<nRowCuts<<" MIR cuts"<<std::endl;
      assert(cs.sizeRowCuts() > 0);
      OsiSolverInterface::ApplyCutsReturnCode rc = siP->applyCuts(cs);

      siP->resolve();

      double lpRelaxAfter= siP->getObjValue();
      std::cout<<"Initial LP value: "<<lpRelax<<std::endl;
      std::cout<<"LP value with cuts: "<<lpRelaxAfter<<std::endl;
      assert( lpRelax < lpRelaxAfter );
      assert(lpRelaxAfter < 964);
    }
    delete siP;
  }

}

//--------------------------------------------------------------------------
// test EKKsolution methods.
void
CglProbingUnitTest(
  const OsiSolverInterface * baseSiP,
  const std::string mpsDir )
{
# ifdef CGL_DEBUG
  int i ;	// define just once
# endif
  CoinRelFltEq eq(0.000001);

  // Test default constructor
  {
    CglProbing aGenerator;
  }
  
  // Test copy & assignment
  {
    CglProbing rhs;
    {
      CglProbing bGenerator;
      CglProbing cGenerator(bGenerator);
      rhs=bGenerator;
    }
  }

  {
    OsiCuts osicuts;
    CglProbing test1;
    OsiSolverInterface  * siP = baseSiP->clone();
    int nColCuts;
    int nRowCuts;
    
    std::string fn = mpsDir+"p0033";
    siP->readMps(fn.c_str(),"mps");
    siP->initialSolve();
    // just unsatisfied variables
    test1.generateCuts(*siP,osicuts);
    nColCuts = osicuts.sizeColCuts();
    nRowCuts = osicuts.sizeRowCuts();
    std::cout<<"There are "<<nRowCuts<<" probing cuts"<<std::endl;
    {
      std::cout<<"there are "<<nColCuts<<" probing column cuts"<<std::endl;

#ifdef CGL_DEBUG
      const double * lo = siP->getColLower();
      const double * up = siP->getColUpper();
      for (i=0; i<nColCuts; i++){
	OsiColCut ccut;
	CoinPackedVector cpv;
	ccut = osicuts.colCut(i);
	cpv = ccut.lbs();
	int n = cpv.getNumElements();
        int j;
	const int * indices = cpv.getIndices();
	double* elements = cpv.getElements();
	for (j=0;j<n;j++) {
	  int icol=indices[j];
	  if (elements[j]>lo[icol])
	    std::cout<<"Can increase lb on "<<icol<<" from "<<lo[icol]<<
	      " to "<<elements[j]<<std::endl;
	}
	cpv = ccut.ubs();
	n = cpv.getNumElements();
	indices = cpv.getIndices();
	elements = cpv.getElements();

	for (j=0;j<n;j++) {
	  int icol=indices[j];
	  if (elements[j]<up[icol])
	    std::cout<<"Can decrease ub on "<<icol<<" from "<<up[icol]<<
	      " to "<<elements[j]<<std::endl;
	}
      }

#endif

    }

#ifdef CGL_DEBUG
    for (i=0; i<nRowCuts; i++){
      OsiRowCut rcut;
      CoinPackedVector rpv;
      const double * colsol = siP->getColSolution();
      rcut = osicuts.rowCut(i);
      rpv = rcut.row();
      const int n = rpv.getNumElements();
      const int * indices = rpv.getIndices();
      double* elements = rpv.getElements();
      double sum2=0.0;
      int k=0;
      double lb=rcut.lb();
      double ub=rcut.ub();
      for (k=0; k<n; k++){
	int column=indices[k];
	sum2 += colsol[column]*elements[k];
      }
      if (sum2 >ub + 1.0e-7 ||sum2 < lb - 1.0e-7) {
	std::cout<<"Cut "<<i<<" lb "<<lb<<" solution "<<sum2<<" ub "<<ub<<std::endl;
//.........这里部分代码省略.........

//.........这里部分代码省略.........
	  std::cout<<"(col="<<column<<",el="<<elements[k]<<",sol="<<
	    colsol[column]<<") ";
	}
	std::cout <<std::endl;
      }
#endif

      if (i-nOldCuts==testCut) {
	assert( eq(rhs,ub));
	assert(n==1);
	for (k=0; k<n; k++){
	  int column=indices[k];
	  assert (eq(cut[column],elements[k]));
	}
	// add cut
	rowLower[3]=-1.0e100;
	rowUpper[3]=ub;
	matrix.appendRow(rpv);
      }
    }
    nOldCuts=nRowCuts;
    // basis 2
    int rowBasis2[4]={1,1,-1,-1};
    int colBasis2[2]={1,1};
    warm.setSize(2,4);
    for (i=0;i<4;i++) {
      if (rowBasis2[i]<0) {
	warm.setArtifStatus(i,CoinWarmStartBasis::atLowerBound);
      } else {
	warm.setArtifStatus(i,CoinWarmStartBasis::basic);
      }
    }
    for (i=0;i<2;i++) {
      if (colBasis2[i]<0) {
	warm.setStructStatus(i,CoinWarmStartBasis::atLowerBound);
      } else {
	warm.setStructStatus(i,CoinWarmStartBasis::basic);
      }
    }

    // solution 2
    double colsol2[2]={2.0,0.5};
    test1.generateCuts(NULL, osicuts, matrix,
		 /*objective,*/ colsol2,
		 colLower, colUpper,
		 rowLower, rowUpper, intVar, &warm);
    nRowCuts = osicuts.sizeRowCuts();
    std::cout<<"There are "<<nRowCuts<<" gomory cuts"<<std::endl;
    assert (nRowCuts==nOldCuts);
    
  }

  // Miplib3 problem p0033
  if (1) {
    // Setup
    OsiSolverInterface  * siP = baseSiP->clone();
    std::string fn(mpsDir+"p0033");
    siP->readMps(fn.c_str(),"mps");
    siP->activateRowCutDebugger("p0033");
    CglGomory test;

    // Solve the LP relaxation of the model and
    // print out ofv for sake of comparison 
    siP->initialSolve();
    double lpRelaxBefore=siP->getObjValue();
    std::cout<<"Initial LP value: "<<lpRelaxBefore<<std::endl;
    assert( eq(lpRelaxBefore, 2520.5717391304347) );

    // Fails with OsiCpx, OsiXpr:
    /**********
    double mycs[] = {0, 1, 0, 0, -2.0837010502455788e-19, 1, 0, 0, 1,
		       0.021739130434782594, 0.35652173913043478, 
		       -6.7220534694101275e-18, 5.3125906451789717e-18, 
		       1, 0, 1.9298798670241979e-17, 0, 0, 0,
		       7.8875708048320448e-18, 0.5, 0, 
		       0.85999999999999999, 1, 1, 0.57999999999999996,
		       1, 0, 1, 0, 0.25, 0, 0.67500000000000004};
    siP->setColSolution(mycs);
    ****/

    OsiCuts cuts;    
    
    // Test generateCuts method
    test.generateCuts(*siP,cuts);
    int nRowCuts = cuts.sizeRowCuts();
    std::cout<<"There are "<<nRowCuts<<" Gomory cuts"<<std::endl;
    assert(cuts.sizeRowCuts() > 0);
    OsiSolverInterface::ApplyCutsReturnCode rc = siP->applyCuts(cuts);
    
    siP->resolve();
    double lpRelaxAfter=siP->getObjValue(); 
    std::cout<<"LP value with cuts: "<<lpRelaxAfter<<std::endl;
    //assert( eq(lpRelaxAfter, 2592.1908295194507) );
    assert( lpRelaxAfter> 2550.0 );
    assert( lpRelaxBefore < lpRelaxAfter );
    assert(lpRelaxAfter < 3089.1);
    
    delete siP;
  } 
}

void
CglRedSplitUnitTest(const OsiSolverInterface *baseSiP,
		    const std::string mpsDir)
{
  // Test default constructor
  {
    CglRedSplit aGenerator;
  }
  
  // Test copy & assignment
  {
    CglRedSplit rhs;
    {
      CglRedSplit bGenerator;
      CglRedSplit cGenerator(bGenerator);
      rhs=bGenerator;
    }
  }

  // Test get/set methods
  {
    CglRedSplit getset;
    CglRedSplitParam gsparam = getset.getParam();
    
    double geps = 10 * gsparam.getEPS();
    gsparam.setEPS(geps);
    double geps2 = gsparam.getEPS();
    assert(geps == geps2);

    double gepse = 10 * gsparam.getEPS_ELIM();
    gsparam.setEPS_ELIM(gepse);
    double gepse2 = gsparam.getEPS_ELIM();
    assert(gepse == gepse2);

    double gmv = 10 * gsparam.getMINVIOL();
    gsparam.setMINVIOL(gmv);
    double gmv2 = gsparam.getMINVIOL();
    assert(gmv == gmv2);

    int gucg = gsparam.getUSE_CG2();
    gucg = 1 - gucg;
    gsparam.setUSE_CG2(gucg);
    int gucg2 = gsparam.getUSE_CG2();
    assert(gucg == gucg2);
  }

  // Test generateCuts
  {
    CglRedSplit gct;
    OsiSolverInterface  *siP = baseSiP->clone();
    std::string fn = mpsDir+"p0033";
    std::string fn2 = mpsDir+"p0033.mps";
    FILE *in_f = fopen(fn2.c_str(), "r");
    if(in_f == NULL) {
      std::cout<<"Can not open file "<<fn2<<std::endl<<"Skip test of CglRedSplit::generateCuts()"<<std::endl;
    }
    else {
      fclose(in_f);
      siP->readMps(fn.c_str(),"mps");
 
      siP->initialSolve();
      double lpRelax = siP->getObjValue();
      
      OsiCuts cs;
      gct.getParam().setMAX_SUPPORT(34);
      gct.getParam().setUSE_CG2(1);
      //      gct.getParam().setUSE_CG2(1);
      gct.generateCuts(*siP, cs);
      int nRowCuts = cs.sizeRowCuts();
      std::cout<<"There are "<<nRowCuts<<" Reduce-and-Split cuts"<<std::endl;
      assert(cs.sizeRowCuts() > 0);
      OsiSolverInterface::ApplyCutsReturnCode rc = siP->applyCuts(cs);
      
      siP->resolve();
      
      double lpRelaxAfter= siP->getObjValue(); 
      std::cout<<"Initial LP value: "<<lpRelax<<std::endl;
      std::cout<<"LP value with cuts: "<<lpRelaxAfter<<std::endl;
      assert( lpRelax < lpRelaxAfter );
      assert(lpRelaxAfter < 3089.1);
    }
    delete siP;
  }

}

//.........这里部分代码省略.........
            {
                OsiCuts cuts;
                test.generateCuts(*siP,cuts);
                cuts.printCuts();
                assert(cuts.sizeRowCuts()==1);
                OsiRowCut aCut = cuts.rowCut(0);
                CoinPackedVector row = aCut.row();
                if (row.getNumElements() == 1)
                {
                    assert(row.getIndices()[0]==1);
                    assert(eq(row.getElements()[0], -1));
                }
                else if (row.getNumElements() == 2)
                {
                    assert(row.getIndices()[0]==0);
                    assert(eq(row.getElements()[0], 0.));
                    assert(row.getIndices()[1]==1);
                    assert(eq(row.getElements()[1], -1));
                }
                assert(eq(aCut.lb(), 0.));
                OsiSolverInterface::ApplyCutsReturnCode rc = siP->applyCuts(cuts);

                siP->resolve();
            }
            delete siP;
        }
    }

    if (1)  //Test on p0033
    {
        // Setup
        OsiSolverInterface  * siP = si->clone();
        std::string fn(mpsDir+"p0033");
        siP->readMps(fn.c_str(),"mps");
        siP->activateRowCutDebugger("p0033");
        CglLandP test;

        // Solve the LP relaxation of the model and
        // print out ofv for sake of comparison
        siP->initialSolve();
        double lpRelaxBefore=siP->getObjValue();
        assert( eq(lpRelaxBefore, 2520.5717391304347) );
#ifdef CGL_DEBUG
        printf("\n\nOrig LP min=%f\n",lpRelaxBefore);
#endif

        OsiCuts cuts;

        // Test generateCuts method
        test.generateCuts(*siP,cuts);
        OsiSolverInterface::ApplyCutsReturnCode rc = siP->applyCuts(cuts);

        siP->resolve();
        double lpRelaxAfter=siP->getObjValue();
        //assert( eq(lpRelaxAfter, 2592.1908295194507) );

        std::cout<<"Relaxation after "<<lpRelaxAfter<<std::endl;
        assert( lpRelaxAfter> 2840. );
#ifdef CGL_DEBUG
        printf("\n\nOrig LP min=%f\n",lpRelaxBefore);
        printf("\n\nFinal LP min=%f\n",lpRelaxAfter);
#endif
        assert( lpRelaxBefore < lpRelaxAfter );

        delete siP;
    }

//.........这里部分代码省略.........

  // internal consistency
  {
    const int ne = 1;
    int inx[ne] = { -3 };
    double el[ne] = { 1.2 };
    OsiColCut r;
    r.setLbs(ne, inx, el);
    OSIUNITTEST_ASSERT_ERROR(!r.consistent(), {}, "osicolcut", "consistent");
  }
  {
    const int ne = 1;
    int inx[ne] = { -3 };
    double el[ne] = { 1.2 };
    OsiColCut r;
    r.setUbs(ne, inx, el);
    OSIUNITTEST_ASSERT_ERROR(!r.consistent(), {}, "osicolcut", "consistent");
  }
  {
    const int ne = 1;
    int inx[ne] = { 100 };
    double el[ne] = { 1.2 };
    const int ne1 = 2;
    int inx1[ne1] = { 50, 100 };
    double el1[ne1] = { 100., 100. };
    OsiColCut r;
    r.setUbs(ne, inx, el);
    r.setLbs(ne1, inx1, el1);
    OSIUNITTEST_ASSERT_ERROR(r.consistent(), {}, "osicolcut", "consistent");

    OsiSolverInterface *imP = baseSiP->clone();
    assert(imP != NULL);
    std::string fn = mpsDir + "exmip1";
    imP->readMps(fn.c_str(), "mps");
    OSIUNITTEST_ASSERT_ERROR(!r.consistent(*imP), {}, "osicolcut", "consistent");
    delete imP;
  }
  {
    const int ne = 1;
    int inx[ne] = { 100 };
    double el[ne] = { 1.2 };
    const int ne1 = 2;
    int inx1[ne1] = { 50, 100 };
    double el1[ne1] = { 100., 1. };
    OsiColCut r;
    r.setUbs(ne, inx, el);
    r.setLbs(ne1, inx1, el1);
    OSIUNITTEST_ASSERT_ERROR(r.consistent(), {}, "osicolcut", "consistent");
  }
  {
    // Test consistent(IntegerModel) method.
    OsiSolverInterface *imP = baseSiP->clone();
    assert(imP != NULL);
    std::string fn = mpsDir + "exmip1";
    imP->readMps(fn.c_str(), "mps");

    OsiColCut cut;
    const int ne = 1;
    int inx[ne] = { 20 };
    double el[ne] = { 0.25 };
    cut.setLbs(ne, inx, el);
    OSIUNITTEST_ASSERT_ERROR(!cut.consistent(*imP), {}, "osicolcut", "consistent(IntegerModel)");

    cut.setLbs(0, NULL, NULL);
    cut.setUbs(ne, inx, el);
    OSIUNITTEST_ASSERT_ERROR(!cut.consistent(*imP), {}, "osicolcut", "consistent(IntegerModel)");

//Solver function
int sci_rmps(char *fname) 
{
    //creating a problem pointer using base class of OsiSolverInterface and
    //instantiate the object using derived class of ClpSolverInterface
    OsiSolverInterface* si = new OsiClpSolverInterface();

    // Error management variable
	SciErr sciErr;

	//data declarations
	int *piAddressVarOne = NULL;                 //pointer used to access argument of the function
	char* ptr;                              	 //pointer to point to address of file name
    double* options_;                            //options to set maximum iterations 
	CheckInputArgument(pvApiCtx, 2,2 );          //Check we have exactly two arguments as input or not
	CheckOutputArgument(pvApiCtx, 6, 6);         //Check we have exactly six arguments on output side or not
    //Getting the input arguments from Scilab
    //Getting the MPS file path
	//Reading mps file
	getStringFromScilab(1,&ptr);

 	std::cout<<ptr;
	
    //get options from Scilab
    if(getFixedSizeDoubleMatrixInList(2 , 2 , 1 , 1 , &options_))
	{
		return 1;
	}

    //Read the MPS file
    si->readMps(ptr);

    //setting options for maximum iterations
    si->setIntParam(OsiMaxNumIteration,options_[0]);

    //Solve the problem
    si->initialSolve();
  
    //Quering about the problem
    //get number of variables
    double numVars_;
    numVars_ = si->getNumCols();
  
    //get number of constraint equations
    double numCons_;
    numCons_ = si->getNumRows();
   
    //Output the solution to Scilab
    //get solution for x
    const double* xValue = si->getColSolution();
   
    //get objective value
    double objValue = si->getObjValue();

    //get Status value
    double status;
    if(si->isProvenOptimal())
    	status=0;
    else if(si->isProvenPrimalInfeasible())
    	status=1;
    else if(si->isProvenDualInfeasible())
        status=2;
    else if(si->isIterationLimitReached())
        status=3;
   	else if(si->isAbandoned())
        status=4;
   	else if(si->isPrimalObjectiveLimitReached())
        status=5;
   	else if(si->isDualObjectiveLimitReached())
        status=6;

    //get number of iterations
    double iterations = si->getIterationCount();

    //get reduced cost 
    const double* reducedCost = si->getReducedCost();
   
    //get dual vector
    const double* dual = si->getRowPrice();
  
    returnDoubleMatrixToScilab(1 , 1 , numVars_ , xValue);
    returnDoubleMatrixToScilab(2 , 1 , 1 , &objValue);
    returnDoubleMatrixToScilab(3 , 1 , 1 , &status);
    returnDoubleMatrixToScilab(4 , 1 , 1 , &iterations);
    returnDoubleMatrixToScilab(5 , 1 , numVars_ , reducedCost);
    returnDoubleMatrixToScilab(6 , 1 , numCons_ , dual);
	
	free(xValue);
	free(dual);
	free(reducedCost);
}

//--------------------------------------------------------------------------
void
CglKnapsackCoverUnitTest(
  const OsiSolverInterface * baseSiP,
  const std::string mpsDir )
{
  int i;
  CoinRelFltEq eq(0.000001);

  // Test default constructor
  {
    CglKnapsackCover kccGenerator;
  }
  
  // Test copy & assignment
  {
    CglKnapsackCover rhs;
    {
      CglKnapsackCover kccGenerator;
      CglKnapsackCover cgC(kccGenerator);
      rhs=kccGenerator;
    }
  }


  // test exactSolveKnapsack
  {  
    CglKnapsackCover kccg;
    const int n=7;
    double c=50;
    double p[n] = {70,20,39,37,7,5,10};
    double w[n] = {31, 10, 20, 19, 4, 3, 6};
    double z;
    int x[n];
    int exactsol = kccg.exactSolveKnapsack(n, c, p, w, z, x);
    assert(exactsol==1);
    assert (z == 107);
    assert (x[0]==1);
    assert (x[1]==0);
    assert (x[2]==0);
    assert (x[3]==1);
    assert (x[4]==0);
    assert (x[5]==0);
    assert (x[6]==0);
  }

  /*
  // Testcase /u/rlh/osl2/mps/scOneInt.mps
  // Model has 3 continous, 2 binary, and 1 general
  // integer variable.
  {
    OsiSolverInterface  * siP = baseSiP->clone();
    int * complement=NULL;
    double * xstar=NULL;

    siP->readMps("../Mps/scOneInt","mps");
    CglKnapsackCover kccg;
    int nCols=siP->getNumCols();
    
    // Test the siP methods for detecting
    // variable type
    int numCont=0, numBinary=0, numIntNonBinary=0, numInt=0;
    for (int thisCol=0; thisCol<nCols; thisCol++) {
      if ( siP->isContinuous(thisCol) ) numCont++;
      if ( siP->isBinary(thisCol) ) numBinary++;
      if ( siP->isIntegerNonBinary(thisCol) ) numIntNonBinary++;
      if ( siP->isInteger(thisCol) ) numInt++;
    }
    assert(numCont==3);
    assert(numBinary==2);
    assert(numIntNonBinary==1);
    assert(numInt==3);
    
    
    // Test initializeCutGenerator
    siP->initialSolve();
    assert(xstar !=NULL);
    for (i=0; i<nCols; i++){
      assert(complement[i]==0);
    }
    int nRows=siP->getNumRows();
    for (i=0; i<nRows; i++){
    int vectorsize = siP->getMatrixByRow()->vectorSize(i);
    assert(vectorsize==2);
    }
    
    kccg.cleanUpCutGenerator(complement,xstar);
    delete siP;
  }
  */  
  
  // Testcase /u/rlh/osl2/mps/tp3.mps
  // Models has 3 cols, 3 rows
  // Row 0 yields a knapsack, others do not.
  {
    // setup
    OsiSolverInterface  * siP = baseSiP->clone();
    std::string fn(mpsDir+"tp3");
    siP->readMps(fn.c_str(),"mps");     
    // All integer variables should be binary.
//.........这里部分代码省略.........