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

//#############################################################################
bool
MibSHeuristic::checkLowerFeasibility(OsiSolverInterface * si, 
				     double * solution)
{

  MibSModel * model = MibSModel_;
  OsiSolverInterface * lSolver = model->bS_->setUpModel(si, true, solution);

  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();

  if(lSolver->isProvenOptimal())
    return true;
  else
    return false;
}

//#############################################################################
bfSol*
MibSHeuristic::getBilevelSolution(const double * sol, double origLower)
{

  /* 
     Find a bilevel feasible solution by solving the LL problem
     for a fixed UL solution, given by the UL portion of sol
  */

  MibSModel * model = MibSModel_;
  OsiSolverInterface * oSolver = model->getSolver();
  OsiSolverInterface * lSolver = model->bS_->setUpModel(oSolver, true, sol);
  
  //double uObjSense(model->getSolver()->getObjSense());
  int lCols(model->getLowerDim());
  int uCols(model->getUpperDim());
  int * lColIndices = model->getLowerColInd();
  int * uColIndices = model->getUpperColInd();
  double etol(etol_);

  int tCols(uCols + lCols); 
  int i(0);

  if(0){
     lSolver->writeLp("bilevelsolver");
     std::cout 
       << "Original Lower-level Solution Value: " << origLower << std::endl;
     for(i = 0; i < lCols; i++){
       std::cout << "lowsol[" << i << "]: " << sol[lColIndices[i]] << std::endl;
     }
  }
  
  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();

  if(lSolver->isProvenOptimal()){

    double objVal(0.0);
    double lowerObj(lSolver->getObjValue());
    double * colsol = new double[tCols];

    for(i = 0; i < uCols; i++){
      colsol[uColIndices[i]] = sol[uColIndices[i]];
    }

    if(0){
      std::cout << "candidate lower solution value: " << origLower << std::endl;
      std::cout << "actual lower solution value: " << lowerObj << std::endl;
    }

    if(fabs(origLower - lowerObj) < etol){
      //original solution was bilevel feasible
      if(0)
	std::cout << "Original solution was bilevel feasible:" << std::endl;
      for(i = 0; i < lCols; i++){
	if(0){
	  std::cout << "lowerportion[" 
		    << i << "]: " << sol[lColIndices[i]] << std::endl; 
	}	
	colsol[lColIndices[i]] = sol[lColIndices[i]];
      }
    }
    else{
      if(0){
	std::cout << "Not bilevel feasible." << std::endl;
      }
      for(i = 0; i < lCols; i++){
	if(0){
	  std::cout << "newportion[" 
		    << i << "]: " << lSolver->getColSolution()[i] << std::endl; 
	}
	colsol[lColIndices[i]] = lSolver->getColSolution()[i];
      }
    }


    for(i = 0; i < tCols; i++)
      objVal += colsol[i] * oSolver->getObjCoefficients()[i];

    bfSol * bfsol = 
      new bfSol(objVal, colsol);
    delete lSolver;
    return bfsol;
  }
  else{
//.........这里部分代码省略.........

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

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

  MibSModel * model = MibSModel_;

  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 * lColIndices = model->getLowerColInd();
  int * uColIndices = model->getUpperColInd();
  double * lObjCoeffs = model->getLowerObjCoeffs();
  
  //int tCols(lCols + uCols);
  int tCols(oSolver->getNumCols());
  //assert(tCols == oSolver->getNumCols());

  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 * nObjCoeffs = new double[tCols];
  int i(0), index(0);
  
  CoinZeroN(nObjCoeffs, tCols);

  for(i = 0; i < lCols; i++){
    index = lColIndices[i];
    nObjCoeffs[index] = lObjCoeffs[i];
  }

  //MibS objective sense is the opposite of OSI's!
  hSolver->setObjSense(objSense);

  hSolver->setObjective(nObjCoeffs);
 
  //double cutoff(model->getCutoff());
  double cutoff(model->getKnowledgeBroker()->getIncumbentValue());

  if(model->getNumSolutions()){
  
    CoinPackedVector objCon;
    //double rhs(cutoff * objSense);
    //double smlTol(1.0);
    double rhs(cutoff);
       
    for(i = 0; i < tCols; i++){
      objCon.insert(i, oSolver->getObjCoefficients()[i] 
		    * oSolver->getObjSense());
    }
    
    hSolver->addRow(objCon, - hSolver->getInfinity(), rhs);
  }
  
  if(0)
     hSolver->writeLp("lobjheurstic");

  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(hSolver->isProvenOptimal()){

    double upperObjVal(0.0);

    /*****************NEW ******************/

    MibSSolution *mibSol = NULL;

    OsiSolverInterface * lSolver = model->bS_->setUpModel(hSolver, true);
//.........这里部分代码省略.........

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

  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);
  }
  else{
    dynamic_cast<OsiSymSolverInterface *> 
      (sSolver)->setSymParam("prep_level", -1);
    
    dynamic_cast<OsiSymSolverInterface *> 
      (sSolver)->setSymParam("verbosity", -2);

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

  //dynamic_cast<OsiSymSolverInterface *> (sSolver)->branchAndBound();

  sSolver->branchAndBound();

  if(sSolver->isProvenOptimal()){

    if(0){
      std::cout << "writing lp file." << std::endl;
      sSolver->writeLp("afterbeta");
      //sSolver->writeMps("afterbeta");
    }
    
    double upperObjVal(0.0);
    double lowerObjVal(0.0);
    

    for(i = 0; i < tCols; i++){
      upperObjVal += 
	sSolver->getColSolution()[i] * oSolver->getObjCoefficients()[i];
      if(0){
	std::cout << "sSolver->getColSolution()[" << i << "] :"
		  << sSolver->getColSolution()[i] << std::endl;
      }
    }
    lowerObjVal = getLowerObj(sSolver->getColSolution(), lObjSense);
    
    if(beta == 1.0){
      
      /*
	fix upper-level objective to current value and 
	reoptimize wrt to lower-level objective
      */
      
      //OsiSolverInterface * nSolver = new OsiCbcSolverInterface();
      OsiSolverInterface * nSolver = new OsiSymSolverInterface();

//#############################################################################
bfSol*
MibSHeuristic::getBilevelSolution1(const double * sol)
{

  /* 
     Find a bilevel feasible solution by solving the LL problem
     for a fixed UL solution, given by the UL portion of sol
  */


  MibSModel * model = MibSModel_;
  OsiSolverInterface * oSolver = model->getSolver();
  OsiSolverInterface * lSolver = new OsiCbcSolverInterface(oSolver);
  
  //double uObjSense(model->getSolver()->getObjSense());
  double lObjSense(model->getLowerObjSense());  
  int lCols(model->getLowerDim());
  int uCols(model->getUpperDim());
  int * lColIndices = model->getLowerColInd();
  int uRowNum = model->getUpperRowNum();
  int lRowNum = model->getLowerRowNum();
  int * uRowIndices = model->getUpperRowInd();
  int * lRowIndices = model->getLowerRowInd();
  int * uColIndices = model->getUpperColInd();
  double * lObjCoeffs = model->getLowerObjCoeffs();

  int tCols(uCols + lCols); 
  int i(0), index(0);

  /* delete the UL rows */
  lSolver->deleteRows(uRowNum, uRowIndices);

  /* Fix the UL variables to their current value in sol */

  for(i = 0; i < uCols; i++){
    index = uColIndices[i];
    lSolver->setColLower(index, sol[index]);
    lSolver->setColUpper(index, sol[index]);
  }

  /* Set the objective to the LL objective coefficients */

  double * nObjCoeffs = new double[tCols];
  CoinZeroN(nObjCoeffs, tCols);
      
  for(i = 0; i < lCols; i++){
    index = lColIndices[i];
    nObjCoeffs[index] = lObjCoeffs[i] * lObjSense;
  }
  
  lSolver->setObjective(nObjCoeffs);

  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();

  if(lSolver->isProvenOptimal()){

    double objVal(0.0);

    for(i = 0; i < tCols; i++)
      objVal += lSolver->getColSolution()[i] * oSolver->getObjCoefficients()[i];
    
    double * colsol = new double[tCols];
    CoinCopyN(lSolver->getColSolution(), tCols, colsol);
 
    bfSol * bfsol = 
      new bfSol(objVal, colsol);
    delete lSolver;
    return bfsol;
  }
  else{
    delete lSolver;
    return NULL;
  }

}

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

//.........这里部分代码省略.........
     if (whichCutsLL == 0){
	sym_set_int_param(env, "generate_cgl_cuts", FALSE);
     }else{
	sym_set_int_param(env, "generate_cgl_gomory_cuts", GENERATE_DEFAULT);
     }
     if (whichCutsLL == 1){
	sym_set_int_param(env, "generate_cgl_knapsack_cuts", 
			  DO_NOT_GENERATE);
	sym_set_int_param(env, "generate_cgl_probing_cuts", 
			  DO_NOT_GENERATE);
	sym_set_int_param(env, "generate_cgl_clique_cuts", 
			  DO_NOT_GENERATE);
	sym_set_int_param(env, "generate_cgl_twomir_cuts", 
			  DO_NOT_GENERATE);
	sym_set_int_param(env, "generate_cgl_flowcover_cuts", 
			  DO_NOT_GENERATE);
     }
  }else if (feasCheckSolver == "CPLEX"){
#ifdef USE_CPLEX
     lSolver->setHintParam(OsiDoReducePrint);
     lSolver->messageHandler()->setLogLevel(0);
     CPXENVptr cpxEnv = 
	dynamic_cast<OsiCpxSolverInterface*>(lSolver)->getEnvironmentPtr();
     assert(cpxEnv);
     CPXsetintparam(cpxEnv, CPX_PARAM_SCRIND, CPX_OFF);
     CPXsetintparam(cpxEnv, CPX_PARAM_THREADS, maxThreadsLL);
#endif
  }
  
  if (warmStartLL && feasCheckSolver == "SYMPHONY"){
     lSolver->resolve();
     setWarmStart(lSolver->getWarmStart());
  }else{
     lSolver->branchAndBound();
  }

  const double * sol = model_->solver()->getColSolution();
  double objVal(lSolver->getObjValue() * model_->getLowerObjSense());
  
  MibSTreeNode * node = static_cast<MibSTreeNode *>(model_->activeNode_);
  MibSTreeNode * parent = 
    static_cast<MibSTreeNode *>(model_->activeNode_->getParent());

  if((!node->isBoundSet()) 
     && (node->getIndex() != 0)){
    double parentBound = parent->getLowerUB();
    node->setLowerUB(parentBound);
    node->setIsBoundSet(true);
  }
  
  if(objVal > node->getLowerUB()){
    
    node->setLowerUB(objVal);
    node->setIsBoundSet(true);
    
  }

  double etol(model_->etol_);
  double lowerObj = getLowerObj(sol, model_->getLowerObjSense());  

  int lN(model_->lowerDim_); // lower-level dimension
  int uN(model_->upperDim_); // lower-level dimension
  if(!optLowerSolution_)
    optLowerSolution_ = new double[lN];

  if(!optLowerSolutionOrd_)

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

            // Now add the routes from this solution to the IP
            for(k=1;k<=V.get_total_number_of_routes();k++)
            {
                // Clean up the route by running INTRA_ROUTE optimizations only
                // using the route_search method of the different local search
                // heuristics, accepting improving moves only (VRPH_DOWNHILL)
                OnePointMove OPM;
                TwoOpt TO;
                ThreeOpt ThO;
                while(OPM.route_search(&V,k,k,VRPH_DOWNHILL|VRPH_INTRA_ROUTE_ONLY )){}
                while(TO.route_search(&V,k,k,VRPH_DOWNHILL|VRPH_INTRA_ROUTE_ONLY  )){};
                while(ThO.route_search(&V,k,VRPH_DOWNHILL|VRPH_INTRA_ROUTE_ONLY )){};

                // Copy route k from the solution to the VRPRoute R
                V.update_route(k,&route);
                route.create_name();
                // Add it to the "route warehouse" - this uses a hash table to keep track
                // of duplicate columns
                status=V.route_wh->add_route(&route);

                if(status!=DUPLICATE_ROUTE)
                {
                    // This route is not currently in the WH and so it cannot be in the
                    // set partitioning problem
                    //OSI_add_route(si,&V,&route);
                    OSI_add_route(si,&V,&route,route_id,orderings);
                    route_id++;
                }
            }
            
            // Set the row RHS's if we need to
            if(first_sol)
            {
                first_sol=false;
                for(int rownum=0;rownum<n;rownum++)
                    si->setRowBounds(rownum,1,1);
                // Note that changing this to >= would be a set covering problem
                // where each customer can be visited by more than one route
            }
        }

        // Now erase all the solutions from the WH
        V.solution_wh->liquidate();

        if(verbose)
        {
            printf("Attempt %02d:  Solving IP with %d columns\n",i,si->getNumCols());
            printf("%d routes in the WH\n",V.route_wh->num_unique_routes);
        }

        // Solve the current set partitioning problem using the MIP solver
        start=clock();
        si->branchAndBound();
        stop=clock();
        mip_time += (stop-start);

        double opt=si->getObjValue();
        x=si->getColSolution();
        last_num_cols=si->getNumCols();
        if(verbose)
            printf("Optimal solution (%d columns) is %f\n",last_num_cols,opt);

        
        // Now recover the solution from the IP solution
        OSI_recover_solution(si, orderings, IP_sol_buff);
       
        if(verbose)
            printf("IP solution has obj. function value: %5.2f\n"
            "Best heuristic  obj. function value: %5.2f\n",
            si->getObjValue(),best_heur_sol);     
    }

    if(verbose)
        printf(
        "\nResults\n"
        "--------\n"
        "After %d runs\n"
        "IP solution has obj. function value: %5.2f\n"
        "Best heuristic  obj. function value: %5.2f\n",
        num_attempts,si->getObjValue(),best_heur_sol);

    // print to stderr since GLPK prints "conflict graph" to stdout (a known bug...)
    // best_heur_sol best_mip_sol heur_time mip_time
    fprintf(stderr,"%5.3f %5.3f %5.3f %5.3f\n", best_heur_sol, si->getObjValue(),
        (double)(heur_time)/CLOCKS_PER_SEC, (double)(mip_time)/CLOCKS_PER_SEC);

    delete V.route_wh;
    delete fresh_solution;
    delete [] sol_buff;
    delete [] IP_sol_buff;
    delete si;

    for(i=0;i<MAX_ROUTES;i++)
        if(orderings[i])
            delete [] orderings[i];

     
    return 0;
}

/*
 * Class:     thebeast_osi_OsiSolverJNI
 * Method:    branchAndBound
 * Signature: (I)V
 */
JNIEXPORT void JNICALL Java_thebeast_osi_OsiSolverJNI_branchAndBound
  (JNIEnv *, jobject, jint ptr){
  OsiSolverInterface* solver = (OsiSolverInterface*) ptr;
  solver->branchAndBound();
}