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

// Create result
void OsiSolverResult::createResult(const OsiSolverInterface &solver, const double *lowerBefore,
  const double *upperBefore)
{
  delete[] primalSolution_;
  delete[] dualSolution_;
  if (solver.isProvenOptimal() && !solver.isDualObjectiveLimitReached()) {
    objectiveValue_ = solver.getObjValue() * solver.getObjSense();
    CoinWarmStartBasis *basis = dynamic_cast< CoinWarmStartBasis * >(solver.getWarmStart());
    assert(basis);
    basis_ = *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());
  } else {
    // infeasible
    objectiveValue_ = COIN_DBL_MAX;
    basis_ = CoinWarmStartBasis();
    ;
    primalSolution_ = NULL;
    dualSolution_ = NULL;
  }
}

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


}

//.........这里部分代码省略.........
    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.
    // Assert that this is true.
    for ( i = 0;  i < siP->getNumCols();  i++ )
      if ( siP->isInteger(i) ) 
        assert(siP->getColUpper()[i]==1.0 && siP->isBinary(i));  
    OsiCuts cs;
    CoinPackedVector krow;
    double b=0;
    int nCols=siP->getNumCols();
    int * complement=new int [nCols];
    double * xstar=new double [nCols];

    CglKnapsackCover kccg;

    // solve LP relaxation
    // a "must" before calling initialization
    siP->initialSolve();
    double lpRelaxBefore=siP->getObjValue();
    std::cout<<"Initial LP value: "<<lpRelaxBefore<<std::endl;
    assert( eq(siP->getObjValue(), 97.185) );
    double mycs[] = {.627, .667558333333, .038};
    siP->setColSolution(mycs);
    const double *colsol = siP->getColSolution(); 
    int k;
    for (k=0; k<nCols; k++){
      xstar[k]=colsol[k];
      complement[k]=0;
    }
    
    // test deriveAKnapsack
    int rind = ( siP->getRowSense()[0] == 'N' ) ? 1 : 0;
    const CoinShallowPackedVector reqdBySunCC = siP->getMatrixByRow()->getVector(rind) ;
    int deriveaknap = kccg.deriveAKnapsack(*siP, cs, krow,b,complement,xstar,rind,reqdBySunCC);
    assert(deriveaknap ==1);
    assert(complement[0]==0);
    assert(complement[1]==1);
    assert(complement[2]==1);
    int inx[3] = {0,1,2};
    double el[3] = {161, 120, 68};
    CoinPackedVector r;
    r.setVector(3,inx,el);
    assert (krow == r);
    //assert (b == 183.0); ????? but x1 and x2 at 1 is valid 
    
    // test findGreedyCover 
    CoinPackedVector cover,remainder;
#if 0
    int findgreedy =  kccg.findGreedyCover( 0, krow, b, xstar, cover, remainder );
    assert( findgreedy == 1 );
    int coveri = cover.getNumElements();

//.........这里部分代码省略.........
  if(hSolver->isProvenOptimal()){

    double upperObjVal(0.0);

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

    MibSSolution *mibSol = NULL;

    OsiSolverInterface * lSolver = model->bS_->setUpModel(hSolver, true);

    if(0){
       lSolver->writeLp("tmp");
    }

    if(0){
       dynamic_cast<OsiCbcSolverInterface *> 
	  (lSolver)->getModelPtr()->messageHandler()->setLogLevel(0);
    }    
    else{
       dynamic_cast<OsiSymSolverInterface *> 
	  (lSolver)->setSymParam("prep_level", -1);
       
       dynamic_cast<OsiSymSolverInterface *> 
	  (lSolver)->setSymParam("verbosity", -2);
       
       dynamic_cast<OsiSymSolverInterface *> 
	  (lSolver)->setSymParam("max_active_nodes", 1);
    }

    lSolver->branchAndBound();

    if (lSolver->isProvenOptimal()){
       const double * sol = hSolver->getColSolution();
       double objVal(lSolver->getObjValue() * objSense);
       double etol(etol_);
       double lowerObj = getLowerObj(sol, objSense);  
       
       double * optUpperSolutionOrd = new double[uCols];
       double * optLowerSolutionOrd = new double[lCols];
       
       CoinZeroN(optUpperSolutionOrd, uCols);
       CoinZeroN(optLowerSolutionOrd, lCols);
       
       if(fabs(objVal - lowerObj) < etol){
	  
	  /** Current solution is bilevel feasible **/
	  
	  for(i = 0; i < tCols; i++)
	     upperObjVal += 
		hSolver->getColSolution()[i] * oSolver->getObjCoefficients()[i];
	  
	  mibSol = new MibSSolution(hSolver->getNumCols(),
				    hSolver->getColSolution(),
				    upperObjVal,
				    model);
	  
	  model->storeSolution(BlisSolutionTypeHeuristic, mibSol);
	  mibSol = NULL;
       }
       else{
	  
	  /* solution is not bilevel feasible, create one that is */
	  
	  const double * uSol = hSolver->getColSolution();
	  const double * lSol = lSolver->getColSolution();
	  int numElements(hSolver->getNumCols());

//.........这里部分代码省略.........
      }
      
      nSolver->setObjective(nObjCoeffs);
      
      CoinPackedVector objCon;
      
      for(i = 0; i < tCols; i++){
	objCon.insert(i, uObjCoeffs[i] * uObjSense);
      }
      
      nSolver->addRow(objCon, upperObjVal, upperObjVal);
      nSolver->writeLp("beta1");
      if(0){
	dynamic_cast<OsiCbcSolverInterface *> 
	  (nSolver)->getModelPtr()->messageHandler()->setLogLevel(0);
      }
      else{
	 dynamic_cast<OsiSymSolverInterface *> 
	    (nSolver)->setSymParam("prep_level", -1);

	 dynamic_cast<OsiSymSolverInterface *> 
	    (nSolver)->setSymParam("verbosity", -2);

	 dynamic_cast<OsiSymSolverInterface *> 
	    (nSolver)->setSymParam("max_active_nodes", 1);
      }

      nSolver->branchAndBound();
     

      double * colsol = new double[tCols];

      if(nSolver->isProvenOptimal()){
	lowerObjVal = nSolver->getObjValue();
	CoinCopyN(nSolver->getColSolution(), tCols, colsol);
      }
      else{
	//just take the current solution
	lowerObjVal = sSolver->getObjValue();
	CoinCopyN(sSolver->getColSolution(), tCols, colsol);
      }

      delete[] nObjCoeffs;
      nObjCoeffs = 0;
      delete sSolver;
      delete nSolver;
      return mcSol(std::make_pair(upperObjVal, lowerObjVal), colsol);
    }
    else if(beta == 0.0){
      
      /*
	fix lower-level objective to current value and 
	reoptimize wrt to upper-level objective
      */
      
      //OsiSolverInterface * nSolver = new OsiCbcSolverInterface();
      OsiSolverInterface * nSolver = new OsiSymSolverInterface();
      nSolver->loadProblem(*oSolver->getMatrixByCol(),
			   oSolver->getColLower(), oSolver->getColUpper(),
			   oSolver->getObjCoefficients(),
			   oSolver->getRowLower(), oSolver->getRowUpper());
      for(j = 0; j < tCols; j++){
	if(oSolver->isInteger(j))
	  nSolver->setInteger(j);
      }
      

/** Create a set of candidate branching objects. */
int 
BlisBranchStrategyPseudo::createCandBranchObjects(int numPassesLeft,
						  double ub)
{
    int bStatus = 0;
    int i, pass, colInd;

    int preferDir, saveLimit;
    int numFirsts  = 0;
    int numInfs = 0;
    int minCount = 0;
    int numLowerTightens = 0;
    int numUpperTightens = 0;
    double lpX, score, infeasibility, downDeg, upDeg, sumDeg = 0.0; 
    
    bool roundAgain, downKeep, downGood, upKeep, upGood;


    int *lbInd = NULL;
    int *ubInd = NULL;
    double *newLB = NULL;
    double *newUB = NULL;

    double *saveUpper = NULL;
    double *saveLower = NULL;
    double *saveSolution = NULL;

    BlisModel *model = dynamic_cast<BlisModel *>(model_);
    OsiSolverInterface *solver = model->solver();
    
    int numCols = model->getNumCols();
    int numObjects = model->numObjects();
    int aveIterations = model->getAveIterations();


    //std::cout <<  "aveIterations = " <<  aveIterations << std::endl;

     //------------------------------------------------------
    // Check if max time is reached or no pass is left.
    //------------------------------------------------------
    
    double timeLimit = model->AlpsPar()->entry(AlpsParams::timeLimit);
    AlpsKnowledgeBroker *broker = model->getKnowledgeBroker();
    bool maxTimeReached = (broker->timer().getTime() > timeLimit);
    bool selectNow = false;
    
    if (maxTimeReached || !numPassesLeft) {
        selectNow = true;
#ifdef BLIS_DEBUG
        printf("PSEUDO: CREATE: maxTimeReached %d, numPassesLeft %d\n", 
               maxTimeReached, numPassesLeft);
#endif
    }
    
    // Store first time objects.
    std::vector<BlisObjectInt *> firstObjects;

    // Store infeasible objects.
    std::vector<BlisObjectInt *> infObjects;

    // TODO: check if sorting is expensive.
    std::multimap<double, BcpsBranchObject*, BlisPseuoGreater> candObjects;

    double objValue = solver->getObjSense() * solver->getObjValue();

    const double * lower = solver->getColLower();
    const double * upper = solver->getColUpper();
    saveSolution = new double[numCols];
    memcpy(saveSolution, solver->getColSolution(), numCols*sizeof(double));

    //--------------------------------------------------
    // Find the infeasible objects.
    // NOTE: we might go round this loop twice if we are feed in
    //       a "feasible" solution.
    //--------------------------------------------------
    
    for (pass = 0; pass < 2; ++pass) {
	
        numInfs = 0;

        BcpsObject * object = NULL;
        BlisObjectInt * intObject = NULL;
            
        infObjects.clear();
        firstObjects.clear();
        
        for (i = 0; i < numObjects; ++i) {
                
            object = model->objects(i);
            infeasibility = object->infeasibility(model, preferDir);
            
            if (infeasibility) {
                
                ++numInfs;
                intObject = dynamic_cast<BlisObjectInt *>(object);
                
                if (intObject) {
                    infObjects.push_back(intObject);
                    
//.........这里部分代码省略.........

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

  }


}

//.........这里部分代码省略.........
     }
     if (whichCutsLL == 1){
	sym_set_int_param(env, "generate_cgl_knapsack_cuts", 
			  DO_NOT_GENERATE);
	sym_set_int_param(env, "generate_cgl_probing_cuts", 
			  DO_NOT_GENERATE);
	sym_set_int_param(env, "generate_cgl_clique_cuts", 
			  DO_NOT_GENERATE);
	sym_set_int_param(env, "generate_cgl_twomir_cuts", 
			  DO_NOT_GENERATE);
	sym_set_int_param(env, "generate_cgl_flowcover_cuts", 
			  DO_NOT_GENERATE);
     }
  }else if (feasCheckSolver == "CPLEX"){
#ifdef USE_CPLEX
     lSolver->setHintParam(OsiDoReducePrint);
     lSolver->messageHandler()->setLogLevel(0);
     CPXENVptr cpxEnv = 
	dynamic_cast<OsiCpxSolverInterface*>(lSolver)->getEnvironmentPtr();
     assert(cpxEnv);
     CPXsetintparam(cpxEnv, CPX_PARAM_SCRIND, CPX_OFF);
     CPXsetintparam(cpxEnv, CPX_PARAM_THREADS, maxThreadsLL);
#endif
  }
  
  if (warmStartLL && feasCheckSolver == "SYMPHONY"){
     lSolver->resolve();
     setWarmStart(lSolver->getWarmStart());
  }else{
     lSolver->branchAndBound();
  }

  const double * sol = model_->solver()->getColSolution();
  double objVal(lSolver->getObjValue() * model_->getLowerObjSense());
  
  MibSTreeNode * node = static_cast<MibSTreeNode *>(model_->activeNode_);
  MibSTreeNode * parent = 
    static_cast<MibSTreeNode *>(model_->activeNode_->getParent());

  if((!node->isBoundSet()) 
     && (node->getIndex() != 0)){
    double parentBound = parent->getLowerUB();
    node->setLowerUB(parentBound);
    node->setIsBoundSet(true);
  }
  
  if(objVal > node->getLowerUB()){
    
    node->setLowerUB(objVal);
    node->setIsBoundSet(true);
    
  }

  double etol(model_->etol_);
  double lowerObj = getLowerObj(sol, model_->getLowerObjSense());  

  int lN(model_->lowerDim_); // lower-level dimension
  int uN(model_->upperDim_); // lower-level dimension
  if(!optLowerSolution_)
    optLowerSolution_ = new double[lN];

  if(!optLowerSolutionOrd_)
    optLowerSolutionOrd_ = new double[lN];
  
  CoinZeroN(optLowerSolution_, lN);
  CoinZeroN(optLowerSolutionOrd_, lN);

//.........这里部分代码省略.........
	  std::cout<<"(col="<<column<<",el="<<elements[k]<<",sol="<<
	    colsol[column]<<") ";
	}
	std::cout <<std::endl;
      }
#endif

      if (i-nOldCuts==testCut) {
	assert( eq(rhs,ub));
	assert(n==1);
	for (k=0; k<n; k++){
	  int column=indices[k];
	  assert (eq(cut[column],elements[k]));
	}
	// add cut
	rowLower[3]=-1.0e100;
	rowUpper[3]=ub;
	matrix.appendRow(rpv);
      }
    }
    nOldCuts=nRowCuts;
    // basis 2
    int rowBasis2[4]={1,1,-1,-1};
    int colBasis2[2]={1,1};
    warm.setSize(2,4);
    for (i=0;i<4;i++) {
      if (rowBasis2[i]<0) {
	warm.setArtifStatus(i,CoinWarmStartBasis::atLowerBound);
      } else {
	warm.setArtifStatus(i,CoinWarmStartBasis::basic);
      }
    }
    for (i=0;i<2;i++) {
      if (colBasis2[i]<0) {
	warm.setStructStatus(i,CoinWarmStartBasis::atLowerBound);
      } else {
	warm.setStructStatus(i,CoinWarmStartBasis::basic);
      }
    }

    // solution 2
    double colsol2[2]={2.0,0.5};
    test1.generateCuts(NULL, osicuts, matrix,
		 /*objective,*/ colsol2,
		 colLower, colUpper,
		 rowLower, rowUpper, intVar, &warm);
    nRowCuts = osicuts.sizeRowCuts();
    std::cout<<"There are "<<nRowCuts<<" gomory cuts"<<std::endl;
    assert (nRowCuts==nOldCuts);
    
  }

  // Miplib3 problem p0033
  if (1) {
    // Setup
    OsiSolverInterface  * siP = baseSiP->clone();
    std::string fn(mpsDir+"p0033");
    siP->readMps(fn.c_str(),"mps");
    siP->activateRowCutDebugger("p0033");
    CglGomory test;

    // Solve the LP relaxation of the model and
    // print out ofv for sake of comparison 
    siP->initialSolve();
    double lpRelaxBefore=siP->getObjValue();
    std::cout<<"Initial LP value: "<<lpRelaxBefore<<std::endl;
    assert( eq(lpRelaxBefore, 2520.5717391304347) );

    // Fails with OsiCpx, OsiXpr:
    /**********
    double mycs[] = {0, 1, 0, 0, -2.0837010502455788e-19, 1, 0, 0, 1,
		       0.021739130434782594, 0.35652173913043478, 
		       -6.7220534694101275e-18, 5.3125906451789717e-18, 
		       1, 0, 1.9298798670241979e-17, 0, 0, 0,
		       7.8875708048320448e-18, 0.5, 0, 
		       0.85999999999999999, 1, 1, 0.57999999999999996,
		       1, 0, 1, 0, 0.25, 0, 0.67500000000000004};
    siP->setColSolution(mycs);
    ****/

    OsiCuts cuts;    
    
    // Test generateCuts method
    test.generateCuts(*siP,cuts);
    int nRowCuts = cuts.sizeRowCuts();
    std::cout<<"There are "<<nRowCuts<<" Gomory cuts"<<std::endl;
    assert(cuts.sizeRowCuts() > 0);
    OsiSolverInterface::ApplyCutsReturnCode rc = siP->applyCuts(cuts);
    
    siP->resolve();
    double lpRelaxAfter=siP->getObjValue(); 
    std::cout<<"LP value with cuts: "<<lpRelaxAfter<<std::endl;
    //assert( eq(lpRelaxAfter, 2592.1908295194507) );
    assert( lpRelaxAfter> 2550.0 );
    assert( lpRelaxBefore < lpRelaxAfter );
    assert(lpRelaxAfter < 3089.1);
    
    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;
  }

}

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

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

BlisReturnStatus
BlisStrongBranch(BlisModel *model, double objValue, int colInd, double x,
		 const double *saveLower, const double *saveUpper,
		 bool &downKeep, bool &downFinished, double &downDeg,
		 bool &upKeep, bool &upFinished, double &upDeg)
{
    BlisReturnStatus status = BlisReturnStatusOk;
    int lpStatus = 0;

    int j, numIntInfDown, numObjInfDown;

    double newObjValue;
    
    OsiSolverInterface * solver = model->solver();
    
    int numCols = solver->getNumCols();
    const double * lower = solver->getColLower();
    const double * upper = solver->getColUpper();

    // Restore bounds
    int numDiff = 0;

    BlisSolution* ksol = NULL;

    int ind = model->getIntObjIndices()[colInd];
    BlisObjectInt *intObj = dynamic_cast<BlisObjectInt *>(model->objects(ind));
    
#ifdef BLIS_DEBUG_MORE
    for (j = 0; j < numCols; ++j) {
	if (saveLower[j] != lower[j]) {
	    //solver->setColLower(j, saveLower[j]);
            ++numDiff;
	}
	if (saveUpper[j] != upper[j]) {
	    //solver->setColUpper(j, saveUpper[j]);
            ++numDiff;
	}
    }
    std::cout << "BEFORE: numDiff = " << numDiff << std::endl;
#endif	 
   
    //------------------------------------------------------
    // Branching down.
    //------------------------------------------------------

    solver->setColUpper(colInd, floor(x));
    solver->solveFromHotStart();
    
    newObjValue = solver->getObjSense() * solver->getObjValue();
    downDeg = newObjValue - objValue;
    
    if (solver->isProvenOptimal()) {
	lpStatus = 0; // optimal
#ifdef BLIS_DEBUG_MORE
        printf("STRONG: COL[%d]: downDeg=%g, x=%g\n", colInd, downDeg, x);
#endif
        // Update pseudocost
        intObj->pseudocost().update(-1, downDeg, x);
        model->setSharedObjectMark(ind);        

        // Check if ip feasible
        ksol = model->feasibleSolution(numIntInfDown, numObjInfDown);
        if (ksol) {
#ifdef BLIS_DEBUG_MORE
            printf("STRONG:Down:found a feasible solution\n");
#endif
            
            model->storeSolution(BlisSolutionTypeStrong, ksol);
	    downKeep = false;
        }
	else {
	    downKeep = true;
	}
	downFinished = true;
    }
    else if (solver->isIterationLimitReached() && 
	     !solver->isDualObjectiveLimitReached()) {
	lpStatus = 2;      // unknown 
	downKeep = true;
	downFinished = false;
    }
    else {
        downDeg = 1.0e20;
	lpStatus = 1; // infeasible
	downKeep = false;
	downFinished = false;
    }       
            
#ifdef BLIS_DEBUG_MORE
    std::cout << "Down: lpStatus = " << lpStatus << std::endl;
#endif
    
    // restore bounds
    numDiff = 0;
    for (j = 0; j < numCols; ++j) {
	if (saveLower[j] != lower[j]) {
	    solver->setColLower(j, saveLower[j]);
            ++numDiff;
	}
	if (saveUpper[j] != upper[j]) {
//.........这里部分代码省略.........

int main (int argc, const char *argv[])
{

  OsiClpSolverInterface solver1;
  //#define USE_OSI_NAMES
#ifdef USE_OSI_NAMES
  // Say we are keeping names (a bit slower this way)
  solver1.setIntParam(OsiNameDiscipline,1);
#endif
  // Read in model using argv[1]
  // and assert that it is a clean model
  std::string mpsFileName;
#if defined(SAMPLEDIR)
  mpsFileName = SAMPLEDIR "/p0033.mps";
#else
  if (argc < 2) {
    fprintf(stderr, "Do not know where to find sample MPS files.\n");
    exit(1);
  }
#endif
  if (argc>=2) mpsFileName = argv[1];
  int numMpsReadErrors = solver1.readMps(mpsFileName.c_str(),"");
  assert(numMpsReadErrors==0);

  // Strip off integer information and save
  int numberColumns = solver1.getNumCols();
  char * integer = new char[numberColumns];
  int i;
  for (i=0;i<numberColumns;i++) {
    if (solver1.isInteger(i)) {
      integer[i]=1;
      solver1.setContinuous(i);
    } else {
      integer[i]=0;
    }
  }
  // Pass to Cbc initialize defaults 
  CbcModel model(solver1);    
  CbcMain0(model);

  // Solve just to show there are no integers
  model.branchAndBound();
  // Set cutoff etc back in model and solver
  model.resetModel();
  // Solver was cloned so get it
  OsiSolverInterface * solver = model.solver();
  // Put back integers.  Here the user could do anything really
#define ADD_DIRECTLY
#ifndef ADD_DIRECTLY
  for (i=0;i<numberColumns;i++) {
    if (integer[i])
      solver->setInteger(i);
  }
#else
  CbcObject ** objects = new CbcObject * [ numberColumns];
  int n=0;
  for (i=0;i<numberColumns;i++) {
    if (integer[i]) {
      CbcSimpleIntegerDynamicPseudoCost * newObject =
	new CbcSimpleIntegerDynamicPseudoCost(&model,i);
      objects[n++]=newObject;
    }
  }
  model.addObjects(n,objects);
  for (i=0;i<n;i++)
    delete objects[i];
  delete [] objects;
#endif
  delete [] integer;
  /* Now go into code for standalone solver
     Could copy arguments and add -quit at end to be safe
     but this will do
  */
  if (argc>2) {
    CbcMain1(argc-1,argv+1,model);
  } else {
    const char * argv2[]={"driver3","-solve","-quit"};
    CbcMain1(3,argv2,model);
  }

  // Print solution if finished (could get from model.bestSolution() as well

  if (solver->getObjValue()*solver->getObjSense()<1.0e50) {
    
    const double * solution = solver->getColSolution();
    
    int iColumn;
    std::cout<<std::setiosflags(std::ios::fixed|std::ios::showpoint)<<std::setw(14);
    
    std::cout<<"--------------------------------------"<<std::endl;
#ifdef USE_OSI_NAMES
    
    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<<" "<<std::setw(8)<<setiosflags(std::ios::left)<<solver->getColName(iColumn)
		 <<resetiosflags(std::ios::adjustfield)<<std::setw(14)<<" "<<value<<std::endl;
    }
#else
    // names may not be in current solver - use original
//.........这里部分代码省略.........