C++ OsiRowCut类代码示例

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

// Add a row cut from elements
void
CglStored::addCut(double lb, double ub, int size, const int * colIndices, const double * elements)
{
  OsiRowCut rc;
  rc.setRow(size,colIndices,elements,false);
  rc.setLb(lb);
  rc.setUb(ub);
  cuts_.insert(rc);
}

//-------------------------------------------------------------------
// Constructor from file
//-------------------------------------------------------------------
CglStored::CglStored (const char * fileName) :
  CglCutGenerator(),
  requiredViolation_(1.0e-5),
  probingInfo_(NULL),
  numberColumns_(0),
  bestSolution_(NULL),
  bounds_(NULL)
{  
  FILE * fp = fopen(fileName,"rb");
  if (fp) {
#ifndef NDEBUG
    size_t numberRead;
#endif
    int maxInCut=0;
    int * index = NULL;
    double * coefficient = NULL;
    double rhs[2];
    int n=0;
    while (n>=0) {
#ifndef NDEBUG
      numberRead = fread(&n,sizeof(int),1,fp);
      assert (numberRead==1);
#else
      fread(&n,sizeof(int),1,fp);
#endif
      if (n<0)
	break;
      if (n>maxInCut) {
	maxInCut=n;
	delete [] index;
	delete [] coefficient;
	index = new int [maxInCut];
	coefficient = new double [maxInCut];
      }
#ifndef NDEBUG
      numberRead = fread(rhs,sizeof(double),2,fp);
      assert (numberRead==2);
#else
      fread(rhs,sizeof(double),2,fp);
#endif
      fread(index,sizeof(int),n,fp);
      fread(coefficient,sizeof(double),n,fp);
      OsiRowCut rc;
      rc.setRow(n,index,coefficient,false);
      rc.setLb(rhs[0]);
      rc.setUb(rhs[1]);   
      cuts_.insert(rc);
    }
    delete [] index;
    delete [] coefficient;
    fclose(fp);
  }
}

void
CglClique::recordClique(const int len, int* indices, OsiCuts& cs)
{
   /* transform relative indices into user indices and order them */
   for (int j = len - 1; j >= 0; j--)
      indices[j] = sp_orig_col_ind[indices[j]];
   std::sort(indices, indices + len);
   OsiRowCut rowcut;
   double* coef = new double[len];
   std::fill(coef, coef + len, 1.0);
   rowcut.setRow(len, indices, coef);
   rowcut.setUb(1.0);
   CoinAbsFltEq equal(1.0e-12);
   cs.insertIfNotDuplicate(rowcut,equal);
   delete[] coef;
}

QuadRow::QuadRow(const OsiRowCut &cut):
  c_(0),
  a_(cut.row()),
  Q_()
  {
    initialize(); 
  }

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

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

OsiRowCut *
BlisConstraint::createOsiRowCut()
{
    double lower = CoinMax(getLbHard(), getLbSoft());
    double upper = CoinMin(getUbHard(), getUbSoft());

    OsiRowCut * cut = new OsiRowCut;
    if (!cut) {
        /* Out of memory. */
        throw CoinError("Out of Memory", "Blis_constraintToOsiCut", "NONE");
    }

    assert(size_ > 0);

    cut->setLb(lower);
    cut->setUb(upper);
    cut->setRow(size_, indices_, values_);

    return cut;
}

// approximate the cone around given point.
// If given point is in interior do nothing.
// if it is on the boundry add support
// We do not expect it to be infeasible for now. This may change in future
// in case of infeasibility we will just call separate routine.
// feas of given solution:
// Note that the numerical accuracy of the input sol
// depends on the underlying solver that created sol, and
// its accuracy level set.
// current strategy:
// the point we are given is not very close to the boundry
// do a computationally cheap projection of the point to te boundry
// and generate cuts that support cone.
// For now we just increase leadin var and project the point to the boundry.
void LorentzCone::approximate(double const * sol, OsiCuts * cuts) {
  // todo(aykut) improve accuracy of given solution
  // we do nothing for this for now.
  // coef array, [2x1, -2x2, -2x3, ... -2xn]
  double * coef = new double[size_];
  double * p = new double[size_];
  for (int i=0; i<size_; ++i) {
    p[i] = sol[members_[i]];
  }
  // 2. compute coefficient from given solution
  if (type()==LORENTZ) {
    // cone is in canonical form
    for (int i=1; i<size_; ++i) {
      coef[i] = 2.0*p[i];
    }
    coef[0] = -2.0*p[0];
  }
  else if (type()==RLORENTZ) {
    // map solution from RLORENTZ space to LORENTZ space, find the projection
    // on LORENTZ, project this point to RLORENTZ and generate cut
    // cone is a rotated cone
    // generate cut from xbar
    coef[0] = -2.0*p[1];
    coef[1] = -2.0*p[0];
    for (int i=2; i<size_; ++i) {
      coef[i] = 2.0*p[i];
    }
  }
  // rhs is allways 0.0
  OsiRowCut * cut = new OsiRowCut();
  cut->setRow(size_, members_, coef);
  cuts->insert(*cut);
  delete[] coef;
  delete[] p;
  delete cut;
}

CbcBranchingObject *
CbcBranchAllDifferent::createCbcBranch(OsiSolverInterface * /*solver*/
                                       , const OsiBranchingInformation * /*info*/,
                                       int /*way*/)
{
    // by default way must be -1
    //assert (way==-1);
    const double * solution = model_->testSolution();
    double * values = new double[numberInSet_];
    int * which = new int[numberInSet_];
    int i;
    for (i = 0; i < numberInSet_; i++) {
        int iColumn = which_[i];
        values[i] = solution[iColumn];
        which[i] = iColumn;
    }
    CoinSort_2(values, values + numberInSet_, which);
    double last = -1.0;
    double closest = 1.0;
    int worst = -1;
    for (i = 0; i < numberInSet_; i++) {
        if (values[i] - last < closest) {
            closest = values[i] - last;
            worst = i - 1;
        }
        last = values[i];
    }
    assert (closest <= 0.99999);
    OsiRowCut down;
    down.setLb(-COIN_DBL_MAX);
    down.setUb(-1.0);
    int pair[2];
    double elements[] = {1.0, -1.0};
    pair[0] = which[worst];
    pair[1] = which[worst+1];
    delete [] values;
    delete [] which;
    down.setRow(2, pair, elements);
    // up is same - just with rhs changed
    OsiRowCut up = down;
    up.setLb(1.0);
    up.setUb(COIN_DBL_MAX);
    // Say is not a fix type branch
    CbcCutBranchingObject * newObject =
        new CbcCutBranchingObject(model_, down, up, false);
    if (model_->messageHandler()->logLevel() > 1)
        printf("creating cut in CbcBranchCut\n");
    return newObject;
}

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

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

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

//-------------------------------------------------------------------
//  Given the model data, a row of the model, and a LP solution,
//  this function tries to generate a violated lifted simple generalized
//  flow cover.
//-------------------------------------------------------------------
bool
CglFlowCover::generateOneFlowCut( const OsiSolverInterface & si,
				  const int rowLen,
				  int* ind,
				  double* coef,
				  char sense,
				  double rhs,
				  OsiRowCut& flowCut,
				  double& violation ) const
{
  bool generated       = false;
  const double* xlp    = si.getColSolution();
  const int numCols    = si.getNumCols();

  double* up           = new double [rowLen];
  double* x            = new double [rowLen];
  double* y            = new double [rowLen];
  CglFlowColType* sign = new CglFlowColType [rowLen];

  int i, j;
  double value, LB, UB;

  CglFlowVLB VLB;
  CglFlowVUB VUB;
  static int count=0;
  ++count;
  CGLFLOW_DEBUG=false;
  doLift=true;
  // Get integer types
  const char * columnType = si.getColType ();
  for (i = 0; i < rowLen; ++i) {
    if ( xlp[ind[i]] - floor(xlp[ind[i]]) > EPSILON_ && ceil(xlp[ind[i]]) - xlp[ind[i]] > EPSILON_ )
      break;
  }

  if (i == rowLen)  {
    delete [] sign;
    delete [] up;
    delete [] x;
    delete [] y;
    return generated;
  }

  //-------------------------------------------------------------------------

  if (sense == 'G') flipRow(rowLen, coef, rhs); // flips everything,
  // but the sense


  if(CGLFLOW_DEBUG) {
    std::cout << "***************************" << std::endl;
    std::cout << "Generate Flow cover -- initial constraint, converted to L sense..." << std::endl;
    std::cout << "Rhs = " << rhs << std::endl;
    std::cout << "coef [var_index]" << " -- " <<  "xlp[var_index]" << '\t' << "vub_coef[vub_index] vub_lp_value OR var_index_col_ub" << std::endl;

    for(int iD = 0; iD < rowLen; ++iD) {
      VUB = getVubs(ind[iD]);

      std::cout << std::setw(5) << coef[iD] << "["  << std::setw(5)  << ind[iD] << "] -- "
		<< std::setw(20) << xlp[ind[iD]] << '\t';
      if (VUB.getVar() != UNDEFINED_) {
	std::cout << std::setw(20) << VUB.getVal() << "[" << std::setw(5) << VUB.getVar() << "]"
		  << std::setw(20) << xlp[VUB.getVar()] << std::endl;
      }
      else
	std::cout << std::setw(20) << si.getColUpper()[ind[iD]] << "       " << std::setw(20) << 1.0 << std::endl;

    }
  }

  //-------------------------------------------------------------------------
  // Generate conservation inequality and capacity equalities from
  // the given row.

  for (i = 0; i < rowLen; ++i) {

    VLB = getVlbs(ind[i]);
    LB = ( VLB.getVar() != UNDEFINED_ ) ?
      VLB.getVal() : si.getColLower()[ind[i]];

    VUB = getVubs(ind[i]);
    UB = ( VUB.getVar() != UNDEFINED_ ) ?
      VUB.getVal() : si.getColUpper()[ind[i]];

    if (LB < -EPSILON_) {   // Only consider rows whose variables are all
      delete [] sign;       // non-negative (LB>= 0).
      delete [] up;
      delete [] x;
      delete [] y;
      return generated;
    }

    if ( columnType[ind[i]]==1 ) {   // Binary variable
      value = coef[i];
      if (value > EPSILON_)
//.........这里部分代码省略.........