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

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

/** scale the cut passed as argument using provided normalization factor*/
void scale(OsiRowCut &cut, double norma)
{
    assert(norma >0.);
    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]/norma);
    }
    cut.setLb(cut.lb()/norma);
    cut.setRow(row);
}

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

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

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

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

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

//.........这里部分代码省略.........
#endif
	    good=false;
	  } else {
	    if (static_cast<double> (ii) * minimumViolationPer_>violation||
		ii>maximumEntries_) {
#ifdef CGL_DEBUG
	      printf("why no cut\n");
#endif
	      if (packed) {
		// sort and see if we can get down to length
		// relax by taking out ones with solution 0.0
		nincut=ii;
		for (k=0;k<nincut;k++) {
		  int jcol=candidate[k];
		  double value = fabs(dj[jcol]);
		  if (solution[jcol])
		  value = -solution[jcol];
		  sortit[k].dj=value;
		  sortit[k].element=element[k];
		  sortit[k].sequence=jcol;
		}
		// sort 
		std::sort(sortit,sortit+nincut,double_double_int_triple_compare());
		nincut = CoinMin(nincut,maximumEntries_);
		sum=0.0;
		for (k=0;k<nincut;k++) {
		  int jcol=sortit[k].sequence;
		  candidate[k]=jcol;
		  element[k]=sortit[k].element;
		  sum+=solution[jcol]*element[k];
		}
		violation = sum-rhs;
		ii=nincut;
		if (violation<minimumViolation_) {
		  good=false;
		}
	      } else { 
		good=false;
	      }
	    }
	  }
	  if (good) {
	    //this assumes not many cuts
	    int j;
#if 0
	    double value=0.0;
	    for (j=0;j<ii;j++) {
	      int icol=candidate[j];
	      value += check[icol]*element[j];
	    }
#else
            CoinPackedVector candidatePv(ii,candidate,element);
            candidatePv.sortIncrIndex();
            double value = candidatePv.dotProduct(check);
#endif

	    for (j=0;j<ncuts;j++) {
	      if (value==hash[j]) {
		//could check equality - quicker just to assume
		break;
	      }
	    }
	    if (j==ncuts) {
	      //new
	      if (ncuts==maxcuts) {
		maxcuts *= 2;
		hash = reinterpret_cast<double *> (realloc(hash,maxcuts*sizeof(double)));
	      }
	      hash[ncuts++]=value;
	      rc.setRow(ii,candidate,element);
#ifdef CGL_DEBUG
	      printf("sum %g rhs %g %d\n",sum,rhs,ii);
	      if (debugger) 
		assert(!debugger->invalidCut(rc)); 
#endif
	      cs.insert(rc);
	    }
	  }
	}
	/* end of adding cut */
      }
    }
  }
  delete [] countcol;
  delete [] element;
  delete [] candidate;
  delete [] sortit;
  delete [] clean;
  delete [] path;
  delete [] stack;
  free(hash);
  delete [] check;
  delete [] mark;
  delete [] starts;
  delete [] lookup;
  delete [] mrow;
  free(rowfound);
  free(to);
  free(cost);
}

//.........这里部分代码省略.........
	  ilo= static_cast<int> (ceil(lo));
	  double up = columnUpper[iColumn];
	  if (up>COIN_INT_MAX)
	    up=COIN_INT_MAX;
	  iup= static_cast<int> (floor(up));
	  vlb_[iColumn]=ilo;
	  vub_[iColumn]=iup;
	}
      }
    }
    if (true) {
    cutInfo_.sep_012_cut(mr_,mc_,mnz_,
				 mtbeg_,mtcnt_, mtind_, mtval_,
				 vlb_, vub_,
				 mrhs_, msense_,
				 solution,
				 info.inTree ? false : true,
				 &cnum,&cnzcnt,
				 &cbeg,&ccnt,&cind,&cval,&crhs,&csense);
    } else {
      int k = 4*mr_+2*mnz_;
      int * temp = new int[k];
      int * mtbeg = temp;
      int * mtcnt = mtbeg + mr_;
      int * mtind = mtcnt+mr_;
      int * mtval = mtind+mnz_;
      int * mrhs = mtval+mnz_;
      char * msense = reinterpret_cast<char*> (mrhs+mr_);
      int i;
      k=0;
      int kel=0;
      for (i=0;i<mr_;i++) {
	int kel2=kel;
	int rhs = mrhs_[i];
	for (int j=mtbeg_[i];j<mtbeg_[i]+mtcnt_[i];j++) {
	  int iColumn=mtind_[j];
	  int value=mtval_[j];
	  if (vlb_[iColumn]<vub_[iColumn]) {
	    mtind[kel]=mtind_[j];
	    mtval[kel++]=mtval_[j];
	  } else {
	    rhs -= vlb_[iColumn]*value;
	  }
	}
	if (kel>kel2) {
	  mtcnt[k]=kel-kel2;
	  mtbeg[k]=kel2;
	  mrhs[k]=rhs;
	  msense[k++]=msense_[i];
	}
      }
      if (kel) {
	cutInfo_.sep_012_cut(k,mc_,kel,
				 mtbeg,mtcnt, mtind, mtval,
				 vlb_, vub_,
				 mrhs, msense,
				 solution,
				 info.inTree ? false : true,
				 &cnum,&cnzcnt,
				 &cbeg,&ccnt,&cind,&cval,&crhs,&csense);
      }
      delete [] temp;
    }
    if (cnum) {
      // add cuts
      double * element = new double[mc_];
      for (int i=0;i<cnum;i++) {
	int n = ccnt[i];
	int start = cbeg[i];
	for (int j=0;j<n;j++) 
	  element[j]=cval[start+j];
	OsiRowCut rc;
	if (csense[i]=='L') {
	  rc.setLb(-COIN_DBL_MAX);
	  rc.setUb(crhs[i]);
	} else if (csense[i]=='G') {
	  rc.setLb(crhs[i]);
	  rc.setUb(COIN_DBL_MAX);
	} else {
	  abort();
	}
	rc.setRow(n,cind+start,element,false);
	if ((flags_&1)!=0)
	  rc.setGloballyValid();
	//double violation = rc.violated(solution);
	//if (violation>1.0e-6)
	  cs.insert(rc);
	  //else
	  //printf("violation of %g\n",violation);
      }
      delete [] element;
      free(cbeg); 
      free(ccnt);
      free(cind);
      free(cval);
      free(crhs);
      free(csense);
    }
  }
}

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

//.........这里部分代码省略.........
                        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;
            }
            nSort = CoinMin(n, numberClean_ / 1000000);
        }
    } else {
#define FIX_IF_LESS -0.1
        // 3 in same row and sum <FIX_IF_LESS?
        int numberRows = matrixByRow_.getNumRows();
        const double * solution = model_->testSolution();
        const int * column = matrixByRow_.getIndices();
        const CoinBigIndex * rowStart = matrixByRow_.getVectorStarts();
        const int * rowLength = matrixByRow_.getVectorLengths();
        double bestSum = 1.0;
        int nBest = -1;
        int kRow = -1;
        OsiSolverInterface * solver = model_->solver();
        for (int i = 0; i < numberRows; i++) {
            int numberUnsatisfied = 0;
            double sum = 0.0;
            for (int j = rowStart[i]; j < rowStart[i] + rowLength[i]; j++) {
                int iColumn = column[j];
                if (solver->isInteger(iColumn)) {
                    double solValue = solution[iColumn];
                    if (solValue > 1.0e-5 && solValue < FIX_IF_LESS) {
                        numberUnsatisfied++;
                        sum += solValue;
                    }
                }
            }
            if (numberUnsatisfied >= 3 && sum < FIX_IF_LESS) {
                // possible
                if (numberUnsatisfied > nBest ||
                        (numberUnsatisfied == nBest && sum < bestSum)) {
                    nBest = numberUnsatisfied;
                    bestSum = sum;
                    kRow = i;
                }
            }
        }
        assert (nBest > 0);
        for (int j = rowStart[kRow]; j < rowStart[kRow] + rowLength[kRow]; j++) {
            int iColumn = column[j];
            if (solver->isInteger(iColumn)) {
                double solValue = solution[iColumn];
                if (solValue > 1.0e-5 && solValue < FIX_IF_LESS) {
                    sort[nSort++] = iColumn;
                }
            }
        }
    }
    OsiRowCut down;
    down.setLb(-COIN_DBL_MAX);
    double rhs = 0.0;
    for (i = 0; i < nSort; i++) {
        int iColumn = sort[i];
        double distanceDown = solution[iColumn] - lower[iColumn];
        double distanceUp = upper[iColumn] - solution[iColumn];
        if (distanceDown < distanceUp) {
            rhs += lower[iColumn];
            dsort[i] = 1.0;
        } else {
            rhs -= upper[iColumn];
            dsort[i] = -1.0;
        }
    }
    down.setUb(rhs);
    down.setRow(nSort, sort, dsort);
    down.setEffectiveness(COIN_DBL_MAX); // so will persist
    delete [] sort;
    delete [] dsort;
    // up is same - just with rhs changed
    OsiRowCut up = down;
    up.setLb(rhs + 1.0);
    up.setUb(COIN_DBL_MAX);
    // Say can fix one way
    CbcCutBranchingObject * newObject =
        new CbcCutBranchingObject(model_, down, up, true);
    if (model_->messageHandler()->logLevel() > 1)
        printf("creating cut in CbcBranchCut\n");
    return newObject;
}

/** Get an inner-approximation constraint obtained by drawing a chord linking the two given points x and x2. 
 * This only applies to nonlinear constraints featuring univariate functions (f(x) <= y).**/
bool
HeuristicInnerApproximation::getMyInnerApproximation(Bonmin::OsiTMINLPInterface &si, OsiCuts &cs, int ind,
    const double * x, const double * x2) {

  int n, m, nnz_jac_g, nnz_h_lag;
  Ipopt::TNLP::IndexStyleEnum index_style;
  Bonmin::TMINLP2TNLP * problem = si.problem(); 
  problem->get_nlp_info(n, m, nnz_jac_g, nnz_h_lag, index_style);
  double infty = si.getInfinity();


  CoinPackedVector cut;
  double lb = -infty;
  double ub = 0;


  double g = 0;
  double g2 = 0;
  double diff = 0;
  double a = 0;
  problem->eval_gi(n, x, 1, ind, g);
  problem->eval_gi(n, x2, 1, ind, g2);
  Bonmin::vector<int> jCol(n);
  int nnz;
  problem->eval_grad_gi(n, x2, 0, ind, nnz, jCol(), NULL);
  Bonmin::vector<double> jValues(nnz);
  problem->eval_grad_gi(n, x2, 0, ind, nnz, NULL, jValues());
  bool add = false;
  //printf("const %i nnz %i\n", ind, nnz);
  for (int i = 0; i < nnz; i++) {
    const int &colIdx = jCol[i];
    if(index_style == Ipopt::TNLP::FORTRAN_STYLE) jCol[i]--;
    diff = x[colIdx] - x2[colIdx];

    if (fabs(diff) >= 1e-8) {
       a = (g - g2) / diff;
       cut.insert(colIdx, a);
       ub = (a * x[colIdx] - g) - fabs(a * x[colIdx] - g)*1e-6;
       //printf("const %i col %i p[col] %g pp[col] %g g %g g2 %g diff %g\n",ind, colIdx, x[colIdx], x2[colIdx], g, g2, diff);
       add = true;
    } else {
       cut.insert(colIdx, jValues[i]);
       //printf("const %i col %i val %g\n",ind, colIdx, jValues[i]);
    }
  }

  if (add) {

    OsiRowCut newCut;
    newCut.setGloballyValidAsInteger(1);
    newCut.setLb(lb);
    
      //********* Perspective Extension ********//
    int binary_id = 0; // index corresponding to the binary variable activating the corresponding constraint
    const int* ids = problem->get_const_xtra_id(); // vector of indices corresponding to the binary variable activating the corresponding constraint
    // Get the index of the corresponding indicator binary variable
    binary_id = (ids == NULL) ? -1 : ids[ind];
    if(binary_id>0) {// If this hyperplane is a linearization of a disjunctive constraint, we link its righthand side to the corresponding indicator binary variable
        cut.insert(binary_id, -ub); // ∂x ≤ ub => ∂x - ub*z ≤ 0
        newCut.setUb(0);
    }
    else
        newCut.setUb(ub);
    //********* Perspective Extension ********//

    newCut.setRow(cut);
    cs.insert(newCut);
    //newCut.print();
    return true;
  }
  return false;
}

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

//.........这里部分代码省略.........
    assert( cover.getElements()[0] == 161.0);
    assert( cover.getElements()[1] == 120.0);
    assert( remainder.getNumElements() == 1);
    assert( remainder.getIndices()[0] == 2);
    assert( remainder.getElements()[0] == 68.0);

    // test liftCoverCut
    CoinPackedVector cut;
    double * rowupper = ekk_rowupper(model);
    double cutRhs = cover.getNumElements() - 1.0;
    kccg.liftCoverCut(b, krow.getNumElements(),
      cover, remainder,
      cut);
    assert ( cut.getNumElements() == 3 );
    assert ( cut.getIndices()[0] == 0 );
    assert ( cut.getIndices()[1] == 1 );
    assert ( cut.getIndices()[2] == 2 );
    assert( cut.getElements()[0] == 1 );
    assert( cut.getElements()[1] == 1 );
    assert( eq(cut.getElements()[2], 0.087719) );
    
    // test liftAndUncomplementAndAdd
    OsiCuts cuts;    
    kccg.liftAndUncomplementAndAdd(*siP.getRowUpper()[0],krow,b,complement,0,
      cover,remainder,cuts);   
    int sizerowcuts = cuts.sizeRowCuts();
    assert ( sizerowcuts== 1 );
    OsiRowCut testRowCut = cuts.rowCut(0);
    CoinPackedVector testRowPV = testRowCut.row(); 
    OsiRowCut sampleRowCut;
    const int sampleSize = 3;
    int sampleCols[sampleSize]={0,1,2};
    double sampleElems[sampleSize]={1.0,-1.0,-0.087719};
    sampleRowCut.setRow(sampleSize,sampleCols,sampleElems);
    sampleRowCut.setLb(-DBL_MAX);
    sampleRowCut.setUb(-0.087719);
    bool equiv =  testRowPV.equivalent(sampleRowCut.row(),CoinRelFltEq(1.0e-05) );
    assert ( equiv );
#endif
    
    // test find PseudoJohnAndEllisCover
    cover.setVector(0,NULL, NULL);
    remainder.setVector(0,NULL,NULL);

    rind = ( siP->getRowSense()[0] == 'N' ) ? 1 : 0;
    int findPJE =  kccg.findPseudoJohnAndEllisCover( rind, krow, 
						     b, xstar, cover, remainder );
    assert( findPJE == 1 );
    assert ( cover.getIndices()[0] == 0 );
    assert ( cover.getIndices()[1] == 2 );
    assert ( cover.getElements()[0] == 161 );    
    assert ( cover.getElements()[1] == 68 );    
    assert ( remainder.getIndices()[0] == 1 );
    assert ( remainder.getElements()[0] == 120 );    
    OsiCuts cuts;    
    kccg.liftAndUncomplementAndAdd((*siP).getRowUpper()[rind],krow,b, complement, rind,
      cover,remainder,cuts);   
    assert (cuts.sizeRowCuts() == 1 );

    OsiRowCut testRowCut = cuts.rowCut(0);
    CoinPackedVector testRowPV = testRowCut.row();


    const int sampleSize = 3;
    int sampleCols[sampleSize]={0,1,2};
    double sampleElems[sampleSize]={1.0, -1.0, -1.0};