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

void
CglClique::generateCuts(const OsiSolverInterface& si, OsiCuts & cs,
			const CglTreeInfo info)
{
   int i;
   bool has_petol_set = petol != -1.0;

   if (! has_petol_set)
      si.getDblParam(OsiPrimalTolerance, petol);
   int numberOriginalRows = si.getNumRows();
   if (info.inTree&&justOriginalRows_)
     numberOriginalRows = info.formulation_rows;
   int numberRowCutsBefore = cs.sizeRowCuts();
   // First select which rows/columns we are interested in.
   if (!setPacking_) {
      selectFractionalBinaries(si);
      if (!sp_orig_row_ind) {
	 selectRowCliques(si,numberOriginalRows);
      }
   } else {
      selectFractionals(si);
      delete[] sp_orig_row_ind;
      sp_numrows = numberOriginalRows;
      //sp_numcols = si.getNumCols();
      sp_orig_row_ind = new int[sp_numrows];
      for (i = 0; i < sp_numrows; ++i)
	 sp_orig_row_ind[i] = i;
   }
   // Just original rows
   if (justOriginalRows_&&info.inTree) 
     sp_numrows = CoinMin(info.formulation_rows,sp_numrows);
     

   createSetPackingSubMatrix(si);
   fgraph.edgenum = createNodeNode();
   createFractionalGraph();

   cl_indices = new int[sp_numcols];
   cl_del_indices = new int[sp_numcols];

   if (do_row_clique)
      find_rcl(cs);
   if (do_star_clique)
      find_scl(cs);
   if (!info.inTree&&((info.options&4)==4||((info.options&8)&&!info.pass))) {
     int numberRowCutsAfter = cs.sizeRowCuts();
     for (int i=numberRowCutsBefore;i<numberRowCutsAfter;i++)
       cs.rowCutPtr(i)->setGloballyValid();
   }

   delete[] cl_indices;     cl_indices = 0;
   delete[] cl_del_indices; cl_del_indices = 0;

   deleteFractionalGraph();
   delete[] node_node;      node_node = 0;
   deleteSetPackingSubMatrix();

   if (! has_petol_set)
      petol = -1;
}

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

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

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

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

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

/** Standard cut generation methods. */
void
OaDecompositionBase::generateCuts(const OsiSolverInterface &si,  OsiCuts & cs,
    const CglTreeInfo info) {
  if (nlp_ == NULL) {
    throw CoinError("Error in cut generator for outer approximation no NLP ipopt assigned", "generateCuts", "OaDecompositionBase");
  }

  // babInfo is used to communicate with the b-and-b solver (Cbc or Bcp).
  BabInfo * babInfo = dynamic_cast<BabInfo *> (si.getAuxiliaryInfo());
  assert(babInfo);
  assert(babInfo->babPtr());
  numSols_ = babInfo->babPtr()->model().getSolutionCount ();
  CglTreeInfo info_copy = info;
  const CbcNode * node = babInfo->babPtr()->model().currentNode();
  info_copy.level = (node == NULL) ? 0 : babInfo->babPtr()->model().currentNode()->depth();
  if(babInfo->hasSolution()) numSols_ ++;
  if (babInfo)
    if (!babInfo->mipFeasible())
      return;

  //Get the continuous solution
  const double *colsol = si.getColSolution();


  vector<double> savedColLower(nlp_->getNumCols());
  CoinCopyN(nlp_->getColLower(), nlp_->getNumCols(), savedColLower());
  vector<double> savedColUpper(nlp_->getNumCols());
  CoinCopyN(nlp_->getColUpper(), nlp_->getNumCols(), savedColUpper());


  OsiBranchingInformation brInfo(nlp_, false);
  brInfo.solution_ = colsol;
  //Check integer infeasibility
  bool isInteger = integerFeasible(*nlp_, brInfo, parameters_.cbcIntegerTolerance_,
                              objects_, nObjects_);


  //Check nodeNumber if it did not change scan savedCuts_ if one is violated force it and exit
  int nodeNumber = babInfo->babPtr()->model().getNodeCount();
  if(nodeNumber == currentNodeNumber_){
#ifdef OA_DEBUG
    printf("OA decomposition recalled from the same node!\n");
#endif
    int numCuts = savedCuts_.sizeRowCuts();
    for(int i = 0 ; i < numCuts ; i++){
       //Check if cuts off solution
       if(savedCuts_.rowCut(i).violated(colsol) > 0.){
#ifdef OA_DEBUG
         printf("A violated saved cut has been found\n");
#endif
         savedCuts_.rowCut(i).setEffectiveness(9.99e99);
         cs.insert(savedCuts_.rowCut(i));
         savedCuts_.eraseRowCut(i);
         return;
         i--; numCuts--;
       }
    }
  }
  else {
    currentNodeNumber_ = nodeNumber;
    savedCuts_.dumpCuts();
  } 
         
  if (!isInteger) {
    if (!doLocalSearch(babInfo))//create sub mip solver.
      return;
  }

  //get the current cutoff
  double cutoff;
  si.getDblParam(OsiDualObjectiveLimit, cutoff);

  // Save solvers state if needed

  solverManip * lpManip = NULL;
  if (lp_ != NULL) {
      assert(lp_ == &si);
      lpManip = new solverManip(lp_, true, leaveSiUnchanged_, true, true);
  }
  else {
    lpManip = new solverManip(si);
  }
  lpManip->setObjects(objects_, nObjects_);

  double milpBound = performOa(cs, *lpManip, babInfo, cutoff, info_copy);

  if(babInfo->hasSolution()){
     babInfo->babPtr()->model().setSolutionCount (numSols_ - 1);
  }

  //Transmit the bound found by the milp
  {
    if (milpBound>-1e100)
    {
      // Also store into solver
      if (babInfo)
        babInfo->setMipBound(milpBound);
    }
  }  //Clean everything :
//.........这里部分代码省略.........

/*
  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 */
int
CbcHeuristicVND::solution(double & solutionValue,
                          double * betterSolution)
{
    numCouldRun_++;
    int returnCode = 0;
    const double * bestSolution = model_->bestSolution();
    if (!bestSolution)
        return 0; // No solution found yet
#ifdef HEURISTIC_INFORM
    printf("Entering heuristic %s - nRuns %d numCould %d when %d\n",
	   heuristicName(),numRuns_,numCouldRun_,when_);
#endif
    if (numberSolutions_ < model_->getSolutionCount()) {
        // new solution - add info
        numberSolutions_ = model_->getSolutionCount();

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

        int i;
        for (i = 0; i < numberIntegers; i++) {
            int iColumn = integerVariable[i];
            const OsiObject * object = model_->object(i);
            // get original bounds
            double originalLower;
            double originalUpper;
            getIntegerInformation( object, originalLower, originalUpper);
            double value = bestSolution[iColumn];
            if (value < originalLower) {
                value = originalLower;
            } else if (value > originalUpper) {
                value = originalUpper;
            }
        }
    }
    int numberNodes = model_->getNodeCount();
    if (howOften_ == 100) {
        if (numberNodes < lastNode_ + 12)
            return 0;
        // Do at 50 and 100
        if ((numberNodes > 40 && numberNodes <= 50) || (numberNodes > 90 && numberNodes < 100))
            numberNodes = howOften_;
    }
    if ((numberNodes % howOften_) == 0 && (model_->getCurrentPassNumber() <= 1 ||
                                           model_->getCurrentPassNumber() == 999999)) {
        lastNode_ = model_->getNodeCount();
        OsiSolverInterface * solver = model_->solver();

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

        const double * currentSolution = solver->getColSolution();
        OsiSolverInterface * newSolver = cloneBut(3); // was model_->continuousSolver()->clone();
        //const double * colLower = newSolver->getColLower();
        //const double * colUpper = newSolver->getColUpper();

        double primalTolerance;
        solver->getDblParam(OsiPrimalTolerance, primalTolerance);

        // Sort on distance
        double * distance = new double [numberIntegers];
        int * which = new int [numberIntegers];

        int i;
        int nFix = 0;
        double tolerance = 10.0 * primalTolerance;
        for (i = 0; i < numberIntegers; i++) {
            int iColumn = integerVariable[i];
            const OsiObject * object = model_->object(i);
            // get original bounds
            double originalLower;
            double originalUpper;
            getIntegerInformation( object, originalLower, originalUpper);
            double valueInt = bestSolution[iColumn];
            if (valueInt < originalLower) {
                valueInt = originalLower;
            } else if (valueInt > originalUpper) {
                valueInt = originalUpper;
            }
            baseSolution_[iColumn] = currentSolution[iColumn];
            distance[i] = fabs(currentSolution[iColumn] - valueInt);
            which[i] = i;
            if (fabs(currentSolution[iColumn] - valueInt) < tolerance)
                nFix++;
        }
        CoinSort_2(distance, distance + numberIntegers, which);
        nDifferent_ = numberIntegers - nFix;
        stepSize_ = nDifferent_ / 10;
        k_ = stepSize_;
        //nFix = numberIntegers-stepSize_;
        for (i = 0; i < nFix; i++) {
            int j = which[i];
            int iColumn = integerVariable[j];
            const OsiObject * object = model_->object(i);
            // get original bounds
//.........这里部分代码省略.........

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

/*
  Randomized Rounding Heuristic
  Returns 1 if solution, 0 if not
*/
int
CbcHeuristicRandRound::solution(double & solutionValue,
                                double * betterSolution)
{
    // rlh: Todo: Memory Cleanup

    //  std::cout << "Entering the Randomized Rounding Heuristic" << std::endl;

    setWhen(1);  // setWhen(1) didn't have the effect I expected (e.g., run once).

    // Run only once.
    //
    //    See if at root node
    bool atRoot = model_->getNodeCount() == 0;
    int passNumber = model_->getCurrentPassNumber();
    //    Just do once
    if (!atRoot || passNumber > 1) {
        // std::cout << "Leaving the Randomized Rounding Heuristic" << std::endl;
        return 0;
    }

    std::cout << "Entering the Randomized Rounding Heuristic" << std::endl;
    typedef struct {
        int numberSolutions;
        int maximumSolutions;
        int numberColumns;
        double ** solution;
        int * numberUnsatisfied;
    } clpSolution;

    double start = CoinCpuTime();
    numCouldRun_++; //
#ifdef HEURISTIC_INFORM
    printf("Entering heuristic %s - nRuns %d numCould %d when %d\n",
	   heuristicName(),numRuns_,numCouldRun_,when_);
#endif
    // Todo: Ask JJHF what "number of times
    // the heuristic could run" means.

    OsiSolverInterface * solver = model_->solver()->clone();
    double primalTolerance ;
    solver->getDblParam(OsiPrimalTolerance, primalTolerance) ;
    OsiClpSolverInterface * clpSolver = dynamic_cast<OsiClpSolverInterface *> (solver);
    assert (clpSolver);
    ClpSimplex * simplex = clpSolver->getModelPtr();

    // Initialize the structure holding the solutions for the Simplex iterations
    clpSolution solutions;
    // Set typeStruct field of ClpTrustedData struct to 1 to indicate
    // desired behavior for  RandRound heuristic (which is what?)
    ClpTrustedData trustedSolutions;
    trustedSolutions.typeStruct = 1;
    trustedSolutions.data = &solutions;
    solutions.numberSolutions = 0;
    solutions.maximumSolutions = 0;
    solutions.numberColumns = simplex->numberColumns();
    solutions.solution = NULL;
    solutions.numberUnsatisfied = NULL;
    simplex->setTrustedUserPointer(&trustedSolutions);

    // Solve from all slack to get some points
    simplex->allSlackBasis();

    // Calling primal() invalidates pointers to some rim vectors,
    // like...row sense (!)
    simplex->primal();

    // 1. Okay - so a workaround would be to copy the data I want BEFORE
    // calling primal.
    // 2. Another approach is to ask the simplex solvers NOT to mess up my
    // rims.
    // 3. See freeCachedResults() for what is getting
    // deleted. Everything else points into the structure.
    // ...or use collower and colupper rather than rowsense.
    // ..store address of where one of these

    // Store the basic problem information
    // -Get the number of columns, rows and rhs vector
    int numCols = clpSolver->getNumCols();
    int numRows = clpSolver->getNumRows();

    // Find the integer variables (use columnType(?))
    // One if not continuous, that is binary or general integer)
    // columnType() = 0 continuous
    //              = 1 binary
    //              = 2 general integer
    bool * varClassInt = new bool[numCols];
    const char* columnType = clpSolver->columnType();
    int numGenInt = 0;
    for (int i = 0; i < numCols; i++) {
        if (clpSolver->isContinuous(i))
            varClassInt[i] = 0;
        else
            varClassInt[i] = 1;
        if (columnType[i] == 2) numGenInt++;
    }
//.........这里部分代码省略.........

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

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