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

//.........这里部分代码省略.........
    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
  {
    CglSimpleRounding cg;
    
    OsiSolverInterface * siP = baseSiP->clone();
    std::string fn = mpsDir+"p0033";
    siP->readMps(fn.c_str(),"mps");
    OsiCuts cuts;
    cg.generateCuts(*siP,cuts);

    // p0033 is the optimal solution to p0033
    int objIndices[14] = { 
       0,  6,  7,  9, 13, 17, 18,
      22, 24, 25, 26, 27, 28, 29 };
    CoinPackedVector p0033(14,objIndices,1.0);

    // test that none of the generated cuts
    // chops off the optimal solution
    int nRowCuts = cuts.sizeRowCuts();
    OsiRowCut rcut;
    CoinPackedVector rpv;
    int i;
    for (i=0; i<nRowCuts; i++){
      rcut = cuts.rowCut(i);
      rpv = rcut.row();
      double p0033Sum = (rpv*p0033).sum();
      double rcutub = rcut.ub();
      assert (p0033Sum <= rcutub);
    }

    // test that the cuts improve the 
    // lp objective function value
    siP->initialSolve();
    double lpRelaxBefore=siP->getObjValue();
    OsiSolverInterface::ApplyCutsReturnCode rc = siP->applyCuts(cuts);
    siP->resolve();
    double lpRelaxAfter=siP->getObjValue(); 
#ifdef CGL_DEBUG
    printf("\n\nOrig LP min=%f\n",lpRelaxBefore);
    printf("Final LP min=%f\n\n",lpRelaxAfter);
#endif
    assert( lpRelaxBefore < lpRelaxAfter );

    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;
  }

}

//.........这里部分代码省略.........
   mat.deleteRows(numNonLinear, nonLinear());
   
   int numcols = nlp.getNumCols();
   Bonmin::vector<double> obj(numcols);
   for (int i = 0; i < numcols; i++)
     obj[i] = 0.;
   
   si.loadProblem(mat, nlp.getColLower(), nlp.getColUpper(), 
                  obj(), rowLow(), rowUp());
   const Bonmin::TMINLP::VariableType* variableType = problem->var_types();
   for (int i = 0; i < n; i++) {
     if ((variableType[i] == Bonmin::TMINLP::BINARY) || (variableType[i] == Bonmin::TMINLP::INTEGER))
       si.setInteger(i);
   }
   if (getObj) {
     bool addObjVar = false;
     if (problem->hasLinearObjective()) {
       double zero;
       Bonmin::vector<double> x0(n, 0.);
       problem->eval_f(n, x0(), 1, zero);
       si.setDblParam(OsiObjOffset, -zero);
       //Copy the linear objective and don't create a dummy variable.
       problem->eval_grad_f(n, x, 1, obj());
       si.setObjective(obj());
     } else {
       addObjVar = true;
     }
   
     if (addObjVar) {
       nlp.addObjectiveFunction(si, x);
     }
   }
   
   // Hassan IA initial description
   int InnerDesc = 1;
   if (InnerDesc == 1) {
     OsiCuts cs;
   
     double * p = CoinCopyOfArray(colLower, n);
     double * pp = CoinCopyOfArray(colLower, n);
     double * up = CoinCopyOfArray(colUpper, n);
   
     for (int i = 0; i < n; i++) {
       if (p[i] < -1e3){
          p[i] = pp[i] = -1e3;
       }
       if (up[i] > 1e2){
          up[i] = 1e2;
       }
     } 

     const int& nbAp = nbAp_;
     printf("Generating approximation with %i points.\n", nbAp);
   
     std::vector<double> step(n);
     int n_lin = 0;
   
     for (int i = 0; i < n; i++) {
       //if ((variableType[i] == Bonmin::TMINLP::BINARY) || (variableType[i] == Bonmin::TMINLP::INTEGER)) {
       if (varTypes[i] == Ipopt::TNLP::LINEAR) {
         n_lin ++;
         step[i] = 0;
         p[i] = pp[i] = up[i] = 0;
       }
       else {
         step[i] = (up[i] - p[i]) / (nbAp);
       }
     }
     printf("Number of linears %i\n", n_lin);
   
     for (int j = 1; j < nbAp; j++) {
   
       for (int i = 0; i < n; i++) {
         pp[i] += step[i];
       }
   
       for (int i = 0; (i < m ); i++) {
         if (constTypes[i] == Ipopt::TNLP::LINEAR) continue;
         bool status = getMyInnerApproximation(nlp, cs, i, p, pp);// Generate a chord connecting the two points
         if(status == false){
           printf("Error in generating inner approximation\n");
           exit(1);
         }
       }
       std::copy(pp, pp+n, p);
      
     }
   
     for(int i = 0; (i< m); i++) {
         if (constTypes[i] == Ipopt::TNLP::LINEAR) continue;
         getMyInnerApproximation(nlp, cs, i, p, up);// Generate a chord connecting the two points
     }

        delete [] p; 
        delete [] pp;
        delete [] up; 
     si.applyCuts(cs);
   }
   printf("************  Done extracting inner approx ********");
  }

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
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;
  }

}

//.........这里部分代码省略.........
   vector<double> obj(numcols);
   for (int i = 0; i < numcols; i++)
     obj[i] = 0.;
   
   si.loadProblem(mat, nlp.getColLower(), nlp.getColUpper(), 
                  obj(), rowLow(), rowUp());
   const Bonmin::TMINLP::VariableType* variableType = problem->var_types();
   for (int i = 0; i < n; i++) {
     if ((variableType[i] == TMINLP::BINARY) || (variableType[i]
         == TMINLP::INTEGER))
       si.setInteger(i);
   }
   if (getObj) {
     bool addObjVar = false;
     if (problem->hasLinearObjective()) {
       double zero;
       vector<double> x0(n, 0.);
       problem->eval_f(n, x0(), 1, zero);
       si.setDblParam(OsiObjOffset, -zero);
       //Copy the linear objective and don't create a dummy variable.
       problem->eval_grad_f(n, x, 1, obj());
       si.setObjective(obj());
     } else {
       addObjVar = true;
     }
   
     if (addObjVar) {
       nlp.addObjectiveFunction(si, x);
     }
   }
   
   // Hassan IA initial description
   int InnerDesc = 1;
   if (InnerDesc == 1) {
     OsiCuts cs;
   
     double * p = CoinCopyOfArray(colLower, n);
     double * pp = CoinCopyOfArray(colLower, n);
     double * up = CoinCopyOfArray(colUpper, n);
   
     const int& nbAp = nbAp_;
     std::vector<int> nbG(m, 0);// Number of generated points for each nonlinear constraint
   
     std::vector<double> step(n);
   
     for (int i = 0; i < n; i++) {
   
       if (colUpper[i] > 1e08) {
         up[i] = 0;
       }
   
       if (colUpper[i] > 1e08 || colLower[i] < -1e08 || (variableType[i]
           == TMINLP::BINARY) || (variableType[i] == TMINLP::INTEGER)) {
         step[i] = 0;
       } else
         step[i] = (up[i] - colLower[i]) / (nbAp);
   
       if (colLower[i] < -1e08) {
         p[i] = 0;
         pp[i] = 0;
       }
   
     }
     vector<double> g_p(m);
     vector<double> g_pp(m);
   
     for (int j = 1; j <= nbAp; j++) {
   
       for (int i = 0; i < n; i++) {
         pp[i] += step[i];
       }
   
       problem->eval_g(n, p, 1, m, g_p());
       problem->eval_g(n, pp, 1, m, g_pp());
       double diff = 0;
       int varInd = 0;
       for (int i = 0; (i < m && constTypes[i] == Ipopt::TNLP::NON_LINEAR); i++) {
         if (varInd == n - 1)
           varInd = 0;
         diff = std::abs(g_p[i] - g_pp[i]);
         if (nbG[i] < nbAp - 1) {
           getMyInnerApproximation(nlp, cs, i, p, pp);// Generate a chord connecting the two points
           p[varInd] = pp[varInd];
           nbG[i]++;
         }
         varInd++;
       }
     }
   
     for(int i = 0; (i< m && constTypes[i] == Ipopt::TNLP::NON_LINEAR); i++) {
      //  getConstraintOuterApproximation(cs, i, colUpper, NULL, true);// Generate Tangents at current point
         getMyInnerApproximation(nlp, cs, i, p, up);// Generate a chord connecting the two points
     }

        delete [] p; 
        delete [] pp;
        delete [] up; 
     si.applyCuts(cs);
   }
  }

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;
  }

}

//.........这里部分代码省略.........
	  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;
  }

}

//.........这里部分代码省略.........
            double colLower[2]={0.0,0.0};
            double colUpper[2]={1.0,1.0};
            double obj[2]={-1,-1};
            int intVar[2]={0,1};

            OsiSolverInterface  * siP = si->clone();
            siP->loadProblem(columnCopy, colLower, colUpper, obj, rowLower, rowUpper);
            siP->setInteger(intVar,2);
            CglLandP test;
            test.setLogLevel(2);
            test.parameter().sepSpace = CglLandP::Full;
            siP->resolve();
            // Test generateCuts method
            {
                OsiCuts cuts;
                test.generateCuts(*siP,cuts);
                cuts.printCuts();
                assert(cuts.sizeRowCuts()==1);
                OsiRowCut aCut = cuts.rowCut(0);
                assert(eq(aCut.lb(), -.0714286));
                CoinPackedVector row = aCut.row();
                if (row.getNumElements() == 1)
                {
                    assert(row.getIndices()[0]==1);
                    assert(eq(row.getElements()[0], -4*.0714286));
                }
                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));
                }
                OsiSolverInterface::ApplyCutsReturnCode rc = siP->applyCuts(cuts);

                siP->resolve();
            }
            if (0)
            {
                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;
        }
    }

//.........这里部分代码省略.........
    kccg.ekk_validateIntType(si);     
    
    // Solve the LP relaxation of the model and
    // print out ofv for sake of comparison 
    ekk_allSlackBasis(model);
    ekk_crash(model,1); 
    ekk_primalSimplex(model,1);
    double lpRelaxBefore=ekk_getRobjvalue(model);
#ifdef CGL_DEBUG
    printf("\n\nOrig LP min=%f\n",lpRelaxBefore);
#endif
    
    // Determine if lp sol is ip optimal
    // Note: no ekk_function to do this
    int nCols=ekk_getInumcols(model);
    double * optLpSol = ekk_colsol(model);
    int ipOpt = 1;
    i=0;
    while (i++<nCols && ipOpt){
      if(optLpSol[i] < 1.0-1.0e-08 && optLpSol[i]> 1.0e-08) ipOpt = 0;
    }
    
    if (ipOpt){
#ifdef CGL_DEBUG
      printf("Lp solution is within ip optimality tolerance\n");
#endif
    }    
    else {
      OsiSolverInterface iModel(model);
      OsiCuts cuts;    
      
      // Test generateCuts method
      kccg.generateCuts(iModel,cuts);
      OsiSolverInterface::ApplyCutsReturnCode rc = iModel.applyCuts(cuts);
      
      ekk_mergeBlocks(model,1);         
      ekk_dualSimplex(model);
      double lpRelaxAfter=ekk_getRobjvalue(model); 
#ifdef CGL_DEBUG
      printf("\n\nFinal LP min=%f\n",lpRelaxAfter);
#endif
      assert( lpRelaxBefore < lpRelaxAfter );
      
      // This may need to be updated as other 
      // minimal cover finders are added
      assert( cuts.sizeRowCuts() == 1 );
      OsiRowCut testRowCut = cuts.rowCut(0);
      CoinPackedVector testRowPV = testRowCut.row();
      
      OsiRowCut sampleRowCut;
      const int sampleSize = 6;
      int sampleCols[sampleSize]={0,1,2,3,4,5};
      double sampleElems[sampleSize]={1.0,1.0,1.0,1.0,0.5, 2.0};
      sampleRowCut.setRow(sampleSize,sampleCols,sampleElems);
      sampleRowCut.setLb(-DBL_MAX);
      sampleRowCut.setUb(3.0);
      bool equiv = testRowPV.equivalent(sampleRowCut.row(),CoinRelFltEq(1.0e-05) );
      assert( testRowPV.equivalent(sampleRowCut.row(),CoinRelFltEq(1.0e-05) ) );
    }
    
    // Exit out of OSL
    ekk_deleteModel(model);
    ekk_endContext(env);
    
  }
#endif

//.........这里部分代码省略.........
      const double * duals = nlp->getRowPrice() + 2 *n;
      vector<double> grad(n, 0); 
      vector<int> indices(n, 0);
      tminlp->eval_grad_f(n, x_sol, false, grad());
      for(int i = 0 ; i < m ; i++){
        if(c_lin[i] == Ipopt::TNLP::LINEAR) continue;
        int nnz;
        tminlp->eval_grad_gi(n, x_sol, false, i, nnz, indices(), NULL);  
        tminlp->eval_grad_gi(n, x_sol, false, i, nnz, NULL, grad());
        for(int k = 0 ; k < nnz ; k++){
          objective[indices[k]] += alpha *duals[i] * grad[k];
        } 
      }
      for(int i = 0 ; i < n ; i++){
         if(variableType[i] != Bonmin::TMINLP::CONTINUOUS)
         objective[i] += alpha * grad[i];
         //if(fabs(objective[i]) < 1e-4) objective[i] = 0;
         else objective[i] = 0;
      }
      std::copy(objective.begin(), objective.end(), std::ostream_iterator<double>(std::cout, " "));
      std::cout<<std::endl;
#endif

      // load the problem to OSI
      OsiSolverInterface *si = mip_->solver();
      assert(si != NULL);
      CoinMessageHandler * handler = model_->messageHandler()->clone();
      si->passInMessageHandler(handler);
      si->messageHandler()->setLogLevel(1);

      si->loadProblem(matrix, model_->solver()->getColLower(), model_->solver()->getColUpper(), objective(), 
                      newRowLower(), newRowUpper());
      si->setInteger(idxIntegers(), static_cast<int>(idxIntegers.size()));
      si->applyCuts(noGoods);

      bool hasFractionnal = true;
      while(hasFractionnal){
        mip_->optimize(DBL_MAX, 0, 60);
        hasFractionnal = false;
#if 0
        bool feasible = false;
        if(mip_->getLastSolution()) {
  	const double* solution = mip_->getLastSolution();
          std::copy(solution, solution + n, newSolution.begin());
  	feasible = true;
  
        }

    if(feasible) {
      // fix the integer variables and solve the NLP
      // also add no good cut
      CoinPackedVector v;
      double lb = 1;
      for (int iColumn=0;iColumn<n;iColumn++) {
	if (variableType[iColumn] != Bonmin::TMINLP::CONTINUOUS) {
	  double value=newSolution[iColumn];
	  if (fabs(floor(value+0.5)-value)>integerTolerance) {
#ifdef DEBUG_BON_HEURISTIC_DIVE_MIP
	    cout<<"It should be infeasible because: "<<endl;
	    cout<<"variable "<<iColumn<<" is not integer"<<endl;
#endif
	    feasible = false;
	    break;
	  }
	  else {
	    value=floor(newSolution[iColumn]+0.5);