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

//--------------------------------------------------------------------------
// test cut debugger methods.
void
OsiRowCutDebuggerUnitTest(const OsiSolverInterface * baseSiP, const std::string & mpsDir)
{
  
  CoinRelFltEq eq;
  
  // Test default constructor
  {
    OsiRowCutDebugger r;
    OSIUNITTEST_ASSERT_ERROR(r.integerVariable_ == NULL, {}, "osirowcutdebugger", "default constructor");
    OSIUNITTEST_ASSERT_ERROR(r.knownSolution_   == NULL, {}, "osirowcutdebugger", "default constructor");
    OSIUNITTEST_ASSERT_ERROR(r.numberColumns_   == 0,    {}, "osirowcutdebugger", "default constructor");
  }
  
  {
    // Get non trivial instance
    OsiSolverInterface * imP = baseSiP->clone();
    std::string fn = mpsDir+"exmip1";
    imP->readMps(fn.c_str(),"mps");
    OSIUNITTEST_ASSERT_ERROR(imP->getNumRows() == 5, {}, "osirowcutdebugger", "read exmip1");
    
    /*
      Activate the debugger. The garbled name here is deliberate; the
      debugger should isolate the portion of the string between '/' and
      '.' (in normal use, this would be the base file name, stripped of
      the prefix and extension).
    */
    imP->activateRowCutDebugger("ab cd /x/ /exmip1.asc");
    
    int i;
    
    // return debugger
    const OsiRowCutDebugger * debugger = imP->getRowCutDebugger();
    OSIUNITTEST_ASSERT_ERROR(debugger != NULL, {}, "osirowcutdebugger", "return debugger");
    OSIUNITTEST_ASSERT_ERROR(debugger->numberColumns_ == 8, {}, "osirowcutdebugger", "return debugger");
    
    const bool type[]={0,0,1,1,0,0,0,0};
    const double values[]= {2.5, 0, 1, 1, 0.5, 3, 0, 0.26315789473684253};
    CoinPackedVector objCoefs(8,imP->getObjCoefficients());
   
    bool type_ok = true;
#if 0
    for (i=0;i<8;i++)
      type_ok &= type[i] == debugger->integerVariable_[i];
    OSIUNITTEST_ASSERT_ERROR(type_ok, {}, "osirowcutdebugger", "???");
#endif
    
    double objValue = objCoefs.dotProduct(values);
    double debuggerObjValue = objCoefs.dotProduct(debugger->knownSolution_);
    OSIUNITTEST_ASSERT_ERROR(eq(objValue, debuggerObjValue), {}, "osirowcutdebugger", "objective value");
    
    OsiRowCutDebugger rhs;
    {
      OsiRowCutDebugger rC1(*debugger);

      OSIUNITTEST_ASSERT_ERROR(rC1.numberColumns_ == 8, {}, "osirowcutdebugger", "copy constructor");
      type_ok = true;
      for (i=0;i<8;i++)
      	type_ok &= type[i] == rC1.integerVariable_[i];
      OSIUNITTEST_ASSERT_ERROR(type_ok, {}, "osirowcutdebugger", "copy constructor");
      OSIUNITTEST_ASSERT_ERROR(eq(objValue,objCoefs.dotProduct(rC1.knownSolution_)), {}, "osirowcutdebugger", "copy constructor");
      
      rhs = rC1;
      OSIUNITTEST_ASSERT_ERROR(rhs.numberColumns_ == 8, {}, "osirowcutdebugger", "assignment operator");
      type_ok = true;
      for (i=0;i<8;i++)
      	type_ok &= type[i] == rhs.integerVariable_[i];
      OSIUNITTEST_ASSERT_ERROR(type_ok, {}, "osirowcutdebugger", "assignment operator");
      OSIUNITTEST_ASSERT_ERROR(eq(objValue,objCoefs.dotProduct(rhs.knownSolution_)), {}, "osirowcutdebugger", "assignment operator");
    }
    // Test that rhs has correct values even though lhs has gone out of scope
    OSIUNITTEST_ASSERT_ERROR(rhs.numberColumns_ == 8, {}, "osirowcutdebugger", "assignment operator");
    type_ok = true;
    for (i=0;i<8;i++)
    	type_ok &= type[i] == rhs.integerVariable_[i];
    OSIUNITTEST_ASSERT_ERROR(type_ok, {}, "osirowcutdebugger", "assignment operator");
    OSIUNITTEST_ASSERT_ERROR(eq(objValue,objCoefs.dotProduct(rhs.knownSolution_)), {}, "osirowcutdebugger", "assignment operator");

    OsiRowCut cut[2];
    
    const int ne = 3;
    int inx[ne] = { 0, 2, 3 };
    double el[ne] = { 1., 1., 1. };
    cut[0].setRow(ne,inx,el);
    cut[0].setUb(5.);
    
    el[1]=5;
    cut[1].setRow(ne,inx,el);
    cut[1].setUb(5);
    OsiCuts cs; 
    cs.insert(cut[0]);
    cs.insert(cut[1]);
    OSIUNITTEST_ASSERT_ERROR(!debugger->invalidCut(cut[0]), {}, "osirowcutdebugger", "recognize (in)valid cut");
    OSIUNITTEST_ASSERT_ERROR( debugger->invalidCut(cut[1]), {}, "osirowcutdebugger", "recognize (in)valid cut");
    OSIUNITTEST_ASSERT_ERROR(debugger->validateCuts(cs,0,2) == 1, {}, "osirowcutdebugger", "recognize (in)valid cut");
    OSIUNITTEST_ASSERT_ERROR(debugger->validateCuts(cs,0,1) == 0, {}, "osirowcutdebugger", "recognize (in)valid cut");
    delete imP;
  }
//.........这里部分代码省略.........

// inner part of dive
int 
CbcHeuristicDive::solution(double & solutionValue, int & numberNodes,
			   int & numberCuts, OsiRowCut ** cuts,
			   CbcSubProblem ** & nodes,
			   double * newSolution)
{
#ifdef DIVE_DEBUG
    int nRoundInfeasible = 0;
    int nRoundFeasible = 0;
#endif
    int reasonToStop = 0;
    double time1 = CoinCpuTime();
    int numberSimplexIterations = 0;
    int maxSimplexIterations = (model_->getNodeCount()) ? maxSimplexIterations_
                               : maxSimplexIterationsAtRoot_;
    // but can't be exactly coin_int_max
    maxSimplexIterations = CoinMin(maxSimplexIterations,COIN_INT_MAX>>3);
    OsiSolverInterface * solver = cloneBut(6); // was model_->solver()->clone();
# ifdef COIN_HAS_CLP
    OsiClpSolverInterface * clpSolver
    = dynamic_cast<OsiClpSolverInterface *> (solver);
    if (clpSolver) {
      ClpSimplex * clpSimplex = clpSolver->getModelPtr();
      int oneSolveIts = clpSimplex->maximumIterations();
      oneSolveIts = CoinMin(1000+2*(clpSimplex->numberRows()+clpSimplex->numberColumns()),oneSolveIts);
      clpSimplex->setMaximumIterations(oneSolveIts);
      if (!nodes) {
        // say give up easily
        clpSimplex->setMoreSpecialOptions(clpSimplex->moreSpecialOptions() | 64);
      } else {
	// get ray
	int specialOptions = clpSimplex->specialOptions();
	specialOptions &= ~0x3100000;
	specialOptions |= 32;
        clpSimplex->setSpecialOptions(specialOptions);
        clpSolver->setSpecialOptions(clpSolver->specialOptions() | 1048576);
	if ((model_->moreSpecialOptions()&16777216)!=0) {
	  // cutoff is constraint
	  clpSolver->setDblParam(OsiDualObjectiveLimit, COIN_DBL_MAX);
	}
      }
    }
# endif
    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 = model_->getDblParam(CbcModel::CbcIntegerTolerance);
    double primalTolerance;
    solver->getDblParam(OsiPrimalTolerance, primalTolerance);

    int numberRows = matrix_.getNumRows();
    assert (numberRows <= solver->getNumRows());
    int numberIntegers = model_->numberIntegers();
    const int * integerVariable = model_->integerVariable();
    double direction = solver->getObjSense(); // 1 for min, -1 for max
    double newSolutionValue = direction * solver->getObjValue();
    int returnCode = 0;
    // Column copy
    const double * element = matrix_.getElements();
    const int * row = matrix_.getIndices();
    const CoinBigIndex * columnStart = matrix_.getVectorStarts();
    const int * columnLength = matrix_.getVectorLengths();
#ifdef DIVE_FIX_BINARY_VARIABLES
    // Row copy
    const double * elementByRow = matrixByRow_.getElements();
    const int * column = matrixByRow_.getIndices();
    const CoinBigIndex * rowStart = matrixByRow_.getVectorStarts();
    const int * rowLength = matrixByRow_.getVectorLengths();
#endif

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

    // vectors to store the latest variables fixed at their bounds
    int* columnFixed = new int [numberIntegers];
    double* originalBound = new double [numberIntegers+2*numberColumns];
    double * lowerBefore = originalBound+numberIntegers;
    double * upperBefore = lowerBefore+numberColumns;
    memcpy(lowerBefore,lower,numberColumns*sizeof(double));
    memcpy(upperBefore,upper,numberColumns*sizeof(double));
    double * lastDjs=newSolution+numberColumns;
    bool * fixedAtLowerBound = new bool [numberIntegers];
    PseudoReducedCost * candidate = new PseudoReducedCost [numberIntegers];
    double * random = new double [numberIntegers];

    int maxNumberAtBoundToFix = static_cast<int> (floor(percentageToFix_ * numberIntegers));
    assert (!maxNumberAtBoundToFix||!nodes);

    // count how many fractional variables
    int numberFractionalVariables = 0;
    for (int i = 0; i < numberIntegers; i++) {
        random[i] = randomNumberGenerator_.randomDouble() + 0.3;
        int iColumn = integerVariable[i];
        double value = newSolution[iColumn];
        if (fabs(floor(value + 0.5) - value) > integerTolerance) {
//.........这里部分代码省略.........

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

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

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

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

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

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


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

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

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

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

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

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

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

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

  lSolver->branchAndBound();

  if(lSolver->isProvenOptimal()){

    double objVal(0.0);

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

}

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

//.........这里部分代码省略.........
		       -6.7220534694101275e-18, 5.3125906451789717e-18, 
		       1, 0, 1.9298798670241979e-17, 0, 0, 0,
		       7.8875708048320448e-18, 0.5, 0, 
		       0.85999999999999999, 1, 1, 0.57999999999999996,
		       1, 0, 1, 0, 0.25, 0, 0.67500000000000004};
    siP->setColSolution(mycs);
#ifdef CGL_DEBUG
    printf("\n\nOrig LP min=%f\n",lpRelaxBefore);
#endif
    
    OsiCuts cuts;    
    
    // Test generateCuts method
    kccg.generateCuts(*siP,cuts);
    OsiSolverInterface::ApplyCutsReturnCode rc = siP->applyCuts(cuts);
    
    siP->resolve();
    double lpRelaxAfter=siP->getObjValue(); 
    assert( eq(lpRelaxAfter, 2829.0597826086955) );
#ifdef CGL_DEBUG
    printf("\n\nOrig LP min=%f\n",lpRelaxBefore);
    printf("\n\nFinal LP min=%f\n",lpRelaxAfter);
#endif
    assert( lpRelaxBefore < lpRelaxAfter );
    
    // the CoinPackedVector p0033 is the optimal
    // IP solution to the miplib problem p0033
    int objIndices[14] = { 
       0,  6,  7,  9, 13, 17, 18,
      22, 24, 25, 26, 27, 28, 29 };
    CoinPackedVector p0033(14,objIndices,1.0);

    // Sanity check
    const double *  objective=siP->getObjCoefficients();
    double ofv =0 ;
    int r;
    for (r=0; r<nCols; r++){
      ofv=ofv + p0033[r]*objective[r];
    }
    CoinRelFltEq eq;
    assert( eq(ofv,3089.0) );

    int nRowCuts = cuts.sizeRowCuts();
    OsiRowCut rcut;
    CoinPackedVector rpv;
    for (i=0; i<nRowCuts; i++){
      rcut = cuts.rowCut(i);
      rpv = rcut.row();
      double p0033Sum = (rpv*p0033).sum();
      assert (p0033Sum <= rcut.ub() );
    }
  
    delete siP;
  } 

  // if a debug file is there then look at it
  {
    FILE * fp = fopen("knapsack.debug","r");
    if (fp) {
      int ncol,nel;
      double up;
      int x = fscanf(fp,"%d %d %lg",&ncol,&nel,&up);
      if (x<=0)
	throw("bad fscanf");
      printf("%d columns, %d elements, upper %g\n",ncol,nel,up);
      double * sol1 = new double[nel];

int
CbcHeuristicNaive::solution(double & solutionValue,
                            double * betterSolution)
{
    numCouldRun_++;
    // See if to do
    bool atRoot = model_->getNodeCount() == 0;
    int passNumber = model_->getCurrentPassNumber();
    if (!when() || (when() == 1 && model_->phase() != 1) || !atRoot || passNumber != 1)
        return 0; // switched off
    // Don't do if it was this heuristic which found solution!
    if (this == model_->lastHeuristic())
        return 0;
    numRuns_++;
    double cutoff;
    model_->solver()->getDblParam(OsiDualObjectiveLimit, cutoff);
    double direction = model_->solver()->getObjSense();
    cutoff *= direction;
    cutoff = CoinMin(cutoff, solutionValue);
    OsiSolverInterface * solver = model_->continuousSolver();
    if (!solver)
        solver = model_->solver();
    const double * colLower = solver->getColLower();
    const double * colUpper = solver->getColUpper();
    const double * objective = solver->getObjCoefficients();

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

    int i;
    bool solutionFound = false;
    CoinWarmStartBasis saveBasis;
    CoinWarmStartBasis * basis =
        dynamic_cast<CoinWarmStartBasis *>(solver->getWarmStart()) ;
    if (basis) {
        saveBasis = * basis;
        delete basis;
    }
    // First just fix all integers as close to zero as possible
    OsiSolverInterface * newSolver = cloneBut(7); // wassolver->clone();
    for (i = 0; i < numberIntegers; i++) {
        int iColumn = integerVariable[i];
        double lower = colLower[iColumn];
        double upper = colUpper[iColumn];
        double value;
        if (lower > 0.0)
            value = lower;
        else if (upper < 0.0)
            value = upper;
        else
            value = 0.0;
        newSolver->setColLower(iColumn, value);
        newSolver->setColUpper(iColumn, value);
    }
    newSolver->initialSolve();
    if (newSolver->isProvenOptimal()) {
        double solValue = newSolver->getObjValue() * direction ;
        if (solValue < cutoff) {
            // we have a solution
            solutionFound = true;
            solutionValue = solValue;
            memcpy(betterSolution, newSolver->getColSolution(),
                   numberColumns*sizeof(double));
            COIN_DETAIL_PRINT(printf("Naive fixing close to zero gave solution of %g\n", solutionValue));
            cutoff = solValue - model_->getCutoffIncrement();
        }
    }
    // Now fix all integers as close to zero if zero or large cost
    int nFix = 0;
    for (i = 0; i < numberIntegers; i++) {
        int iColumn = integerVariable[i];
        double lower = colLower[iColumn];
        double upper = colUpper[iColumn];
        double value;
        if (fabs(objective[i]) > 0.0 && fabs(objective[i]) < large_) {
            nFix++;
            if (lower > 0.0)
                value = lower;
            else if (upper < 0.0)
                value = upper;
            else
                value = 0.0;
            newSolver->setColLower(iColumn, value);
            newSolver->setColUpper(iColumn, value);
        } else {
            // set back to original
            newSolver->setColLower(iColumn, lower);
            newSolver->setColUpper(iColumn, upper);
        }
    }
    const double * solution = solver->getColSolution();
    if (nFix) {
        newSolver->setWarmStart(&saveBasis);
        newSolver->setColSolution(solution);
        newSolver->initialSolve();
        if (newSolver->isProvenOptimal()) {
            double solValue = newSolver->getObjValue() * direction ;
            if (solValue < cutoff) {
                // try branch and bound
//.........这里部分代码省略.........

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

//-------------------------------------------------------------------
// Generate cuts
//-------------------------------------------------------------------
void CglAllDifferent::generateCuts(const OsiSolverInterface & si, OsiCuts & cs,
			      const CglTreeInfo ) const
{
#ifndef NDEBUG
  int nCols=si.getNumCols();
#endif
  int i;
  const double * lower = si.getColLower();
  const double * upper = si.getColUpper();
#ifdef CGL_DEBUG
  const OsiRowCutDebugger * debugger = si.getRowCutDebugger();
  if (debugger&&debugger->onOptimalPath(si)) {
    printf("On optimal path %d\n",nPath);
    nPath++;
    int nCols=si.getNumCols();
    const double * solution = si.getColSolution();
    const double * optimal = debugger->optimalSolution();
    const double * objective = si.getObjCoefficients();
    double objval1=0.0,objval2=0.0;
    for (i=0;i<nCols;i++) {
#if CGL_DEBUG>1
      printf("%d %g %g %g %g\n",i,lower[i],solution[i],upper[i],optimal[i]);
#endif
      objval1 += solution[i]*objective[i];
      objval2 += optimal[i]*objective[i];
      assert(optimal[i]>=lower[i]&&optimal[i]<=upper[i]);
    }
    printf("current obj %g, integer %g\n",objval1,objval2);
  }
#endif
  int * lo = new int[numberDifferent_];
  int * up = new int[numberDifferent_];
  for (i=0;i<numberDifferent_;i++) {
    int iColumn = originalWhich_[i];
    assert (iColumn<nCols);
    lo[i]  = static_cast<int> (lower[iColumn]);
    assert (floor(lower[iColumn]+0.5)==lower[iColumn]);
    up[i]  = static_cast<int> (upper[iColumn]);
    assert (floor(upper[iColumn]+0.5)==upper[iColumn]);
    assert (up[i]>=lo[i]);
  }
  // We are going to assume we can just have one big 2d array!
  // Could save by going to bits
  // also could skip sets where all are fixed
  // could do some of above by separate first pass
  // once a variable fixed - can take out of list
  // so need to redo complete stuff (including temp which_) every big pass
  int offset = COIN_INT_MAX;
  int maxValue = -COIN_INT_MAX;
  int numberLook=0;
  // copies
  //int * which = new int [numberTotal];
  //int * start = new int [numberSets_+1];
  for (i=0;i<numberSets_;i++) {
    for (int j=start_[i];j<start_[i+1];j++) {
      int k=which_[j];
      offset = CoinMin(offset,lo[k]);
      maxValue = CoinMax(maxValue,up[k]);
    }
    numberLook++;
    int gap = maxValue-offset+1;
    double size = static_cast<double> (gap) * numberDifferent_;
    if (size>1.0e7) {
      if (logLevel_)
        printf("Only looking at %d sets\n",numberLook);
      break;
    }
  }
  // Which sets a variable is in
  int * back = new int [start_[numberSets_]];
  int * backStart = new int[numberDifferent_+1];
  memset(backStart,0,(numberDifferent_+1)*sizeof(int));
  int numberTotal = start_[numberLook];
  for (i=0;i<numberTotal;i++) {
    int k=which_[i];
    // note +1
    backStart[k+1]++;
  }
  int n=0;
  for (i=0;i<numberDifferent_;i++) {
    int nThis = backStart[i+1];
    backStart[i+1]=n;
    n+= nThis;
  }
  // at end all backStart correct!
  for (i=0;i<numberLook;i++) {
    for (int j=start_[i];j<start_[i+1];j++) {
      int k=which_[j];
      // note +1
      int iPut = backStart[k+1];
      back[iPut]=i;
      backStart[k+1]=iPut+1;
    }
  }
  // value is possible for variable k if possible[k*gap+value] is nonzero
  int gap = maxValue-offset+1;
  char * possible = new char[gap*numberDifferent_];
//.........这里部分代码省略.........