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

/*===========================================================================*
  Scan through the variables and select those that are binary and are at a
  fractional level.
 *===========================================================================*/
void
CglClique::selectFractionalBinaries(const OsiSolverInterface& si)
{
   // extract the primal tolerance from the solver
   double lclPetol = 0.0;
   si.getDblParam(OsiPrimalTolerance, lclPetol);

   const int numcols = si.getNumCols();
   if (petol<0.0) {
     // do all if not too many
     int n=0;
     for (int i = 0; i < numcols; ++i) {
       if (si.isBinary(i))
	 n++;
     }
     if (n<5000)
       lclPetol=-1.0e-5;
   }
   const double* x = si.getColSolution();
   std::vector<int> fracind;
   int i;
   for (i = 0; i < numcols; ++i) {
      if (si.isBinary(i) && x[i] > lclPetol && x[i] < 1-petol)
	 fracind.push_back(i);
   }
   sp_numcols = static_cast<int>(fracind.size());
   sp_orig_col_ind = new int[sp_numcols];
   sp_colsol = new double[sp_numcols];
   for (i = 0; i < sp_numcols; ++i) {
      sp_orig_col_ind[i] = fracind[i];
      sp_colsol[i] = x[fracind[i]];
   }
}

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

//-------------------------------------------------------------------
// Generate Stored cuts
//------------------------------------------------------------------- 
void 
CglStoredUser::generateCuts(const OsiSolverInterface & si, OsiCuts & cs,
			     const CglTreeInfo info) const
{
  // Get basic problem information
  const double * solution = si.getColSolution();
  if (info.inTree&&info.pass>numberPasses_) {
    // only continue if integer feasible
    int numberColumns=si.getNumCols(); 
    int i;
    const double * colUpper = si.getColUpper();
    const double * colLower = si.getColLower();
    int numberAway=0;
    for (i=0;i<numberColumns;i++) {
      double value = solution[i];
      // In case slightly away from bounds
      value = CoinMax(colLower[i],value);
      value = CoinMin(colUpper[i],value);
      if (si.isInteger(i)&&fabs(value-fabs(value+0.5))>1.0e-5) 
	numberAway++;
    }
    if (numberAway)
      return; // let code branch
  }
  int numberRowCuts = cuts_.sizeRowCuts();
  for (int i=0;i<numberRowCuts;i++) {
    const OsiRowCut * rowCutPointer = cuts_.rowCutPtr(i);
    double violation = rowCutPointer->violated(solution);
    if (violation>=requiredViolation_)
      cs.insert(*rowCutPointer);
  }
}

// Create result
void OsiSolverResult::createResult(const OsiSolverInterface &solver, const double *lowerBefore,
  const double *upperBefore)
{
  delete[] primalSolution_;
  delete[] dualSolution_;
  if (solver.isProvenOptimal() && !solver.isDualObjectiveLimitReached()) {
    objectiveValue_ = solver.getObjValue() * solver.getObjSense();
    CoinWarmStartBasis *basis = dynamic_cast< CoinWarmStartBasis * >(solver.getWarmStart());
    assert(basis);
    basis_ = *basis;
    int numberRows = basis_.getNumArtificial();
    int numberColumns = basis_.getNumStructural();
    assert(numberColumns == solver.getNumCols());
    assert(numberRows == solver.getNumRows());
    primalSolution_ = CoinCopyOfArray(solver.getColSolution(), numberColumns);
    dualSolution_ = CoinCopyOfArray(solver.getRowPrice(), numberRows);
    fixed_.addBranch(-1, numberColumns, lowerBefore, solver.getColLower(),
      upperBefore, solver.getColUpper());
  } else {
    // infeasible
    objectiveValue_ = COIN_DBL_MAX;
    basis_ = CoinWarmStartBasis();
    ;
    primalSolution_ = NULL;
    dualSolution_ = NULL;
  }
}

/* Generate cuts for the model data contained in si.
   The generated cuts are inserted into and returned in the
   collection of cuts cs.
*/
bool
AbcCutGenerator::generateCuts( OsiCuts & cs , bool fullScan)
{
    int howOften = whenCutGenerator_;
    if (howOften == -100)
	return false;
    if (howOften > 0)
	howOften = howOften % 1000000;
    else 
	howOften = 1;
    if (!howOften)
	howOften = 1;
    bool returnCode = false;
    OsiSolverInterface * solver = model_->solver();

#if defined(ABC_DEBUG_MORE)
    std::cout << "model_->getNodeCount() = " << model_->getNodeCount()
	      << std::endl;
#endif


    if (fullScan || (model_->getNodeCount() % howOften) == 0 ) {
	CglProbing* generator =
	    dynamic_cast<CglProbing*>(generator_);
	if (!generator) {
	    generator_->generateCuts(*solver,cs);
	} else {
	    // Probing - return tight column bounds
	   CglTreeInfo info;
	    generator->generateCutsAndModify(*solver,cs,&info);
	    const double * tightLower = generator->tightLower();
	    const double * lower = solver->getColLower();
	    const double * tightUpper = generator->tightUpper();
	    const double * upper = solver->getColUpper();
	    const double * solution = solver->getColSolution();
	    int j;
	    int numberColumns = solver->getNumCols();
	    double primalTolerance = 1.0e-8;
	    for (j=0; j<numberColumns; j++) {
		if (tightUpper[j] == tightLower[j] &&
		    upper[j] > lower[j]) {
		    // fix
		    solver->setColLower(j, tightLower[j]);
		    solver->setColUpper(j, tightUpper[j]);
		    if (tightLower[j] > solution[j] + primalTolerance ||
			tightUpper[j] < solution[j] - primalTolerance)
			returnCode = true;
		}
	    }
	}
    }
    return returnCode;
}

// Creates a branching object from this infeasible object.
BcpsBranchObject * 
BlisObjectInt::createBranchObject(BcpsModel *m, int direction) const
{
    BlisModel *model = dynamic_cast<BlisModel* >(m);
    OsiSolverInterface * solver = model->solver();
    
    double integerTolerance = model->BlisPar()->entry(BlisParams::integerTol);
    
    const double * solution = solver->getColSolution();
    const double * lower = solver->getColLower();
    const double * upper = solver->getColUpper();
    
    double value = solution[columnIndex_];
    //std::cout << "COL"<< columnIndex_ << ": x = " << value << std::endl;
    
    // Force value in bounds.
    value = CoinMax(value, lower[columnIndex_]);
    value = CoinMin(value, upper[columnIndex_]);
    
    double nearest = floor(value + 0.5);
    
    assert (upper[columnIndex_] > lower[columnIndex_]);
    
    int hotstartStrategy = model->getHotstartStrategy();
    
    if (hotstartStrategy <= 0) {
        if (fabs(value - nearest) < integerTolerance) {
            // Already integeral.
            std::cout << "ERROR: COL" << columnIndex_ << ": x=" << value 
                      << ", nearest=" << nearest 
                      << ", intTol=" << integerTolerance << std::endl;
            assert(0);
        }
    } 
    else {
	const double * incumbent = model->incumbent();
	double targetValue = incumbent[columnIndex_];
	if (direction > 0) {
	    value = targetValue - 0.1;
	}
	else {
	    value = targetValue + 0.1;
	}
    }
    
    return new BlisBranchObjectInt(model, objectIndex_, direction, value);
}

  void 
  LinearCutsGenerator::generateCuts(const OsiSolverInterface &solver, OsiCuts &cs,
                     const CglTreeInfo info) const {

    //const OsiTMINLPInterface * tmp = dynamic_cast<const OsiTMINLPInterface *>(&solver);
    OsiTMINLPInterface * nlp = dynamic_cast<OsiTMINLPInterface *>(solver.clone());//const_cast<OsiTMINLPInterface *>(tmp);
    assert(nlp);
    OuterApprox oa;
    //si.writeMps("toto");
    int numberRows = nlp->getNumRows();
    for(int i = 0 ; i < 5 ; i++){
      nlp->resolve();
      OsiClpSolverInterface si;
      oa(*nlp, &si, solver.getColSolution(), true); 
      si.resolve();
      OsiCuts cuts;
      for(std::list<Coin::SmartPtr<CuttingMethod> >::const_iterator i = methods_.begin() ;
          i != methods_.end() ; i++){
         (*i)->cgl->generateCuts(si, cuts, info);
      }
      std::vector<OsiRowCut *> mycuts(cuts.sizeRowCuts());
      for(int i = 0 ; i < cuts.sizeRowCuts() ; i++){
        mycuts[i] = cuts.rowCutPtr(i);
        cs.insert(*mycuts[i]);
      }
      nlp->applyRowCuts(mycuts.size(), const_cast<const OsiRowCut **> (&mycuts[0]));
    }

    // Take off slack cuts
    std::vector<int> kept;
    int numberRowsNow = nlp->getNumRows();
    int * del = new int [numberRowsNow-numberRows];
    nlp->resolve();
    
    const double * activity = nlp->getRowActivity();
    const double * lb = nlp->getRowLower();
    const double * ub = nlp->getRowUpper();
    CoinRelFltEq eq(1e-06);
    //int nDelete=0;
    for (int i=numberRowsNow -1;i>=numberRows;i--) {
      if ( !(eq(activity[i], lb[i]) || eq(activity[i], ub[i])) )
        cs.eraseRowCut(i - numberRows);
    }
    delete [] del;
    delete nlp;
  }

bool
OaDecompositionBase::OaDebug::checkInteger(const OsiSolverInterface &nlp, 
                                           std::ostream & os) const {
   const double * colsol = nlp.getColSolution();
   int numcols = nlp.getNumCols();
  for (int i = 0 ; i < numcols ; i++) {
    if (nlp.isInteger(i)) {
      if (fabs(colsol[i]) - floor(colsol[i] + 0.5) >
          1e-07) {
        std::cerr<<"Integer infeasible point (should not be), integer infeasibility for variable "<<i
        <<" is, "<<fabs(colsol[i] - floor(colsol[i] + 0.5))<<std::endl;
      }
    }
    return true;
  }

}

// Compute the infeasibility based on currently relax solution.
double 
BlisObjectInt::infeasibility(BcpsModel *m, int & preferredWay) const
{
    BlisModel * model = dynamic_cast<BlisModel *>(m);
    OsiSolverInterface * solver = model->solver();
    
    const double * solution =  solver->getColSolution();

    const double * lower = solver->getColLower();
    const double * upper = solver->getColUpper();

    double value = solution[columnIndex_];
    
    value = std::max(value, lower[columnIndex_]);
    value = std::min(value, upper[columnIndex_]);
    
    //printf("%d %g %g %g %g\n",columnIndex_,value,lower[columnIndex_],
    //   solution[columnIndex_],upper[columnIndex_]);

    double nearest = floor(value + (1.0 - breakEven_));
    double integerTolerance = model->BlisPar()->entry(BlisParams::integerTol);
    if (nearest > value) {
	preferredWay = 1;
    }
    else {
	preferredWay = -1;
    }
    
    double weight = fabs(value - nearest);

    // normalize so weight is 0.5 at break even
    if (nearest < value) {
	weight = (0.5/breakEven_) * weight;
    }
    else {
	weight = (0.5/(1.0 - breakEven_)) * weight;
    }
    
    if (fabs(value - nearest) <= integerTolerance) {
	return 0.0;
    }
    else {
	return weight;
    }
}

// Returns true if current solution satsifies one side of branch
bool 
OsiSolverBranch::feasibleOneWay(const OsiSolverInterface & solver) const
{
  bool feasible = false;
  int numberColumns = solver.getNumCols();
  const double * columnLower = solver.getColLower();
  const double * columnUpper = solver.getColUpper();
  const double * columnSolution = solver.getColSolution();
  double primalTolerance;
  solver.getDblParam(OsiPrimalTolerance,primalTolerance);
  for (int base = 0; base<4; base +=2) {
    feasible=true;
    int i;
    for (i=start_[base];i<start_[base+1];i++) {
      int iColumn = indices_[i];
      if (iColumn<numberColumns) {
        double value = CoinMax(bound_[i],columnLower[iColumn]);
        if (columnSolution[iColumn]<value-primalTolerance) {
          feasible=false;
          break;
        }
      } else {
        abort(); // do later (other stuff messed up anyway - e.g. CBC)
      }
    }
    if (!feasible)
      break;
    for (i=start_[base+1];i<start_[base+2];i++) {
      int iColumn = indices_[i];
      if (iColumn<numberColumns) {
        double value = CoinMin(bound_[i],columnUpper[iColumn]);
        if (columnSolution[iColumn]>value+primalTolerance) {
          feasible=false;
          break;
        }
      } else {
        abort(); // do later (other stuff messed up anyway - e.g. CBC)
      }
    }
    if (feasible)
      break; // OK this way
  }
  return feasible;
}

// Generate cuts
void
CglFakeClique::generateCuts(const OsiSolverInterface& si, OsiCuts & cs,
			const CglTreeInfo info) const
{
  if (fakeSolver_) {
    assert (si.getNumCols()==fakeSolver_->getNumCols());
    fakeSolver_->setColLower(si.getColLower());
    fakeSolver_->setColSolution(si.getColSolution());
    fakeSolver_->setColUpper(si.getColUpper());
    CglClique::generateCuts(*fakeSolver_,cs,info);
    if (probing_) {
      // get and set branch and bound cutoff
      double cutoff;
      si.getDblParam(OsiDualObjectiveLimit,cutoff);
      fakeSolver_->setDblParam(OsiDualObjectiveLimit,cutoff);
      probing_->generateCuts(*fakeSolver_,cs,info);
    }
  } else {
    // just use real solver
    CglClique::generateCuts(si,cs,info);
  }
}

// Force this object within exiting bounds, then fix the bounds at the
// the nearest integer value. Assume solution value is within tolerance of
// the nearest integer.
void 
BlisObjectInt::feasibleRegion(BcpsModel *m)
{

    BlisModel *model = dynamic_cast<BlisModel* >(m);
    OsiSolverInterface * solver = model->solver();
  
    const double * solution =  solver->getColSolution();
    const double * lower = solver->getColLower();
    const double * upper = solver->getColUpper();

    double value = solution[columnIndex_];

    // 1) Force value to be in bounds.
    value = CoinMax(value, lower[columnIndex_]);
    value = CoinMin(value, upper[columnIndex_]);
    
    double nearest = floor(value + 0.5);

    // 2) Fix variable at the nearest integer
    assert (fabs(value - nearest) <= 0.01);
    solver->setColLower(columnIndex_, nearest);
    solver->setColUpper(columnIndex_, nearest);
}

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

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

  if(hSolver->isProvenOptimal()){

    double upperObjVal(0.0);

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

    MibSSolution *mibSol = NULL;

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

    if(0){
       lSolver->writeLp("tmp");
    }

    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()){
       const double * sol = hSolver->getColSolution();
       double objVal(lSolver->getObjValue() * objSense);
       double etol(etol_);
       double lowerObj = getLowerObj(sol, objSense);  
       
       double * optUpperSolutionOrd = new double[uCols];
       double * optLowerSolutionOrd = new double[lCols];
       
       CoinZeroN(optUpperSolutionOrd, uCols);
       CoinZeroN(optLowerSolutionOrd, lCols);
       
       if(fabs(objVal - lowerObj) < etol){
	  
	  /** Current solution is bilevel feasible **/
	  
	  for(i = 0; i < tCols; i++)
	     upperObjVal += 
		hSolver->getColSolution()[i] * oSolver->getObjCoefficients()[i];
	  
	  mibSol = new MibSSolution(hSolver->getNumCols(),
				    hSolver->getColSolution(),
				    upperObjVal,
				    model);
	  
	  model->storeSolution(BlisSolutionTypeHeuristic, mibSol);
	  mibSol = NULL;
       }
       else{
	  
	  /* solution is not bilevel feasible, create one that is */
	  
	  const double * uSol = hSolver->getColSolution();
	  const double * lSol = lSolver->getColSolution();

//.........这里部分代码省略.........
  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();
      nSolver->loadProblem(*oSolver->getMatrixByCol(),
			   oSolver->getColLower(), oSolver->getColUpper(),
			   oSolver->getObjCoefficients(),
			   oSolver->getRowLower(), oSolver->getRowUpper());
      for(j = 0; j < tCols; j++){
	if(oSolver->isInteger(j))
	  nSolver->setInteger(j);
      }
      

      CoinZeroN(nObjCoeffs, tCols);
      
      for(i = 0; i < lCols; i++){
	index = lColIndices[i];
	nObjCoeffs[index] = lObjCoeffs[i] * lObjSense;
      }