C++ CoinPackedVector类代码示例

//.........这里部分代码省略.........

					if (infeasible)
						*result = SCIP_CUTOFF;
					else //if (*result != SCIP_CUTOFF)
						*result = SCIP_SEPARATED;
				}
				else
					*result = SCIP_INFEASIBLE;
				addedPoolCut = true;
				break;
			}
		}
		if (addedPoolCut)
		{
			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 */

//===========================================================================//
DecompConstraintSet * ATM_DecompApp::createModelRelax2(const int d){
   
   UtilPrintFuncBegin(m_osLog, m_classTag,
		      "createModelRelax2()", m_appParam.LogLevel, 2);
   
   int                a, t, nRows, pairIndex, nRowsMax;
   double             rhs, coefA, coefC;
   pair<int,int>      adP;
   vector<int>::const_iterator vi;
   
   const int           nSteps    = m_appParam.NumSteps;
   const int           nAtms     = m_instance.getNAtms();
   const vector<int> & pairsAD   = m_instance.getPairsAD();   

   const double      * a_ad      = m_instance.get_a_ad(); //dense storage
   const double      * b_ad      = m_instance.get_b_ad(); //dense storage
   const double      * c_ad      = m_instance.get_c_ad(); //dense storage
   const double      * d_ad      = m_instance.get_d_ad(); //dense storage
   const double      * e_ad      = m_instance.get_e_ad(); //dense storage
   const double      * w_ad      = m_instance.get_w_ad(); //dense storage
   const double      * B_d       = m_instance.get_B_d();

   const int           nCols     = numCoreCols();
   if(m_appParam.UseTightModel)
      nRowsMax  = 3*nAtms + 1 + getNAtmsSteps();
   else
      nRowsMax  = 2*nAtms + 1 + (3*getNAtmsSteps());

   //---
   //--- A'[d] for d in D (independent blocks)
   //---    for a in A
   //---       f+[a,d] - f-[a,d] =  
   //---         a[a,d] sum{t in T} (t/n) x1[a,t] +
   //---         b[a,d] x2[a]                     +     
   //---         c[a,d] sum{t in T} (t/n) z[a,t]  + 
   //---         d[a,d] x3[a]                     +
   //---         e[a,d]
   //---       f-[a,d] <= w[a,d] * v[a,d]
   //---    sum{a in A} (f+[a,d] - f-[a,d]) <= B[d]
   //---    for a in A, t in T:
   //---       z[a,t] <= x1[a,t],            
   //---       z[a,t] <= x2[a],              
   //---       z[a,t] >= x1[a,t] + x2[a] - 1.
   //---
   
   DecompConstraintSet * model = new DecompConstraintSet();
   CoinAssertHint(model, "Error: Out of Memory");

   model->M = new CoinPackedMatrix(false, 0.0, 0.0);
   CoinAssertHint(model->M, "Error: Out of Memory");
   model->M->setDimensions(0, nCols);
   model->reserve(nRowsMax, nCols);

   CoinPackedVector rowBudget;
   nRows     = 0;
   pairIndex = 0;
   for(vi = pairsAD.begin(); vi != pairsAD.end(); vi++){
      adP = m_instance.getIndexADInv(*vi);
      if(d != adP.second){
	 pairIndex++;
	 continue;      
      }
      a = adP.first;

      //---
      //--- sum{a in A} (f+[a,d] - f-[a,d]) <= B[d]
      //---
      rowBudget.insert(getColOffset_fp() + pairIndex,  1.0);
      rowBudget.insert(getColOffset_fm() + pairIndex, -1.0);


      //---
      //---       f+[a,d] - f-[a,d] =  
      //---         a[a,d] sum{t in T} (t/n) x1[a,t] +
      //---         b[a,d] x2[a]                     +     
      //---         c[a,d] sum{t in T} (t/n) z[a,t]  + 
      //---         d[a,d] x3[a]                     +
      //---         e[a,d]
      //---
      CoinPackedVector row;
      string rowName = "demand_def(a_" 
	 + m_instance.getAtmName(a) + ",d_" 
	 + m_instance.getDateName(d) + ")";
      row.insert(colIndex_fp(pairIndex),  1.0);
      row.insert(colIndex_fm(pairIndex), -1.0);
      coefA = -a_ad[*vi] / nSteps;
      coefC = -c_ad[*vi] / nSteps;
      if(m_appParam.UseTightModel){
	 for(t = 0; t <= nSteps; t++){ 
	    row.insert(colIndex_x1(a,t), t * coefA);
	    row.insert(colIndex_z (a,t), t * coefC);
	 }
      }
      else{
	 for(t = 0; t < nSteps; t++){ 
	    row.insert(colIndex_x1(a,t), (t+1) * coefA);
	    row.insert(colIndex_z (a,t), (t+1) * coefC);
	 }
      }
//.........这里部分代码省略.........

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

//===========================================================================//
int ATM_DecompApp::createConZtoX(DecompConstraintSet * model,
				 const int             atmIndex){

   //---
   //---  for a in A, t in T:
   //---     z[a,t] = x1[a,t] * x2[a] 
   //---         <==>
   //--- OLD:
   //---  for a in A:
   //---     sum{t in T}  x1[a,t] <= 1 (where, T = 1..n)
   //---  for a in A, t in T:
   //---     z[a,t] >= 0 <= 1,                  
   //---     z[a,t] <= x1[a,t],            
   //---     z[a,t] <= x2[a],              
   //---     z[a,t] >= x1[a,t] + x2[a] - 1.
   //--- NEW: 
   //---  for a in A:
   //---     sum{t in T}  x1[a,t] = 1 (where, T = 0..n)
   //---     sum{t in T}   z[a,t] = x2[a]
   //---  for a in A, t in T:
   //---     z[a,t] >= 0 <= 1,                  
   //---     z[a,t] <= x1[a,t].
   //---
   int       t;
   int       nRows = 0;
   const int nSteps = m_appParam.NumSteps;
   if(m_appParam.UseTightModel){
      for(t = 0; t <= nSteps; t++){
	 CoinPackedVector row;
	 string strAT = "(a_" + m_instance.getAtmName(atmIndex) 
	    + ",t_" + UtilIntToStr(t) + ")";
	 string rowName = "ztox1" + strAT;	 
	 row.insert(colIndex_z (atmIndex,t),  1.0);
	 row.insert(colIndex_x1(atmIndex,t), -1.0);
	 model->appendRow(row, -DecompInf, 0.0, rowName);
	 nRows++;
      }
      CoinPackedVector row2;
      string rowName2 = "ztox2(a_" 
	 + m_instance.getAtmName(atmIndex) + ")"; 
      row2.insert(colIndex_x2(atmIndex), -1.0);
      for(t = 0; t <= nSteps; t++){
	 row2.insert(colIndex_z (atmIndex,t),  1.0);
      }
      model->appendRow(row2, 0.0, 0.0, rowName2);	 
      nRows++;         
   }
   else{
      for(t = 0; t < nSteps; t++){
	 CoinPackedVector row1, row2, row3;
	 string strAT = "(a_" + m_instance.getAtmName(atmIndex) 
	    + ",t_" + UtilIntToStr(t+1) + ")";
	 string rowName1 = "ztox1" + strAT;
	 string rowName2 = "ztox2" + strAT;
	 string rowName3 = "ztox3" + strAT;
	 
	 row1.insert(colIndex_z (atmIndex,t),  1.0);
	 row1.insert(colIndex_x1(atmIndex,t), -1.0);
	 model->appendRow(row1, -DecompInf, 0.0, rowName1);
	 
	 row2.insert(colIndex_z (atmIndex,t),  1.0);
	 row2.insert(colIndex_x2(atmIndex)  , -1.0);
	 model->appendRow(row2, -DecompInf, 0.0, rowName2);
	 
	 row3.insert(colIndex_z (atmIndex,t),  1.0);
	 row3.insert(colIndex_x1(atmIndex,t), -1.0);
	 row3.insert(colIndex_x2(atmIndex)  , -1.0);
	 model->appendRow(row3, -1.0, DecompInf, rowName3);
	 
	 nRows+=3;         
      }
   }
   return nRows;
}

// --------------------------------------------------------------------- //
int GAP_DecompApp::createModelPartAP(DecompConstraintSet* model)
{
   int             i, j, colIndex;
   int             status     = GAPStatusOk;
   int             nTasks     = m_instance.getNTasks();    //n
   int             nMachines  = m_instance.getNMachines(); //m
   int             nCols      = nTasks * nMachines;
   int             nRows      = nTasks;
   UtilPrintFuncBegin(m_osLog, m_classTag,
                      "createModelPartAP()", m_appParam.LogLevel, 2);
   //---
   //--- Build the core model constraints (AP = assignment problem).
   //---
   //--- m is number of machines (index i)
   //--- n is number of tasks    (index j)
   //---
   //--- sum{i in 1..m} x[i,j] = 1, j in 1..n
   //---
   //--- Example structure: m=3, n=4
   //---  x   x   x     = 1 [j=1]
   //---   x   x   x    = 1 [j=2]
   //---    x   x   x   = 1 [j=3]
   //---     x   x   x  = 1 [j=4]
   //---
   //---
   //--- Allocate an empty row-ordered CoinPackedMatrix. Since we plan
   //---   to add rows, set the column dimension and let the row dimension
   //---   be set dynamically.
   //---
   model->M = new CoinPackedMatrix(false, 0.0, 0.0);
   CoinAssertHint(model->M, "Error: Out of Memory");
   model->M->setDimensions(0, nCols);
   //---
   //--- we know the sizes needed, so reserve space for them (for efficiency)
   //---
   model->reserve(nRows, nCols);

   //---
   //--- create one row per task
   //---   rowNames are not needed, they are used for debugging
   //---
   for (j = 0; j < nTasks; j++) {
      CoinPackedVector row;
      string           rowName = "a(j_" + UtilIntToStr(j) + ")";

      for (i = 0; i < nMachines; i++) {
         colIndex = getIndexIJ(i, j);
         row.insert(colIndex, 1.0);
      }

      model->appendRow(row, 1.0, 1.0, rowName);
   }

   //---
   //--- set the col upper and lower bounds (all in [0,1])
   //---
   UtilFillN(model->colLB, nCols,  0.0);
   UtilFillN(model->colUB, nCols,  1.0);
   //---
   //--- set the indices of the integer variables of model
   //---   (all vars are binary)
   //---
   UtilIotaN(model->integerVars, nCols, 0);
   //---
   //--- set column names for debugging
   //---
   colIndex = 0;

   for (i = 0; i < nMachines; i++) {
      for (j = 0; j < nTasks; j++) {
         string colName = "x("
                          + UtilIntToStr(colIndex) + "_"
                          + UtilIntToStr(i) + "," + UtilIntToStr(j) + ")";
         model->colNames.push_back(colName);
         colIndex++;
      }
   }

   UtilPrintFuncEnd(m_osLog, m_classTag,
                    "createModelPartAP()", m_appParam.LogLevel, 2);
   return status;
}

bool OSI_X_SolverWrapper::setup(const LP_Constraints & cstraints) //cstraints <-> constraints
{
  if ( si == nullptr )
  {
    return false;
  }
  assert(nbParams_ == cstraints.nbParams_);

  const unsigned int NUMVAR = cstraints.constraint_mat_.cols();
  std::vector<double>
    col_lb(NUMVAR), // the column lower bounds
    col_ub(NUMVAR); // the column upper bounds

  this->nbParams_ = NUMVAR;

  si->setObjSense( ((cstraints.bminimize_) ? 1 : -1) );

  const Mat & A = cstraints.constraint_mat_;

  //Equality constraint will be done by two constraints due to the API limitation ( >= & <=).
  const size_t nbLine = A.rows() +
    std::count(cstraints.vec_sign_.begin(), cstraints.vec_sign_.end(), LP_Constraints::LP_EQUAL);

  // Define default lower and upper bound to [-inf, inf]
  std::vector<double>
    row_lb(nbLine, -si->getInfinity()), // the row lower bounds
    row_ub(nbLine, si->getInfinity()); // the row upper bounds

  std::unique_ptr<CoinPackedMatrix> matrix(new CoinPackedMatrix(false,0,0));
  matrix->setDimensions(0, NUMVAR);

  //-- Add row-wise constraint
  size_t indexRow = 0;
  for (int i=0; i < A.rows(); ++i)
  {
    const Vec temp = A.row(i);
    CoinPackedVector row;

    if (cstraints.vec_sign_[i] == LP_Constraints::LP_EQUAL ||
        cstraints.vec_sign_[i] == LP_Constraints::LP_LESS_OR_EQUAL)
    {
      const int coef = 1;
      for ( int j = 0; j < A.cols(); ++j )
      {
        row.insert(j, coef * temp.data()[j]);
      }
      row_ub[indexRow] = coef * cstraints.constraint_objective_(i);
      matrix->appendRow(row);
      ++indexRow;
    }

    if (cstraints.vec_sign_[i] == LP_Constraints::LP_EQUAL ||
        cstraints.vec_sign_[i] == LP_Constraints::LP_GREATER_OR_EQUAL)
    {
      const int coef = -1;
      for ( int j = 0; j < A.cols(); ++j )
      {
        row.insert(j, coef * temp.data()[j]);
      }
      row_ub[indexRow] = coef * cstraints.constraint_objective_(i);
      matrix->appendRow(row);
      ++indexRow;
    }
  }

  //-- Setup bounds for all the parameters
  if (cstraints.vec_bounds_.size() == 1)
  {
    // Setup the same bound for all the parameters
    std::fill(col_lb.begin(), col_lb.end(), cstraints.vec_bounds_[0].first);
    std::fill(col_ub.begin(), col_ub.end(), cstraints.vec_bounds_[0].second);
  }
  else // each parameter have its own bounds
  {
    for (int i=0; i < this->nbParams_; ++i)
    {
      col_lb[i] = cstraints.vec_bounds_[i].first;
      col_ub[i] = cstraints.vec_bounds_[i].second;
    }
  }

  si->loadProblem(*matrix, &col_lb[0], &col_ub[0], cstraints.vec_cost_.empty() ? nullptr : &cstraints.vec_cost_[0], &row_lb[0], &row_ub[0] );

  return true;
}

//.........这里部分代码省略.........
			"setUpModel", "MibsBilevel");
#endif
     }else{
	throw CoinError("Unknown solver chosen",
			"setUpModel", "MibsBilevel");
     }
	
     int * integerVars = new int[lCols];
     double * objCoeffs = new double[lCols];
     
     CoinFillN(integerVars, lCols, 0);
     //CoinZeroN(objCoeffs, lCols);
     
     int intCnt(0);
  
     /** Fill in array of lower-level integer variables **/
     
     for(i = 0; i < lCols; i++){
	index1 = lColIndices[i];
	if(oSolver->isInteger(index1)){
	   integerVars[intCnt] = i;
	   intCnt++;
	}
     }
     
     CoinDisjointCopyN(lObjCoeffs, lCols, objCoeffs);
     
     CoinPackedMatrix * newMat = new CoinPackedMatrix(false, 0, 0);
     newMat->setDimensions(0, lCols);
     double tmp(0.0);
     
     /*
       for(i = 0; i < lRows; i++){
       CoinPackedVector row;
       index1 = lRowIndices[i];
       start = matStarts[index1];
       end = start + matrix->getVectorSize(index1);
       for(j = start; j < end; j++){
       index2 = matIndices[j];
       //tmp = findIndex(index, lCols, lColIndices);
       tmp = binarySearch(0, lCols - 1, index2, lColIndices);
       if(tmp > -1)
       row.insert(tmp, matElements[j]);
       }
       newMat->appendRow(row);
       }
     */
     for(i = 0; i < lRows; i++){
	CoinPackedVector row;
	index1 = lRowIndices[i];
	for(j = 0; j < lCols; j++){
	   index2 = lColIndices[j];
	   tmp = matrix->getCoefficient(index1, index2);
	   row.insert(j, tmp);
	}
	newMat->appendRow(row);
     }

     /*
       nSolver->assignProblem(newMat, colLb, colUb,
       objCoeffs, rowLb, rowUb);
     */
     
     nSolver->loadProblem(*newMat, colLb, colUb,
			  objCoeffs, rowLb, rowUb);
     

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

//.........这里部分代码省略.........
    model.setDblParam(CbcModel::CbcMaximumSeconds,60.0*minutes);
  }
  // Switch off most output
  if (model.getNumCols()<3000) {
    model.messageHandler()->setLogLevel(1);
    //model.solver()->messageHandler()->setLogLevel(0);
  } else {
    model.messageHandler()->setLogLevel(2);
    model.solver()->messageHandler()->setLogLevel(1);
  }
  // Default strategy will leave cut generators as they exist already
  // so cutsOnlyAtRoot (1) ignored
  // numberStrong (2) is 5 (default)
  // numberBeforeTrust (3) is 5 (default is 0)
  // printLevel (4) defaults (0)
  CbcStrategyDefault strategy(true,5,5);
  // Set up pre-processing to find sos if wanted
  if (preProcess)
    strategy.setupPreProcessing(2);
  model.setStrategy(strategy);

  // Go round adding cuts to cutoff last solution
  // Stop after finding 20 best solutions
  for (int iPass=0;iPass<20;iPass++) {
    time1 = CoinCpuTime();
    // Do complete search
    model.branchAndBound();
    
    std::cout<<mpsFileName<<" took "<<CoinCpuTime()-time1<<" seconds, "
             <<model.getNodeCount()<<" nodes with objective "
             <<model.getObjValue()
             <<(!model.status() ? " Finished" : " Not finished")
             <<std::endl;
    // Stop if infeasible
    if (model.isProvenInfeasible())
      break;
    // Print solution if finished - we can't get names from Osi! - so get from OsiClp
    
    assert (model.getMinimizationObjValue()<1.0e50);
    OsiSolverInterface * solver = model.solver();
    int numberColumns = solver->getNumCols();
    
    const double * solution = model.bestSolution();
    //const double * lower = solver->getColLower();
    //const double * upper = solver->getColUpper();

    // Get names from solver1 (as OsiSolverInterface may lose)
    std::vector<std::string> columnNames = *solver1.getModelPtr()->columnNames();
    
    int iColumn;
    std::cout<<std::setiosflags(std::ios::fixed|std::ios::showpoint)<<std::setw(14);
    
    std::cout<<"--------------------------------------"<<std::endl;
    for (iColumn=0;iColumn<numberColumns;iColumn++) {
      double value=solution[iColumn];
      if (fabs(value)>1.0e-7&&solver->isInteger(iColumn)) 
	std::cout<<std::setw(6)<<iColumn<<" "
                 <<columnNames[iColumn]<<" "
                 <<value
          //<<" "<<lower[iColumn]<<" "<<upper[iColumn]
                 <<std::endl;
    }
    std::cout<<"--------------------------------------"<<std::endl;
  
    std::cout<<std::resetiosflags(std::ios::fixed|std::ios::showpoint|std::ios::scientific);
    /* Now add cut to reference copy.
       resetting to reference copy also gets rid of best solution so we
       should either save best solution, reset, add cut OR
       add cut to reference copy then reset - this is doing latter
    */
    OsiSolverInterface * refSolver = model.referenceSolver();
    const double * bestSolution = model.bestSolution();
    const double * originalLower = refSolver->getColLower();
    const double * originalUpper = refSolver->getColUpper();
    CoinPackedVector cut;
    double rhs = 1.0;
    for (iColumn=0;iColumn<numberColumns;iColumn++) {
      double value=bestSolution[iColumn];
      if (solver->isInteger(iColumn)) {
        // only works for 0-1 variables
        assert (originalLower[iColumn]==0.0&&
                originalUpper[iColumn]==1.0);
        // double check integer
        assert (fabs(floor(value+0.5)-value)<1.0e-5);
        if (value>0.5) {
          // at 1.0
          cut.insert(iColumn,-1.0);
          rhs -= 1.0;
        } else {
          // at 0.0
          cut.insert(iColumn,1.0);
        }
      }
    }
    // now add cut
    refSolver->addRow(cut,rhs,COIN_DBL_MAX);
    model.resetToReferenceSolver();
  }
  return 0;
}    

//--------------------------------------------------------------------------
// ** At present this does not use any solver
void
CglGomoryUnitTest(
  const OsiSolverInterface * baseSiP,
  const std::string mpsDir )
{
  CoinRelFltEq eq(0.000001);

  // Test default constructor
  {
    CglGomory aGenerator;
    assert (aGenerator.getLimit()==50);
    assert (aGenerator.getAway()==0.05);
  }
  
  // Test copy & assignment etc
  {
    CglGomory rhs;
    {
      CglGomory bGenerator;
      bGenerator.setLimit(99);
      bGenerator.setAway(0.2);
      CglGomory cGenerator(bGenerator);
      rhs=bGenerator;
      assert (rhs.getLimit()==99);
      assert (rhs.getAway()==0.2);
    }
  }

  // Test explicit form - all integer (pg 125 Wolsey)
  if (1) {
    OsiCuts osicuts;
    CglGomory test1;
    int i;
    int nOldCuts=0,nRowCuts;
 
    // matrix data
    //deliberate hiccup of 2 between 0 and 1
    CoinBigIndex start[5]={0,4,7,8,9};
    int length[5]={2,3,1,1,1};
    int rows[11]={0,2,-1,-1,0,1,2,0,1,2};
    double elements[11]={7.0,2.0,1.0e10,1.0e10,-2.0,1.0,-2.0,1,1,1};
    CoinPackedMatrix matrix(true,3,5,8,elements,rows,start,length);
    
    // rim data (objective not used just yet)
    double rowLower[5]={14.0,3.0,3.0,1.0e10,1.0e10};
    double rowUpper[5]={14.0,3.0,3.0,-1.0e10,-1.0e10};
    double colLower[7]={0.0,0.0,0.0,0.0,0.0,0.0,0.0};
    double colUpper[7]={100.0,100.0,100.0,100.0,100.0,100.0,100.0};
  
    // integer
    char intVar[7]={2,2,2,2,2,2,2};

    // basis 1
    int rowBasis1[3]={-1,-1,-1};
    int colBasis1[5]={1,1,-1,-1,1};
    CoinWarmStartBasis warm;
    warm.setSize(5,3);
    for (i=0;i<3;i++) {
      if (rowBasis1[i]<0) {
	warm.setArtifStatus(i,CoinWarmStartBasis::atLowerBound);
      } else {
	warm.setArtifStatus(i,CoinWarmStartBasis::basic);
      }
    }
    for (i=0;i<5;i++) {
      if (colBasis1[i]<0) {
	warm.setStructStatus(i,CoinWarmStartBasis::atLowerBound);
      } else {
	warm.setStructStatus(i,CoinWarmStartBasis::basic);
      }
    }

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

//#############################################################################
void 
MibSHeuristic::lowerObjHeuristic()
{

  /* 
     optimize wrt to lower-level objective 
     over current feasible lp feasible region 
  */

  MibSModel * model = MibSModel_;

  OsiSolverInterface * oSolver = model->getSolver();
#ifndef COIN_HAS_SYMPHONY
  OsiSolverInterface * hSolver = new OsiCbcSolverInterface();
#else
  OsiSolverInterface* hSolver = new OsiSymSolverInterface();
#endif
  
  double objSense(model->getLowerObjSense());  
  int lCols(model->getLowerDim());
  int uCols(model->getUpperDim());
  int * lColIndices = model->getLowerColInd();
  int * uColIndices = model->getUpperColInd();
  double * lObjCoeffs = model->getLowerObjCoeffs();
  
  //int tCols(lCols + uCols);
  int tCols(oSolver->getNumCols());
  //assert(tCols == oSolver->getNumCols());

  hSolver->loadProblem(*oSolver->getMatrixByCol(),
		       oSolver->getColLower(), oSolver->getColUpper(),
		       oSolver->getObjCoefficients(),
		       oSolver->getRowLower(), oSolver->getRowUpper());

  int j(0);
  for(j = 0; j < tCols; j++){
    if(oSolver->isInteger(j))
      hSolver->setInteger(j);
  }

  double * nObjCoeffs = new double[tCols];
  int i(0), index(0);
  
  CoinZeroN(nObjCoeffs, tCols);

  for(i = 0; i < lCols; i++){
    index = lColIndices[i];
    nObjCoeffs[index] = lObjCoeffs[i];
  }

  //MibS objective sense is the opposite of OSI's!
  hSolver->setObjSense(objSense);

  hSolver->setObjective(nObjCoeffs);
 
  //double cutoff(model->getCutoff());
  double cutoff(model->getKnowledgeBroker()->getIncumbentValue());

  if(model->getNumSolutions()){
  
    CoinPackedVector objCon;
    //double rhs(cutoff * objSense);
    //double smlTol(1.0);
    double rhs(cutoff);
       
    for(i = 0; i < tCols; i++){
      objCon.insert(i, oSolver->getObjCoefficients()[i] 
		    * oSolver->getObjSense());
    }
    
    hSolver->addRow(objCon, - hSolver->getInfinity(), rhs);
  }
  
  if(0)
     hSolver->writeLp("lobjheurstic");

#ifndef COIN_HAS_SYMPHONY
  dynamic_cast<OsiCbcSolverInterface *> 
     (hSolver)->getModelPtr()->messageHandler()->setLogLevel(0);
#else
  dynamic_cast<OsiSymSolverInterface *> 
     (hSolver)->setSymParam("prep_level", -1);
  
  dynamic_cast<OsiSymSolverInterface *> 
     (hSolver)->setSymParam("verbosity", -2);
  
  dynamic_cast<OsiSymSolverInterface *> 
     (hSolver)->setSymParam("max_active_nodes", 1);
#endif

  hSolver->branchAndBound();

  if(hSolver->isProvenOptimal()){

    double upperObjVal(0.0);

    /*****************NEW ******************/

    MibSSolution *mibSol = NULL;
//.........这里部分代码省略.........

//===========================================================================//
void OSDipApp::createModels() {

	UtilPrintFuncBegin(m_osLog, m_classTag, "createModels()",
			m_appParam.LogLevel, 2);
	int i;
	int j;
	const int nCols = m_osInterface.getVariableNumber();
	const int nRows = m_osInterface.getConstraintNumber();
	
	try{

		//First the define the objective function over the entire variable space
		//Create the memory for the objective  function
		m_objective = new double[nCols];
		for (i = 0; i < nCols; i++) {
			m_objective[i] = m_osInterface.getObjectiveFunctionCoeff()[i];
			//std::cout << "obj coeff = " << m_objective[i] << std::endl;
		}
		setModelObjective( m_objective);
	
		//---
		//--- Construct the core matrix.
		//---
		int nRowsRelax, nRowsCore;
		nRowsRelax = 0;
		nRowsCore = 0;
  
	

		std::vector<OtherConstraintOption*> otherConstraintOptions;
		std::vector<OtherConstraintOption*>::iterator vit;
		
		

		//
		// Now construct the block matrices
		//
		int *rowsRelax;
		int whichBlock;
		DecompConstraintSet *modelRelax = NULL;

		std::set<int> blockVars; //variables indexes in the specific block
		std::set<int> blockVarsAll; //all variable indexes that appear in a block
		std::set<int> blockConAll; //all constraint indexes that appear in a block
		std::set<int>::iterator sit;
		CoinPackedVector *row;
		int *rowVars;
		int rowSize;
		if (m_osInterface.m_osoption != NULL
				&& m_osInterface.m_osoption->getNumberOfOtherConstraintOptions()
						> 0) {
	
			otherConstraintOptions
					= m_osInterface.m_osoption->getOtherConstraintOptions("Dip");
			
			//iterate over the vector of contraint options
			for (vit = otherConstraintOptions.begin(); vit
					!= otherConstraintOptions.end(); vit++) {
	
				// see if we have a Core Constraint Set
	
				if( ( (*vit)->name.compare("constraintSet") == 0)
						&& ( (*vit)->type.compare("Block") == 0)) {
					
					
					//get the block number
	
					//ch = new char[(*vit)->value.size() + 1];
					//ch[(*vit)->value.size()] = 0;
					//memcpy(ch, (*vit)->value.c_str(), (*vit)->value.size());
					//whichBlock = atoi(ch);
					//delete ch;
					whichBlock = atoi( (*vit)->value.c_str() );
			
	
					
					// first get the number of constraints in this block
					nRowsRelax = (*vit)->numberOfCon;
					rowsRelax = new int[nRowsRelax];
	
					//now get the variable indexes for just this block
					//first clear indexes from a previous block
					blockVars.clear();
					for (i = 0; i < nRowsRelax; i++) {
						rowsRelax[i] = (*vit)->con[i]->idx;
						
						if( (*vit)->con[i]->idx >= nRows) throw ErrorClass( "found an invalid row index in OSoL file");
						
						m_blocks[ whichBlock].push_back( rowsRelax[i] );
						
						//also add to the set of all rows
						
						if (blockConAll.find(  (*vit)->con[i]->idx ) == blockConAll.end()) {
							blockConAll.insert(  (*vit)->con[i]->idx );	
						}					
						
						//add the variables for this row to the set blockVars
						row = m_osInterface.getRow(rowsRelax[i]);
						rowSize = row->getNumElements();
//.........这里部分代码省略.........

//.........这里部分代码省略.........

      CouNumber 
	xi    = colsol [ind],
	lb    = lower [ind],
	ub    = upper [ind],
	delta = ub-lb;

      if (fabs (delta) > COUENNE_EPS) {

	CouNumber normlambda = lambda / (delta*delta),
	  coeff = normlambda * (lb + ub - 2. * xi);

	a0 += normlambda * (lb*ub - xi*xi);

	//a0 += coeff * xi - normlambda * (xi - lb) * (ub - xi);
	//a0 += normlambda * lb * ub;
	Qxs [ind] += coeff;
	//Qxs [ind] += normlambda * (lb + ub);
      }// else coeff = 0.;

      //      a0 += lambda [k] * lower  [ind] * upper  [ind];
      //      a0 -= lambda [k] * colsol [ind] * colsol [ind];

      //Qxs [ind] -= lambda [k] * (colsol [ind]) * 2;
    }

  // Count the number of non-zeroes
  int nnz = 0;
  for (int i=0; i < numcols ; i++)
    if (fabs (Qxs [i]) > COUENNE_EPS)
      nnz++;

#ifdef DEBUG
  printf ("2Qx = (");for(int i=0;i<numcols;i++)printf("%g ",Qxs[i]);printf (")[%g], %d nz\n",a0,nnz);
#endif

  // Pack the vector into a CoinPackedVector and generate the cut.
  CoinPackedVector a (false);
  a.reserve (nnz);

#ifdef DEBUG
  CouNumber lhs = 0, lhsc = 0,
    *optimum = cg -> Problem () -> bestSol (),
    *current = cg -> Problem () -> X ();
#endif

  for (int i=0; i < numcols; i++)

    if (fabs (Qxs [i]) > 1.0e-21) { // why 1.0e-21? Look at CoinPackedMatrix.cpp:2188

      // compute violation
#ifdef DEBUG
      if (optimum) {
	printf ("%+g * %g ", Qxs [i], optimum [i]);
	lhs  += Qxs [i] * optimum [i];
      }
      lhsc += Qxs [i] * current [i];
#endif
      a.insert (i, Qxs [i]);
    }

  OsiRowCut cut;
  cut.setRow (a);

  delete [] Qxs;
    
  if (varVal < exprVal) { //(lambda == dCoeffLo_) {

     cut.setUb (a0);

#ifdef DEBUG
     if (optimum && (lhs - a0 > COUENNE_EPS)) {
       printf ("cut violates optimal solution: %g > %g\n", lhs, a0);
       cut.print ();
     }
     if (lhsc < a0 + COUENNE_EPS){
       printf ("cut (+) is not cutting: ");
       cut.print ();
     }
#endif
     //     cut.setLb(-COUENNE_INFINITY);
  }
  else {

    cut.setLb (a0);
#ifdef DEBUG
    if (optimum && (lhs - a0 < -COUENNE_EPS)) {
       printf ("cut violates optimal solution: %g < %g\n", lhs, a0);
       cut.print ();
    }
    if (lhsc > a0 - COUENNE_EPS){
       printf ("cut (-) is not cutting: ");
       cut.print ();
     }
#endif
    //    cut.setUb(COUENNE_INFINITY);
  }

  cs.insert (cut);
}

//===========================================================================//
void MCF_DecompApp::createModelCore(DecompConstraintSet* model)
{
   //---
   //--- MASTER (A''):
   //---      sum{k in K} x[k,i,j] >= l[i,j],       for all (i,j) in A
   //---      sum{k in K} x[k,i,j] <= u[i,j],       for all (i,j) in A
   //---      x[k,i,j] integer >= l[i,j] <= u[i,j], for all (i,j) in A
   //---
   int   k, a, colIndex;
   int   numCommodities = m_instance.m_numCommodities;
   int   numArcs        = m_instance.m_numArcs;
   int   numCols        = numCommodities * numArcs;
   int   numRows        = 2              * numArcs;
   MCF_Instance::arc* arcs = m_instance.m_arcs;
   UtilPrintFuncBegin(m_osLog, m_classTag,
                      "createModelCore()", m_appParam.LogLevel, 2);
   //---
   //--- create space for the model matrix (row-majored)
   //---
   model->M = new CoinPackedMatrix(false, 0.0, 0.0);

   if (!model->M) {
      throw UtilExceptionMemory("createModelCore", "MCF_DecompApp");
   }

   model->M->setDimensions(0, numCols);
   model->reserve(numRows, numCols);
   //---
   //--- create the rows and set the col/row bounds
   //---
   UtilFillN(model->colLB, numCols,  0.0);
   UtilFillN(model->colUB, numCols,  m_infinity);

   for (a = 0; a < numArcs; a++) {
      CoinPackedVector row;
      double           arcLB = arcs[a].lb;
      double           arcUB = arcs[a].ub;

      for (k = 0; k < numCommodities; k++) {
         colIndex = k * numArcs + a;
         model->colLB[colIndex] = arcLB;
         model->colUB[colIndex] = arcUB;
         row.insert(colIndex, 1.0);
      }

      //TODO: any issue with range constraints?
      model->appendRow(row, -m_infinity, arcUB);
      string rowNameUB = "capUB(" +
                         UtilIntToStr(a)            + "_" +
                         UtilIntToStr(arcs[a].tail) + "," +
                         UtilIntToStr(arcs[a].head) + ")";
      model->rowNames.push_back(rowNameUB);
      model->appendRow(row, arcLB,  m_infinity);
      string rowNameLB = "capLB(" +
                         UtilIntToStr(a)            + "_" +
                         UtilIntToStr(arcs[a].tail) + "," +
                         UtilIntToStr(arcs[a].head) + ")";
      model->rowNames.push_back(rowNameLB);
   }

   //---
   //--- create column names (helps with debugging)
   //---
   for (k = 0; k < numCommodities; k++) {
      for (a = 0; a < numArcs; a++) {
         string colName = "x(comm_" + UtilIntToStr(k) + "," +
                          UtilIntToStr(a)            + "_" +
                          UtilIntToStr(arcs[a].tail) + "," +
                          UtilIntToStr(arcs[a].head) + ")";
         model->colNames.push_back(colName);
      }
   }

   //---
   //--- set the indices of the integer variables of model
   //---
   UtilIotaN(model->integerVars, numCols, 0);
   UtilPrintFuncEnd(m_osLog, m_classTag,
                    "createModelCore()", m_appParam.LogLevel, 2);
}

void
CoinShallowPackedVectorUnitTest()
{
  CoinRelFltEq eq;
  int i;
  // Test default constructor
  {
    CoinShallowPackedVector r;
    assert( r.indices_==NULL );
    assert( r.elements_==NULL );
    assert( r.nElements_==0 );
  }

  // Test set and get methods
  const int ne = 4;
  int inx[ne] = { 1, 3, 4, 7 };
  double el[ne] = { 1.2, 3.4, 5.6, 7.8 };
  {
    CoinShallowPackedVector r;    
    assert( r.getNumElements()==0 );
    
    // Test setting/getting elements with int* & double* vectors
    r.setVector( ne, inx, el );
    assert( r.getNumElements()==ne );
    for ( i=0; i<ne; i++ ) {
      assert( r.getIndices()[i]  == inx[i] );
      assert( r.getElements()[i] == el[i]  );
    }
    assert ( r.getMaxIndex()==7 );
    assert ( r.getMinIndex()==1 );

    // try to clear it
    r.clear();
    assert( r.indices_==NULL );
    assert( r.elements_==NULL );
    assert( r.nElements_==0 );

    // Test setting/getting elements with indices out of order  
    const int ne2 = 5;
    int inx2[ne2] = { 2, 4, 8, 14, 3 };
    double el2[ne2] = { 2.2, 4.4, 6.6, 8.8, 3.3 };
 
    r.setVector(ne2,inx2,el2);
    
    assert( r.getNumElements()==ne2 );    
    for (i = 0; i < ne2; ++i) {
       assert( r.getIndices()[i]==inx2[i] );
       assert( r.getElements()[i]==el2[i] );
    }
    
    assert ( r.getMaxIndex()==14 );
    assert ( r.getMinIndex()==2 );
    // try to call it once more
    assert ( r.getMaxIndex()==14 );
    assert ( r.getMinIndex()==2 );

    CoinShallowPackedVector r1(ne2,inx2,el2);
    assert( r == r1 );

    // assignment operator
    r1.clear();
    r1 = r;
    assert( r == r1 );

    // assignment from packed vector
    CoinPackedVector pv1(ne2,inx2,el2);
    r1 = pv1;
    assert( r == r1 );

    // construction
    CoinShallowPackedVector r2(r1);
    assert( r2 == r );
    
    // construction from packed vector
    CoinShallowPackedVector r3(pv1);
    assert( r3 == r );

    // test duplicate indices
    {
      const int ne3 = 4;
      int inx3[ne3] = { 2, 4, 2, 3 };
      double el3[ne3] = { 2.2, 4.4, 8.8, 6.6 };
      r.setVector(ne3,inx3,el3, false);
      assert(r.testForDuplicateIndex() == false);
      bool errorThrown = false;
      try {
        r.setTestForDuplicateIndex(true);
      }
      catch (CoinError& e) {
        errorThrown = true;
      }
      assert( errorThrown );

      r.clear();
      errorThrown = false;
      try {
	 r.setVector(ne3,inx3,el3);
      }
      catch (CoinError& e) {
        errorThrown = true;
//.........这里部分代码省略.........