C++ OsiSolverInterface::getNumRows方法代码示例

//-------------------------------------------------------------------
// Constructor from solver
//-------------------------------------------------------------------
OsiSolverResult::OsiSolverResult(const OsiSolverInterface &solver, const double *lowerBefore,
  const double *upperBefore)
  : objectiveValue_(COIN_DBL_MAX)
  , primalSolution_(NULL)
  , dualSolution_(NULL)
{
  if (solver.isProvenOptimal() && !solver.isDualObjectiveLimitReached()) {
    objectiveValue_ = solver.getObjValue() * solver.getObjSense();
    CoinWarmStartBasis *basis = dynamic_cast< CoinWarmStartBasis * >(solver.getWarmStart());
    assert(basis);
    basis_ = *basis;
    delete basis;
    int numberRows = basis_.getNumArtificial();
    int numberColumns = basis_.getNumStructural();
    assert(numberColumns == solver.getNumCols());
    assert(numberRows == solver.getNumRows());
    primalSolution_ = CoinCopyOfArray(solver.getColSolution(), numberColumns);
    dualSolution_ = CoinCopyOfArray(solver.getRowPrice(), numberRows);
    fixed_.addBranch(-1, numberColumns, lowerBefore, solver.getColLower(),
      upperBefore, solver.getColUpper());
  }
}

// Generate cuts
void
CglFakeClique::generateCuts(const OsiSolverInterface& si, OsiCuts & cs,
			const CglTreeInfo info)
{
  if (fakeSolver_) {
    assert (si.getNumCols()==fakeSolver_->getNumCols());
    fakeSolver_->setColLower(si.getColLower());
    const double * solution = si.getColSolution();
    fakeSolver_->setColSolution(solution);
    fakeSolver_->setColUpper(si.getColUpper());
    // get and set branch and bound cutoff
    double cutoff;
    si.getDblParam(OsiDualObjectiveLimit,cutoff);
    fakeSolver_->setDblParam(OsiDualObjectiveLimit,COIN_DBL_MAX);
#ifdef COIN_HAS_CLP
    OsiClpSolverInterface * clpSolver
      = dynamic_cast<OsiClpSolverInterface *> (fakeSolver_);
    if (clpSolver) {
      // fix up fake solver
      const ClpSimplex * siSimplex = clpSolver->getModelPtr();
      // need to set djs
      memcpy(siSimplex->primalColumnSolution(),
	     si.getReducedCost(),si.getNumCols()*sizeof(double));
      fakeSolver_->setDblParam(OsiDualObjectiveLimit,cutoff);
    }
#endif
    const CoinPackedMatrix * matrixByRow = si.getMatrixByRow();
    const double * elementByRow = matrixByRow->getElements();
    const int * column = matrixByRow->getIndices();
    const CoinBigIndex * rowStart = matrixByRow->getVectorStarts();
    const int * rowLength = matrixByRow->getVectorLengths();
    const double * rowUpper = si.getRowUpper();
    const double * rowLower = si.getRowLower();
    
    // Scan all rows looking for possibles
    int numberRows = si.getNumRows();
    double tolerance = 1.0e-3;
    for (int iRow=0;iRow<numberRows;iRow++) {
      CoinBigIndex start = rowStart[iRow];
      CoinBigIndex end = start + rowLength[iRow];
      double upRhs = rowUpper[iRow]; 
      double loRhs = rowLower[iRow]; 
      double sum = 0.0;
      for (CoinBigIndex j=start;j<end;j++) {
	int iColumn=column[j];
	double value = elementByRow[j];
	sum += solution[iColumn]*value;
      }
      if (sum<loRhs-tolerance||sum>upRhs+tolerance) {
	// add as cut
	OsiRowCut rc;
	rc.setLb(loRhs);
	rc.setUb(upRhs);
	rc.setRow(end-start,column+start,elementByRow+start,false);
	CoinAbsFltEq equal(1.0e-12);
	cs.insertIfNotDuplicate(rc,equal);
      }
    }
    CglClique::generateCuts(*fakeSolver_,cs,info);
    if (probing_) {
      probing_->generateCuts(*fakeSolver_,cs,info);
    }
  } else {
    // just use real solver
    CglClique::generateCuts(si,cs,info);
  }
}

// Create a list of rows which might yield cuts
// The possible parameter is a list to cut down search
void CglOddHole::createRowList( const OsiSolverInterface & si,
		      const int * possible)
{
  // Get basic problem information
  int nRows=si.getNumRows(); 
  
  const CoinPackedMatrix * rowCopy = si.getMatrixByRow();

  const int * column = rowCopy->getIndices();
  const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
  const int * rowLength = rowCopy->getVectorLengths(); 
  
  int rowIndex;
  delete [] suitableRows_;
  numberRows_=nRows;

  const double * rowElements = rowCopy->getElements();
  const double * rowupper = si.getRowUpper();
  const double * rowlower = si.getRowLower();
  const double * collower = si.getColLower();
  const double * colupper = si.getColUpper();

  suitableRows_=new int[nRows];
  if (possible) {
    memcpy(suitableRows_,possible,nRows*sizeof(int));
  } else {
    int i;
    for (i=0;i<nRows;i++) {
      suitableRows_[i]=1;
    }
  }
  for (rowIndex=0; rowIndex<nRows; rowIndex++){
    double rhs1=rowupper[rowIndex];
    double rhs2=rowlower[rowIndex];
    if (suitableRows_[rowIndex]) {
      int i;
      bool goodRow=true;
      for (i=rowStart[rowIndex];
	   i<rowStart[rowIndex]+rowLength[rowIndex];i++) {
	int thisCol=column[i];
	if (colupper[thisCol]-collower[thisCol]>epsilon_) {
	  // could allow general integer variables but unlikely
	  if (!si.isBinary(thisCol) ) {
	    goodRow=false;
	    break;
	  }
	  if (fabs(rowElements[i]-1.0)>epsilon_) {
	    goodRow=false;
	    break;
	  }
	} else {
	  rhs1 -= collower[thisCol]*rowElements[i];
	  rhs2 -= collower[thisCol]*rowElements[i];
	}
      }
      if (fabs(rhs1-1.0)>epsilon_&&fabs(rhs2-1.0)>epsilon_) {
	goodRow=false;
      }
      if (goodRow) {
	suitableRows_[rowIndex]=1;
      } else {
	suitableRows_[rowIndex]=0;
      }
    }
  }
}

//-------------------------------------------------------------------------------
// Generate three cycle cuts
//------------------------------------------------------------------- 
void CglOddHole::generateCuts(const OsiSolverInterface & si, OsiCuts & cs,
			      const CglTreeInfo info)
{
  // Get basic problem information
  int nRows=si.getNumRows(); 
  int nCols=si.getNumCols(); 
  
  const CoinPackedMatrix * rowCopy = si.getMatrixByRow();

  // Could do cliques and extra OSL cliques
  // For moment just easy ones
  
  // If no information exists then get a list of suitable rows
  // If it does then suitable rows are subset of information
  
  CglOddHole temp;
  int * checkRow = new int[nRows];
  int i;
  if (!suitableRows_) {
    for (i=0;i<nRows;i++) {
      checkRow[i]=1;
    }
  } else {
    // initialize and extend rows to current size
    memset(checkRow,0,nRows*sizeof(int));
    memcpy(checkRow,suitableRows_,CoinMin(nRows,numberRows_)*sizeof(int));
  }
  temp.createRowList(si,checkRow);
  // now cut down further by only allowing rows with fractional solution
  double * solution = new double[nCols];
  memcpy(solution,si.getColSolution(),nCols*sizeof(double));
  const int * column = rowCopy->getIndices();
  const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
  const int * rowLength = rowCopy->getVectorLengths(); 
  const double * collower = si.getColLower();
  const double * colupper = si.getColUpper();
  int * suitable = temp.suitableRows_;

  // At present I am using new and delete as easier to see arrays in debugger
  int * fixed = new int[nCols]; // mark fixed columns 

  for (i=0;i<nCols;i++) {
    if (si.isBinary(i) ) {
      fixed[i]=0;
      if (colupper[i]-collower[i]<epsilon_) {
	solution[i]=0.0;
	fixed[i]=2;
      } else if (solution[i]<epsilon_) {
	solution[i]=0.0;
	fixed[i]=-1;
      } else if (solution[i]>onetol_) {
	solution[i]=1.0;
	fixed[i]=+1;
      }
    } else {
      //mark as fixed even if not (can not intersect any interesting rows)
      solution[i]=0.0;
      fixed[i]=3;
    }
  }
  // first do packed
  const double * rowlower = si.getRowLower();
  const double * rowupper = si.getRowUpper();
  for (i=0;i<nRows;i++) {
    if (suitable[i]) {
      int k;
      double sum=0.0;
      if (rowupper[i]>1.001) suitable[i]=-1;
      for (k=rowStart[i]; k<rowStart[i]+rowLength[i];k++) {
	int icol=column[k];
	if (!fixed[icol]) sum += solution[icol];
      }
      if (sum<0.9) suitable[i]=-1; //say no good
    }
  }
#ifdef CGL_DEBUG
  const OsiRowCutDebugger * debugger = si.getRowCutDebugger();
  if (debugger&&!debugger->onOptimalPath(si))
    debugger = NULL;
#else
  const OsiRowCutDebugger * debugger = NULL;
#endif
  temp.generateCuts(debugger, *rowCopy,solution,
		    si.getReducedCost(),cs,suitable,fixed,info,true);
  // now cover
  //if no >= then skip
  bool doCover=false;
  int nsuitable=0;
  for (i=0;i<nRows;i++) {
    suitable[i]=abs(suitable[i]);
    if (suitable[i]) {
      int k;
      double sum=0.0;
      if (rowlower[i]<0.999) sum=2.0;
      if (rowupper[i]>1.001) doCover=true;
      for (k=rowStart[i]; k<rowStart[i]+rowLength[i];k++) {
	int icol=column[k];
//.........这里部分代码省略.........

//.........这里部分代码省略.........
	SmiCoreData *smiCore = new SmiCoreData(&ocsi,nstg,cstg,rstg);

	cout << "ModelDiscrete: generated Core data" << endl;

	// Create discrete distribution
	SmiDiscreteDistribution *smiDD = new SmiDiscreteDistribution(smiCore);

	cout << "ModelDiscrete: adding random variables" << endl;

	int index=0;
	for (jj=0;jj<nindp;jj++)
	{
		SmiDiscreteRV *smiRV = new SmiDiscreteRV(1);
		for (ii=0;ii<nsamp[jj];ii++)
		{
			CoinPackedVector empty_vec;
			CoinPackedMatrix empty_mat;
			CoinPackedVector cpv_rlo ;
			CoinPackedVector cpv_rup ;
			cpv_rlo.insert(n_first_stg_rows + jj, demand[index+ii]);
			cpv_rup.insert(n_first_stg_rows + jj, demand[index+ii]);
			smiRV->addEvent(empty_mat,empty_vec,empty_vec,empty_vec,cpv_rlo,cpv_rup,dprobs[index+ii]);
			cpv_rlo.clear();
			cpv_rup.clear();
		}
		smiDD->addDiscreteRV(smiRV);
		index+=nsamp[jj];
	}

	assert(smiDD->getNumRV() == nindp);
	cout << "ModelDiscrete: added " << nindp << " random variables" << endl;



	cout << "ModelDiscrete: processing into scenarios" << endl;

	smiModel->processDiscreteDistributionIntoScenarios(smiDD);

	// load problem data into OsiSolver
	smiModel->loadOsiSolverData();
	// get Osi pointer
	OsiSolverInterface *smiOsi = smiModel->getOsiSolverInterface();
	// set some parameters
	smiOsi->setHintParam(OsiDoPresolveInInitial,true);
	smiOsi->setHintParam(OsiDoScale,true);
	smiOsi->setHintParam(OsiDoCrash,true);
	// solve using Osi Solver
	smiOsi->initialSolve();
	// test optimal value
    	assert(fabs(smiOsi->getObjValue()-1566.042)<0.01);

	// test solutions
	const double *dsoln = smiOsi->getColSolution();
	double objSum = 0.0;

	/* The canonical way to traverse the tree:
	   For each scenario, get the leaf node.
	   Then get the parent.  Repeat until parent is NULL.
	   (Only the root node has a NULL parent.)
	 */

	for(int is=0; is<smiModel->getNumScenarios(); ++is)
	{
		/* this loop calculates the scenario objective value */
		double scenSum = 0.0;

		// start with leaf node
		SmiScnNode *node = smiModel->getLeafNode(is);

		// leaf node probability is the scenario probability
		double scenprob = node->getModelProb();

		while (node != NULL)
		{
			// getColStart returns the starting index of node in OSI model
			for(int j=node->getColStart(); j<node->getColStart()+node->getNumCols(); ++j)
			{
				// getCoreColIndex returns the corresponding Core index
				// in the original (user's) ordering
				scenSum += dobj[node->getCoreColIndex(j)]*dsoln[j];


			}

			// get parent of node
			node = node->getParent();
		}
		objSum += scenSum*scenprob;
	}

	assert(fabs(smiOsi->getObjValue()-objSum) < 0.01);

		// print results
		printf("Solved stochastic program %s\n", name);
		printf("Number of rows: %d\n",smiOsi->getNumRows());
		printf("Number of cols: %d\n",smiOsi->getNumCols());
		printf("Optimal value: %g\n",smiOsi->getObjValue());


}

int * analyze(OsiClpSolverInterface * solverMod, int & numberChanged,
		     double & increment, bool changeInt,
		     CoinMessageHandler * generalMessageHandler, bool noPrinting)
{
    bool noPrinting_ = noPrinting;
    OsiSolverInterface * solver = solverMod->clone();
    char generalPrint[200];
    if (0) {
        // just get increment
        CbcModel model(*solver);
        model.analyzeObjective();
        double increment2 = model.getCutoffIncrement();
        printf("initial cutoff increment %g\n", increment2);
    }
    const double *objective = solver->getObjCoefficients() ;
    const double *lower = solver->getColLower() ;
    const double *upper = solver->getColUpper() ;
    int numberColumns = solver->getNumCols() ;
    int numberRows = solver->getNumRows();
    double direction = solver->getObjSense();
    int iRow, iColumn;

    // Row copy
    CoinPackedMatrix matrixByRow(*solver->getMatrixByRow());
    const double * elementByRow = matrixByRow.getElements();
    const int * column = matrixByRow.getIndices();
    const CoinBigIndex * rowStart = matrixByRow.getVectorStarts();
    const int * rowLength = matrixByRow.getVectorLengths();

    // Column copy
    CoinPackedMatrix  matrixByCol(*solver->getMatrixByCol());
    const double * element = matrixByCol.getElements();
    const int * row = matrixByCol.getIndices();
    const CoinBigIndex * columnStart = matrixByCol.getVectorStarts();
    const int * columnLength = matrixByCol.getVectorLengths();

    const double * rowLower = solver->getRowLower();
    const double * rowUpper = solver->getRowUpper();

    char * ignore = new char [numberRows];
    int * changed = new int[numberColumns];
    int * which = new int[numberRows];
    double * changeRhs = new double[numberRows];
    memset(changeRhs, 0, numberRows*sizeof(double));
    memset(ignore, 0, numberRows);
    numberChanged = 0;
    int numberInteger = 0;
    for (iColumn = 0; iColumn < numberColumns; iColumn++) {
        if (upper[iColumn] > lower[iColumn] + 1.0e-8 && solver->isInteger(iColumn))
            numberInteger++;
    }
    bool finished = false;
    while (!finished) {
        int saveNumberChanged = numberChanged;
        for (iRow = 0; iRow < numberRows; iRow++) {
            int numberContinuous = 0;
            double value1 = 0.0, value2 = 0.0;
            bool allIntegerCoeff = true;
            double sumFixed = 0.0;
            int jColumn1 = -1, jColumn2 = -1;
            for (CoinBigIndex j = rowStart[iRow]; j < rowStart[iRow] + rowLength[iRow]; j++) {
                int jColumn = column[j];
                double value = elementByRow[j];
                if (upper[jColumn] > lower[jColumn] + 1.0e-8) {
                    if (!solver->isInteger(jColumn)) {
                        if (numberContinuous == 0) {
                            jColumn1 = jColumn;
                            value1 = value;
                        } else {
                            jColumn2 = jColumn;
                            value2 = value;
                        }
                        numberContinuous++;
                    } else {
                        if (fabs(value - floor(value + 0.5)) > 1.0e-12)
                            allIntegerCoeff = false;
                    }
                } else {
                    sumFixed += lower[jColumn] * value;
                }
            }
            double low = rowLower[iRow];
            if (low > -1.0e20) {
                low -= sumFixed;
                if (fabs(low - floor(low + 0.5)) > 1.0e-12)
                    allIntegerCoeff = false;
            }
            double up = rowUpper[iRow];
            if (up < 1.0e20) {
                up -= sumFixed;
                if (fabs(up - floor(up + 0.5)) > 1.0e-12)
                    allIntegerCoeff = false;
            }
            if (!allIntegerCoeff)
                continue; // can't do
            if (numberContinuous == 1) {
                // see if really integer
                // This does not allow for complicated cases
                if (low == up) {
                    if (fabs(value1) > 1.0e-3) {
//.........这里部分代码省略.........

void
CglLandP::CachedData::getData(const OsiSolverInterface &si)
{
    int nBasics = si.getNumRows();
    int nNonBasics = si.getNumCols();
    if (basis_ != NULL)
        delete basis_;
    basis_ = dynamic_cast<CoinWarmStartBasis *> (si.getWarmStart());
    if (!basis_)
        throw NoBasisError();

    if (nBasics_ > 0 || nBasics != nBasics_)
    {
        delete [] basics_;
        basics_ = NULL;
    }
    if (basics_ == NULL)
    {
        basics_ = new int[nBasics];
        nBasics_ = nBasics;
    }

    if (nNonBasics_ > 0 || nNonBasics != nNonBasics_)
    {
        delete [] nonBasics_;
        nonBasics_ = NULL;
    }
    if (nonBasics_ == NULL)
    {
        nonBasics_ = new int[nNonBasics];
        nNonBasics_ = nNonBasics;
    }
    int n = nBasics + nNonBasics;
    if ( nBasics_ + nNonBasics_ > 0 || nBasics_ + nNonBasics_ != n)
    {
        delete [] colsol_;
        delete [] integers_;
        integers_ = NULL;
        colsol_ = NULL;
        slacks_ = NULL;
    }
    if (colsol_ == NULL)
    {
        colsol_ = new double[n];
        slacks_ = &colsol_[nNonBasics];
    }

    if (integers_ == NULL)
    {
        integers_ = new bool[n];
    }

    const double * rowLower = si.getRowLower();
    const double * rowUpper = si.getRowUpper();
    //determine which slacks are integer
    const CoinPackedMatrix * m = si.getMatrixByCol();
    const double * elems = m->getElements();
    const int * inds = m->getIndices();
    const CoinBigIndex * starts = m->getVectorStarts();
    const int * lengths = m->getVectorLengths();
    //    int numElems = m->getNumElements();
    int numCols = m->getNumCols();
    assert(numCols == nNonBasics_);
    //   int numRows = m->getNumRows();
    CoinFillN(integers_ ,n, true);
    for (int i = 0 ;  i < numCols ; i++)
    {
        if (si.isContinuous(i))
            integers_[i] = false;
    }
    bool * integerSlacks = integers_ + numCols;
    for (int i = 0 ; i < nBasics ; i++)
    {
        if (rowLower[i] > -1e50 && INT_INFEAS(rowLower[i]) > 1e-15)
            integerSlacks[i] = false;
        if (rowUpper[i] < 1e50 && INT_INFEAS(rowUpper[i]) > 1e-15)
            integerSlacks[i] = false;
    }
    for (int i = 0 ;  i < numCols ; i++)
    {
        CoinBigIndex end = starts[i] + lengths[i];
        if (integers_[i])
        {
            for (CoinBigIndex k=starts[i] ; k < end; k++)
            {
                if (integerSlacks[inds[k]] && INT_INFEAS(elems[k])>1e-15 )
                    integerSlacks[inds[k]] = false;
            }
        }
        else
        {
            for (CoinBigIndex k=starts[i] ; k < end; k++)
            {
                if (integerSlacks[inds[k]])
                    integerSlacks[inds[k]] = false;
            }
        }
    }

    CoinCopyN(si.getColSolution(), si.getNumCols(), colsol_);
//.........这里部分代码省略.........

void
CglLandP::generateCuts(const OsiSolverInterface & si, OsiCuts & cs,
                       const CglTreeInfo info )
{
    if ((info.pass == 0) && !info.inTree)
    {
        numrows_ = si.getNumRows();
    }
// scanExtraCuts(cs, si.getColSolution());
    Parameters params = params_;
    params.rhsWeight = numrows_ + 2;

    handler_->message(CUT_GAP, messages_)<<info.pass<<si.getObjValue() <<CoinMessageEol;

    if (info.inTree)   //put lower pivot limit
    {
        params.pivotLimit = std::min(params.pivotLimit, params.pivotLimitInTree);
        params.countMistakenRc = true;
    }
    if (params.timeLimit < 0)
    {
        params.pivotLimit = 0;
    }

    assert(si.basisIsAvailable());


#ifdef APPEND_ROW
    OsiSolverInterface * t_si = si.clone();
    if (params.modularize)
    {
        int new_idx = si.getNumCols();
        int v_idx[1] = {new_idx};
        double v_val[1] = {-1};
        CoinPackedVector v(1, v_idx, v_val, false);
        t_si->addCol(CoinPackedVector(), 0, 1, 0);
        t_si->setInteger(new_idx);
        t_si->addRow(v,0, 0);
        t_si->resolve();
    }
#else
    const OsiSolverInterface * t_si = &si;
#endif

    cached_.getData(*t_si);
    CglLandPSimplex landpSi(*t_si, cached_, params, validator_);
    if (params.generateExtraCuts == CglLandP::AllViolatedMigs)
    {
        landpSi.genThisBasisMigs(cached_, params);
    }
    landpSi.setLogLevel(handler_->logLevel());
    int nCut = 0;

    std::vector<int> indices;
    getSortedFractionalIndices(indices,cached_, params);

#ifndef NDEBUG
    int numrows = si.getNumRows();
#endif

#ifdef DO_STAT
    //Get informations on current optimum
    {
        OsiSolverInterface * gapTester = si.clone();
        gapTester->resolve();

        roundsStats_.analyseOptimalBasis(gapTester,info.pass, numrows_);
        delete gapTester;
    }
#endif

    params_.timeLimit += CoinCpuTime();
    CoinRelFltEq eq(1e-04);

    for (unsigned int i = 0; i < indices.size() && nCut < params.maxCutPerRound &&
            nCut < cached_.nBasics_ ; i++)
    {

        //Check for time limit
        int iRow = indices[i];
        assert(iRow < numrows);
        OsiRowCut cut;
        int code=1;
        OsiSolverInterface * ncSi = NULL;

        if (params.pivotLimit != 0)
        {
            ncSi = t_si->clone();
            landpSi.setSi(ncSi);
            ncSi->setDblParam(OsiDualObjectiveLimit, COIN_DBL_MAX);
            ncSi->messageHandler()->setLogLevel(0);
        }

        int generated = 0;
        if (params.pivotLimit == 0)
        {
            generated = landpSi.generateMig(iRow, cut, params);
        }
        else
        {
//.........这里部分代码省略.........

//--------------------------------------------------------------------------
// test EKKsolution methods.
void
CglOddHoleUnitTest(
  const OsiSolverInterface * baseSiP,
  const std::string mpsDir )
{
  CoinRelFltEq eq(0.000001);

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


  // test on simple case
  {  
    const int nRows=3;
    const int nCols=3;
    const int nEls=6;
    const double elem[]={1.0,1.0,1.0,1.0,1.0,1.0};
    const int row[]={0,1,0,2,1,2};
    const CoinBigIndex start[]={0,2,4};
    const int len[]={2,2,2};
    CoinPackedMatrix matrix(true,nRows,nCols,nEls,elem,row,start,len);
    const double sol[]={0.5,0.5,0.5};
    const double dj[]={0,0,0};
    const int which[]={1,1,1};
    const int fixed[]={0,0,0};
    OsiCuts cs;
    CglOddHole test1;
    CglTreeInfo info;
    info.randomNumberGenerator=NULL;
    test1.generateCuts(NULL,matrix,sol,dj,cs,which,fixed,info,true);
    CoinPackedVector check;
    int index[] = {0,1,2};
    double el[] = {1,1,1};
    check.setVector(3,index,el);
    //assert (cs.sizeRowCuts()==2);
    assert (cs.sizeRowCuts()==1);
    // sort Elements in increasing order
    CoinPackedVector rpv=cs.rowCut(0).row();
    rpv.sortIncrIndex();
    assert (check==rpv);
  }
  
  // Testcase /u/rlh/osl2/mps/scOneInt.mps
  // Model has 3 continous, 2 binary, and 1 general
  // integer variable.
  {
    OsiSolverInterface  * siP = baseSiP->clone();
    
    std::string fn = mpsDir+"scOneInt";
    siP->readMps(fn.c_str(),"mps");
#if 0
    CglOddHole cg;
    int nCols=siP->getNumCols();
    
    // Test the siP methods for detecting
    // variable type
    int numCont=0, numBinary=0, numIntNonBinary=0, numInt=0;
    for (int thisCol=0; thisCol<nCols; thisCol++) {
      if ( siP->isContinuous(thisCol) ) numCont++;
      if ( siP->isBinary(thisCol) ) numBinary++;
      if ( siP->isIntegerNonBinary(thisCol) ) numIntNonBinary++;
      if ( siP->isInteger(thisCol) ) numInt++;
    }
    assert(numCont==3);
    assert(numBinary==2);
    assert(numIntNonBinary==1);
    assert(numInt==3);
    
    
    // Test initializeCutGenerator
    siP->initialSolve();
    assert(xstar !=NULL);
    for (i=0; i<nCols; i++){
      assert(complement[i]==0);
    }
    int nRows=siP->getNumRows();
    for (i=0; i<nRows; i++){
    int vectorsize = siP->getMatrixByRow()->getVectorSize(i);
    assert(vectorsize==2);
    }
    
    kccg.cleanUpCutGenerator(complement,xstar);
#endif  
    delete siP;
  }

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

//-------------------------------------------------------------------
// Determine row types. Find the VUBS and VLBS.
//-------------------------------------------------------------------
void
CglFlowCover::flowPreprocess(const OsiSolverInterface& si) const
{
  CoinPackedMatrix matrixByRow(*si.getMatrixByRow());

  int numRows = si.getNumRows();
  int numCols = si.getNumCols();

  const char* sense        = si.getRowSense();
  const double* RHS        = si.getRightHandSide();

  const double* coefByRow  = matrixByRow.getElements();
  const int* colInds       = matrixByRow.getIndices();
  const int* rowStarts     = matrixByRow.getVectorStarts();
  const int* rowLengths    = matrixByRow.getVectorLengths();
  int iRow      = -1;
  int iCol      = -1;

  numCols_ = numCols;     // Record col and row numbers for copy constructor
  numRows_ = numRows;

  if (rowTypes_ != 0) {
    delete [] rowTypes_; rowTypes_ = 0;
  }
  rowTypes_ = new CglFlowRowType [numRows];// Destructor will free memory
  // Get integer types
  const char * columnType = si.getColType (true);

  // Summarize the row type infomation.
  int numUNDEFINED   = 0;
  int numVARUB       = 0;
  int numVARLB       = 0;
  int numVAREQ       = 0;
  int numMIXUB       = 0;
  int numMIXEQ       = 0;
  int numNOBINUB     = 0;
  int numNOBINEQ     = 0;
  int numSUMVARUB    = 0;
  int numSUMVAREQ    = 0;
  int numUNINTERSTED = 0;

  int* ind     = new int [numCols];
  double* coef = new double [numCols];
  for (iRow = 0; iRow < numRows; ++iRow) {
    int rowLen   = rowLengths[iRow];
    char sen     = sense[iRow];
    double rhs   = RHS[iRow];

    CoinDisjointCopyN(colInds + rowStarts[iRow], rowLen, ind);
    CoinDisjointCopyN(coefByRow + rowStarts[iRow], rowLen, coef);

    CglFlowRowType rowType = determineOneRowType(si, rowLen, ind, coef,
						 sen, rhs);

    rowTypes_[iRow] = rowType;

    switch(rowType) {
    case  CGLFLOW_ROW_UNDEFINED:
      ++numUNDEFINED;
      break;
    case  CGLFLOW_ROW_VARUB:
      ++numVARUB;
      break;
    case  CGLFLOW_ROW_VARLB:
      ++numVARLB;
      break;
    case  CGLFLOW_ROW_VAREQ:
      ++numVAREQ;
      break;
    case  CGLFLOW_ROW_MIXUB:
      ++numMIXUB;
      break;
    case  CGLFLOW_ROW_MIXEQ:
      ++numMIXEQ;
      break;
    case  CGLFLOW_ROW_NOBINUB:
      ++numNOBINUB;
      break;
    case  CGLFLOW_ROW_NOBINEQ:
      ++numNOBINEQ;
      break;
    case  CGLFLOW_ROW_SUMVARUB:
      ++numSUMVARUB;
      break;
    case  CGLFLOW_ROW_SUMVAREQ:
      ++numSUMVAREQ;
      break;
    case  CGLFLOW_ROW_UNINTERSTED:
      ++numUNINTERSTED;
      break;
    default:
      throw CoinError("Unknown row type", "flowPreprocess",
		      "CglFlowCover");
    }

  }
  delete [] ind;  ind  = NULL;
//.........这里部分代码省略.........

/*===========================================================================*
  Create the set packing submatrix
 *===========================================================================*/
void
CglClique::createSetPackingSubMatrix(const OsiSolverInterface& si)
{
   sp_col_start = new int[sp_numcols+1];
   sp_row_start = new int[sp_numrows+1];
   std::fill(sp_col_start, sp_col_start + (sp_numcols+1), 0);
   std::fill(sp_row_start, sp_row_start + (sp_numrows+1), 0);

   int i, j;

   const CoinPackedMatrix& mcol = *si.getMatrixByCol();
   const int numrows = si.getNumRows();
   int* clique = new int[numrows];
   std::fill(clique, clique+numrows, -1);
   for (i = 0; i < sp_numrows; ++i)
      clique[sp_orig_row_ind[i]] = i;

   for (j = 0; j < sp_numcols; ++j) {
      const CoinShallowPackedVector& vec = mcol.getVector(sp_orig_col_ind[j]);
      const int* ind = vec.getIndices();
      for (i = vec.getNumElements() - 1; i >= 0; --i) {
	 if (clique[ind[i]] >= 0) {
	    ++sp_col_start[j];
	    ++sp_row_start[clique[ind[i]]];
	 }
      }
   }

   std::partial_sum(sp_col_start, sp_col_start+sp_numcols, sp_col_start);
   std::rotate(sp_col_start, sp_col_start+sp_numcols,
	       sp_col_start + (sp_numcols+1));
   std::partial_sum(sp_row_start, sp_row_start+sp_numrows, sp_row_start);
   std::rotate(sp_row_start, sp_row_start+sp_numrows,
	       sp_row_start + (sp_numrows+1));
   const int nzcnt = sp_col_start[sp_numcols];
   assert(nzcnt == sp_row_start[sp_numrows]);
/*
  Now create the vectors with row indices for each column (sp_col_ind) and
  column indices for each row (sp_row_ind). It turns out that
  CoinIsOrthogonal assumes that the row indices for a given column are listed
  in ascending order. This is *not* a solver-independent assumption! At best,
  one can hope that the underlying solver will produce an index vector that's
  either ascending or descending. Under that assumption, compare the first
  and last entries and proceed accordingly. Eventually some solver will come
  along that hands back an index vector in random order, and CoinIsOrthogonal
  will break.  Until then, try and avoid the cost of a sort.
*/
   sp_col_ind = new int[nzcnt];
   sp_row_ind = new int[nzcnt];
   int last=0;
   for (j = 0; j < sp_numcols; ++j) {
      const CoinShallowPackedVector& vec = mcol.getVector(sp_orig_col_ind[j]);
      const int len = vec.getNumElements();
      const int* ind = vec.getIndices();
      if (ind[0] < ind[len-1]) {
	for (i = 0; i < len; ++i) {
	   const int sp_row = clique[ind[i]];
	   if (sp_row >= 0) {
	      sp_col_ind[sp_col_start[j]++] = sp_row;
	      sp_row_ind[sp_row_start[sp_row]++] = j;
	   }
	}
      }
      else {
	for (i = len-1; i >= 0; --i) {
	   const int sp_row = clique[ind[i]];
	   if (sp_row >= 0) {
	      sp_col_ind[sp_col_start[j]++] = sp_row;
	      sp_row_ind[sp_row_start[sp_row]++] = j;
	   }
	}
      }
      // sort
      std::sort(sp_col_ind+last,sp_col_ind+sp_col_start[j]);
      last=sp_col_start[j];
   }
   std::rotate(sp_col_start, sp_col_start+sp_numcols,
	       sp_col_start + (sp_numcols+1));
   sp_col_start[0] = 0;
   std::rotate(sp_row_start, sp_row_start+sp_numrows,
	       sp_row_start + (sp_numrows+1));
   sp_row_start[0] = 0;

   delete[] clique;
}

//Solver function
int sci_rmps(char *fname) 
{
    //creating a problem pointer using base class of OsiSolverInterface and
    //instantiate the object using derived class of ClpSolverInterface
    OsiSolverInterface* si = new OsiClpSolverInterface();

    // Error management variable
	SciErr sciErr;

	//data declarations
	int *piAddressVarOne = NULL;                 //pointer used to access argument of the function
	char* ptr;                              	 //pointer to point to address of file name
    double* options_;                            //options to set maximum iterations 
	CheckInputArgument(pvApiCtx, 2,2 );          //Check we have exactly two arguments as input or not
	CheckOutputArgument(pvApiCtx, 6, 6);         //Check we have exactly six arguments on output side or not
    //Getting the input arguments from Scilab
    //Getting the MPS file path
	//Reading mps file
	getStringFromScilab(1,&ptr);

 	std::cout<<ptr;
	
    //get options from Scilab
    if(getFixedSizeDoubleMatrixInList(2 , 2 , 1 , 1 , &options_))
	{
		return 1;
	}

    //Read the MPS file
    si->readMps(ptr);

    //setting options for maximum iterations
    si->setIntParam(OsiMaxNumIteration,options_[0]);

    //Solve the problem
    si->initialSolve();
  
    //Quering about the problem
    //get number of variables
    double numVars_;
    numVars_ = si->getNumCols();
  
    //get number of constraint equations
    double numCons_;
    numCons_ = si->getNumRows();
   
    //Output the solution to Scilab
    //get solution for x
    const double* xValue = si->getColSolution();
   
    //get objective value
    double objValue = si->getObjValue();

    //get Status value
    double status;
    if(si->isProvenOptimal())
    	status=0;
    else if(si->isProvenPrimalInfeasible())
    	status=1;
    else if(si->isProvenDualInfeasible())
        status=2;
    else if(si->isIterationLimitReached())
        status=3;
   	else if(si->isAbandoned())
        status=4;
   	else if(si->isPrimalObjectiveLimitReached())
        status=5;
   	else if(si->isDualObjectiveLimitReached())
        status=6;

    //get number of iterations
    double iterations = si->getIterationCount();

    //get reduced cost 
    const double* reducedCost = si->getReducedCost();
   
    //get dual vector
    const double* dual = si->getRowPrice();
  
    returnDoubleMatrixToScilab(1 , 1 , numVars_ , xValue);
    returnDoubleMatrixToScilab(2 , 1 , 1 , &objValue);
    returnDoubleMatrixToScilab(3 , 1 , 1 , &status);
    returnDoubleMatrixToScilab(4 , 1 , 1 , &iterations);
    returnDoubleMatrixToScilab(5 , 1 , numVars_ , reducedCost);
    returnDoubleMatrixToScilab(6 , 1 , numCons_ , dual);
	
	free(xValue);
	free(dual);
	free(reducedCost);
}

//.........这里部分代码省略.........
        Increment the jj'th random variable by incr[jj]
		and generate new sample data.
    ***** */

    /* sample space increment initialized to 1 */
    incr = (int *) malloc( nindp*sizeof(int) );
    for (jj=0;jj<nindp;jj++) incr[jj] = 1;

    for (int iss=1;iss<ns;iss++) {
		iii=iss; jj=0;
		while ( !(iii%nsamp[jj]) ) {
			iii /= nsamp[jj];
			incr[jj] = -incr[jj];
			jj++;
		}
		dp /= dprobs[ indx[jj] ];
		indx[jj] += incr[jj];
		dp *= dprobs[ indx[jj] ];

		// set data
		drlo[n_first_stg_rows + jj] = demand[ indx[jj] ];
		drup[n_first_stg_rows + jj] = demand[ indx[jj] ];

		cpv_rlo.setElement(cpv_rlo.findIndex(n_first_stg_rows + jj),demand[ indx[jj] ]);
		cpv_rup.setElement(cpv_rup.findIndex(n_first_stg_rows + jj),demand[ indx[jj] ]);

		// genScenario
		is = smiModel->generateScenario(osiCore,NULL,NULL,NULL,NULL,
			&cpv_rlo,&cpv_rup,branch,anc,dp);


	}

	assert(ns==smiModel->getNumScenarios());
	cout << "ModelScenario: Finished adding scenarios" << endl;



	// load problem data into OsiSolver
	smiModel->loadOsiSolverData();
	// get Osi pointer
	OsiSolverInterface *smiOsi = smiModel->getOsiSolverInterface();
	// set some parameters
	smiOsi->setHintParam(OsiDoPresolveInInitial,true);
	smiOsi->setHintParam(OsiDoScale,true);
	smiOsi->setHintParam(OsiDoCrash,true);
	// solve using Osi Solver
	smiOsi->initialSolve();
	// test optimal value
    	assert(fabs(smiOsi->getObjValue()-1566.042)<0.01);

	// test solutions
	const double *dsoln = smiOsi->getColSolution();
	double objSum = 0.0;

	/* The canonical way to traverse the tree:
	   For each scenario, get the leaf node.
	   Then get the parent.  Repeat until parent is NULL.
	   (Only the root node has a NULL parent.)
	 */


	for(is=0; is<ns; ++is)
	{
		/* this loop calculates the scenario objective value */
		double scenSum = 0.0;

		// start with leaf node
		SmiScnNode *node = smiModel->getLeafNode(is);

		// leaf node probability is the scenario probability
		double scenprob = node->getModelProb();

		while (node != NULL)
		{

			// getColStart returns the starting index of node in OSI model
			for(int j=node->getColStart(); j<node->getColStart()+node->getNumCols(); ++j)
			{
				// getCoreColIndex returns the corresponding Core index
				// in the original (user's) ordering
				scenSum += dobj[node->getCoreColIndex(j)]*dsoln[j];


			}
			// get parent of node
			node = node->getParent();
		}
		objSum += scenSum*scenprob;
	}

	assert(fabs(smiOsi->getObjValue()-objSum) < 0.01);

		// print results
		printf("Solved stochastic program %s\n", name);
		printf("Number of rows: %d\n",smiOsi->getNumRows());
		printf("Number of cols: %d\n",smiOsi->getNumCols());
		printf("Optimal value: %g\n",smiOsi->getObjValue());

}

//--------------------------------------------------------------------------
// test cut debugger methods.
void
OsiRowCutDebuggerUnitTest(const OsiSolverInterface * baseSiP, const std::string & mpsDir)
{
  
  CoinRelFltEq eq;
  
  // Test default constructor
  {
    OsiRowCutDebugger r;
    OSIUNITTEST_ASSERT_ERROR(r.integerVariable_ == NULL, {}, "osirowcutdebugger", "default constructor");
    OSIUNITTEST_ASSERT_ERROR(r.knownSolution_   == NULL, {}, "osirowcutdebugger", "default constructor");
    OSIUNITTEST_ASSERT_ERROR(r.numberColumns_   == 0,    {}, "osirowcutdebugger", "default constructor");
  }
  
  {
    // Get non trivial instance
    OsiSolverInterface * imP = baseSiP->clone();
    std::string fn = mpsDir+"exmip1";
    imP->readMps(fn.c_str(),"mps");
    OSIUNITTEST_ASSERT_ERROR(imP->getNumRows() == 5, {}, "osirowcutdebugger", "read exmip1");
    
    /*
      Activate the debugger. The garbled name here is deliberate; the
      debugger should isolate the portion of the string between '/' and
      '.' (in normal use, this would be the base file name, stripped of
      the prefix and extension).
    */
    imP->activateRowCutDebugger("ab cd /x/ /exmip1.asc");
    
    int i;
    
    // return debugger
    const OsiRowCutDebugger * debugger = imP->getRowCutDebugger();
    OSIUNITTEST_ASSERT_ERROR(debugger != NULL, {}, "osirowcutdebugger", "return debugger");
    OSIUNITTEST_ASSERT_ERROR(debugger->numberColumns_ == 8, {}, "osirowcutdebugger", "return debugger");
    
    const bool type[]={0,0,1,1,0,0,0,0};
    const double values[]= {2.5, 0, 1, 1, 0.5, 3, 0, 0.26315789473684253};
    CoinPackedVector objCoefs(8,imP->getObjCoefficients());
   
    bool type_ok = true;
#if 0
    for (i=0;i<8;i++)
      type_ok &= type[i] == debugger->integerVariable_[i];
    OSIUNITTEST_ASSERT_ERROR(type_ok, {}, "osirowcutdebugger", "???");
#endif
    
    double objValue = objCoefs.dotProduct(values);
    double debuggerObjValue = objCoefs.dotProduct(debugger->knownSolution_);
    OSIUNITTEST_ASSERT_ERROR(eq(objValue, debuggerObjValue), {}, "osirowcutdebugger", "objective value");
    
    OsiRowCutDebugger rhs;
    {
      OsiRowCutDebugger rC1(*debugger);

      OSIUNITTEST_ASSERT_ERROR(rC1.numberColumns_ == 8, {}, "osirowcutdebugger", "copy constructor");
      type_ok = true;
      for (i=0;i<8;i++)
      	type_ok &= type[i] == rC1.integerVariable_[i];
      OSIUNITTEST_ASSERT_ERROR(type_ok, {}, "osirowcutdebugger", "copy constructor");
      OSIUNITTEST_ASSERT_ERROR(eq(objValue,objCoefs.dotProduct(rC1.knownSolution_)), {}, "osirowcutdebugger", "copy constructor");
      
      rhs = rC1;
      OSIUNITTEST_ASSERT_ERROR(rhs.numberColumns_ == 8, {}, "osirowcutdebugger", "assignment operator");
      type_ok = true;
      for (i=0;i<8;i++)
      	type_ok &= type[i] == rhs.integerVariable_[i];
      OSIUNITTEST_ASSERT_ERROR(type_ok, {}, "osirowcutdebugger", "assignment operator");
      OSIUNITTEST_ASSERT_ERROR(eq(objValue,objCoefs.dotProduct(rhs.knownSolution_)), {}, "osirowcutdebugger", "assignment operator");
    }
    // Test that rhs has correct values even though lhs has gone out of scope
    OSIUNITTEST_ASSERT_ERROR(rhs.numberColumns_ == 8, {}, "osirowcutdebugger", "assignment operator");
    type_ok = true;
    for (i=0;i<8;i++)
    	type_ok &= type[i] == rhs.integerVariable_[i];
    OSIUNITTEST_ASSERT_ERROR(type_ok, {}, "osirowcutdebugger", "assignment operator");
    OSIUNITTEST_ASSERT_ERROR(eq(objValue,objCoefs.dotProduct(rhs.knownSolution_)), {}, "osirowcutdebugger", "assignment operator");

    OsiRowCut cut[2];
    
    const int ne = 3;
    int inx[ne] = { 0, 2, 3 };
    double el[ne] = { 1., 1., 1. };
    cut[0].setRow(ne,inx,el);
    cut[0].setUb(5.);
    
    el[1]=5;
    cut[1].setRow(ne,inx,el);
    cut[1].setUb(5);
    OsiCuts cs; 
    cs.insert(cut[0]);
    cs.insert(cut[1]);
    OSIUNITTEST_ASSERT_ERROR(!debugger->invalidCut(cut[0]), {}, "osirowcutdebugger", "recognize (in)valid cut");
    OSIUNITTEST_ASSERT_ERROR( debugger->invalidCut(cut[1]), {}, "osirowcutdebugger", "recognize (in)valid cut");
    OSIUNITTEST_ASSERT_ERROR(debugger->validateCuts(cs,0,2) == 1, {}, "osirowcutdebugger", "recognize (in)valid cut");
    OSIUNITTEST_ASSERT_ERROR(debugger->validateCuts(cs,0,1) == 0, {}, "osirowcutdebugger", "recognize (in)valid cut");
    delete imP;
  }
//.........这里部分代码省略.........

//-----------------------------------------------------------------------------
// Generate Lift-and-Project cuts
//------------------------------------------------------------------- 
void CglLiftAndProject::generateCuts(const OsiSolverInterface& si, OsiCuts& cs,
				     const CglTreeInfo /*info*/)
{
  // Assumes the mixed 0-1 problem 
  //
  //   min {cx: <Atilde,x> >= btilde} 
  //
  // is in canonical form with all bounds,
  // including x_t>=0, -x_t>=-1 for x_t binary,
  // explicitly stated in the constraint matrix. 
  // See ~/COIN/Examples/Cgl2/cgl2.cpp 
  // for a general purpose "convert" function. 

  // Reference [BCC]: Balas, Ceria, and Corneujols,
  // "A lift-and-project cutting plane algorithm
  // for mixed 0-1 program", Math Prog 58, (1993) 
  // 295-324.

  // This implementation uses Normalization 1.

  // Given canonical problem and
  // the lp-relaxation solution, x,
  // the LAP cut generator attempts to construct
  // a cut for every x_j such that 0<x_j<1
  // [BCC:307]
 

  // x_j is the strictly fractional binary variable
  // the cut is generated from
  int j = 0; 

  // Get basic problem information
  // let Atilde be an m by n matrix
  const int m = si.getNumRows(); 
  const int n = si.getNumCols(); 
  const double * x = si.getColSolution();

  // Remember - Atildes may have gaps..
  const CoinPackedMatrix * Atilde = si.getMatrixByRow();
  const double * AtildeElements =  Atilde->getElements();
  const int * AtildeIndices =  Atilde->getIndices();
  const CoinBigIndex * AtildeStarts = Atilde->getVectorStarts();
  const int * AtildeLengths = Atilde->getVectorLengths();  
  const int AtildeFullSize = AtildeStarts[m];
  const double * btilde = si.getRowLower();

  // Set up memory for system (10) [BCC:307]
  // (the problem over the norm intersected 
  //  with the polar cone)
  // 
  // min <<x^T,Atilde^T>,u> + x_ju_0
  // s.t.
  //     <B,w> = (0,...,0,beta_,beta)^T
  //        w  is nonneg for all but the
  //           last two entries, which are free.
  // where 
  // w = (u,v,v_0,u_0)in BCC notation 
  //      u and v are m-vectors; u,v >=0
  //      v_0 and u_0 are free-scalars, and
  //  
  // B = Atilde^T  -Atilde^T  -e_j e_j
  //     btilde^T   e_0^T      0   0
  //     e_0^T      btilde^T   1   0

  // ^T indicates Transpose
  // e_0 is a (AtildeNCols x 1) vector of all zeros 
  // e_j is e_0 with a 1 in the jth position

  // Storing B in column order. B is a (n+2 x 2m+2) matrix 
  // But need to allow for possible gaps in Atilde.
  // At each iteration, only need to change 2 cols and objfunc
  // Sane design of OsiSolverInterface does not permit mucking
  // with matrix.
  // Because we must delete and add cols to alter matrix,
  // and we can only add columns on the end of the matrix
  // put the v_0 and u_0 columns on the end.
  // rather than as described in [BCC]
 
  // Initially allocating B with space for v_0 and u_O cols
  // but not populating, for efficiency.

  // B without u_0 and v_0 is a (n+2 x 2m) size matrix.

  int twoM = 2*m;
  int BNumRows = n+2;
  int BNumCols = twoM+2;
  int BFullSize = 2*AtildeFullSize+twoM+3;
  double * BElements = new double[BFullSize];
  int * BIndices = new int[BFullSize];
  CoinBigIndex * BStarts = new CoinBigIndex [BNumCols+1];
  int * BLengths = new int[BNumCols];


  int i, ij, k=0;
  int nPlus1=n+1;
  int offset = AtildeStarts[m]+m;
  for (i=0; i<m; i++){
//.........这里部分代码省略.........