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

// Apply bounds
void OsiSolverBranch::applyBounds(OsiSolverInterface &solver, int way) const
{
  int base = way + 1;
  assert(way == -1 || way == 1);
  int numberColumns = solver.getNumCols();
  const double *columnLower = solver.getColLower();
  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]);
      solver.setColLower(iColumn, value);
    } else {
      int iRow = iColumn - numberColumns;
      const double *rowLower = solver.getRowLower();
      double value = CoinMax(bound_[i], rowLower[iRow]);
      solver.setRowLower(iRow, value);
    }
  }
  const double *columnUpper = solver.getColUpper();
  for (i = start_[base + 1]; i < start_[base + 2]; i++) {
    int iColumn = indices_[i];
    if (iColumn < numberColumns) {
      double value = CoinMin(bound_[i], columnUpper[iColumn]);
      solver.setColUpper(iColumn, value);
    } else {
      int iRow = iColumn - numberColumns;
      const double *rowUpper = solver.getRowUpper();
      double value = CoinMin(bound_[i], rowUpper[iRow]);
      solver.setRowUpper(iRow, value);
    }
  }
}

// Useful constructor
CbcFollowOn::CbcFollowOn (CbcModel * model)
        : CbcObject(model)
{
    assert (model);
    OsiSolverInterface * solver = model_->solver();
    matrix_ = *solver->getMatrixByCol();
    matrix_.removeGaps();
    matrix_.setExtraGap(0.0);
    matrixByRow_ = *solver->getMatrixByRow();
    int numberRows = matrix_.getNumRows();

    rhs_ = new int[numberRows];
    int i;
    const double * rowLower = solver->getRowLower();
    const double * rowUpper = solver->getRowUpper();
    // Row copy
    const double * elementByRow = matrixByRow_.getElements();
    const int * column = matrixByRow_.getIndices();
    const CoinBigIndex * rowStart = matrixByRow_.getVectorStarts();
    const int * rowLength = matrixByRow_.getVectorLengths();
    for (i = 0; i < numberRows; i++) {
        rhs_[i] = 0;
        double value = rowLower[i];
        if (value == rowUpper[i]) {
            if (floor(value) == value && value >= 1.0 && value < 10.0) {
                // check elements
                bool good = true;
                for (int j = rowStart[i]; j < rowStart[i] + rowLength[i]; j++) {
                    int iColumn = column[j];
                    if (!solver->isBinary(iColumn))
                        good = false;
                    double elValue = elementByRow[j];
                    if (floor(elValue) != elValue || value < 1.0)
                        good = false;
                }
                if (good)
                    rhs_[i] = static_cast<int> (value);
            }
        }
    }
}

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

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

// Generate cuts
void
CglFakeClique::generateCuts(const OsiSolverInterface& si, OsiCuts & cs,
			const CglTreeInfo info)
{
  if (fakeSolver_) {
    assert (si.getNumCols()==fakeSolver_->getNumCols());
    fakeSolver_->setColLower(si.getColLower());
    const double * solution = si.getColSolution();
    fakeSolver_->setColSolution(solution);
    fakeSolver_->setColUpper(si.getColUpper());
    // get and set branch and bound cutoff
    double cutoff;
    si.getDblParam(OsiDualObjectiveLimit,cutoff);
    fakeSolver_->setDblParam(OsiDualObjectiveLimit,COIN_DBL_MAX);
#ifdef COIN_HAS_CLP
    OsiClpSolverInterface * clpSolver
      = dynamic_cast<OsiClpSolverInterface *> (fakeSolver_);
    if (clpSolver) {
      // fix up fake solver
      const ClpSimplex * siSimplex = clpSolver->getModelPtr();
      // need to set djs
      memcpy(siSimplex->primalColumnSolution(),
	     si.getReducedCost(),si.getNumCols()*sizeof(double));
      fakeSolver_->setDblParam(OsiDualObjectiveLimit,cutoff);
    }
#endif
    const CoinPackedMatrix * matrixByRow = si.getMatrixByRow();
    const double * elementByRow = matrixByRow->getElements();
    const int * column = matrixByRow->getIndices();
    const CoinBigIndex * rowStart = matrixByRow->getVectorStarts();
    const int * rowLength = matrixByRow->getVectorLengths();
    const double * rowUpper = si.getRowUpper();
    const double * rowLower = si.getRowLower();
    
    // Scan all rows looking for possibles
    int numberRows = si.getNumRows();
    double tolerance = 1.0e-3;
    for (int iRow=0;iRow<numberRows;iRow++) {
      CoinBigIndex start = rowStart[iRow];
      CoinBigIndex end = start + rowLength[iRow];
      double upRhs = rowUpper[iRow]; 
      double loRhs = rowLower[iRow]; 
      double sum = 0.0;
      for (CoinBigIndex j=start;j<end;j++) {
	int iColumn=column[j];
	double value = elementByRow[j];
	sum += solution[iColumn]*value;
      }
      if (sum<loRhs-tolerance||sum>upRhs+tolerance) {
	// add as cut
	OsiRowCut rc;
	rc.setLb(loRhs);
	rc.setUb(upRhs);
	rc.setRow(end-start,column+start,elementByRow+start,false);
	CoinAbsFltEq equal(1.0e-12);
	cs.insertIfNotDuplicate(rc,equal);
      }
    }
    CglClique::generateCuts(*fakeSolver_,cs,info);
    if (probing_) {
      probing_->generateCuts(*fakeSolver_,cs,info);
    }
  } else {
    // just use real solver
    CglClique::generateCuts(si,cs,info);
  }
}

// Create a list of rows which might yield cuts
// The possible parameter is a list to cut down search
void CglOddHole::createRowList( const OsiSolverInterface & si,
		      const int * possible)
{
  // Get basic problem information
  int nRows=si.getNumRows(); 
  
  const CoinPackedMatrix * rowCopy = si.getMatrixByRow();

  const int * column = rowCopy->getIndices();
  const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
  const int * rowLength = rowCopy->getVectorLengths(); 
  
  int rowIndex;
  delete [] suitableRows_;
  numberRows_=nRows;

  const double * rowElements = rowCopy->getElements();
  const double * rowupper = si.getRowUpper();
  const double * rowlower = si.getRowLower();
  const double * collower = si.getColLower();
  const double * colupper = si.getColUpper();

  suitableRows_=new int[nRows];
  if (possible) {
    memcpy(suitableRows_,possible,nRows*sizeof(int));
  } else {
    int i;
    for (i=0;i<nRows;i++) {
      suitableRows_[i]=1;
    }
  }
  for (rowIndex=0; rowIndex<nRows; rowIndex++){
    double rhs1=rowupper[rowIndex];
    double rhs2=rowlower[rowIndex];
    if (suitableRows_[rowIndex]) {
      int i;
      bool goodRow=true;
      for (i=rowStart[rowIndex];
	   i<rowStart[rowIndex]+rowLength[rowIndex];i++) {
	int thisCol=column[i];
	if (colupper[thisCol]-collower[thisCol]>epsilon_) {
	  // could allow general integer variables but unlikely
	  if (!si.isBinary(thisCol) ) {
	    goodRow=false;
	    break;
	  }
	  if (fabs(rowElements[i]-1.0)>epsilon_) {
	    goodRow=false;
	    break;
	  }
	} else {
	  rhs1 -= collower[thisCol]*rowElements[i];
	  rhs2 -= collower[thisCol]*rowElements[i];
	}
      }
      if (fabs(rhs1-1.0)>epsilon_&&fabs(rhs2-1.0)>epsilon_) {
	goodRow=false;
      }
      if (goodRow) {
	suitableRows_[rowIndex]=1;
      } else {
	suitableRows_[rowIndex]=0;
      }
    }
  }
}

//-------------------------------------------------------------------------------
// Generate three cycle cuts
//------------------------------------------------------------------- 
void CglOddHole::generateCuts(const OsiSolverInterface & si, OsiCuts & cs,
			      const CglTreeInfo info)
{
  // Get basic problem information
  int nRows=si.getNumRows(); 
  int nCols=si.getNumCols(); 
  
  const CoinPackedMatrix * rowCopy = si.getMatrixByRow();

  // Could do cliques and extra OSL cliques
  // For moment just easy ones
  
  // If no information exists then get a list of suitable rows
  // If it does then suitable rows are subset of information
  
  CglOddHole temp;
  int * checkRow = new int[nRows];
  int i;
  if (!suitableRows_) {
    for (i=0;i<nRows;i++) {
      checkRow[i]=1;
    }
  } else {
    // initialize and extend rows to current size
    memset(checkRow,0,nRows*sizeof(int));
    memcpy(checkRow,suitableRows_,CoinMin(nRows,numberRows_)*sizeof(int));
  }
  temp.createRowList(si,checkRow);
  // now cut down further by only allowing rows with fractional solution
  double * solution = new double[nCols];
  memcpy(solution,si.getColSolution(),nCols*sizeof(double));
  const int * column = rowCopy->getIndices();
  const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
  const int * rowLength = rowCopy->getVectorLengths(); 
  const double * collower = si.getColLower();
  const double * colupper = si.getColUpper();
  int * suitable = temp.suitableRows_;

  // At present I am using new and delete as easier to see arrays in debugger
  int * fixed = new int[nCols]; // mark fixed columns 

  for (i=0;i<nCols;i++) {
    if (si.isBinary(i) ) {
      fixed[i]=0;
      if (colupper[i]-collower[i]<epsilon_) {
	solution[i]=0.0;
	fixed[i]=2;
      } else if (solution[i]<epsilon_) {
	solution[i]=0.0;
	fixed[i]=-1;
      } else if (solution[i]>onetol_) {
	solution[i]=1.0;
	fixed[i]=+1;
      }
    } else {
      //mark as fixed even if not (can not intersect any interesting rows)
      solution[i]=0.0;
      fixed[i]=3;
    }
  }
  // first do packed
  const double * rowlower = si.getRowLower();
  const double * rowupper = si.getRowUpper();
  for (i=0;i<nRows;i++) {
    if (suitable[i]) {
      int k;
      double sum=0.0;
      if (rowupper[i]>1.001) suitable[i]=-1;
      for (k=rowStart[i]; k<rowStart[i]+rowLength[i];k++) {
	int icol=column[k];
	if (!fixed[icol]) sum += solution[icol];
      }
      if (sum<0.9) suitable[i]=-1; //say no good
    }
  }
#ifdef CGL_DEBUG
  const OsiRowCutDebugger * debugger = si.getRowCutDebugger();
  if (debugger&&!debugger->onOptimalPath(si))
    debugger = NULL;
#else
  const OsiRowCutDebugger * debugger = NULL;
#endif
  temp.generateCuts(debugger, *rowCopy,solution,
		    si.getReducedCost(),cs,suitable,fixed,info,true);
  // now cover
  //if no >= then skip
  bool doCover=false;
  int nsuitable=0;
  for (i=0;i<nRows;i++) {
    suitable[i]=abs(suitable[i]);
    if (suitable[i]) {
      int k;
      double sum=0.0;
      if (rowlower[i]<0.999) sum=2.0;
      if (rowupper[i]>1.001) doCover=true;
      for (k=rowStart[i]; k<rowStart[i]+rowLength[i];k++) {
	int icol=column[k];
//.........这里部分代码省略.........

OsiSolverInterface *
expandKnapsack(CoinModel & model, int * whichColumn, int * knapsackStart,
               int * knapsackRow, int &numberKnapsack,
               CglStored & stored, int logLevel,
               int fixedPriority, int SOSPriority, CoinModel & tightenedModel)
{
    int maxTotal = numberKnapsack;
    // load from coin model
    OsiSolverLink *si = new OsiSolverLink();
    OsiSolverInterface * finalModel = NULL;
    si->setDefaultMeshSize(0.001);
    // need some relative granularity
    si->setDefaultBound(100.0);
    si->setDefaultMeshSize(0.01);
    si->setDefaultBound(100000.0);
    si->setIntegerPriority(1000);
    si->setBiLinearPriority(10000);
    si->load(model, true, logLevel);
    // get priorities
    const int * priorities = model.priorities();
    int numberColumns = model.numberColumns();
    if (priorities) {
        OsiObject ** objects = si->objects();
        int numberObjects = si->numberObjects();
        for (int iObj = 0; iObj < numberObjects; iObj++) {
            int iColumn = objects[iObj]->columnNumber();
            if (iColumn >= 0 && iColumn < numberColumns) {
#ifndef NDEBUG
                OsiSimpleInteger * obj =
                    dynamic_cast <OsiSimpleInteger *>(objects[iObj]) ;
#endif
                assert (obj);
                int iPriority = priorities[iColumn];
                if (iPriority > 0)
                    objects[iObj]->setPriority(iPriority);
            }
        }
        if (fixedPriority > 0) {
            si->setFixedPriority(fixedPriority);
        }
        if (SOSPriority < 0)
            SOSPriority = 100000;
    }
    CoinModel coinModel = *si->coinModel();
    assert(coinModel.numberRows() > 0);
    tightenedModel = coinModel;
    int numberRows = coinModel.numberRows();
    // Mark variables
    int * whichKnapsack = new int [numberColumns];
    int iRow, iColumn;
    for (iColumn = 0; iColumn < numberColumns; iColumn++)
        whichKnapsack[iColumn] = -1;
    int kRow;
    bool badModel = false;
    // analyze
    if (logLevel > 1) {
        for (iRow = 0; iRow < numberRows; iRow++) {
            /* Just obvious one at first
            positive non unit coefficients
            all integer
            positive rowUpper
            for now - linear (but further down in code may use nonlinear)
            column bounds should be tight
            */
            //double lower = coinModel.getRowLower(iRow);
            double upper = coinModel.getRowUpper(iRow);
            if (upper < 1.0e10) {
                CoinModelLink triple = coinModel.firstInRow(iRow);
                bool possible = true;
                int n = 0;
                int n1 = 0;
                while (triple.column() >= 0) {
                    int iColumn = triple.column();
                    const char *  el = coinModel.getElementAsString(iRow, iColumn);
                    if (!strcmp("Numeric", el)) {
                        if (coinModel.columnLower(iColumn) == coinModel.columnUpper(iColumn)) {
                            triple = coinModel.next(triple);
                            continue; // fixed
                        }
                        double value = coinModel.getElement(iRow, iColumn);
                        if (value < 0.0) {
                            possible = false;
                        } else {
                            n++;
                            if (value == 1.0)
                                n1++;
                            if (coinModel.columnLower(iColumn) < 0.0)
                                possible = false;
                            if (!coinModel.isInteger(iColumn))
                                possible = false;
                            if (whichKnapsack[iColumn] >= 0)
                                possible = false;
                        }
                    } else {
                        possible = false; // non linear
                    }
                    triple = coinModel.next(triple);
                }
                if (n - n1 > 1 && possible) {
                    double lower = coinModel.getRowLower(iRow);
//.........这里部分代码省略.........

//-----------------------------------------------------------------------------
// Generate Lift-and-Project cuts
//------------------------------------------------------------------- 
void CglLiftAndProject::generateCuts(const OsiSolverInterface& si, OsiCuts& cs,
				     const CglTreeInfo /*info*/)
{
  // Assumes the mixed 0-1 problem 
  //
  //   min {cx: <Atilde,x> >= btilde} 
  //
  // is in canonical form with all bounds,
  // including x_t>=0, -x_t>=-1 for x_t binary,
  // explicitly stated in the constraint matrix. 
  // See ~/COIN/Examples/Cgl2/cgl2.cpp 
  // for a general purpose "convert" function. 

  // Reference [BCC]: Balas, Ceria, and Corneujols,
  // "A lift-and-project cutting plane algorithm
  // for mixed 0-1 program", Math Prog 58, (1993) 
  // 295-324.

  // This implementation uses Normalization 1.

  // Given canonical problem and
  // the lp-relaxation solution, x,
  // the LAP cut generator attempts to construct
  // a cut for every x_j such that 0<x_j<1
  // [BCC:307]
 

  // x_j is the strictly fractional binary variable
  // the cut is generated from
  int j = 0; 

  // Get basic problem information
  // let Atilde be an m by n matrix
  const int m = si.getNumRows(); 
  const int n = si.getNumCols(); 
  const double * x = si.getColSolution();

  // Remember - Atildes may have gaps..
  const CoinPackedMatrix * Atilde = si.getMatrixByRow();
  const double * AtildeElements =  Atilde->getElements();
  const int * AtildeIndices =  Atilde->getIndices();
  const CoinBigIndex * AtildeStarts = Atilde->getVectorStarts();
  const int * AtildeLengths = Atilde->getVectorLengths();  
  const int AtildeFullSize = AtildeStarts[m];
  const double * btilde = si.getRowLower();

  // Set up memory for system (10) [BCC:307]
  // (the problem over the norm intersected 
  //  with the polar cone)
  // 
  // min <<x^T,Atilde^T>,u> + x_ju_0
  // s.t.
  //     <B,w> = (0,...,0,beta_,beta)^T
  //        w  is nonneg for all but the
  //           last two entries, which are free.
  // where 
  // w = (u,v,v_0,u_0)in BCC notation 
  //      u and v are m-vectors; u,v >=0
  //      v_0 and u_0 are free-scalars, and
  //  
  // B = Atilde^T  -Atilde^T  -e_j e_j
  //     btilde^T   e_0^T      0   0
  //     e_0^T      btilde^T   1   0

  // ^T indicates Transpose
  // e_0 is a (AtildeNCols x 1) vector of all zeros 
  // e_j is e_0 with a 1 in the jth position

  // Storing B in column order. B is a (n+2 x 2m+2) matrix 
  // But need to allow for possible gaps in Atilde.
  // At each iteration, only need to change 2 cols and objfunc
  // Sane design of OsiSolverInterface does not permit mucking
  // with matrix.
  // Because we must delete and add cols to alter matrix,
  // and we can only add columns on the end of the matrix
  // put the v_0 and u_0 columns on the end.
  // rather than as described in [BCC]
 
  // Initially allocating B with space for v_0 and u_O cols
  // but not populating, for efficiency.

  // B without u_0 and v_0 is a (n+2 x 2m) size matrix.

  int twoM = 2*m;
  int BNumRows = n+2;
  int BNumCols = twoM+2;
  int BFullSize = 2*AtildeFullSize+twoM+3;
  double * BElements = new double[BFullSize];
  int * BIndices = new int[BFullSize];
  CoinBigIndex * BStarts = new CoinBigIndex [BNumCols+1];
  int * BLengths = new int[BNumCols];


  int i, ij, k=0;
  int nPlus1=n+1;
  int offset = AtildeStarts[m]+m;
  for (i=0; i<m; i++){
//.........这里部分代码省略.........

void
CglLandP::CachedData::getData(const OsiSolverInterface &si)
{
    int nBasics = si.getNumRows();
    int nNonBasics = si.getNumCols();
    if (basis_ != NULL)
        delete basis_;
    basis_ = dynamic_cast<CoinWarmStartBasis *> (si.getWarmStart());
    if (!basis_)
        throw NoBasisError();

    if (nBasics_ > 0 || nBasics != nBasics_)
    {
        delete [] basics_;
        basics_ = NULL;
    }
    if (basics_ == NULL)
    {
        basics_ = new int[nBasics];
        nBasics_ = nBasics;
    }

    if (nNonBasics_ > 0 || nNonBasics != nNonBasics_)
    {
        delete [] nonBasics_;
        nonBasics_ = NULL;
    }
    if (nonBasics_ == NULL)
    {
        nonBasics_ = new int[nNonBasics];
        nNonBasics_ = nNonBasics;
    }
    int n = nBasics + nNonBasics;
    if ( nBasics_ + nNonBasics_ > 0 || nBasics_ + nNonBasics_ != n)
    {
        delete [] colsol_;
        delete [] integers_;
        integers_ = NULL;
        colsol_ = NULL;
        slacks_ = NULL;
    }
    if (colsol_ == NULL)
    {
        colsol_ = new double[n];
        slacks_ = &colsol_[nNonBasics];
    }

    if (integers_ == NULL)
    {
        integers_ = new bool[n];
    }

    const double * rowLower = si.getRowLower();
    const double * rowUpper = si.getRowUpper();
    //determine which slacks are integer
    const CoinPackedMatrix * m = si.getMatrixByCol();
    const double * elems = m->getElements();
    const int * inds = m->getIndices();
    const CoinBigIndex * starts = m->getVectorStarts();
    const int * lengths = m->getVectorLengths();
    //    int numElems = m->getNumElements();
    int numCols = m->getNumCols();
    assert(numCols == nNonBasics_);
    //   int numRows = m->getNumRows();
    CoinFillN(integers_ ,n, true);
    for (int i = 0 ;  i < numCols ; i++)
    {
        if (si.isContinuous(i))
            integers_[i] = false;
    }
    bool * integerSlacks = integers_ + numCols;
    for (int i = 0 ; i < nBasics ; i++)
    {
        if (rowLower[i] > -1e50 && INT_INFEAS(rowLower[i]) > 1e-15)
            integerSlacks[i] = false;
        if (rowUpper[i] < 1e50 && INT_INFEAS(rowUpper[i]) > 1e-15)
            integerSlacks[i] = false;
    }
    for (int i = 0 ;  i < numCols ; i++)
    {
        CoinBigIndex end = starts[i] + lengths[i];
        if (integers_[i])
        {
            for (CoinBigIndex k=starts[i] ; k < end; k++)
            {
                if (integerSlacks[inds[k]] && INT_INFEAS(elems[k])>1e-15 )
                    integerSlacks[inds[k]] = false;
            }
        }
        else
        {
            for (CoinBigIndex k=starts[i] ; k < end; k++)
            {
                if (integerSlacks[inds[k]])
                    integerSlacks[inds[k]] = false;
            }
        }
    }

    CoinCopyN(si.getColSolution(), si.getNumCols(), colsol_);
//.........这里部分代码省略.........

// See if rounding will give solution
// Sets value of solution
// Assumes rhs for original matrix still okay
// At present only works with integers 
// Fix values if asked for
// Returns 1 if solution, 0 if not
int
AbcRounding::solution(double & solutionValue,
		      double * betterSolution)
{
    // Get a copy of original matrix (and by row for rounding);
    matrix_ = *(model_->solver()->getMatrixByCol());
    matrixByRow_ = *(model_->solver()->getMatrixByRow());
    seed_=1;

    OsiSolverInterface * solver = model_->solver();
    const double * lower = solver->getColLower();
    const double * upper = solver->getColUpper();
    const double * rowLower = solver->getRowLower();
    const double * rowUpper = solver->getRowUpper();
    const double * solution = solver->getColSolution();
    const double * objective = solver->getObjCoefficients();
    double integerTolerance = 1.0e-5;
    //model_->getDblParam(AbcModel::AbcIntegerTolerance);
    double primalTolerance;
    solver->getDblParam(OsiPrimalTolerance, primalTolerance);

    int numberRows = matrix_.getNumRows();

    int numberIntegers = model_->numberIntegers();
    const int * integerVariable = model_->integerVariable();
    int i;
    double direction = solver->getObjSense();
    double newSolutionValue = direction * solver->getObjValue();
    int returnCode = 0;

    // Column copy
    const double * element = matrix_.getElements();
    const int * row = matrix_.getIndices();
    const int * columnStart = matrix_.getVectorStarts();
    const int * columnLength = matrix_.getVectorLengths();
    // Row copy
    const double * elementByRow = matrixByRow_.getElements();
    const int * column = matrixByRow_.getIndices();
    const int * rowStart = matrixByRow_.getVectorStarts();
    const int * rowLength = matrixByRow_.getVectorLengths();

    // Get solution array for heuristic solution
    int numberColumns = solver->getNumCols();
    double * newSolution = new double [numberColumns];
    memcpy(newSolution, solution, numberColumns * sizeof(double));

    double * rowActivity = new double[numberRows];
    memset(rowActivity, 0, numberRows*sizeof(double));
    for (i = 0; i < numberColumns; i++) {
	int j;
	double value = newSolution[i];
	if (value) {
	    for (j = columnStart[i];
		 j < columnStart[i] + columnLength[i]; j++) {
		int iRow = row[j];
		rowActivity[iRow] += value*element[j];
	    }
	}
    }
    // check was feasible - if not adjust (cleaning may move)
    for (i = 0; i < numberRows; i++) {
	if(rowActivity[i] < rowLower[i]) {
	    //assert (rowActivity[i]>rowLower[i]-1000.0*primalTolerance);
	    rowActivity[i] = rowLower[i];
	} else if(rowActivity[i] > rowUpper[i]) {
	    //assert (rowActivity[i]<rowUpper[i]+1000.0*primalTolerance);
	    rowActivity[i] = rowUpper[i];
	}
    }
    for (i = 0; i < numberIntegers; i++) {
	int iColumn = integerVariable[i];
	double value = newSolution[iColumn];
	if (fabs(floor(value + 0.5) - value) > integerTolerance) {
	    double below = floor(value);
	    double newValue = newSolution[iColumn];
	    double cost = direction * objective[iColumn];
	    double move;
	    if (cost > 0.0) {
		// try up
		move = 1.0 - (value - below);
	    } else if (cost < 0.0) {
		// try down
		move = below - value;
	    } else {
		// won't be able to move unless we can grab another variable
		// just for now go down
		move = below-value;
	    }
	    newValue += move;
	    newSolution[iColumn] = newValue;
	    newSolutionValue += move * cost;
	    int j;
	    for (j = columnStart[iColumn];
		 j < columnStart[iColumn] + columnLength[iColumn]; j++) {
//.........这里部分代码省略.........

// This looks at solution and sets bounds to contain solution
void 
CbcLink::feasibleRegion()
{
  int j;
  int firstNonZero=-1;
  int lastNonZero = -1;
  OsiSolverInterface * solver = model_->solver();
  const double * solution = model_->testSolution();
  const double * upper = solver->getColUpper();
  double integerTolerance = 
    model_->getDblParam(CbcModel::CbcIntegerTolerance);
  double weight = 0.0;
  double sum =0.0;

  int base=0;
  for (j=0;j<numberMembers_;j++) {
    for (int k=0;k<numberLinks_;k++) {
      int iColumn = which_[base+k];
      double value = CoinMax(0.0,solution[iColumn]);
      sum += value;
      if (value>integerTolerance&&upper[iColumn]) {
        weight += weights_[j]*value;
        if (firstNonZero<0)
          firstNonZero=j;
        lastNonZero=j;
      }
    }
    base += numberLinks_;
  }
#ifdef DISTANCE
  if (lastNonZero-firstNonZero>sosType_-1) {
    /* may still be satisfied.
       For LOS type 2 we might wish to move coding around
       and keep initial info in model_ for speed
    */
    int iWhere;
    bool possible=false;
    for (iWhere=firstNonZero;iWhere<=lastNonZero;iWhere++) {
      if (fabs(weight-weights_[iWhere])<1.0e-8) {
	possible=true;
	break;
      }
    }
    if (possible) {
      // One could move some of this (+ arrays) into model_
      const CoinPackedMatrix * matrix = solver->getMatrixByCol();
      const double * element = matrix->getMutableElements();
      const int * row = matrix->getIndices();
      const CoinBigIndex * columnStart = matrix->getVectorStarts();
      const int * columnLength = matrix->getVectorLengths();
      const double * rowSolution = solver->getRowActivity();
      const double * rowLower = solver->getRowLower();
      const double * rowUpper = solver->getRowUpper();
      int numberRows = matrix->getNumRows();
      double * array = new double [numberRows];
      CoinZeroN(array,numberRows);
      int * which = new int [numberRows];
      int n=0;
      int base=numberLinks_*firstNonZero;
      for (j=firstNonZero;j<=lastNonZero;j++) {
	for (int k=0;k<numberLinks_;k++) {
	  int iColumn = which_[base+k];
	  double value = CoinMax(0.0,solution[iColumn]);
	  if (value>integerTolerance&&upper[iColumn]) {
	    value = CoinMin(value,upper[iColumn]);
	    for (int j=columnStart[iColumn];j<columnStart[iColumn]+columnLength[iColumn];j++) {
	      int iRow = row[j];
	      double a = array[iRow];
	      if (a) {
		a += value*element[j];
		if (!a)
		  a = 1.0e-100;
	      } else {
		which[n++]=iRow;
		a=value*element[j];
		assert (a);
	      }
	      array[iRow]=a;
	    }
	  }
	}
	base += numberLinks_;
      }
      base=numberLinks_*iWhere;
      for (int k=0;k<numberLinks_;k++) {
	int iColumn = which_[base+k];
	const double value = 1.0;
	for (int j=columnStart[iColumn];j<columnStart[iColumn]+columnLength[iColumn];j++) {
	  int iRow = row[j];
	  double a = array[iRow];
	  if (a) {
	    a -= value*element[j];
	    if (!a)
	      a = 1.0e-100;
	  } else {
	    which[n++]=iRow;
	    a=-value*element[j];
	    assert (a);
	  }
//.........这里部分代码省略.........

// Infeasibility - large is 0.5
double 
CbcLink::infeasibility(int & preferredWay) const
{
  int j;
  int firstNonZero=-1;
  int lastNonZero = -1;
  OsiSolverInterface * solver = model_->solver();
  const double * solution = model_->testSolution();
  //const double * lower = solver->getColLower();
  const double * upper = solver->getColUpper();
  double integerTolerance = 
    model_->getDblParam(CbcModel::CbcIntegerTolerance);
  double weight = 0.0;
  double sum =0.0;

  // check bounds etc
  double lastWeight=-1.0e100;
  int base=0;
  for (j=0;j<numberMembers_;j++) {
    for (int k=0;k<numberLinks_;k++) {
      int iColumn = which_[base+k];
      //if (lower[iColumn])
      //throw CoinError("Non zero lower bound in CBCLink","infeasibility","CbcLink");
      if (lastWeight>=weights_[j]-1.0e-7)
        throw CoinError("Weights too close together in CBCLink","infeasibility","CbcLink");
      double value = CoinMax(0.0,solution[iColumn]);
      sum += value;
      if (value>integerTolerance&&upper[iColumn]) {
        // Possibly due to scaling a fixed variable might slip through
        if (value>upper[iColumn]+1.0e-8) {
          // Could change to #ifdef CBC_DEBUG
#ifndef NDEBUG
          if (model_->messageHandler()->logLevel()>1)
            printf("** Variable %d (%d) has value %g and upper bound of %g\n",
                   iColumn,j,value,upper[iColumn]);
#endif
        } 
	value = CoinMin(value,upper[iColumn]);
        weight += weights_[j]*value;
        if (firstNonZero<0)
          firstNonZero=j;
        lastNonZero=j;
      }
    }
    base += numberLinks_;
  }
  double valueInfeasibility;
  preferredWay=1;
  if (lastNonZero-firstNonZero>=sosType_) {
    // find where to branch
    assert (sum>0.0);
    weight /= sum;
    valueInfeasibility = lastNonZero-firstNonZero+1;
    valueInfeasibility *= 0.5/((double) numberMembers_);
    //#define DISTANCE
#ifdef DISTANCE
    assert (sosType_==1); // code up
    /* may still be satisfied.
       For LOS type 2 we might wish to move coding around
       and keep initial info in model_ for speed
    */
    int iWhere;
    bool possible=false;
    for (iWhere=firstNonZero;iWhere<=lastNonZero;iWhere++) {
      if (fabs(weight-weights_[iWhere])<1.0e-8) {
	possible=true;
	break;
      }
    }
    if (possible) {
      // One could move some of this (+ arrays) into model_
      const CoinPackedMatrix * matrix = solver->getMatrixByCol();
      const double * element = matrix->getMutableElements();
      const int * row = matrix->getIndices();
      const CoinBigIndex * columnStart = matrix->getVectorStarts();
      const int * columnLength = matrix->getVectorLengths();
      const double * rowSolution = solver->getRowActivity();
      const double * rowLower = solver->getRowLower();
      const double * rowUpper = solver->getRowUpper();
      int numberRows = matrix->getNumRows();
      double * array = new double [numberRows];
      CoinZeroN(array,numberRows);
      int * which = new int [numberRows];
      int n=0;
      int base=numberLinks_*firstNonZero;
      for (j=firstNonZero;j<=lastNonZero;j++) {
	for (int k=0;k<numberLinks_;k++) {
	  int iColumn = which_[base+k];
	  double value = CoinMax(0.0,solution[iColumn]);
	  if (value>integerTolerance&&upper[iColumn]) {
	    value = CoinMin(value,upper[iColumn]);
	    for (int j=columnStart[iColumn];j<columnStart[iColumn]+columnLength[iColumn];j++) {
	      int iRow = row[j];
	      double a = array[iRow];
	      if (a) {
		a += value*element[j];
		if (!a)
		  a = 1.0e-100;
	      } else {
//.........这里部分代码省略.........

/*
  First tries setting a variable to better value.  If feasible then
  tries setting others.  If not feasible then tries swaps
  Returns 1 if solution, 0 if not 
  The main body of this routine implements an O((q^2)/2) brute force search
  around the current solution, for q = number of integer variables. Call this
  the inc/dec heuristic.  For each integer variable x<i>, first decrement the
  value. Then, for integer variables x<i+1>, ..., x<q-1>, try increment and
  decrement. If one of these permutations produces a better solution,
  remember it.  Then repeat, with x<i> incremented. If we find a better
  solution, update our notion of current solution and continue.

  The net effect is a greedy walk: As each improving pair is found, the
  current solution is updated and the search continues from this updated
  solution.

  Way down at the end, we call solutionFix, which will create a drastically
  restricted problem based on variables marked as used, then do mini-BaC on
  the restricted problem. This can occur even if we don't try the inc/dec
  heuristic. This would be more obvious if the inc/dec heuristic were broken
  out as a separate routine and solutionFix had a name that reflected where
  it was headed.

  The return code of 0 is grossly overloaded, because it maps to a return
  code of 0 from solutionFix, which is itself grossly overloaded. See
  comments in solutionFix and in CbcHeuristic::smallBranchAndBound.
  */
int
CbcHeuristicLocal::solution(double & solutionValue,
                            double * betterSolution)
{
/*
  Execute only if a new solution has been discovered since the last time we
  were called.
*/

    numCouldRun_++;
    // See if frequency kills off idea
    int swap = swap_%100;
    int skip = swap_/100;
    int nodeCount = model_->getNodeCount();
    if (nodeCount<lastRunDeep_+skip && nodeCount != lastRunDeep_+1) 
      return 0;
    if (numberSolutions_ == model_->getSolutionCount() &&
	(numberSolutions_ == howOftenShallow_ ||
	 nodeCount < lastRunDeep_+2*skip))
        return 0;
    howOftenShallow_ = numberSolutions_;
    numberSolutions_ = model_->getSolutionCount();
    if (nodeCount<lastRunDeep_+skip ) 
      return 0;
    lastRunDeep_ = nodeCount;
    howOftenShallow_ = numberSolutions_;

    if ((swap%10) == 2) {
        // try merge
        return solutionFix( solutionValue, betterSolution, NULL);
    }
/*
  Exclude long (column), thin (row) systems.

  Given the n^2 nature of the search, more than 100,000 columns could get
  expensive. But I don't yet see the rationale for the second part of the
  condition (cols > 10*rows). And cost is proportional to number of integer
  variables --- shouldn't we use that?

  Why wait until we have more than one solution?
*/
    if ((model_->getNumCols() > 100000 && model_->getNumCols() >
            10*model_->getNumRows()) || numberSolutions_ <= 1)
        return 0; // probably not worth it
    // worth trying

    OsiSolverInterface * solver = model_->solver();
    const double * rowLower = solver->getRowLower();
    const double * rowUpper = solver->getRowUpper();
    const double * solution = model_->bestSolution();
/*
  Shouldn't this test be redundant if we've already checked that
  numberSolutions_ > 1? Stronger: shouldn't this be an assertion?
*/
    if (!solution)
        return 0; // No solution found yet
    const double * objective = solver->getObjCoefficients();
    double primalTolerance;
    solver->getDblParam(OsiPrimalTolerance, primalTolerance);

    int numberRows = matrix_.getNumRows();

    int numberIntegers = model_->numberIntegers();
    const int * integerVariable = model_->integerVariable();

    int i;
    double direction = solver->getObjSense();
    double newSolutionValue = model_->getObjValue() * direction;
    int returnCode = 0;
    numRuns_++;
    // Column copy
    const double * element = matrix_.getElements();
    const int * row = matrix_.getIndices();
//.........这里部分代码省略.........

int * analyze(OsiClpSolverInterface * solverMod, int & numberChanged,
		     double & increment, bool changeInt,
		     CoinMessageHandler * generalMessageHandler, bool noPrinting)
{
    bool noPrinting_ = noPrinting;
    OsiSolverInterface * solver = solverMod->clone();
    char generalPrint[200];
    if (0) {
        // just get increment
        CbcModel model(*solver);
        model.analyzeObjective();
        double increment2 = model.getCutoffIncrement();
        printf("initial cutoff increment %g\n", increment2);
    }
    const double *objective = solver->getObjCoefficients() ;
    const double *lower = solver->getColLower() ;
    const double *upper = solver->getColUpper() ;
    int numberColumns = solver->getNumCols() ;
    int numberRows = solver->getNumRows();
    double direction = solver->getObjSense();
    int iRow, iColumn;

    // Row copy
    CoinPackedMatrix matrixByRow(*solver->getMatrixByRow());
    const double * elementByRow = matrixByRow.getElements();
    const int * column = matrixByRow.getIndices();
    const CoinBigIndex * rowStart = matrixByRow.getVectorStarts();
    const int * rowLength = matrixByRow.getVectorLengths();

    // Column copy
    CoinPackedMatrix  matrixByCol(*solver->getMatrixByCol());
    const double * element = matrixByCol.getElements();
    const int * row = matrixByCol.getIndices();
    const CoinBigIndex * columnStart = matrixByCol.getVectorStarts();
    const int * columnLength = matrixByCol.getVectorLengths();

    const double * rowLower = solver->getRowLower();
    const double * rowUpper = solver->getRowUpper();

    char * ignore = new char [numberRows];
    int * changed = new int[numberColumns];
    int * which = new int[numberRows];
    double * changeRhs = new double[numberRows];
    memset(changeRhs, 0, numberRows*sizeof(double));
    memset(ignore, 0, numberRows);
    numberChanged = 0;
    int numberInteger = 0;
    for (iColumn = 0; iColumn < numberColumns; iColumn++) {
        if (upper[iColumn] > lower[iColumn] + 1.0e-8 && solver->isInteger(iColumn))
            numberInteger++;
    }
    bool finished = false;
    while (!finished) {
        int saveNumberChanged = numberChanged;
        for (iRow = 0; iRow < numberRows; iRow++) {
            int numberContinuous = 0;
            double value1 = 0.0, value2 = 0.0;
            bool allIntegerCoeff = true;
            double sumFixed = 0.0;
            int jColumn1 = -1, jColumn2 = -1;
            for (CoinBigIndex j = rowStart[iRow]; j < rowStart[iRow] + rowLength[iRow]; j++) {
                int jColumn = column[j];
                double value = elementByRow[j];
                if (upper[jColumn] > lower[jColumn] + 1.0e-8) {
                    if (!solver->isInteger(jColumn)) {
                        if (numberContinuous == 0) {
                            jColumn1 = jColumn;
                            value1 = value;
                        } else {
                            jColumn2 = jColumn;
                            value2 = value;
                        }
                        numberContinuous++;
                    } else {
                        if (fabs(value - floor(value + 0.5)) > 1.0e-12)
                            allIntegerCoeff = false;
                    }
                } else {
                    sumFixed += lower[jColumn] * value;
                }
            }
            double low = rowLower[iRow];
            if (low > -1.0e20) {
                low -= sumFixed;
                if (fabs(low - floor(low + 0.5)) > 1.0e-12)
                    allIntegerCoeff = false;
            }
            double up = rowUpper[iRow];
            if (up < 1.0e20) {
                up -= sumFixed;
                if (fabs(up - floor(up + 0.5)) > 1.0e-12)
                    allIntegerCoeff = false;
            }
            if (!allIntegerCoeff)
                continue; // can't do
            if (numberContinuous == 1) {
                // see if really integer
                // This does not allow for complicated cases
                if (low == up) {
                    if (fabs(value1) > 1.0e-3) {
//.........这里部分代码省略.........