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

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

//.........这里部分代码省略.........
                        std::cout << "x[15] = " << saveSolution[15]
                                  << std::endl;
                    }
#endif

                    intObject = NULL;
                }
                else {
                    // TODO: currently all are integer objects.
#ifdef BLIS_DEBU
                    assert(0);
#endif
                }

            }
        }

        if (numInfs) {
#ifdef BLIS_DEBUG_MORE
            std::cout << "REL: numInfs = " << numInfs
                      << std::endl;
#endif
            break;
        }
        else if (pass == 0) {
            // The first pass and is IP feasible.

#ifdef BLIS_DEBUG
            std::cout << "REL: given a feasible sol" << std::endl;
#endif

            roundAgain = false;
            CoinWarmStartBasis * ws =
                dynamic_cast<CoinWarmStartBasis*>(solver->getWarmStart());
            if (!ws) break;

            // Force solution values within bounds
            for (i = 0; i < numCols; ++i) {
                lpX = saveSolution[i];
                if (lpX < lower[i]) {
                    saveSolution[i] = lower[i];
                    roundAgain = true;
                    ws->setStructStatus(i, CoinWarmStartBasis::atLowerBound);
                }
                else if (lpX > upper[i]) {
                    saveSolution[i] = upper[i];
                    roundAgain = true;
                    ws->setStructStatus(i, CoinWarmStartBasis::atUpperBound);
                }
            }

            if (roundAgain) {
                // Need resolve and do the second round selection.
                solver->setWarmStart(ws);
                delete ws;

                // Resolve.
                solver->resolve();

                if (!solver->isProvenOptimal()) {
                    // Become infeasible, can do nothing.
                    bStatus = -2;
                    goto TERM_CREATE;
                }
                else {
                    // Save new lp solution.

//.........这里部分代码省略.........
                                  << std::endl;
                    }
#endif

                    intObject = NULL;
                }
                else {
                    // TODO: currently all are integer objects.
#ifdef BLIS_DEBU
                    assert(0);
#endif
                }
                
            }
        }
            
        if (numInfs) {
#if 0
            std::cout << "PSEUDO: numInfs = " << numInfs
                      << std::endl;
#endif
            break;
        }
        else if (pass == 0) {
            // The first pass and is IP feasible.
            
#if 1
            std::cout << "ERROR: PSEUDO: given a integer feasible sol, no fraction variable" << std::endl;
            assert(0);
#endif      
            
            roundAgain = false;
            CoinWarmStartBasis * ws = 
                dynamic_cast<CoinWarmStartBasis*>(solver->getWarmStart());
            if (!ws) break;
            
            // Force solution values within bounds
            for (i = 0; i < numCols; ++i) {
                lpX = saveSolution[i];
                if (lpX < lower[i]) {
                    saveSolution[i] = lower[i];
                    roundAgain = true;
                    ws->setStructStatus(i, CoinWarmStartBasis::atLowerBound);
                } 
                else if (lpX > upper[i]) {
                    saveSolution[i] = upper[i];
                    roundAgain = true;
                    ws->setStructStatus(i, CoinWarmStartBasis::atUpperBound);
                } 
            }
            
            if (roundAgain) {
                // Need resolve and do the second round selection.
                solver->setWarmStart(ws);
                delete ws;
                
                // Resolve.
                solver->resolve();
		
                if (!solver->isProvenOptimal()) {
                    // Become infeasible, can do nothing. 
                    bStatus = -2;
                    goto TERM_CREATE;
                }
                else {
                    // Save new lp solution.

//.........这里部分代码省略.........
  //   delete column for v_0
  // endFor
  // clean up memory
  // return 0;

  int * nVectorIndices = new int[n];
  CoinIotaN(nVectorIndices, n, 0);

  bool haveWarmStart = false;
  bool equalObj1, equalObj2;
  CoinRelFltEq eq;

  double v_0Elements[2] = {-1,1};
  double u_0Elements[1] = {1};

  CoinWarmStart * warmStart = 0;

  double * ustar = new double[m];
  CoinFillN(ustar, m, 0.0);

  double* alpha = new double[n];
  CoinFillN(alpha, n, 0.0);

  for (j=0;j<n;j++){
    if (!si.isBinary(j)) continue; // Better to ask coneSi? No! 
                                   // coneSi has no binInfo.
    equalObj1=eq(x[j],0);
    equalObj2=eq(x[j],1);
    if (equalObj1 || equalObj2) continue;
    // IMPROVEME: if (haveWarmStart) check if j attractive;

    // AskLL:wanted to declare u_0 and v_0 packedVec outside loop
    // and setIndices, but didn't see a method to do that(?)
    // (Could "insert". Seems inefficient)
    int v_0Indices[2]={j,nPlus1};
    int u_0Indices[1]={j};
    // 
    CoinPackedVector  v_0(2,v_0Indices,v_0Elements,false);
    CoinPackedVector  u_0(1,u_0Indices,u_0Elements,false);

#if CGL_DEBUG
    const CoinPackedMatrix *see1 = coneSi->getMatrixByRow();
#endif

    coneSi->addCol(v_0,-solverINFINITY,solverINFINITY,0);
    coneSi->addCol(u_0,-solverINFINITY,solverINFINITY,x[j]);
    if(haveWarmStart) {
      coneSi->setWarmStart(warmStart);
      coneSi->resolve();
    }
    else {

#if CGL_DEBUG
      const CoinPackedMatrix *see2 = coneSi->getMatrixByRow();
#endif

      coneSi->initialSolve();
    }
    if(coneSi->isProvenOptimal()){
      warmStart = coneSi->getWarmStart();
      haveWarmStart=true;
      const double * wstar = coneSi->getColSolution();
      CoinDisjointCopyN(wstar, m, ustar);
      Atilde->transposeTimes(ustar,alpha);
      alpha[j]+=wstar[BNumCols-1]; 
      
#if debug
      int p;
      double sum;
      for(p=0;p<n;p++)sum+=alpha[p]*x[p];
      if (sum<=beta_){
	throw CoinError("Cut not violated",
			"cutGeneration",
			"CglLiftAndProject");
      }
#endif

      // add <alpha^T,x> >= beta_ to cutset
      OsiRowCut rc;
      rc.setRow(n,nVectorIndices,alpha);
      rc.setLb(beta_);
      rc.setUb(solverINFINITY);
      cs.insert(rc);
    }
    // delete col for u_o and v_0
    coneSi->deleteCols(2,delCols);

    // clean up memory
  }
  // clean up
  delete [] alpha;
  delete [] ustar;
  delete [] nVectorIndices;
  // BMatrix, BColLowers,BColUppers, BObjective, BRowLowers, BRowUppers
  // are all freed by OsiSolverInterface destructor (?)
  delete [] BLengths;
  delete [] BStarts;
  delete [] BIndices;
  delete [] BElements;
}

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

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

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

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

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

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

//.........这里部分代码省略.........
     sym_set_int_param(env, "max_sp_size", 100);
     sym_set_int_param(env, "do_reduced_cost_fixing", FALSE);
     if (whichCutsLL == 0){
	sym_set_int_param(env, "generate_cgl_cuts", FALSE);
     }else{
	sym_set_int_param(env, "generate_cgl_gomory_cuts", GENERATE_DEFAULT);
     }
     if (whichCutsLL == 1){
	sym_set_int_param(env, "generate_cgl_knapsack_cuts", 
			  DO_NOT_GENERATE);
	sym_set_int_param(env, "generate_cgl_probing_cuts", 
			  DO_NOT_GENERATE);
	sym_set_int_param(env, "generate_cgl_clique_cuts", 
			  DO_NOT_GENERATE);
	sym_set_int_param(env, "generate_cgl_twomir_cuts", 
			  DO_NOT_GENERATE);
	sym_set_int_param(env, "generate_cgl_flowcover_cuts", 
			  DO_NOT_GENERATE);
     }
  }else if (feasCheckSolver == "CPLEX"){
#ifdef USE_CPLEX
     lSolver->setHintParam(OsiDoReducePrint);
     lSolver->messageHandler()->setLogLevel(0);
     CPXENVptr cpxEnv = 
	dynamic_cast<OsiCpxSolverInterface*>(lSolver)->getEnvironmentPtr();
     assert(cpxEnv);
     CPXsetintparam(cpxEnv, CPX_PARAM_SCRIND, CPX_OFF);
     CPXsetintparam(cpxEnv, CPX_PARAM_THREADS, maxThreadsLL);
#endif
  }
  
  if (warmStartLL && feasCheckSolver == "SYMPHONY"){
     lSolver->resolve();
     setWarmStart(lSolver->getWarmStart());
  }else{
     lSolver->branchAndBound();
  }

  const double * sol = model_->solver()->getColSolution();
  double objVal(lSolver->getObjValue() * model_->getLowerObjSense());
  
  MibSTreeNode * node = static_cast<MibSTreeNode *>(model_->activeNode_);
  MibSTreeNode * parent = 
    static_cast<MibSTreeNode *>(model_->activeNode_->getParent());

  if((!node->isBoundSet()) 
     && (node->getIndex() != 0)){
    double parentBound = parent->getLowerUB();
    node->setLowerUB(parentBound);
    node->setIsBoundSet(true);
  }
  
  if(objVal > node->getLowerUB()){
    
    node->setLowerUB(objVal);
    node->setIsBoundSet(true);
    
  }

  double etol(model_->etol_);
  double lowerObj = getLowerObj(sol, model_->getLowerObjSense());  

  int lN(model_->lowerDim_); // lower-level dimension
  int uN(model_->upperDim_); // lower-level dimension
  if(!optLowerSolution_)
    optLowerSolution_ = new double[lN];

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

//.........这里部分代码省略.........
        }
/*
  End of loop to try increment/decrement of integer variables.

  newSolutionValue does not necessarily match the current newSolution, and
  bestChange simply reflects the best single change. Still, that's sufficient
  to indicate that there's been at least one change. Check that we really do
  have a valid solution.
*/
        if (totalChange + newSolutionValue < solutionValue) {
            // paranoid check
            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];
                    }
                }
            }
            int numberBad = 0;
            double sumBad = 0.0;
            // check was approximately feasible
            for (i = 0; i < numberRows; i++) {
                if (rowActivity[i] < rowLower[i]) {
                    sumBad += rowLower[i] - rowActivity[i];
                    if (rowActivity[i] < rowLower[i] - 10.0*primalTolerance)
                        numberBad++;
                } else if (rowActivity[i] > rowUpper[i]) {
                    sumBad += rowUpper[i] - rowActivity[i];
                    if (rowActivity[i] > rowUpper[i] + 10.0*primalTolerance)
                        numberBad++;
                }
            }
            if (!numberBad) {
                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 = newSolution[iColumn];
                    // if away from lower bound mark that fact
                    if (value > originalLower) {
                        used_[iColumn] = numberSolutions_;
                    }
                }
/*
  Copy the solution to the array returned to the client. Grab a basis from
  the solver (which, if it exists, is almost certainly infeasible, but it
  should be ok for a dual start). The value returned as solutionValue is
  conservative because of handling of newSolutionValue and bestChange, as
  described above.
*/
                // new solution
                memcpy(betterSolution, newSolution, numberColumns*sizeof(double));
                CoinWarmStartBasis * basis =
                    dynamic_cast<CoinWarmStartBasis *>(solver->getWarmStart()) ;
                if (basis) {
                    model_->setBestSolutionBasis(* basis);
                    delete basis;
                }
                returnCode = 1;
                solutionValue = newSolutionValue + bestChange;
            } else {
                // bad solution - should not happen so debug if see message
                COIN_DETAIL_PRINT(printf("Local search got bad solution with %d infeasibilities summing to %g\n",
					 numberBad, sumBad));
            }
        }
    }
/*
  We're done. Clean up.
*/
    delete [] newSolution;
    delete [] rowActivity;
    delete [] way;
    delete [] cost;
    delete [] save;
    delete [] mark;
/*
  Do we want to try swapping values between solutions?
  swap_ is set elsewhere; it's not adjusted during heuristic execution.

  Again, redundant test. We shouldn't be here if numberSolutions_ = 1.
*/
    if (numberSolutions_ > 1 && (swap%10) == 1) {
        // try merge
        int returnCode2 = solutionFix( solutionValue, betterSolution, NULL);
        if (returnCode2)
            returnCode = 1;
    }
    return returnCode;
}