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

/** scale the cut passed as argument*/
void scale(OsiRowCut &cut)
{
    double rhs = fabs(cut.lb());
    CoinPackedVector row;
    row.reserve(cut.row().getNumElements());
    for (int i = 0 ; i < cut.row().getNumElements() ; i++)
    {
        row.insert(cut.row().getIndices()[i], cut.row().getElements()[i]/rhs);
    }
    cut.setLb(cut.lb()/rhs);
    cut.setRow(row);
}

/* Returns
   0 unchanged
   1 strengthened
   2 weakened
   3 deleted
*/
int
CbcCutSubsetModifier::modify(const OsiSolverInterface * /*solver*/,
                             OsiRowCut & cut)
{
    int n = cut.row().getNumElements();
    if (!n)
        return 0;
    const int * column = cut.row().getIndices();
    //const double * element = cut.row().getElements();
    int returnCode = 0;
    for (int i = 0; i < n; i++) {
        if (column[i] >= firstOdd_) {
            returnCode = 3;
            break;
        }
    }
#ifdef COIN_DETAIL
    if (!returnCode) {
        const double * element = cut.row().getElements();
        printf("%g <= ", cut.lb());
        for (int i = 0; i < n; i++) {
            printf("%g*x%d ", element[i], column[i]);
        }
        printf("<= %g\n", cut.ub());
    }
#endif
    //return 3;
    return returnCode;
}

//----------------------------------------------------------------
// == operator 
//-------------------------------------------------------------------
bool
OsiRowCut::operator==(const OsiRowCut& rhs) const
{
  if ( this->OsiCut::operator!=(rhs) ) return false;
  if ( row() != rhs.row() ) return false;
  if ( lb() != rhs.lb() ) return false;
  if ( ub() != rhs.ub() ) return false;
  return true;
}

void
OsiTestSolverInterface::applyRowCut(const OsiRowCut& rc)
{
   const int rownum = getNumRows();
   const double lb = rc.lb();
   const double ub = rc.ub();
   rowRimResize_(rownum + 1);
   rowprice_[rownum] = 0.0;
   rowlower_[rownum] = lb;
   rowupper_[rownum] = ub;
   convertBoundToSense(lb, ub,
		       rowsense_[rownum], rhs_[rownum], rowrange_[rownum]);

   updateRowMatrix_();
   rowMatrix_.appendRow(rc.row());
   colMatrixCurrent_ = false;
}

/* Insert a row cut unless it is a duplicate (CoinRelFltEq)*/
void 
OsiCuts::insertIfNotDuplicate( OsiRowCut & rc , CoinRelFltEq treatAsSame)
{
  double newLb = rc.lb();
  double newUb = rc.ub();
  CoinPackedVector vector = rc.row();
  int numberElements =vector.getNumElements();
  int * newIndices = vector.getIndices();
  double * newElements = vector.getElements();
  CoinSort_2(newIndices,newIndices+numberElements,newElements);
  bool notDuplicate=true;
  int numberRowCuts = sizeRowCuts();
  for ( int i =0; i<numberRowCuts;i++) {
    const OsiRowCut * cutPtr = rowCutPtr(i);
    if (cutPtr->row().getNumElements()!=numberElements)
      continue;
    if (!treatAsSame(cutPtr->lb(),newLb))
      continue;
    if (!treatAsSame(cutPtr->ub(),newUb))
      continue;
    const CoinPackedVector * thisVector = &(cutPtr->row());
    const int * indices = thisVector->getIndices();
    const double * elements = thisVector->getElements();
    int j;
    for(j=0;j<numberElements;j++) {
      if (indices[j]!=newIndices[j])
	break;
      if (!treatAsSame(elements[j],newElements[j]))
	break;
    }
    if (j==numberElements) {
      notDuplicate=false;
      break;
    }
  }
  if (notDuplicate) {
    OsiRowCut * newCutPtr = new OsiRowCut();
    newCutPtr->setLb(newLb);
    newCutPtr->setUb(newUb);
    newCutPtr->setRow(vector);
    rowCutPtrs_.push_back(newCutPtr);
  }
}

/** Clean an OsiCut
\return 1 if min violation is too small
\return 2 if small coefficient can not be removed
\return 3 if dynamic is too big
\return 4 if too many non zero element*/
int
Validator::cleanCut(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();

    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());
    int * indices = vec->getIndices();
    double * elems = vec->getElements();
    int n = vec->getNumElements();

    /** 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-8;
    double smallest = 1e100;
    double biggest = 0;
    for (int i = 0 ; i < n ; i++)
    {
        double val = fabs(elems[i]);
        if (val <= par.getEPS())   //try to remove coef
        {
            if (val>0 && val<1e-20)
            {
                offset++;
                continue;
                throw;
            }
            if (val==0)
            {
                offset++;
                continue;
            }

            int & iCol = indices[i];
            if (elems[i]>0. && colUpper[iCol] < 10000.)
            {
                offset++;
                rhs -= elems[i] * colUpper[iCol];
                elems[i]=0;
            }
            else if (elems[i]<0. && colLower[iCol] > -10000.)
            {
                offset++;
                rhs -= elems[i] * colLower[iCol];
                elems[i]=0.;
            }
            else
            {
#ifdef DEBUG
                std::cout<<"Small coefficient : "<<elems[i]<<" bounds : ["<<colLower[iCol]<<", "<<colUpper[iCol]<<std::endl;
#endif
                numRejected_[SmallCoefficient]++;
                return SmallCoefficient;
            }
        }

        else   //Not a small coefficient keep it
        {
            smallest = std::min(val,smallest);
            biggest = std::max (val,biggest);
            if (biggest > maxRatio_ * smallest)
            {
#ifdef DEBUG
                std::cout<<"Whaooo "<<biggest/smallest<<std::endl;
#endif
                numRejected_[BigDynamic]++;
                return BigDynamic;
            }
            if (offset)   //if offset is zero current values are ok otherwise translate
            {
                int i2 = i - offset;
                indices[i2] = indices[i];
                elems[i2] = elems[i];
            }
        }
    }
    if ((n - offset) > maxNnz)
    {
        numRejected_[DenseCut] ++;
        return DenseCut;
    }
//.........这里部分代码省略.........

//--------------------------------------------------------------------------
// test the simple rounding cut generators methods.
void
CglSimpleRoundingUnitTest(
  const OsiSolverInterface * baseSiP,
  const std::string mpsDir )
{

  // Test default constructor
  {
    CglSimpleRounding cg;
  }

  // Test copy & assignment
  {
    CglSimpleRounding rhs;
    {
      CglSimpleRounding cg;
      CglSimpleRounding cgC(cg);
      rhs=cg;
    }
  }

  // Test gcd and gcdn
  {
    CglSimpleRounding cg;
    int v = cg.gcd(122,356);
    assert(v==2);
    v=cg.gcd(356,122);
    assert(v==2);
    v=cg.gcd(54,67);
    assert(v==1);
    v=cg.gcd(67,54);
    assert(v==1);
    v=cg.gcd(485,485);
    assert(v==485);
    v=cg.gcd(17*13,17*23);
    assert( v==17);
    v=cg.gcd(17*13*5,17*23);
    assert( v==17);
    v=cg.gcd(17*13*23,17*23);
    assert(v==17*23);

    int a[4] = {12, 20, 32, 400};
    v= cg.gcdv(4,a);
    assert(v== 4);
    int b[4] = {782, 4692, 51, 2754};
    v= cg.gcdv(4,b);
    assert(v== 17);
    int c[4] = {50, 40, 30, 10};
    v= cg.gcdv(4,c);
    assert(v== 10);
  }


  // Test generate cuts method on exmip1.5.mps
  {
    CglSimpleRounding cg;
    
    OsiSolverInterface * siP = baseSiP->clone();
    std::string fn = mpsDir+"exmip1.5.mps";
    siP->readMps(fn.c_str(),"");
    OsiCuts cuts;
    cg.generateCuts(*siP,cuts);

    // there should be 3 cuts
    int nRowCuts = cuts.sizeRowCuts();
    assert(nRowCuts==3);

    // get the last "sr"=simple rounding cut that was derived
    OsiRowCut srRowCut2 = cuts.rowCut(2); 
    CoinPackedVector srRowCutPV2 = srRowCut2.row();

    // this is what the last cut should look like: i.e. the "solution"
    const int solSize = 2;
    int solCols[solSize]={2,3};
    double solCoefs[solSize]={5.0, 4.0};
    OsiRowCut solRowCut;
    solRowCut.setRow(solSize,solCols,solCoefs);
    solRowCut.setLb(-COIN_DBL_MAX);
    solRowCut.setUb(2.0);

    // Test for equality between the derived cut and the solution cut

    // Note: testing two OsiRowCuts are equal invokes testing two
    // CoinPackedVectors are equal which invokes testing two doubles
    // are equal.  Usually not a good idea to test that two doubles are equal, 
    // but in this cut the "doubles" represent integer values. Also allow that
    // different solvers have different orderings in packed vectors, which may
    // not match the ordering defined for solRowCut.

    assert(srRowCut2.OsiCut::operator==(solRowCut)) ;
    assert(srRowCut2.row().isEquivalent(solRowCut.row())) ;
    assert(srRowCut2.lb() == solRowCut.lb()) ;
    assert(srRowCut2.ub() == solRowCut.ub()) ;

    delete siP;
  }

  // Test generate cuts method on p0033
//.........这里部分代码省略.........

//--------------------------------------------------------------------------
// test EKKsolution methods.
void
CglProbingUnitTest(
  const OsiSolverInterface * baseSiP,
  const std::string mpsDir )
{
# ifdef CGL_DEBUG
  int i ;	// define just once
# endif
  CoinRelFltEq eq(0.000001);

  // Test default constructor
  {
    CglProbing aGenerator;
  }
  
  // Test copy & assignment
  {
    CglProbing rhs;
    {
      CglProbing bGenerator;
      CglProbing cGenerator(bGenerator);
      rhs=bGenerator;
    }
  }

  {
    OsiCuts osicuts;
    CglProbing test1;
    OsiSolverInterface  * siP = baseSiP->clone();
    int nColCuts;
    int nRowCuts;
    
    std::string fn = mpsDir+"p0033";
    siP->readMps(fn.c_str(),"mps");
    siP->initialSolve();
    // just unsatisfied variables
    test1.generateCuts(*siP,osicuts);
    nColCuts = osicuts.sizeColCuts();
    nRowCuts = osicuts.sizeRowCuts();
    std::cout<<"There are "<<nRowCuts<<" probing cuts"<<std::endl;
    {
      std::cout<<"there are "<<nColCuts<<" probing column cuts"<<std::endl;

#ifdef CGL_DEBUG
      const double * lo = siP->getColLower();
      const double * up = siP->getColUpper();
      for (i=0; i<nColCuts; i++){
	OsiColCut ccut;
	CoinPackedVector cpv;
	ccut = osicuts.colCut(i);
	cpv = ccut.lbs();
	int n = cpv.getNumElements();
        int j;
	const int * indices = cpv.getIndices();
	double* elements = cpv.getElements();
	for (j=0;j<n;j++) {
	  int icol=indices[j];
	  if (elements[j]>lo[icol])
	    std::cout<<"Can increase lb on "<<icol<<" from "<<lo[icol]<<
	      " to "<<elements[j]<<std::endl;
	}
	cpv = ccut.ubs();
	n = cpv.getNumElements();
	indices = cpv.getIndices();
	elements = cpv.getElements();

	for (j=0;j<n;j++) {
	  int icol=indices[j];
	  if (elements[j]<up[icol])
	    std::cout<<"Can decrease ub on "<<icol<<" from "<<up[icol]<<
	      " to "<<elements[j]<<std::endl;
	}
      }

#endif

    }

#ifdef CGL_DEBUG
    for (i=0; i<nRowCuts; i++){
      OsiRowCut rcut;
      CoinPackedVector rpv;
      const double * colsol = siP->getColSolution();
      rcut = osicuts.rowCut(i);
      rpv = rcut.row();
      const int n = rpv.getNumElements();
      const int * indices = rpv.getIndices();
      double* elements = rpv.getElements();
      double sum2=0.0;
      int k=0;
      double lb=rcut.lb();
      double ub=rcut.ub();
      for (k=0; k<n; k++){
	int column=indices[k];
	sum2 += colsol[column]*elements[k];
      }
      if (sum2 >ub + 1.0e-7 ||sum2 < lb - 1.0e-7) {
	std::cout<<"Cut "<<i<<" lb "<<lb<<" solution "<<sum2<<" ub "<<ub<<std::endl;
//.........这里部分代码省略.........

//.........这里部分代码省略.........
    // solution 1
    double colsol1[5]={20.0/7.0,3.0,0.0,0.0,23.0/7.0};
    test1.generateCuts(NULL, osicuts, matrix,
		       /*objective,*/ colsol1,
		 colLower, colUpper,
		 rowLower, rowUpper, intVar, &warm);
    nRowCuts = osicuts.sizeRowCuts();
    std::cout<<"There are "<<nRowCuts<<" gomory cuts"<<std::endl;
    assert (nRowCuts==2);
    // cuts always <=
    int testCut=0; // test first cut as stronger
    double rhs=-6.0;
    double testCut1[5]={0.0,0.0,-1.0,-2.0,0.0};
    double * cut = testCut1;
    double * colsol = colsol1;
    for (i=nOldCuts; i<nRowCuts; i++){
      OsiRowCut rcut;
      CoinPackedVector rpv;
      rcut = osicuts.rowCut(i);
      rpv = rcut.row();
      const int n = rpv.getNumElements();
      const int * indices = rpv.getIndices();
      double* elements = rpv.getElements();
      double sum2=0.0;
      int k=0;
      for (k=0; k<n; k++){
	int column=indices[k];
	sum2 += colsol[column]*elements[k];
      }

      double ub=rcut.ub();

#ifdef CGL_DEBUG
      double lb=rcut.lb();
      if (sum2 >ub + 1.0e-7 ||sum2 < lb - 1.0e-7) {
	std::cout<<"Cut "<<i<<" lb "<<lb<<" solution "<<sum2<<" ub "<<ub<<std::endl;
	for (k=0; k<n; k++){
	  int column=indices[k];
	  std::cout<<"(col="<<column<<",el="<<elements[k]<<",sol="<<
	    colsol[column]<<") ";
	}
	std::cout <<std::endl;
      }
#endif

      if (i-nOldCuts==testCut) {
	assert( eq(rhs,ub));
	assert(n==2);
	for (k=0; k<n; k++){
	  int column=indices[k];
	  assert (eq(cut[column],elements[k]));
	}
	// add cut
	// explicit slack
	matrix.setDimensions(-1,6);
	rpv.insert(5,1.0*7.0); // to get cut in book
	rowLower[3]=ub;
	rowUpper[3]=ub;
	matrix.appendRow(rpv);
      }
    }
    nOldCuts=nRowCuts;
    // basis 2
    int rowBasis2[4]={-1,-1,-1,-1};
    int colBasis2[6]={1,1,1,1,-1,-1};
    warm.setSize(6,4);

SCIP_RETCODE SCIPconshdlrBenders::sepaBenders(
		SCIP * scip,
		SCIP_CONSHDLR * conshdlr,
		SCIP_SOL * sol,
		whereFrom where,
		SCIP_RESULT * result)
{
	OsiCuts cs; /**< Benders cut placeholder */
	SCIP_Real * vals = NULL; /**< current solution */

#if 1
	if (scip_checkpriority_ < 0)
	{
		/** consider incumbent solutions only */
		double primObj = SCIPgetPrimalbound(scip);
		double currObj = SCIPgetSolOrigObj(scip, sol);
		if (SCIPisLT(scip, primObj, currObj))
		{
			DSPdebugMessage(" -> primObj %e currObj %e\n", primObj, currObj);
			return SCIP_OKAY;
		}
	}
#endif

	/** allocate memory */
	SCIP_CALL(SCIPallocMemoryArray(scip, &vals, nvars_));

	/** get current solution */
	SCIP_CALL(SCIPgetSolVals(scip, sol, nvars_, vars_, vals));

	/** TODO The following filter does not work, meaning that it provides suboptimal solution.
	 * I do not know the reason. */
#if 0
	double maxviol = 1.e-10;
	for (int j = 0; j < nvars_ - naux_; ++j)
	{
		SCIP_VARTYPE vartype = SCIPvarGetType(vars_[j]);
		if (vartype == SCIP_VARTYPE_CONTINUOUS) continue;

		double viol = 0.5 - fabs(vals[j] - floor(vals[j]) - 0.5);
		if (viol > maxviol)
			maxviol = viol;
	}
	DSPdebugMessage("maximum violation %e\n", maxviol);

	if (where != from_scip_check &&
		where != from_scip_enfolp &&
		where != from_scip_enfops &&
		maxviol > 1.e-7)
	{
		printf("where %d maxviol %e\n", where, maxviol);
		/** free memory */
		SCIPfreeMemoryArray(scip, &vals);
		return SCIP_OKAY;
	}
#endif

#ifdef DSP_DEBUG2
	double minvals = COIN_DBL_MAX;
	double maxvals = -COIN_DBL_MAX;
	double sumvals = 0.;
	double ssvals  = 0.;
	//printf("nvars_ %d naux_ %d nAuxvars_ %d\n", nvars_, naux_, tss_->nAuxvars_);
	for (int j = 0; j < nvars_ - naux_; ++j)
	{
//		if (vals[j] < 0 || vals[j] > 1)
//			printf("solution %d has value %e.\n", j, vals[j]);
		sumvals += vals[j];
		ssvals  += vals[j] * vals[j];
		minvals = minvals > vals[j] ? vals[j] : minvals;
		maxvals = maxvals < vals[j] ? vals[j] : maxvals;
	}
	DSPdebugMessage("solution: min %e max %e avg %e sum %e two-norm %e\n",
			minvals, maxvals, sumvals / nvars_, sumvals, sqrt(ssvals));
#endif

#define SCAN_GLOBAL_CUT_POOL
#ifdef SCAN_GLOBAL_CUT_POOL
	if (SCIPgetStage(scip) == SCIP_STAGE_SOLVING ||
		SCIPgetStage(scip) == SCIP_STAGE_SOLVED ||
		SCIPgetStage(scip) == SCIP_STAGE_EXITSOLVE)
	{
		bool addedPoolCut = false;
		int numPoolCuts = SCIPgetNPoolCuts(scip);
		int numCutsToScan = 100;
		SCIP_CUT ** poolcuts = SCIPgetPoolCuts(scip);
		for (int i = numPoolCuts - 1; i >= 0; --i)
		{
			if (i < 0) break;
			if (numCutsToScan == 0) break;

			/** retrieve row */
			SCIP_ROW * poolcutrow = SCIPcutGetRow(poolcuts[i]);

			/** benders? */
			if (strcmp(SCIProwGetName(poolcutrow), "benders") != 0)
				continue;

			/** counter */
			numCutsToScan--;
//.........这里部分代码省略.........

//.........这里部分代码省略.........
        // z=              -0.25 s1  -0.25 s2             =  -1
        // Gomory cut from variable x1 is x2 <= 0.5
        // Can be improved by first pivoting s2 in and s4 out, then s1 in and s3 out
        // to x2 <= 0.25
        {
            int start[2] = {0,4};
            int length[2] = {4,4};
            int rows[8] = {0,1,2,3,0,1,2,3};
            double elements[8] = {2.0,-2.0,7.0,-7.0,2.0,2.0,4.0,4.0};
            CoinPackedMatrix  columnCopy(true,4,2,8,elements,rows,start,length);

            double rowLower[4]={-COIN_DBL_MAX,-COIN_DBL_MAX,
                                -COIN_DBL_MAX,-COIN_DBL_MAX};
            double rowUpper[4]={3.,1.,8.,1.};
            double colLower[2]={0.0,0.0};
            double colUpper[2]={1.0,1.0};
            double obj[2]={-1,-1};
            int intVar[2]={0,1};

            OsiSolverInterface  * siP = si->clone();
            siP->loadProblem(columnCopy, colLower, colUpper, obj, rowLower, rowUpper);
            siP->setInteger(intVar,2);
            CglLandP test;
            test.setLogLevel(2);
            test.parameter().sepSpace = CglLandP::Full;
            siP->resolve();
            // Test generateCuts method
            {
                OsiCuts cuts;
                test.generateCuts(*siP,cuts);
                cuts.printCuts();
                assert(cuts.sizeRowCuts()==1);
                OsiRowCut aCut = cuts.rowCut(0);
                assert(eq(aCut.lb(), -.0714286));
                CoinPackedVector row = aCut.row();
                if (row.getNumElements() == 1)
                {
                    assert(row.getIndices()[0]==1);
                    assert(eq(row.getElements()[0], -4*.0714286));
                }
                else if (row.getNumElements() == 2)
                {
                    assert(row.getIndices()[0]==0);
                    assert(eq(row.getElements()[0], 0.));
                    assert(row.getIndices()[1]==1);
                    assert(eq(row.getElements()[1], -1));
                }
                OsiSolverInterface::ApplyCutsReturnCode rc = siP->applyCuts(cuts);

                siP->resolve();
            }
            if (0)
            {
                OsiCuts cuts;
                test.generateCuts(*siP,cuts);
                cuts.printCuts();
                assert(cuts.sizeRowCuts()==1);
                OsiRowCut aCut = cuts.rowCut(0);
                CoinPackedVector row = aCut.row();
                if (row.getNumElements() == 1)
                {
                    assert(row.getIndices()[0]==1);
                    assert(eq(row.getElements()[0], -1));
                }
                else if (row.getNumElements() == 2)
                {

//.........这里部分代码省略.........
	upper[i]=1.0;
	int marked=start[i];
	start[i]=last;
	if (marked)
	  last++;
      }
      siP->loadProblem(ncol,1,start,row,el1,NULL,upper,NULL,&lo,&up);
      // use upper for solution
      memset(upper,0,ncol*sizeof(double));
      for (i=0;i<nel;i++) {
	int icol=col1[i];
	upper[icol]=sol1[i];
	siP->setInteger(icol);
      }
      siP->setColSolution(upper);
      delete [] sol1;
      delete [] el1;
      delete [] col1;
      delete [] start;
      delete [] row;
      delete [] upper;
      CglKnapsackCover kccg;
      
      OsiCuts cuts;    
      
      // Test generateCuts method
      kccg.generateCuts(*siP,cuts);
      // print out and compare to known cuts
      int numberCuts = cuts.sizeRowCuts();
      if (numberCuts) {
	for (i=0;i<numberCuts;i++) {
	  OsiRowCut * thisCut = cuts.rowCutPtr(i);
	  int n=thisCut->row().getNumElements();
	  printf("Cut %d has %d entries, rhs %g %g =>",i,n,thisCut->lb(),
		 thisCut->ub());
	  int j;
	  const int * index = thisCut->row().getIndices();
	  const double * element = thisCut->row().getElements();
	  for (j=0;j<n;j++) {
	    printf(" (%d,%g)",index[j],element[j]);
	  }
	  printf("\n");
	}
      }
      fclose(fp);
    }
  }

  // Testcase /u/rlh/osl2/mps/p0201
  // Miplib3 problem p0282
  // Test that no cuts chop off the optimal ip solution
  {
    // Setup
    OsiSolverInterface  * siP = baseSiP->clone();
    std::string fn(mpsDir+"p0201");
    siP->readMps(fn.c_str(),"mps");
    // All integer variables should be binary.
    // Assert that this is true.
    for ( i = 0;  i < siP->getNumCols();  i++ )
      if ( siP->isInteger(i) ) 
        assert(siP->getColUpper()[i]==1.0 && siP->isBinary(i));    

    const int nCols=siP->getNumCols();
    CglKnapsackCover kccg;
    
    // Solve the LP relaxation of the model and