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

int
main(void)
{
   // Create a problem pointer.  We use the base class here.
   OsiSolverInterface *si;

   // When we instantiate the object, we need a specific derived class.
   si = new OsiClpSolverInterface;

   // Read in an mps file.  This one's from the MIPLIB library.
   si->readMps("../../Data/Sample/p0033");

   // Solve the (relaxation of the) problem
   si->initialSolve();

   // Check the solution
   if ( si->isProvenOptimal() ) { 
      std::cout << "Found optimal solution!" << std::endl; 
      std::cout << "Objective value is " << si->getObjValue() << std::endl;

      int n = si->getNumCols();
      const double *solution;
      solution = si->getColSolution();
      // We could then print the solution or examine it.
   } else {
      std::cout << "Didn't find optimal solution." << std::endl;
      // Could then check other status functions.
   }

   return 0;
}

void SmpsIO(const char * const name )
{
		SmiScnModel smi;

		// read SMPS model from files
		//	<name>.core, <name>.time, and <name>.stoch
		smi.readSmps(name);

		// generate OSI solver object
		// 	here we use OsiClp
		OsiClpSolverInterface *clp = new OsiClpSolverInterface();

		// set solver object for SmiScnModel
		smi.setOsiSolverHandle(*clp);

		// load solver data
		// 	this step generates the deterministic equivalent
		//	and returns an OsiSolver object
		OsiSolverInterface *osiStoch = smi.loadOsiSolverData();

		// set some nice Hints to the OSI solver
		osiStoch->setHintParam(OsiDoPresolveInInitial,true);
		osiStoch->setHintParam(OsiDoScale,true);
		osiStoch->setHintParam(OsiDoCrash,true);

		// solve
		osiStoch->initialSolve();

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

		// print solution to file
		string outfilename(name);
		const string suffix(".out");
		outfilename = outfilename + suffix;
		FILE *fp = fopen(outfilename.c_str(),"w");
		int numScenarios=smi.getNumScenarios();

		for (int i=0 ; i<numScenarios; ++i) {
			double *dsoln=NULL;
			int numCols=0;
			fprintf(fp,"Scenario %d \n",i);
			dsoln = smi.getColSolution(i,&numCols);
			for (int j=0; j<numCols; j++)
				fprintf(fp,"%g \n",dsoln[j]);
			free(dsoln);
		}
		fclose(fp);


}

// Override OsiSolver::solve() because Soplex doesn't have branch & bound
int OsiSpxSolver::solve() {
	OsiSolverInterface* si = OsiSolver::getSolver();
	if(!OsiSolver::prepareSolve()) {
		return UNSAT;
	}

	timer.reset();
	try {
		si->initialSolve();
	} catch (CoinError err) {
		err.print(true);
		return UNSAT;
	}
	time = timer.timeElapsed();

	if (si->isProvenOptimal())
		return SAT;
	else if (si->isProvenPrimalInfeasible())
		return UNSAT;
	else
		return UNKNOWN;
}

void
CglResidualCapacityUnitTest(const OsiSolverInterface *baseSiP,
			    const std::string mpsDir)
{
  // Test default constructor
  {
    CglResidualCapacity aGenerator;
  }
  
  // Test copy & assignment
  {
    CglResidualCapacity rhs;
    {
      CglResidualCapacity bGenerator;
      CglResidualCapacity cGenerator(bGenerator);
      rhs=bGenerator;
    }
  }

  // Test get/set methods
  {
    CglResidualCapacity getset;
    
    double geps = 10 * getset.getEpsilon();
    getset.setEpsilon(geps);
    double geps2 = getset.getEpsilon();
    assert(geps == geps2);

    double gtol = 10 * getset.getTolerance();
    getset.setTolerance(gtol);
    double gtol2 = getset.getTolerance();
    assert(gtol == gtol2);

    int gpre = getset.getDoPreproc();
    gpre = (gpre + 1) % 3 - 1;
    getset.setDoPreproc(gpre);
    int gpre2 = getset.getDoPreproc();
    assert(gpre == gpre2);
  }

  // Test generateCuts
  {
    CglResidualCapacity gct;
    OsiSolverInterface  *siP = baseSiP->clone();
    std::string fn = mpsDir+"capPlan1";
    std::string fn2 = mpsDir+"capPlan1.mps";
    FILE *in_f = fopen(fn2.c_str(), "r");
    if(in_f == NULL) {
      std::cout<<"Can not open file "<<fn2<<std::endl<<"Skip test of CglResidualCapacity::generateCuts()"<<std::endl;
    }
    else {
      fclose(in_f);
      siP->readMps(fn.c_str(),"mps");
 
      siP->initialSolve();
      double lpRelax = siP->getObjValue();
      
      OsiCuts cs;
      gct.setDoPreproc(1); // Needed for DyLP
      gct.generateCuts(*siP, cs);
      int nRowCuts = cs.sizeRowCuts();
      std::cout<<"There are "<<nRowCuts<<" Residual Capacity cuts"<<std::endl;
      assert(cs.sizeRowCuts() > 0);
      OsiSolverInterface::ApplyCutsReturnCode rc = siP->applyCuts(cs);
      
      siP->resolve();
      
      double lpRelaxAfter= siP->getObjValue(); 
      std::cout<<"Initial LP value: "<<lpRelax<<std::endl;
      std::cout<<"LP value with cuts: "<<lpRelaxAfter<<std::endl;
      assert( lpRelax < lpRelaxAfter );
      assert(lpRelaxAfter < 964);
    }
    delete siP;
  }

}

void
CglMixedIntegerRoundingUnitTest(const OsiSolverInterface *baseSiP,
			    const std::string mpsDir)
{
  // Test default constructor
  {
    CglMixedIntegerRounding aGenerator;
  }

  // Test copy & assignment
  {
    CglMixedIntegerRounding rhs;
    {
      CglMixedIntegerRounding bGenerator;
      CglMixedIntegerRounding cGenerator(bGenerator);
      rhs=bGenerator;
    }
  }

  // Test get/set methods
  {
    CglMixedIntegerRounding getset;

    int gagg = 10 * getset.getMAXAGGR_();
    getset.setMAXAGGR_(gagg);
    int gagg2 = getset.getMAXAGGR_();
    assert(gagg == gagg2);

    bool gmult = !getset.getMULTIPLY_();
    getset.setMULTIPLY_(gmult);
    bool gmult2 = getset.getMULTIPLY_();
    assert(gmult == gmult2);

    int gcrit = getset.getCRITERION_();
    gcrit = (gcrit) % 3 + 1;
    getset.setCRITERION_(gcrit);
    int gcrit2 = getset.getCRITERION_();
    assert(gcrit == gcrit2);

    int gpre = getset.getDoPreproc();
    gpre = (gpre + 1) % 3 - 1;
    getset.setDoPreproc(gpre);
    int gpre2 = getset.getDoPreproc();
    assert(gpre == gpre2);
  }

  // Test generateCuts
  {
    CglMixedIntegerRounding gct;
    OsiSolverInterface  *siP = baseSiP->clone();
    std::string fn = mpsDir+"capPlan1";
    std::string fn2 = mpsDir+"capPlan1.mps";
    FILE *in_f = fopen(fn2.c_str(), "r");
    if(in_f == NULL) {
      std::cout<<"Can not open file "<<fn2<<std::endl<<"Skip test of CglMixedIntegerRounding::generateCuts()"<<std::endl;
    }
    else {
      fclose(in_f);
      siP->readMps(fn.c_str(),"mps");

      siP->initialSolve();
      double lpRelax = siP->getObjValue();

      OsiCuts cs;
      gct.generateCuts(*siP, cs);
      int nRowCuts = cs.sizeRowCuts();
      std::cout<<"There are "<<nRowCuts<<" MIR cuts"<<std::endl;
      assert(cs.sizeRowCuts() > 0);
      OsiSolverInterface::ApplyCutsReturnCode rc = siP->applyCuts(cs);

      siP->resolve();

      double lpRelaxAfter= siP->getObjValue();
      std::cout<<"Initial LP value: "<<lpRelax<<std::endl;
      std::cout<<"LP value with cuts: "<<lpRelaxAfter<<std::endl;
      assert( lpRelax < lpRelaxAfter );
      assert(lpRelaxAfter < 964);
    }
    delete siP;
  }

}

// mex function usage:
//  [x,y,status] = mexosi(n_vars,n_cons,A,x_lb,x_ub,c,Ax_lb,Ax_ub,isMIP,isQP,vartype,Q,options)
//                        0      1      2 3    4    5 6     7     8     9    10     11 12
void mexFunction(int nlhs, mxArray *plhs[],int nrhs, const mxArray *prhs[])
{
	// Enable printing in MATLAB
	int loglevel = 0;
	DerivedHandler *mexprinter = new DerivedHandler(); // assumed open	
	mexprinter->setLogLevel(loglevel);		 
	// check that we have the right number of inputs
	if(nrhs < 10) mexErrMsgTxt("At least 10 inputs required in call to mexosi. Bug in osi.m?...");
	// check that we have the right number of outputs
	if(nlhs < 3) mexErrMsgTxt("At least 3 ouptuts required in call to mexosi. Bug in osi.m?...");

    // Get pointers to input values
	const int    n_vars = (int)*mxGetPr(prhs[0]);
	const int    n_cons = (int)*mxGetPr(prhs[1]);
	const mxArray    *A =  prhs[2];
	const double  *x_lb =  mxGetPr(prhs[3]);
	const double  *x_ub =  mxGetPr(prhs[4]);
	const double     *c =  mxGetPr(prhs[5]);
	const double *Ax_lb =  mxGetPr(prhs[6]);
	const double *Ax_ub =  mxGetPr(prhs[7]);
	const bool    isMIP = (bool)*(mxLogical*)mxGetData(prhs[8]);
	const bool     isQP = (bool)*(mxLogical*)mxGetData(prhs[9]);
	mxLogical *isinteger = (mxLogical*)mxGetData(prhs[10]);
	const mxArray*    Q = prhs[11];
	// process the options
	int  returnStatus = 0;
	// extract row/col/value data from A
	const mwIndex * A_col_starts = mxGetJc(A);
	const mwIndex * A_row_index  = mxGetIr(A);
	const double  * A_data       = mxGetPr(A);
    // figure out the number of non-zeros in A
    int nnz = (int)(A_col_starts[n_vars] - A_col_starts[0]); // number of non-zeros
    //mexPrintf("nnz = %d, n_vars = %d, n_cons = %d\n",nnz,n_vars,n_cons);

    // we need to convert these into other types of indices for Coin to use them
    std::vector<CoinBigIndex> A_col_starts_coin(A_col_starts,A_col_starts+n_vars+1);
    std::vector<int>          A_row_index_coin(A_row_index,A_row_index+nnz);
    
    // declare the solver
    OsiSolverInterface* pSolver;
// 	initialize the solver
    if ( isMIP ) {
        pSolver = new OsiCbcSolverInterface;
    } else {
        pSolver = new OsiClpSolverInterface;
    }
   
//	OsiCbcSolverInterface is deprecated and CbcModel should be used instead but don't
//	know how to get that working with loadProblem. 
//     OsiCbcSolverInterface solver1;
//     CbcModel model(solver1);
//     CbcMain0(model);
//     OsiSolverInterface * pSolver = model.solver();
    
	if (nrhs>12) { // get stuff out of the options structure if provided
		// Finish me
	}
    
//     mexPrintf("Setting Log Level to 0.\n");
	// load the problem
	mexPrintf("Loading the problem.\n");
	pSolver->loadProblem( n_vars, n_cons, // problem size
						  &A_col_starts_coin[0], &A_row_index_coin[0], A_data, // the A matrix
						  x_lb,  x_ub, c, // the objective and bounds
						  Ax_lb, Ax_ub ); // the constraint bounds
    
//     pSolver->messageHandler()->setLogLevel(0); // This doesn't seem to work
    pSolver->setHintParam(OsiDoReducePrint,true,OsiHintTry);
    
	// deal with integer inputs
	if ( isMIP ) {
		for(int i=0;i<n_vars;i++) {
			if (isinteger[i]) pSolver->setInteger(i);
		}
	}
	if (isQP) {
		error("QP is not working yet");
		// need to call loadQuadraticObjective here ???
	}
    
//     CbcModel model(pSolver);
//     model.solver()->setHintParam(OsiDoReducePrint,true,OsiHintTry);
    
    
	// solve the problem
	//mexPrintf("Trying to solve the problem.\n");
    if (isMIP) {
        pSolver->branchAndBound();
//         model.branchAndBound();
    } else {
//         model.initialSolve();
        pSolver->initialSolve();
    }
	
	// Allocate memory for return data
    plhs[0] = mxCreateDoubleMatrix(n_vars,1, mxREAL); // for the solution
    plhs[1] = mxCreateDoubleMatrix(n_cons,1, mxREAL); // for the constraint prices
//.........这里部分代码省略.........

int
main(void)
{
   // Create a problem pointer.  We use the base class here.
   OsiSolverInterface *si;

   // When we instantiate the object, we need a specific derived class.
   si = new OSIXXX;

   // Build our own instance from scratch

   /*
    * This section adapted from Matt Galati's example 
    * on the COIN-OR Tutorial website.
    *
    * Problem from Bertsimas, Tsitsiklis page 21
    *  
    *  optimal solution: x* = (1,1)
    *  
    *  minimize -1 x0 - 1 x1
    *  s.t       1 x0 + 2 x1 <= 3
    *            2 x0 + 1 x1 <= 3
    *              x0        >= 0
    *              x1        >= 0
    */

   int n_cols = 2;
   double *objective    = new double[n_cols];//the objective coefficients
   double *col_lb       = new double[n_cols];//the column lower bounds
   double *col_ub       = new double[n_cols];//the column upper bounds

   //Define the objective coefficients.
   //minimize -1 x0 - 1 x1
   objective[0] = -1.0;
   objective[1] = -1.0;

   //Define the variable lower/upper bounds.
   // x0 >= 0   =>  0 <= x0 <= infinity
   // x1 >= 0   =>  0 <= x1 <= infinity
   col_lb[0] = 0.0;
   col_lb[1] = 0.0;
   col_ub[0] = si->getInfinity();
   col_ub[1] = si->getInfinity();
     
   int n_rows = 2;
   double *row_lb = new double[n_rows]; //the row lower bounds
   double *row_ub = new double[n_rows]; //the row upper bounds
     
   //Define the constraint matrix.
   CoinPackedMatrix *matrix =  new CoinPackedMatrix(false,0,0);
   matrix->setDimensions(0, n_cols);

   //1 x0 + 2 x1 <= 3  =>  -infinity <= 1 x0 + 2 x2 <= 3
   CoinPackedVector row1;
   row1.insert(0, 1.0);
   row1.insert(1, 2.0);
   row_lb[0] = -1.0 * si->getInfinity();
   row_ub[0] = 3.0;
   matrix->appendRow(row1);

   //2 x0 + 1 x1 <= 3  =>  -infinity <= 2 x0 + 1 x1 <= 3
   CoinPackedVector row2;
   row2.insert(0, 2.0);
   row2.insert(1, 1.0);
   row_lb[1] = -1.0 * si->getInfinity();
   row_ub[1] = 3.0;
   matrix->appendRow(row2);

   //load the problem to OSI
   si->loadProblem(*matrix, col_lb, col_ub, objective, row_lb, row_ub);

   //write the MPS file to a file called example.mps
   si->writeMps("example");

  

   // Solve the (relaxation of the) problem
   si->initialSolve();

   // Check the solution
   if ( si->isProvenOptimal() ) { 
      std::cout << "Found optimal solution!" << std::endl; 
      std::cout << "Objective value is " << si->getObjValue() << std::endl;

      int n = si->getNumCols();
      const double *solution;
      solution = si->getColSolution();
      // We could then print the solution or examine it.
   } else {
      std::cout << "Didn't find optimal solution." << std::endl;
      // Could then check other status functions.
   }

   return 0;
}

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


}

//.........这里部分代码省略.........
    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
  {
    CglSimpleRounding cg;
    
    OsiSolverInterface * siP = baseSiP->clone();
    std::string fn = mpsDir+"p0033";
    siP->readMps(fn.c_str(),"mps");
    OsiCuts cuts;
    cg.generateCuts(*siP,cuts);

    // p0033 is the optimal solution to p0033
    int objIndices[14] = { 
       0,  6,  7,  9, 13, 17, 18,
      22, 24, 25, 26, 27, 28, 29 };
    CoinPackedVector p0033(14,objIndices,1.0);

    // test that none of the generated cuts
    // chops off the optimal solution
    int nRowCuts = cuts.sizeRowCuts();
    OsiRowCut rcut;
    CoinPackedVector rpv;
    int i;
    for (i=0; i<nRowCuts; i++){
      rcut = cuts.rowCut(i);
      rpv = rcut.row();
      double p0033Sum = (rpv*p0033).sum();
      double rcutub = rcut.ub();
      assert (p0033Sum <= rcutub);
    }

    // test that the cuts improve the 
    // lp objective function value
    siP->initialSolve();
    double lpRelaxBefore=siP->getObjValue();
    OsiSolverInterface::ApplyCutsReturnCode rc = siP->applyCuts(cuts);
    siP->resolve();
    double lpRelaxAfter=siP->getObjValue(); 
#ifdef CGL_DEBUG
    printf("\n\nOrig LP min=%f\n",lpRelaxBefore);
    printf("Final LP min=%f\n\n",lpRelaxAfter);
#endif
    assert( lpRelaxBefore < lpRelaxAfter );

    delete siP;

  }


}

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

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

void
CglRedSplitUnitTest(const OsiSolverInterface *baseSiP,
		    const std::string mpsDir)
{
  // Test default constructor
  {
    CglRedSplit aGenerator;
  }
  
  // Test copy & assignment
  {
    CglRedSplit rhs;
    {
      CglRedSplit bGenerator;
      CglRedSplit cGenerator(bGenerator);
      rhs=bGenerator;
    }
  }

  // Test get/set methods
  {
    CglRedSplit getset;
    CglRedSplitParam gsparam = getset.getParam();
    
    double geps = 10 * gsparam.getEPS();
    gsparam.setEPS(geps);
    double geps2 = gsparam.getEPS();
    assert(geps == geps2);

    double gepse = 10 * gsparam.getEPS_ELIM();
    gsparam.setEPS_ELIM(gepse);
    double gepse2 = gsparam.getEPS_ELIM();
    assert(gepse == gepse2);

    double gmv = 10 * gsparam.getMINVIOL();
    gsparam.setMINVIOL(gmv);
    double gmv2 = gsparam.getMINVIOL();
    assert(gmv == gmv2);

    int gucg = gsparam.getUSE_CG2();
    gucg = 1 - gucg;
    gsparam.setUSE_CG2(gucg);
    int gucg2 = gsparam.getUSE_CG2();
    assert(gucg == gucg2);
  }

  // Test generateCuts
  {
    CglRedSplit gct;
    OsiSolverInterface  *siP = baseSiP->clone();
    std::string fn = mpsDir+"p0033";
    std::string fn2 = mpsDir+"p0033.mps";
    FILE *in_f = fopen(fn2.c_str(), "r");
    if(in_f == NULL) {
      std::cout<<"Can not open file "<<fn2<<std::endl<<"Skip test of CglRedSplit::generateCuts()"<<std::endl;
    }
    else {
      fclose(in_f);
      siP->readMps(fn.c_str(),"mps");
 
      siP->initialSolve();
      double lpRelax = siP->getObjValue();
      
      OsiCuts cs;
      gct.getParam().setMAX_SUPPORT(34);
      gct.getParam().setUSE_CG2(1);
      //      gct.getParam().setUSE_CG2(1);
      gct.generateCuts(*siP, cs);
      int nRowCuts = cs.sizeRowCuts();
      std::cout<<"There are "<<nRowCuts<<" Reduce-and-Split cuts"<<std::endl;
      assert(cs.sizeRowCuts() > 0);
      OsiSolverInterface::ApplyCutsReturnCode rc = siP->applyCuts(cs);
      
      siP->resolve();
      
      double lpRelaxAfter= siP->getObjValue(); 
      std::cout<<"Initial LP value: "<<lpRelax<<std::endl;
      std::cout<<"LP value with cuts: "<<lpRelaxAfter<<std::endl;
      assert( lpRelax < lpRelaxAfter );
      assert(lpRelaxAfter < 3089.1);
    }
    delete siP;
  }

}

//.........这里部分代码省略.........
                CoinPackedVector row = aCut.row();
                if (row.getNumElements() == 1)
                {
                    assert(row.getIndices()[0]==1);
                    assert(eq(row.getElements()[0], -1));
                }
                else if (row.getNumElements() == 2)
                {
                    assert(row.getIndices()[0]==0);
                    assert(eq(row.getElements()[0], 0.));
                    assert(row.getIndices()[1]==1);
                    assert(eq(row.getElements()[1], -1));
                }
                assert(eq(aCut.lb(), 0.));
                OsiSolverInterface::ApplyCutsReturnCode rc = siP->applyCuts(cuts);

                siP->resolve();
            }
            delete siP;
        }
    }

    if (1)  //Test on p0033
    {
        // Setup
        OsiSolverInterface  * siP = si->clone();
        std::string fn(mpsDir+"p0033");
        siP->readMps(fn.c_str(),"mps");
        siP->activateRowCutDebugger("p0033");
        CglLandP test;

        // Solve the LP relaxation of the model and
        // print out ofv for sake of comparison
        siP->initialSolve();
        double lpRelaxBefore=siP->getObjValue();
        assert( eq(lpRelaxBefore, 2520.5717391304347) );
#ifdef CGL_DEBUG
        printf("\n\nOrig LP min=%f\n",lpRelaxBefore);
#endif

        OsiCuts cuts;

        // Test generateCuts method
        test.generateCuts(*siP,cuts);
        OsiSolverInterface::ApplyCutsReturnCode rc = siP->applyCuts(cuts);

        siP->resolve();
        double lpRelaxAfter=siP->getObjValue();
        //assert( eq(lpRelaxAfter, 2592.1908295194507) );

        std::cout<<"Relaxation after "<<lpRelaxAfter<<std::endl;
        assert( lpRelaxAfter> 2840. );
#ifdef CGL_DEBUG
        printf("\n\nOrig LP min=%f\n",lpRelaxBefore);
        printf("\n\nFinal LP min=%f\n",lpRelaxAfter);
#endif
        assert( lpRelaxBefore < lpRelaxAfter );

        delete siP;
    }
    if (1)  //test again with modularization
    {
        // Setup
        OsiSolverInterface  * siP = si->clone();
        std::string fn(mpsDir+"p0033");
        siP->readMps(fn.c_str(),"mps");

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

//--------------------------------------------------------------------------
void
CglKnapsackCoverUnitTest(
  const OsiSolverInterface * baseSiP,
  const std::string mpsDir )
{
  int i;
  CoinRelFltEq eq(0.000001);

  // Test default constructor
  {
    CglKnapsackCover kccGenerator;
  }
  
  // Test copy & assignment
  {
    CglKnapsackCover rhs;
    {
      CglKnapsackCover kccGenerator;
      CglKnapsackCover cgC(kccGenerator);
      rhs=kccGenerator;
    }
  }


  // test exactSolveKnapsack
  {  
    CglKnapsackCover kccg;
    const int n=7;
    double c=50;
    double p[n] = {70,20,39,37,7,5,10};
    double w[n] = {31, 10, 20, 19, 4, 3, 6};
    double z;
    int x[n];
    int exactsol = kccg.exactSolveKnapsack(n, c, p, w, z, x);
    assert(exactsol==1);
    assert (z == 107);
    assert (x[0]==1);
    assert (x[1]==0);
    assert (x[2]==0);
    assert (x[3]==1);
    assert (x[4]==0);
    assert (x[5]==0);
    assert (x[6]==0);
  }

  /*
  // Testcase /u/rlh/osl2/mps/scOneInt.mps
  // Model has 3 continous, 2 binary, and 1 general
  // integer variable.
  {
    OsiSolverInterface  * siP = baseSiP->clone();
    int * complement=NULL;
    double * xstar=NULL;

    siP->readMps("../Mps/scOneInt","mps");
    CglKnapsackCover kccg;
    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()->vectorSize(i);
    assert(vectorsize==2);
    }
    
    kccg.cleanUpCutGenerator(complement,xstar);
    delete siP;
  }
  */  
  
  // Testcase /u/rlh/osl2/mps/tp3.mps
  // Models has 3 cols, 3 rows
  // Row 0 yields a knapsack, others do not.
  {
    // setup
    OsiSolverInterface  * siP = baseSiP->clone();
    std::string fn(mpsDir+"tp3");
    siP->readMps(fn.c_str(),"mps");     
    // All integer variables should be binary.
//.........这里部分代码省略.........

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

}