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

BcpsBranchObject * 
BlisObjectInt::preferredNewFeasible(BcpsModel *m) const
{
    BlisModel *model = dynamic_cast<BlisModel* >(m);
    OsiSolverInterface * solver = model->solver();
    
    double value = solver->getColSolution()[columnIndex_];
    
    double nearest = floor(value + 0.5);
    double integerTolerance = model->BlisPar()->entry(BlisParams::integerTol);

    assert(fabs(value - nearest) <= integerTolerance);

    double dj = solver->getObjSense()*solver->getReducedCost()[columnIndex_];

    BlisBranchObjectInt * object = NULL;

    if (dj >= 0.0) {
	// Better go down
	if (nearest > originalLower_ + 0.5) {
	    // Has room to go down
	    object = new BlisBranchObjectInt(model,
                                             objectIndex_,
                                             -1,
                                             nearest - 1.0,
                                             nearest - 1.0);
	}
    } 
    else {
	// Better go up
	if (nearest < originalUpper_ - 0.5) {
	    // Has room to go up
	    object = new BlisBranchObjectInt(model, 
                                             objectIndex_, 
                                             -1,
                                             nearest + 1.0,
                                             nearest + 1.0);
	}
    }

    return object;
}

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

    const double* reducedCost = NULL;
    // See if not NLP
    if (model_->solverCharacteristics()->reducedCostsAccurate())
        reducedCost = solver->getReducedCost();

    int iteration = 0;
    while (numberFractionalVariables) {
        iteration++;

        // initialize any data
        initializeData();

        // select a fractional variable to bound
        int bestColumn = -1;
        int bestRound; // -1 rounds down, +1 rounds up
        bool canRound = selectVariableToBranch(solver, newSolution,
                                               bestColumn, bestRound);
        // if the solution is not trivially roundable, we don't try to round;
        // if the solution is trivially roundable, we try to round. However,
        // if the rounded solution is worse than the current incumbent,
        // then we don't round and proceed normally. In this case, the
        // bestColumn will be a trivially roundable variable
        if (canRound) {
            // check if by rounding all fractional variables
            // we get a solution with an objective value
            // better than the current best integer solution
            double delta = 0.0;
            for (int i = 0; i < numberIntegers; i++) {
                int iColumn = integerVariable[i];
                double value = newSolution[iColumn];
                if (fabs(floor(value + 0.5) - value) > integerTolerance) {
                    assert(downLocks_[i] == 0 || upLocks_[i] == 0);
                    double obj = objective[iColumn];
                    if (downLocks_[i] == 0 && upLocks_[i] == 0) {
                        if (direction * obj >= 0.0)
                            delta += (floor(value) - value) * obj;

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

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

//Solver function
int sci_rmps(char *fname) 
{
    //creating a problem pointer using base class of OsiSolverInterface and
    //instantiate the object using derived class of ClpSolverInterface
    OsiSolverInterface* si = new OsiClpSolverInterface();

    // Error management variable
	SciErr sciErr;

	//data declarations
	int *piAddressVarOne = NULL;                 //pointer used to access argument of the function
	char* ptr;                              	 //pointer to point to address of file name
    double* options_;                            //options to set maximum iterations 
	CheckInputArgument(pvApiCtx, 2,2 );          //Check we have exactly two arguments as input or not
	CheckOutputArgument(pvApiCtx, 6, 6);         //Check we have exactly six arguments on output side or not
    //Getting the input arguments from Scilab
    //Getting the MPS file path
	//Reading mps file
	getStringFromScilab(1,&ptr);

 	std::cout<<ptr;
	
    //get options from Scilab
    if(getFixedSizeDoubleMatrixInList(2 , 2 , 1 , 1 , &options_))
	{
		return 1;
	}

    //Read the MPS file
    si->readMps(ptr);

    //setting options for maximum iterations
    si->setIntParam(OsiMaxNumIteration,options_[0]);

    //Solve the problem
    si->initialSolve();
  
    //Quering about the problem
    //get number of variables
    double numVars_;
    numVars_ = si->getNumCols();
  
    //get number of constraint equations
    double numCons_;
    numCons_ = si->getNumRows();
   
    //Output the solution to Scilab
    //get solution for x
    const double* xValue = si->getColSolution();
   
    //get objective value
    double objValue = si->getObjValue();

    //get Status value
    double status;
    if(si->isProvenOptimal())
    	status=0;
    else if(si->isProvenPrimalInfeasible())
    	status=1;
    else if(si->isProvenDualInfeasible())
        status=2;
    else if(si->isIterationLimitReached())
        status=3;
   	else if(si->isAbandoned())
        status=4;
   	else if(si->isPrimalObjectiveLimitReached())
        status=5;
   	else if(si->isDualObjectiveLimitReached())
        status=6;

    //get number of iterations
    double iterations = si->getIterationCount();

    //get reduced cost 
    const double* reducedCost = si->getReducedCost();
   
    //get dual vector
    const double* dual = si->getRowPrice();
  
    returnDoubleMatrixToScilab(1 , 1 , numVars_ , xValue);
    returnDoubleMatrixToScilab(2 , 1 , 1 , &objValue);
    returnDoubleMatrixToScilab(3 , 1 , 1 , &status);
    returnDoubleMatrixToScilab(4 , 1 , 1 , &iterations);
    returnDoubleMatrixToScilab(5 , 1 , numVars_ , reducedCost);
    returnDoubleMatrixToScilab(6 , 1 , numCons_ , dual);
	
	free(xValue);
	free(dual);
	free(reducedCost);
}

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

//.........这里部分代码省略.........
        for (i = 0; i < n; i++) {
            int iColumn = integerVariable[i];
            if (used_[iColumn] <= allow) {
                newSolver->setColUpper(iColumn, colLower[iColumn]);
                nFix2++;
            } else {
                break;
            }
        }
        delete [] which;
        nFix += nFix2;
#ifdef CLP_INVESTIGATE2
        printf("Number fixed increased from %d to %d\n",
               nFix - nFix2, nFix);
#endif
    }
    if (nFix*10 > numberIntegers) {
        returnCode = smallBranchAndBound(newSolver, numberNodes_, newSolution, objectiveValue,
                                         objectiveValue, "CbcHeuristicLocal");
 /*
  -2 is return due to user event, and -1 is overloaded with what look to be
  two contradictory meanings.
*/
       if (returnCode < 0) {
            returnCode = 0; // returned on size
            int numberColumns = newSolver->getNumCols();
            int numberContinuous = numberColumns - numberIntegers;
            if (numberContinuous > 2*numberIntegers &&
                    nFix*10 < numberColumns) {
#define LOCAL_FIX_CONTINUOUS
#ifdef LOCAL_FIX_CONTINUOUS
                //const double * colUpper = newSolver->getColUpper();
                const double * colLower = newSolver->getColLower();
                int nAtLb = 0;
                //double sumDj=0.0;
                const double * dj = newSolver->getReducedCost();
                double direction = newSolver->getObjSense();
                for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
                    if (!newSolver->isInteger(iColumn)) {
                        if (!used_[iColumn]) {
                            //double djValue = dj[iColumn]*direction;
                            nAtLb++;
                            //sumDj += djValue;
                        }
                    }
                }
                if (nAtLb) {
                    // fix some continuous
                    double * sort = new double[nAtLb];
                    int * which = new int [nAtLb];
                    //double threshold = CoinMax((0.01*sumDj)/static_cast<double>(nAtLb),1.0e-6);
                    int nFix2 = 0;
                    for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
                        if (!newSolver->isInteger(iColumn)) {
                            if (!used_[iColumn]) {
                                double djValue = dj[iColumn] * direction;
                                if (djValue > 1.0e-6) {
                                    sort[nFix2] = -djValue;
                                    which[nFix2++] = iColumn;
                                }
                            }
                        }
                    }
                    CoinSort_2(sort, sort + nFix2, which);
                    int divisor = 2;
                    nFix2 = CoinMin(nFix2, (numberColumns - nFix) / divisor);
                    for (int i = 0; i < nFix2; i++) {
                        int iColumn = which[i];
                        newSolver->setColUpper(iColumn, colLower[iColumn]);
                    }
                    delete [] sort;
                    delete [] which;
#ifdef CLP_INVESTIGATE2
                    printf("%d integers have zero value, and %d continuous fixed at lb\n",
                           nFix, nFix2);
#endif
                    returnCode = smallBranchAndBound(newSolver,
                                                     numberNodes_, newSolution,
                                                     objectiveValue,
                                                     objectiveValue, "CbcHeuristicLocal");
                    if (returnCode < 0)
                        returnCode = 0; // returned on size
                }
#endif
            }
        }
    }
/*
  If the result is complete exploration with a solution (3) or proven
  infeasibility (2), we could generate a cut (the AI folks would call it a
  nogood) to prevent us from going down this route in the future.
*/
    if ((returnCode&2) != 0) {
        // could add cut
        returnCode &= ~2;
    }

    delete newSolver;
    return returnCode;
}