C++ OsiRowCut::violated方法代码示例

/**Clean cut 2, different algorithm. First check the dynamic of the cut if < maxRatio scale to a biggest coef of 1
   otherwise scale it so that biggest coeff is 1 and try removing tinys ( < 1/maxRatio) either succeed or fail */
int
Validator::cleanCut2(OsiRowCut & aCut, const double * solCut, const OsiSolverInterface &si, const CglParam &/* par */,
                     const double * origColLower, const double * origColUpper) const
{
    /** Compute fill-in in si */
    int numcols = si.getNumCols();
    // int numrows = si.getNumRows();
    const double * colLower = (origColLower) ? origColLower : si.getColLower();
    const double * colUpper = (origColUpper) ? origColUpper : si.getColUpper();

    int maxNnz = static_cast<int> ( maxFillIn_ * static_cast<double> (numcols));

    double rhs = aCut.lb();
    assert (aCut.ub()> 1e50);

    CoinPackedVector *vec = const_cast<CoinPackedVector *>(&aCut.row());
    //  vec->sortIncrIndex();

    int * indices = vec->getIndices();
    double * elems = vec->getElements();
    int n = vec->getNumElements();
    if (n==0)
    {
        numRejected_[EmptyCut]++;
        return EmptyCut;
    }
    /** First compute violation if it is too small exit */
    double violation = aCut.violated(solCut);
    if (violation < minViolation_)
        return 1;

    /** Now relax get dynamic and remove tiny elements */
    int offset = 0;
    rhs -= 1e-10;
    double smallest = fabs(rhs);
    double biggest = smallest;
    double veryTiny = 1e-20;
    for (int i = 0 ; i < n ; i++)
    {
        double val = fabs(elems[i]);
        if (val > veryTiny)   //tiny should be very very small
        {
            smallest = std::min(val,smallest);
            biggest = std::max (val,biggest);
        }
    }

    if (biggest > 1e9)
    {
#ifdef DEBUG
        std::cout<<"Whaooo "<<biggest/smallest<<std::endl;
#endif
        numRejected_[BigDynamic]++;
        return BigDynamic;
    }

    //rescale the cut so that biggest is 1e1.
    double toBeBiggest = rhsScale_;
    rhs *= (toBeBiggest / biggest);
    toBeBiggest /= biggest;
    for (int i = 0 ; i < n ; i++)
    {
        elems[i] *= toBeBiggest;
    }


    if (biggest > maxRatio_ * smallest)   //we have to remove some small coefficients
    {
        double myTiny = biggest * toBeBiggest / maxRatio_;
        veryTiny *= toBeBiggest ;
        for (int i = 0 ; i < n ; i++)
        {
            double val = fabs(elems[i]);
            if (val < myTiny)
            {
                if (val< veryTiny)
                {
                    offset++;
                    continue;
                }
                int & iCol = indices[i];
                if (elems[i]>0. && colUpper[iCol] < 1000.)
                {
                    offset++;
                    rhs -= elems[i] * colUpper[iCol];
                    elems[i]=0;
                }
                else if (elems[i]<0. && colLower[iCol] > -1000.)
                {
                    offset++;
                    rhs -= elems[i] * colLower[iCol];
                    elems[i]=0.;
                }
                else
                {
                    numRejected_[SmallCoefficient]++;
                    return SmallCoefficient;
                }
//.........这里部分代码省略.........

//.........这里部分代码省略.........
			DSPdebugMessage("Added pool cut\n");
			/** free memory */
			SCIPfreeMemoryArray(scip, &vals);
			return SCIP_OKAY;
		}
	}
#endif

	/** generate Benders cuts */
	assert(tss_);
	tss_->generateCuts(nvars_, vals, &cs);

	/** If found Benders cuts */
	for (int i = 0; i < cs.sizeCuts(); ++i)
	{
		/** get cut pointer */
		OsiRowCut * rc = cs.rowCutPtr(i);
		if (!rc) continue;

		const CoinPackedVector cutrow = rc->row();
		if (cutrow.getNumElements() == 0) continue;

		/** is optimality cut? */
		bool isOptimalityCut = false;
		for (int j = nvars_ - naux_; j < nvars_; ++j)
		{
			if (cutrow.getMaxIndex() == j)
			{
				isOptimalityCut = true;
				break;
			}
		}

		double efficacy = rc->violated(vals) / cutrow.twoNorm();
		SCIP_Bool isEfficacious = efficacy > 1.e-6;

#define KK_TEST
#ifdef KK_TEST
		if (SCIPgetStage(scip) == SCIP_STAGE_INITSOLVE ||
			SCIPgetStage(scip) == SCIP_STAGE_SOLVING)
		{
			/** create empty row */
			SCIP_ROW * row = NULL;
			SCIP_CALL(SCIPcreateEmptyRowCons(scip, &row, conshdlr, "benders", rc->lb(), SCIPinfinity(scip),
					FALSE, /**< is row local? */
					FALSE, /**< is row modifiable? */
					FALSE  /**< is row removable? can this be TRUE? */));

			/** cache the row extension and only flush them if the cut gets added */
			SCIP_CALL(SCIPcacheRowExtensions(scip, row));

			/** collect all non-zero coefficients */
			for (int j = 0; j < cutrow.getNumElements(); ++j)
				SCIP_CALL(SCIPaddVarToRow(scip, row, vars_[cutrow.getIndices()[j]], cutrow.getElements()[j]));

			DSPdebugMessage("found Benders (%s) cut: act=%f, lhs=%f, norm=%f, eff=%f, min=%f, max=%f (range=%f)\n",
				isOptimalityCut ? "opti" : "feas",
				SCIPgetRowLPActivity(scip, row), SCIProwGetLhs(row), SCIProwGetNorm(row),
				SCIPgetCutEfficacy(scip, sol, row),
				SCIPgetRowMinCoef(scip, row), SCIPgetRowMaxCoef(scip, row),
				SCIPgetRowMaxCoef(scip, row)/SCIPgetRowMinCoef(scip, row));

			/** flush all changes before adding cut */
			SCIP_CALL(SCIPflushRowExtensions(scip, row));

			DSPdebugMessage("efficacy %e isEfficatious %d\n", efficacy, isEfficacious);