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

//#############################################################################
void 
MibSHeuristic::objCutHeuristic()
{

  /* Solve the LP relaxation with the new constraint d^2 y <= d^y* */

  MibSModel * model = MibSModel_;

  //OsiSolverInterface * oSolver = model->origLpSolver_;
  OsiSolverInterface * oSolver = model->getSolver();
  //OsiSolverInterface * hSolver = new OsiCbcSolverInterface();
  OsiSolverInterface * hSolver = new OsiSymSolverInterface();

  double objSense(model->getLowerObjSense());  
  int lCols(model->getLowerDim());
  int uCols(model->getUpperDim());
  int tCols(lCols + uCols);
  int * lColIndices = model->getLowerColInd();
  int * uColIndices = model->getUpperColInd();
  double * lObjCoeffs = model->getLowerObjCoeffs();

  hSolver->loadProblem(*oSolver->getMatrixByCol(),
		       oSolver->getColLower(), oSolver->getColUpper(),
		       oSolver->getObjCoefficients(),
		       oSolver->getRowLower(), oSolver->getRowUpper());

  int j(0);
  for(j = 0; j < tCols; j++){
    if(oSolver->isInteger(j))
      hSolver->setInteger(j);
  }

  double * optLowerSolutionOrd = model->bS_->optLowerSolutionOrd_;

  CoinPackedVector objCon;
  int i(0), index(0);
  double rhs(0.0);

  for(i = 0; i < lCols; i++){
    index = lColIndices[i];
    objCon.insert(index, lObjCoeffs[i] * objSense);
    //should this be ordered? and should lObjCoeffs by at index?
    //rhs += optLowerSolutionOrd_[i] * lObjCoeffs[i] * objSense;
    rhs += optLowerSolutionOrd[i] * lObjCoeffs[i] * objSense;
  }

  //Hmm, I think this was wrong before...?
  //  hSolver->addRow(objCon, - hSolver->getInfinity(), rhs);
  hSolver->addRow(objCon, rhs, hSolver->getInfinity());

  /* optimize w.r.t. to the UL objective with the new row */
  if(0){
    dynamic_cast<OsiCbcSolverInterface *> 
      (hSolver)->getModelPtr()->messageHandler()->setLogLevel(0);
  }
  else{
    dynamic_cast<OsiSymSolverInterface *> 
      (hSolver)->setSymParam("prep_level", -1);
    
    dynamic_cast<OsiSymSolverInterface *> 
      (hSolver)->setSymParam("verbosity", -2);

    dynamic_cast<OsiSymSolverInterface *> 
      (hSolver)->setSymParam("max_active_nodes", 1);
  }

  hSolver->branchAndBound();

  if(0)
    hSolver->writeLp("objcutheuristic");

  if(hSolver->isProvenOptimal()){

    MibSSolution *mibSol = NULL;

    OsiSolverInterface * lSolver = model->bS_->setUpModel(hSolver, true);

    if(0){
       dynamic_cast<OsiCbcSolverInterface *> 
	  (lSolver)->getModelPtr()->messageHandler()->setLogLevel(0);
    }    
    else{
       dynamic_cast<OsiSymSolverInterface *> 
	  (lSolver)->setSymParam("prep_level", -1);
       
       dynamic_cast<OsiSymSolverInterface *> 
	  (lSolver)->setSymParam("verbosity", -2);
       
       dynamic_cast<OsiSymSolverInterface *> 
	  (lSolver)->setSymParam("max_active_nodes", 1);
    }

    lSolver->branchAndBound();

    const double * sol = hSolver->getColSolution();
    double objVal(lSolver->getObjValue() * objSense);
    double etol(etol_);
    double lowerObj = getLowerObj(sol, objSense);  
    
//.........这里部分代码省略.........

void 
HeuristicInnerApproximation::extractInnerApproximation(Bonmin::OsiTMINLPInterface & nlp, OsiSolverInterface &si,
  const double * x, bool getObj) {
   printf("************  Start extracting inner approx");
   int n;
   int m;
   int nnz_jac_g;
   int nnz_h_lag;
   Ipopt::TNLP::IndexStyleEnum index_style;
   Bonmin::TMINLP2TNLP * problem = nlp.problem(); 
   //Get problem information
   problem->get_nlp_info(n, m, nnz_jac_g, nnz_h_lag, index_style);
   
   Bonmin::vector<int> jRow(nnz_jac_g);
   Bonmin::vector<int> jCol(nnz_jac_g);
   Bonmin::vector<double> jValues(nnz_jac_g);
   problem->eval_jac_g(n, NULL, 0, m, nnz_jac_g, jRow(), jCol(), NULL);
   if(index_style == Ipopt::TNLP::FORTRAN_STYLE)//put C-style
   {
     for(int i = 0 ; i < nnz_jac_g ; i++){
       jRow[i]--;
       jCol[i]--;
     }
   }
   
   //get Jacobian
   problem->eval_jac_g(n, x, 1, m, nnz_jac_g, NULL, NULL,
       jValues());
   
   Bonmin::vector<double> g(m);
   problem->eval_g(n, x, 1, m, g());
   
   Bonmin::vector<int> nonLinear(m);
   //store non linear constraints (which are to be removed from IA)
   int numNonLinear = 0;
   const double * rowLower = nlp.getRowLower();
   const double * rowUpper = nlp.getRowUpper();
   const double * colLower = nlp.getColLower();
   const double * colUpper = nlp.getColUpper();
   assert(m == nlp.getNumRows());
   double infty = si.getInfinity();
   double nlp_infty = nlp.getInfinity();
   Bonmin::vector<Ipopt::TNLP::LinearityType>  constTypes(m);
   Bonmin::vector<Ipopt::TNLP::LinearityType>  varTypes(n);
   problem->get_constraints_linearity(m, constTypes());
   problem->get_variables_linearity(n, varTypes());
   for (int i = 0; i < m; i++) {
     if (constTypes[i] == Ipopt::TNLP::NON_LINEAR) {
       nonLinear[numNonLinear++] = i;
     }
   }
   Bonmin::vector<double> rowLow(m - numNonLinear);
   Bonmin::vector<double> rowUp(m - numNonLinear);
   int ind = 0;
   for (int i = 0; i < m; i++) {
     if (constTypes[i] != Ipopt::TNLP::NON_LINEAR) {
       if (rowLower[i] > -nlp_infty) {
         //   printf("Lower %g ", rowLower[i]);
         rowLow[ind] = (rowLower[i]);
       } else
         rowLow[ind] = -infty;
       if (rowUpper[i] < nlp_infty) {
         //   printf("Upper %g ", rowUpper[i]);
         rowUp[ind] = (rowUpper[i]);
       } else
         rowUp[ind] = infty;
       ind++;
     }
   
   }
   
   CoinPackedMatrix mat(true, jRow(), jCol(), jValues(), nnz_jac_g);
   mat.setDimensions(m, n); // In case matrix was empty, this should be enough
   
   //remove non-linear constraints
   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 {
//.........这里部分代码省略.........

//#############################################################################
mcSol 
MibSHeuristic::solveSubproblem(double beta)
{

  /* 
     optimize wrt to weighted upper-level objective 
     over current feasible lp feasible region 
  */

  MibSModel * model = MibSModel_;
  OsiSolverInterface * oSolver = model->getSolver();
  //OsiSolverInterface * sSolver = new OsiCbcSolverInterface();  
  OsiSolverInterface* sSolver = new OsiSymSolverInterface();
  //sSolver = oSolver->clone();
  //OsiSolverInterface * sSolver = tmpSolver;
  //OsiSolverInterface * tmpSolver = new OsiSolverInterface(oSolver);
  
  double uObjSense(oSolver->getObjSense());
  double lObjSense(model->getLowerObjSense());  
  int lCols(model->getLowerDim());
  int uCols(model->getUpperDim());
  int * lColIndices = model->getLowerColInd();
  int * uColIndices = model->getUpperColInd();
  double * lObjCoeffs = model->getLowerObjCoeffs();
  const double * uObjCoeffs = oSolver->getObjCoefficients();

  double etol(etol_);
  int tCols(uCols + lCols); 

  assert(tCols == oSolver->getNumCols());


  sSolver->loadProblem(*oSolver->getMatrixByCol(),
		       oSolver->getColLower(), oSolver->getColUpper(),
		       oSolver->getObjCoefficients(),
		       oSolver->getRowLower(), oSolver->getRowUpper());

  int j(0);
  for(j = 0; j < tCols; j++){
    if(oSolver->isInteger(j))
      sSolver->setInteger(j);
  }


  double * nObjCoeffs = new double[tCols];
  int i(0), index(0);
  
  CoinZeroN(nObjCoeffs, tCols);
  
  /* Multiply the UL columns of the UL objective by beta */
  for(i = 0; i < uCols; i++){
    index = uColIndices[i];
    if(fabs(uObjCoeffs[index]) > etol)
      nObjCoeffs[index] = beta * uObjCoeffs[index] * uObjSense;
    else 
      nObjCoeffs[index] = 0.0;
  }
    
  /* Multiply the LL columns of the UL objective by beta */
  for(i = 0; i < lCols; i++){
    index = lColIndices[i];
    if(fabs(uObjCoeffs[index]) > etol)
      nObjCoeffs[index] = beta* uObjCoeffs[index] * uObjSense;
    else
      nObjCoeffs[index] = 0.0;
  }
  
  /* Add the LL columns of the LL objective multiplied by (1 - beta) */
  for(i = 0; i < lCols; i++){
    index = lColIndices[i];
    if(fabs(lObjCoeffs[i]) > etol)
      nObjCoeffs[index] += (1 - beta) * lObjCoeffs[i] * lObjSense;
  }
  
  sSolver->setObjective(nObjCoeffs);

  //int i(0);
  if(0){
    for(i = 0; i < sSolver->getNumCols(); i++){
      std::cout << "betaobj " << sSolver->getObjCoefficients()[i] << std::endl;
    }
  }

  if(0){
     sSolver->writeLp("afterbeta");
     //sSolver->writeMps("afterbeta");
  }
  
  if(0){  
    for(i = 0; i < sSolver->getNumCols(); i++){
      std::cout << "obj " << sSolver->getObjCoefficients()[i] << std::endl;
      std::cout << "upper " << sSolver->getColUpper()[i] << std::endl;
      std::cout << "lower " << sSolver->getColLower()[i] << std::endl;
    }
  }

  if(0){
    dynamic_cast<OsiCbcSolverInterface *> 
      (sSolver)->getModelPtr()->messageHandler()->setLogLevel(0);
//.........这里部分代码省略.........

void 
HeuristicInnerApproximation::extractInnerApproximation(OsiTMINLPInterface & nlp, OsiSolverInterface &si,
                                                       const double * x, bool getObj) {
   int n;
   int m;
   int nnz_jac_g;
   int nnz_h_lag;
   Ipopt::TNLP::IndexStyleEnum index_style;
   TMINLP2TNLP * problem = nlp.problem(); 
   //Get problem information
   problem->get_nlp_info(n, m, nnz_jac_g, nnz_h_lag, index_style);
   
   vector<int> jRow(nnz_jac_g);
   vector<int> jCol(nnz_jac_g);
   vector<double> jValues(nnz_jac_g);
   problem->eval_jac_g(n, NULL, 0, m, nnz_jac_g, jRow(), jCol(), NULL);
   if(index_style == Ipopt::TNLP::FORTRAN_STYLE)//put C-style
   {
     for(int i = 0 ; i < nnz_jac_g ; i++){
       jRow[i]--;
       jCol[i]--;
     }
   }
   
   //get Jacobian
   problem->eval_jac_g(n, x, 1, m, nnz_jac_g, NULL, NULL,
       jValues());
   
   vector<double> g(m);
   problem->eval_g(n, x, 1, m, g());
   
   vector<int> nonLinear(m);
   //store non linear constraints (which are to be removed from IA)
   int numNonLinear = 0;
   const double * rowLower = nlp.getRowLower();
   const double * rowUpper = nlp.getRowUpper();
   const double * colLower = nlp.getColLower();
   const double * colUpper = nlp.getColUpper();
   assert(m == nlp.getNumRows());
   double infty = si.getInfinity();
   double nlp_infty = nlp.getInfinity();
   vector<Ipopt::TNLP::LinearityType>  constTypes(m);
   problem->get_constraints_linearity(m, constTypes());
   for (int i = 0; i < m; i++) {
     if (constTypes[i] == Ipopt::TNLP::NON_LINEAR) {
       nonLinear[numNonLinear++] = i;
     }
   }
   vector<double> rowLow(m - numNonLinear);
   vector<double> rowUp(m - numNonLinear);
   int ind = 0;
   for (int i = 0; i < m; i++) {
     if (constTypes[i] != Ipopt::TNLP::NON_LINEAR) {
       if (rowLower[i] > -nlp_infty) {
         //   printf("Lower %g ", rowLower[i]);
         rowLow[ind] = (rowLower[i]);
       } else
         rowLow[ind] = -infty;
       if (rowUpper[i] < nlp_infty) {
         //   printf("Upper %g ", rowUpper[i]);
         rowUp[ind] = (rowUpper[i]);
       } else
         rowUp[ind] = infty;
       ind++;
     }
   
   }
   
   CoinPackedMatrix mat(true, jRow(), jCol(), jValues(), nnz_jac_g);
   mat.setDimensions(m, n); // In case matrix was empty, this should be enough
   
   //remove non-linear constraints
   mat.deleteRows(numNonLinear, nonLinear());
   
   int numcols = nlp.getNumCols();
   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;
     }
//.........这里部分代码省略.........

// mex function usage:
//  [x,y,status] = mexosi(n_vars,n_cons,A,x_lb,x_ub,c,Ax_lb,Ax_ub,isMIP,isQP,vartype,Q,options)
//                        0      1      2 3    4    5 6     7     8     9    10     11 12
void mexFunction(int nlhs, mxArray *plhs[],int nrhs, const mxArray *prhs[])
{
	// Enable printing in MATLAB
	int loglevel = 0;
	DerivedHandler *mexprinter = new DerivedHandler(); // assumed open	
	mexprinter->setLogLevel(loglevel);		 
	// check that we have the right number of inputs
	if(nrhs < 10) mexErrMsgTxt("At least 10 inputs required in call to mexosi. Bug in osi.m?...");
	// check that we have the right number of outputs
	if(nlhs < 3) mexErrMsgTxt("At least 3 ouptuts required in call to mexosi. Bug in osi.m?...");

    // Get pointers to input values
	const int    n_vars = (int)*mxGetPr(prhs[0]);
	const int    n_cons = (int)*mxGetPr(prhs[1]);
	const mxArray    *A =  prhs[2];
	const double  *x_lb =  mxGetPr(prhs[3]);
	const double  *x_ub =  mxGetPr(prhs[4]);
	const double     *c =  mxGetPr(prhs[5]);
	const double *Ax_lb =  mxGetPr(prhs[6]);
	const double *Ax_ub =  mxGetPr(prhs[7]);
	const bool    isMIP = (bool)*(mxLogical*)mxGetData(prhs[8]);
	const bool     isQP = (bool)*(mxLogical*)mxGetData(prhs[9]);
	mxLogical *isinteger = (mxLogical*)mxGetData(prhs[10]);
	const mxArray*    Q = prhs[11];
	// process the options
	int  returnStatus = 0;
	// extract row/col/value data from A
	const mwIndex * A_col_starts = mxGetJc(A);
	const mwIndex * A_row_index  = mxGetIr(A);
	const double  * A_data       = mxGetPr(A);
    // figure out the number of non-zeros in A
    int nnz = (int)(A_col_starts[n_vars] - A_col_starts[0]); // number of non-zeros
    //mexPrintf("nnz = %d, n_vars = %d, n_cons = %d\n",nnz,n_vars,n_cons);

    // we need to convert these into other types of indices for Coin to use them
    std::vector<CoinBigIndex> A_col_starts_coin(A_col_starts,A_col_starts+n_vars+1);
    std::vector<int>          A_row_index_coin(A_row_index,A_row_index+nnz);
    
    // declare the solver
    OsiSolverInterface* pSolver;
// 	initialize the solver
    if ( isMIP ) {
        pSolver = new OsiCbcSolverInterface;
    } else {
        pSolver = new OsiClpSolverInterface;
    }
   
//	OsiCbcSolverInterface is deprecated and CbcModel should be used instead but don't
//	know how to get that working with loadProblem. 
//     OsiCbcSolverInterface solver1;
//     CbcModel model(solver1);
//     CbcMain0(model);
//     OsiSolverInterface * pSolver = model.solver();
    
	if (nrhs>12) { // get stuff out of the options structure if provided
		// Finish me
	}
    
//     mexPrintf("Setting Log Level to 0.\n");
	// load the problem
	mexPrintf("Loading the problem.\n");
	pSolver->loadProblem( n_vars, n_cons, // problem size
						  &A_col_starts_coin[0], &A_row_index_coin[0], A_data, // the A matrix
						  x_lb,  x_ub, c, // the objective and bounds
						  Ax_lb, Ax_ub ); // the constraint bounds
    
//     pSolver->messageHandler()->setLogLevel(0); // This doesn't seem to work
    pSolver->setHintParam(OsiDoReducePrint,true,OsiHintTry);
    
	// deal with integer inputs
	if ( isMIP ) {
		for(int i=0;i<n_vars;i++) {
			if (isinteger[i]) pSolver->setInteger(i);
		}
	}
	if (isQP) {
		error("QP is not working yet");
		// need to call loadQuadraticObjective here ???
	}
    
//     CbcModel model(pSolver);
//     model.solver()->setHintParam(OsiDoReducePrint,true,OsiHintTry);
    
    
	// solve the problem
	//mexPrintf("Trying to solve the problem.\n");
    if (isMIP) {
        pSolver->branchAndBound();
//         model.branchAndBound();
    } else {
//         model.initialSolve();
        pSolver->initialSolve();
    }
	
	// Allocate memory for return data
    plhs[0] = mxCreateDoubleMatrix(n_vars,1, mxREAL); // for the solution
    plhs[1] = mxCreateDoubleMatrix(n_cons,1, mxREAL); // for the constraint prices
//.........这里部分代码省略.........

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

   // Build our own instance from scratch

   /*
    * This section adapted from Matt Galati's example 
    * on the COIN-OR Tutorial website.
    *
    * Problem from Bertsimas, Tsitsiklis page 21
    *  
    *  optimal solution: x* = (1,1)
    *  
    *  minimize -1 x0 - 1 x1
    *  s.t       1 x0 + 2 x1 <= 3
    *            2 x0 + 1 x1 <= 3
    *              x0        >= 0
    *              x1        >= 0
    */

   int n_cols = 2;
   double *objective    = new double[n_cols];//the objective coefficients
   double *col_lb       = new double[n_cols];//the column lower bounds
   double *col_ub       = new double[n_cols];//the column upper bounds

   //Define the objective coefficients.
   //minimize -1 x0 - 1 x1
   objective[0] = -1.0;
   objective[1] = -1.0;

   //Define the variable lower/upper bounds.
   // x0 >= 0   =>  0 <= x0 <= infinity
   // x1 >= 0   =>  0 <= x1 <= infinity
   col_lb[0] = 0.0;
   col_lb[1] = 0.0;
   col_ub[0] = si->getInfinity();
   col_ub[1] = si->getInfinity();
     
   int n_rows = 2;
   double *row_lb = new double[n_rows]; //the row lower bounds
   double *row_ub = new double[n_rows]; //the row upper bounds
     
   //Define the constraint matrix.
   CoinPackedMatrix *matrix =  new CoinPackedMatrix(false,0,0);
   matrix->setDimensions(0, n_cols);

   //1 x0 + 2 x1 <= 3  =>  -infinity <= 1 x0 + 2 x2 <= 3
   CoinPackedVector row1;
   row1.insert(0, 1.0);
   row1.insert(1, 2.0);
   row_lb[0] = -1.0 * si->getInfinity();
   row_ub[0] = 3.0;
   matrix->appendRow(row1);

   //2 x0 + 1 x1 <= 3  =>  -infinity <= 2 x0 + 1 x1 <= 3
   CoinPackedVector row2;
   row2.insert(0, 2.0);
   row2.insert(1, 1.0);
   row_lb[1] = -1.0 * si->getInfinity();
   row_ub[1] = 3.0;
   matrix->appendRow(row2);

   //load the problem to OSI
   si->loadProblem(*matrix, col_lb, col_ub, objective, row_lb, row_ub);

   //write the MPS file to a file called example.mps
   si->writeMps("example");

  

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

//.........这里部分代码省略.........
     
     /*
       for(i = 0; i < lRows; i++){
       CoinPackedVector row;
       index1 = lRowIndices[i];
       start = matStarts[index1];
       end = start + matrix->getVectorSize(index1);
       for(j = start; j < end; j++){
       index2 = matIndices[j];
       //tmp = findIndex(index, lCols, lColIndices);
       tmp = binarySearch(0, lCols - 1, index2, lColIndices);
       if(tmp > -1)
       row.insert(tmp, matElements[j]);
       }
       newMat->appendRow(row);
       }
     */
     for(i = 0; i < lRows; i++){
	CoinPackedVector row;
	index1 = lRowIndices[i];
	for(j = 0; j < lCols; j++){
	   index2 = lColIndices[j];
	   tmp = matrix->getCoefficient(index1, index2);
	   row.insert(j, tmp);
	}
	newMat->appendRow(row);
     }

     /*
       nSolver->assignProblem(newMat, colLb, colUb,
       objCoeffs, rowLb, rowUb);
     */
     
     nSolver->loadProblem(*newMat, colLb, colUb,
			  objCoeffs, rowLb, rowUb);
     
     for(i = 0; i < intCnt; i++){
	nSolver->setInteger(integerVars[i]);
     }
     //nSolver->setInteger(integerVars, intCnt);
     
     nSolver->setObjSense(objSense); //1 min; -1 max
     
     nSolver->setHintParam(OsiDoReducePrint, true, OsiHintDo);

#if 0
     if(0){
	dynamic_cast<OsiCbcSolverInterface *> 
	   (nSolver)->getModelPtr()->messageHandler()->setLogLevel(0);
     }
     else{
	dynamic_cast<OsiSymSolverInterface *> 
	   (nSolver)->setSymParam("prep_level", -1);
	
	dynamic_cast<OsiSymSolverInterface *> 
	   (nSolver)->setSymParam("verbosity", -2);
	
	dynamic_cast<OsiSymSolverInterface *> 
	   (nSolver)->setSymParam("max_active_nodes", 1);
     }
#endif
     delete [] integerVars;

  }else{
     nSolver = solver_;
  }

//.........这里部分代码省略.........
        //    7x1 +  4x2 <= 8
        //   -7x1 +  4x2 <= 1
        //     x1, x2 >= 0 and x1, x2 integer
        // Slacks are s1, s2, s3, s4



        //Test that problem is correct
        // Optimal Basis is x1, x2, s3, s4 with tableau
        //    x1            0.25 s1  -0.25 s2             =  0.5
        //           x2     0.25 s1   0.25 s2             =  1
        //                 -2.75 s1   0.75 s2    s3       =  0.5
        //                  0.75 s1  -2.75 s2        s4   =  0.5
        // z=              -0.25 s1  -0.25 s2             =  -1
        // Gomory cut from variable x1 is x2 <= 0.5
        // Can be improved by first pivoting s2 in and s4 out, then s1 in and s3 out
        // to x2 <= 0.25
        {
            int start[2] = {0,4};
            int length[2] = {4,4};
            int rows[8] = {0,1,2,3,0,1,2,3};
            double elements[8] = {2.0,-2.0,7.0,-7.0,2.0,2.0,4.0,4.0};
            CoinPackedMatrix  columnCopy(true,4,2,8,elements,rows,start,length);

            double rowLower[4]={-COIN_DBL_MAX,-COIN_DBL_MAX,
                                -COIN_DBL_MAX,-COIN_DBL_MAX};
            double rowUpper[4]={3.,1.,8.,1.};
            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)
            {

//.........这里部分代码省略.........
      if (x<=0)
	throw("bad fscanf");
      printf("%d columns, %d elements, upper %g\n",ncol,nel,up);
      double * sol1 = new double[nel];
      double * el1 = new double[nel];
      int * col1 = new int[nel];
      CoinBigIndex * start = new CoinBigIndex [ncol+1];
      memset(start,0,ncol*sizeof(CoinBigIndex ));
      int * row = new int[nel];
      int i;
      for (i=0;i<nel;i++) {
	x=fscanf(fp,"%d %lg %lg",col1+i,el1+i,sol1+i);
	if (x<=0)
	  throw("bad fscanf");
	printf("[%d, e=%g, v=%g] ",col1[i],el1[i],sol1[i]);
	start[col1[i]]=1;
	row[i]=0;
      }
      printf("\n");
      // Setup
      OsiSolverInterface  * siP = baseSiP->clone();
      
      double lo=-1.0e30;
      double * upper = new double[ncol];
      start[ncol]=nel;
      int last=0;
      for (i=0;i<ncol;i++) {
	upper[i]=1.0;
	int marked=start[i];
	start[i]=last;
	if (marked)
	  last++;
      }
      siP->loadProblem(ncol,1,start,row,el1,NULL,upper,NULL,&lo,&up);
      // use upper for solution
      memset(upper,0,ncol*sizeof(double));
      for (i=0;i<nel;i++) {
	int icol=col1[i];
	upper[icol]=sol1[i];
	siP->setInteger(icol);
      }
      siP->setColSolution(upper);
      delete [] sol1;
      delete [] el1;
      delete [] col1;
      delete [] start;
      delete [] row;
      delete [] upper;
      CglKnapsackCover kccg;
      
      OsiCuts cuts;    
      
      // Test generateCuts method
      kccg.generateCuts(*siP,cuts);
      // print out and compare to known cuts
      int numberCuts = cuts.sizeRowCuts();
      if (numberCuts) {
	for (i=0;i<numberCuts;i++) {
	  OsiRowCut * thisCut = cuts.rowCutPtr(i);
	  int n=thisCut->row().getNumElements();
	  printf("Cut %d has %d entries, rhs %g %g =>",i,n,thisCut->lb(),
		 thisCut->ub());
	  int j;
	  const int * index = thisCut->row().getIndices();
	  const double * element = thisCut->row().getElements();
	  for (j=0;j<n;j++) {

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

#if 0
      // Get dual multipliers and build gradient of the lagrangean
      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;