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

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

void
CglClique::selectRowCliques(const OsiSolverInterface& si,int numOriginalRows)
{
   const int numrows = si.getNumRows();
   std::vector<int> clique(numrows, 1);

   int i, j, k;
   
   // First scan through the binary fractional variables and see where do they
   // have a 1 coefficient
   const CoinPackedMatrix& mcol = *si.getMatrixByCol();
   for (j = 0; j < sp_numcols; ++j) {
      const CoinShallowPackedVector& vec = mcol.getVector(sp_orig_col_ind[j]);
      const int* ind = vec.getIndices();
      const double* elem = vec.getElements();
      for (i = vec.getNumElements() - 1; i >= 0; --i) {
	 if (elem[i] != 1.0) {
	    clique[ind[i]] = 0;
	 }
      }
   }

   // Now check the sense and rhs (by checking rowupper) and the rest of the
   // coefficients 
   const CoinPackedMatrix& mrow = *si.getMatrixByRow();
   const double* rub = si.getRowUpper();
   for (i = 0; i < numrows; ++i) {
      if (rub[i] != 1.0||i>=numOriginalRows) {
	 clique[i] = 0;
	 continue;
      }
      if (clique[i] == 1) {
	 const CoinShallowPackedVector& vec = mrow.getVector(i);
	 const double* elem = vec.getElements();
	 for (j = vec.getNumElements() - 1; j >= 0; --j) {
	    if (elem[j] < 0) {
	       clique[i] = 0;
	       break;
	    }
	 }
      }
   }

   // Finally collect the still standing rows into sp_orig_row_ind
   sp_numrows = std::accumulate(clique.begin(), clique.end(), 0);
   sp_orig_row_ind = new int[sp_numrows];
   for (i = 0, k = 0; i < numrows; ++i) {
      if (clique[i] == 1) {
	 sp_orig_row_ind[k++] = i;
      }
   }
}

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

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

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

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

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

//.........这里部分代码省略.........
    // test findGreedyCover 
    CoinPackedVector cover,remainder;
#if 0
    int findgreedy =  kccg.findGreedyCover( 0, krow, b, xstar, cover, remainder );
    assert( findgreedy == 1 );
    int coveri = cover.getNumElements();
    assert( cover.getNumElements() == 2);
    coveri = cover.getIndices()[0];
    assert( cover.getIndices()[0] == 0);
    assert( cover.getIndices()[1] == 1);
    assert( cover.getElements()[0] == 161.0);
    assert( cover.getElements()[1] == 120.0);
    assert( remainder.getNumElements() == 1);
    assert( remainder.getIndices()[0] == 2);
    assert( remainder.getElements()[0] == 68.0);

    // test liftCoverCut
    CoinPackedVector cut;
    double * rowupper = ekk_rowupper(model);
    double cutRhs = cover.getNumElements() - 1.0;
    kccg.liftCoverCut(b, krow.getNumElements(),
      cover, remainder,
      cut);
    assert ( cut.getNumElements() == 3 );
    assert ( cut.getIndices()[0] == 0 );
    assert ( cut.getIndices()[1] == 1 );
    assert ( cut.getIndices()[2] == 2 );
    assert( cut.getElements()[0] == 1 );
    assert( cut.getElements()[1] == 1 );
    assert( eq(cut.getElements()[2], 0.087719) );
    
    // test liftAndUncomplementAndAdd
    OsiCuts cuts;    
    kccg.liftAndUncomplementAndAdd(*siP.getRowUpper()[0],krow,b,complement,0,
      cover,remainder,cuts);   
    int sizerowcuts = cuts.sizeRowCuts();
    assert ( sizerowcuts== 1 );
    OsiRowCut testRowCut = cuts.rowCut(0);
    CoinPackedVector testRowPV = testRowCut.row(); 
    OsiRowCut sampleRowCut;
    const int sampleSize = 3;
    int sampleCols[sampleSize]={0,1,2};
    double sampleElems[sampleSize]={1.0,-1.0,-0.087719};
    sampleRowCut.setRow(sampleSize,sampleCols,sampleElems);
    sampleRowCut.setLb(-DBL_MAX);
    sampleRowCut.setUb(-0.087719);
    bool equiv =  testRowPV.equivalent(sampleRowCut.row(),CoinRelFltEq(1.0e-05) );
    assert ( equiv );
#endif
    
    // test find PseudoJohnAndEllisCover
    cover.setVector(0,NULL, NULL);
    remainder.setVector(0,NULL,NULL);

    rind = ( siP->getRowSense()[0] == 'N' ) ? 1 : 0;
    int findPJE =  kccg.findPseudoJohnAndEllisCover( rind, krow, 
						     b, xstar, cover, remainder );
    assert( findPJE == 1 );
    assert ( cover.getIndices()[0] == 0 );
    assert ( cover.getIndices()[1] == 2 );
    assert ( cover.getElements()[0] == 161 );    
    assert ( cover.getElements()[1] == 68 );    
    assert ( remainder.getIndices()[0] == 1 );
    assert ( remainder.getElements()[0] == 120 );    
    OsiCuts cuts;    
    kccg.liftAndUncomplementAndAdd((*siP).getRowUpper()[rind],krow,b, complement, rind,

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

/* Does a lot of the work,
   Returns 0 if no good, 1 if dj, 2 if clean, 3 if both
   10 if branching on ones away from bound
*/
int
CbcBranchToFixLots::shallWe() const
{
    int returnCode = 0;
    OsiSolverInterface * solver = model_->solver();
    int numberRows = matrixByRow_.getNumRows();
    //if (numberRows!=solver->getNumRows())
    //return 0;
    const double * solution = model_->testSolution();
    const double * lower = solver->getColLower();
    const double * upper = solver->getColUpper();
    const double * dj = solver->getReducedCost();
    int i;
    int numberIntegers = model_->numberIntegers();
    const int * integerVariable = model_->integerVariable();
    if (numberClean_ > 1000000) {
        int wanted = numberClean_ % 1000000;
        int * sort = new int[numberIntegers];
        double * dsort = new double[numberIntegers];
        int nSort = 0;
        for (i = 0; i < numberIntegers; i++) {
            int iColumn = integerVariable[i];
            if (upper[iColumn] > lower[iColumn]) {
                if (!mark_ || !mark_[iColumn]) {
                    double distanceDown = solution[iColumn] - lower[iColumn];
                    double distanceUp = upper[iColumn] - solution[iColumn];
                    double distance = CoinMin(distanceDown, distanceUp);
                    if (distance > 0.001 && distance < 0.5) {
                        dsort[nSort] = distance;
                        sort[nSort++] = iColumn;
                    }
                }
            }
        }
        // sort
        CoinSort_2(dsort, dsort + nSort, sort);
        int n = 0;
        double sum = 0.0;
        for (int k = 0; k < nSort; k++) {
            sum += dsort[k];
            if (sum <= djTolerance_)
                n = k;
            else
                break;
        }
        delete [] sort;
        delete [] dsort;
        return (n >= wanted) ? 10 : 0;
    }
    double integerTolerance =
        model_->getDblParam(CbcModel::CbcIntegerTolerance);
    // make smaller ?
    double tolerance = CoinMin(1.0e-8, integerTolerance);
    // How many fixed are we aiming at
    int wantedFixed = static_cast<int> (static_cast<double>(numberIntegers) * fractionFixed_);
    if (djTolerance_ < 1.0e10) {
        int nSort = 0;
        int numberFixed = 0;
        for (i = 0; i < numberIntegers; i++) {
            int iColumn = integerVariable[i];
            if (upper[iColumn] > lower[iColumn]) {
                if (!mark_ || !mark_[iColumn]) {
                    if (solution[iColumn] < lower[iColumn] + tolerance) {
                        if (dj[iColumn] > djTolerance_) {
                            nSort++;
                        }
                    } else if (solution[iColumn] > upper[iColumn] - tolerance) {
                        if (dj[iColumn] < -djTolerance_) {
                            nSort++;
                        }
                    }
                }
            } else {
                numberFixed++;
            }
        }
        if (numberFixed + nSort < wantedFixed && !alwaysCreate_) {
            returnCode = 0;
        } else if (numberFixed < wantedFixed) {
            returnCode = 1;
        } else {
            returnCode = 0;
        }
    }
    if (numberClean_) {
        // see how many rows clean
        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();
        const double * columnLower = solver->getColLower();
        const double * columnUpper = solver->getColUpper();
//.........这里部分代码省略.........

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

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