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

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


}

/*
  Function for parameters that are enabled/disabled with a hint.
*/
int pushCbcOsiHint (CoinParam *param)

{
    assert (param != 0) ;
    CbcOsiParam *osiParam = dynamic_cast<CbcOsiParam *>(param) ;
    assert (osiParam != 0) ;
    OsiSolverInterface *osi = osiParam->obj() ;
    assert(osi != 0) ;
    /*
      Setup to return nonfatal/fatal error (1/-1) by default, so that all we need
      to do is correct to 0 (no error) if we're successful.
    */
    int retval ;
    if (CoinParamUtils::isInteractive()) {
        retval = 1 ;
    } else {
        retval = -1 ;
    }
    /*
      Set the sense for the hint.
    */
    std::string kwd = param->kwdVal() ;
    bool sense ;
    if (kwd == "off") {
        sense = false ;
    } else {
        sense = true ;
    }
    /*
      Grab the parameter code and translate to an OSI parameter key.
    */
    CbcOsiParam::CbcOsiParamCode code = osiParam->paramCode() ;
    OsiHintParam key ;
    switch (code) {
    case CbcOsiParam::PRESOLVE: {
        key = OsiDoPresolveInInitial ;
        break ;
    }
    case CbcOsiParam::SCALING: {
        key = OsiDoScale ;
        break ;
    }
    default: {
        std::cerr << "pushCbcOsiHint: no equivalent OsiHintParam for "
                  << "parameter code `" << code << "'." << std::endl ;
        retval = -1 ;
        break ;
    }
    }

    bool setOK = osi->setHintParam(key, sense, OsiHintDo) ;

    if (setOK) {
        return (0) ;
    } else {
        return (retval) ;
    }
}

void* problem_convert_to_osi(Problem *p)
{
    int i;
    double rowLb, rowUb;
    OsiSolverInterface *solver = new OsiClpSolverInterface();
    CoinBuild cb;

    solver->setIntParam(OsiNameDiscipline, 2);
    solver->messageHandler()->setLogLevel(0);
    solver->setHintParam(OsiDoReducePrint,true,OsiHintTry);

    for(i = 0; i < p->numCols; i++)
    {
        solver->addCol(0, NULL, NULL, p->colLb[i], p->colUb[i], p->objCoef[i]);
        solver->setColName(i, p->colName[i]);
        
        if(p->colType[i] == CONTINUOUS)
            solver->setContinuous(i);
        else 
            solver->setInteger(i);
    }

    for(i = 0; i < p->numRows; i++)
    {
        switch(p->rowSense[i])
        {
            case 'E':
                rowLb = p->rhs[i];
                rowUb = p->rhs[i];
            break;

            case 'L':
                rowLb = -p->infty;
                rowUb = p->rhs[i];
            break;

            case 'G':
                rowLb = p->rhs[i];
                rowUb = p->infty;
            break;

            default:
                fprintf(stderr, "Error: invalid type of constraint!\n");
                exit(EXIT_FAILURE);
        }

        cb.addRow(p->rowNElements[i], p->idxsByRow[i], p->coefsByRow[i], rowLb, rowUb);
    }

    solver->addRows(cb);
    
    for(i = 0; i < p->numRows; i++)
        solver->setRowName(i, p->rowName[i]);

    return solver;
}

int pushCbcOsiLogLevel (CoinParam *param)

{
    assert (param != 0) ;
    CbcOsiParam *osiParam = dynamic_cast<CbcOsiParam *>(param) ;
    assert (osiParam != 0) ;
    OsiSolverInterface *osi = osiParam->obj() ;
    assert(osi != 0) ;

    int lvl = param->intVal() ;
    /*
      Setup to return nonfatal/fatal error (1/-1) by default, so that all we need
      to do is correct to 0 (no error) if we're successful.
    */
    int retval ;
    if (CoinParamUtils::isInteractive()) {
        retval = 1 ;
    } else {
        retval = -1 ;
    }
    /*
      Now try to do the right thing with a hint. Harder to say -- assume that log
      level 1 is `normal'.
    */
    OsiHintStrength strength ;
    bool sense ;
    if (lvl < 1) {
        strength = OsiHintDo ;
        sense = true ;
    } else if (lvl == 1) {
        strength = OsiHintIgnore ;
        sense = true ;
    } else if (lvl == 2) {
        strength = OsiHintTry ;
        sense = false ;
    } else {
        strength = OsiHintDo ;
        sense = false ;
    }

    bool setOK = osi->setHintParam(OsiDoReducePrint, sense, strength) ;

    /*
      Recover the message handler and set the log level directly.
    */
    CoinMessageHandler *hndl = osi->messageHandler() ;
    assert (hndl != 0) ;
    hndl->setLogLevel(lvl) ;

    if (setOK) {
        return (0) ;
    } else {
        return (retval) ;
    }
}

    ILPSolverCbc::ILPSolverCbc()
    {
        // CbcModel assumes ownership over solver and deletes it in its destructor.
        OsiSolverInterface* solver = new OsiClpSolverInterface();

        // Output should come from CBC, not from CBCs solver.
        solver->messageHandler()->setLogLevel(0);
        solver->setHintParam(OsiDoReducePrint, 1);

        d_model.assignSolver(solver, true);
        set_default_parameters(this);
    }

//.........这里部分代码省略.........
  testingMessage( "Testing OsiGlpkSolverInterface\n" );
  OSIUNITTEST_CATCH_ERROR(OsiGlpkSolverInterfaceUnitTest(mpsDir,netlibDir), {}, "glpk", "osiglpk unittest");
#endif
  
#ifdef COIN_HAS_MSK
  testingMessage( "Testing OsiMskSolverInterface\n" );
  OSIUNITTEST_CATCH_ERROR(OsiMskSolverInterfaceUnitTest(mpsDir,netlibDir), {}, "mosek", "osimsk unittest");
#endif

#ifdef COIN_HAS_GRB
  testingMessage( "Testing OsiGrbSolverInterface\n" );
  OSIUNITTEST_CATCH_ERROR(OsiGrbSolverInterfaceUnitTest(mpsDir,netlibDir), {}, "gurobi", "osigrb unittest");
#endif

#ifdef COIN_HAS_SOPLEX
  testingMessage( "Testing OsiSpxSolverInterface\n" );
  OSIUNITTEST_CATCH_ERROR(OsiSpxSolverInterfaceUnitTest(mpsDir,netlibDir), {}, "soplex", "osispx unittest");
#endif

/*
  Each solver has run its specialised unit test. Check now to see if we need to
  run through the Netlib problems.
*/
  if (parms.find("-testOsiSolverInterface") != parms.end())
  {
    // Create vector of solver interfaces
    std::vector<OsiSolverInterface*> vecSi;
#   if COIN_HAS_XPR
    OsiSolverInterface * xprSi = new OsiXprSolverInterface;
    vecSi.push_back(xprSi);
#   endif
#   if COIN_HAS_CPX
    OsiSolverInterface * cpxSi = new OsiCpxSolverInterface;
    vecSi.push_back(cpxSi);
#   endif
#   if COIN_HAS_GLPK
    OsiSolverInterface * glpkSi = new OsiGlpkSolverInterface;
    glpkSi->setHintParam(OsiDoPresolveInInitial,true,OsiHintTry) ;
    glpkSi->setHintParam(OsiDoReducePrint,true,OsiHintDo) ;
    vecSi.push_back(glpkSi);
#   endif
#   if COIN_HAS_MSK
    OsiSolverInterface * MskSi = new OsiMskSolverInterface;
    vecSi.push_back(MskSi);
#   endif
#   if COIN_HAS_GRB
    OsiSolverInterface * grbSi = new OsiGrbSolverInterface;
    vecSi.push_back(grbSi);
#   endif
#   if COIN_HAS_SOPLEX
    OsiSolverInterface * spxSi = new OsiSpxSolverInterface;
    vecSi.push_back(spxSi);
#   endif
#   ifdef USETESTSOLVER
/*
  The test solver is normally Vol, which can't do any of the netlib problems.
  So let's not try.
*/
    {
      std::string solverName ;
      OsiSolverInterface * testSi = new OsiTestSolverInterface;
      testSi->getStrParam(OsiSolverName,solverName) ;
      if (solverName != "vol")
      { vecSi.push_back(testSi); }
      else
      { testingMessage("Test solver vol cannot do Netlib. Skipping.\n") ; }
    }
#   endif

    if (vecSi.size() > 0)
    { testingMessage( "Testing OsiSolverInterface on Netlib problems.\n" );
      OSIUNITTEST_CATCH_ERROR(OsiSolverInterfaceMpsUnitTest(vecSi,netlibDir), {}, "osi", "netlib unittest");
    }

    unsigned int i;
    for (i=0; i<vecSi.size(); i++)
      delete vecSi[i];
  }
  else {
    testingMessage( "***Skipped Testing of OsiSolverInterface on Netlib problems***\n" );
    testingMessage( "***use -testOsiSolverInterface to run them.***\n" );
  }

/*
  We're done. Report on the results.
*/
  std::cout.flush();
  outcomes.print();

  int nerrors;
  int nerrors_expected;
  outcomes.getCountBySeverity(TestOutcome::ERROR, nerrors, nerrors_expected);

  if (nerrors > nerrors_expected)
    std::cerr << "Tests completed with " << nerrors - nerrors_expected << " unexpected errors." << std::endl ;
  else
    std::cerr << "All tests completed successfully\n";

  return nerrors - nerrors_expected;
}

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

//.........这里部分代码省略.........
       tmp = binarySearch(0, lCols - 1, index2, lColIndices);
       if(tmp > -1)
       row.insert(tmp, matElements[j]);
       }
       newMat->appendRow(row);
       }
     */
     for(i = 0; i < lRows; i++){
	CoinPackedVector row;
	index1 = lRowIndices[i];
	for(j = 0; j < lCols; j++){
	   index2 = lColIndices[j];
	   tmp = matrix->getCoefficient(index1, index2);
	   row.insert(j, tmp);
	}
	newMat->appendRow(row);
     }

     /*
       nSolver->assignProblem(newMat, colLb, colUb,
       objCoeffs, rowLb, rowUb);
     */
     
     nSolver->loadProblem(*newMat, colLb, colUb,
			  objCoeffs, rowLb, rowUb);
     
     for(i = 0; i < intCnt; i++){
	nSolver->setInteger(integerVars[i]);
     }
     //nSolver->setInteger(integerVars, intCnt);
     
     nSolver->setObjSense(objSense); //1 min; -1 max
     
     nSolver->setHintParam(OsiDoReducePrint, true, OsiHintDo);

#if 0
     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);
     }
#endif
     delete [] integerVars;

  }else{
     nSolver = solver_;
  }

#define SYM_VERSION_IS_WS strcmp(SYMPHONY_VERSION, "WS")  

#if SYMPHONY_VERSION_IS_WS
  if (feasCheckSolver == "SYMPHONY" && probType == 1 && warmStartLL &&
      !newOsi && doDualFixing){ //Interdiction

     /** Get upper bound from best known (feasible) lower level solution and try 
	 to fix additional variables by sensitivity analysis **/

//#############################################################################
void 
MibSBilevel::checkBilevelFeasiblity(bool isRoot)
{
  
  int cutStrategy =
    model_->MibSPar_->entry(MibSParams::cutStrategy);

  bool warmStartLL =
    model_->MibSPar_->entry(MibSParams::warmStartLL);

  int maxThreadsLL =
    model_->MibSPar_->entry(MibSParams::maxThreadsLL);

  int whichCutsLL =
    model_->MibSPar_->entry(MibSParams::whichCutsLL);

  int probType =
    model_->MibSPar_->entry(MibSParams::bilevelProblemType);

  std::string feasCheckSolver =
     model_->MibSPar_->entry(MibSParams::feasCheckSolver);

  if (warmStartLL && (feasCheckSolver == "SYMPHONY") && solver_){
     solver_ = setUpModel(model_->getSolver(), false);
  }else{
     if (solver_){
	delete solver_;
     }
     solver_ = setUpModel(model_->getSolver(), true);
  }

  OsiSolverInterface *lSolver = solver_;

  //CoinWarmStart * ws = getWarmStart();
  //if (ws != NULL){
  //   lSolver->setWarmStart(ws);
  //}
  //delete ws;

  if(1)
    lSolver->writeLp("lowerlevel");

  if (feasCheckSolver == "Cbc"){
    dynamic_cast<OsiCbcSolverInterface *> 
      (lSolver)->getModelPtr()->messageHandler()->setLogLevel(0);
  }else if (feasCheckSolver == "SYMPHONY"){
     //dynamic_cast<OsiSymSolverInterface *> 
     // (lSolver)->setSymParam("prep_level", -1);
    
     sym_environment *env = dynamic_cast<OsiSymSolverInterface *> 
	(lSolver)->getSymphonyEnvironment();

     if (warmStartLL){
	sym_set_int_param(env, "keep_warm_start", TRUE);
	if (probType == 1){ //Interdiction
	   sym_set_int_param(env, "should_use_rel_br", FALSE);
	   sym_set_int_param(env, "use_hot_starts", FALSE);
	   sym_set_int_param(env, "should_warmstart_node", TRUE);
	   sym_set_int_param(env, "sensitivity_analysis", TRUE);
	   sym_set_int_param(env, "sensitivity_bounds", TRUE);
	   sym_set_int_param(env, "set_obj_upper_lim", FALSE);
	}
     }
     //Always uncomment for debugging!!
     sym_set_int_param(env, "do_primal_heuristic", FALSE);
     sym_set_int_param(env, "verbosity", -2);
     sym_set_int_param(env, "prep_level", -1);
     sym_set_int_param(env, "max_active_nodes", maxThreadsLL);
     sym_set_int_param(env, "tighten_root_bounds", FALSE);
     sym_set_int_param(env, "max_sp_size", 100);
     sym_set_int_param(env, "do_reduced_cost_fixing", FALSE);
     if (whichCutsLL == 0){
	sym_set_int_param(env, "generate_cgl_cuts", FALSE);
     }else{
	sym_set_int_param(env, "generate_cgl_gomory_cuts", GENERATE_DEFAULT);
     }
     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
  }
//.........这里部分代码省略.........

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

}

//.........这里部分代码省略.........
	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 main(int argc, char **argv)
{
    VRPH_version();

    int i, j, k, n, status, num_attempts, *sol_buff, *IP_sol_buff;
    char in_file[200];
    double lambda, best_heur_sol=VRP_INFINITY;
    bool first_sol=false, bootstrap=false;;
    VRPSolution *fresh_solution;
    OsiSolverInterface *si;
    const double *x;
    int last_num_cols=0, route_id=0;
    time_t start, stop;
    int *orderings[MAX_ROUTES];
    for(i=0;i<MAX_ROUTES;i++)
        orderings[i]=NULL;
    

    // Set timing counters to 0
    heur_time=mip_time=0;

    // Check arguments
    if(argc<5)
    {
        fprintf(stderr,"Usage: %s -f input_file -n num_runs [-v,-b,-c max_columns -d cols_to_delete]\n",
            argv[0]);
        fprintf(stderr,"\t Will solve the problem num_solutions times and add the routes\n");
        fprintf(stderr,"\t to a set partitioning problem.\n");
        fprintf(stderr,"\t Other options:\n");
        fprintf(stderr,"\t -v runs in verbose mode\n");
        fprintf(stderr,"\t -b will use bootstrapping where we send the set partitioning\n"
                       "\t    solution back to the metaheuristic solver\n");
        fprintf(stderr,"\t -c max_columns will allow this many active columns/variables in the IP.\n");
        fprintf(stderr,"\t    Default value is max_columns=500\n");
        fprintf(stderr,"\t -d num_cols_to_delete will delete this many columns once we have too many\n");
        fprintf(stderr,"\t    in the IP. Default value is num_cols_to_delete=100\n");
        exit(-1);
    }

    // Set defaults
    verbose=false;
    max_columns=500;
    num_cols_to_delete=100;

    // Parse command line
    for(i=0;i<argc;i++)
    {
        if(strcmp(argv[i],"-f")==0)
            strcpy(in_file,argv[i+1]);
        if(strcmp(argv[i],"-n")==0)
            num_attempts=atoi(argv[i+1]);
        if(strcmp(argv[i],"-v")==0)
            verbose=true;
         if(strcmp(argv[i],"-b")==0)
            bootstrap=true;
         if(strcmp(argv[i],"-c")==0)
             max_columns=atoi(argv[i+1]);
         if(strcmp(argv[i],"-d")==0)
             num_cols_to_delete=atoi(argv[i+1]);
    }

    // This is the # of non-VRPH_DEPOT nodes
    n=VRPGetDimension(in_file);
    // This will be used to import/export solutions
    fresh_solution = new VRPSolution(n);
    // Create buffers for importing solutions
    sol_buff= new int[n+2];
    IP_sol_buff = new int[n+2];

    // Declare an OSI interface
    si=new OsiGlpkSolverInterface;
    si->setIntParam(OsiNameDiscipline,2);

    for(i=0;i<n;i++)
    {
        si->addRow(0,NULL,NULL,1,1);
    }

    // Declare a VRP of the right size and import the file
    VRP V(n);
    ClarkeWright CW(n);
    VRPRoute route(n);
    V.read_TSPLIB_file(in_file);
    // Set up a "route warehouse" to store the routes to be added to the IP
    V.route_wh=new VRPRouteWarehouse(HASH_TABLE_SIZE);
  
    // Set up a minimization problem
    si->setObjSense(1);
    // Set to error only output
    si->setHintParam(OsiDoReducePrint,true, OsiHintDo);
    // Unfortunately GLPK still prints out something regarding the conflict graph

    for(i=0;i<num_attempts;i++)
    {
        if(i==0 || !bootstrap) 
        {
            lambda=.5+1.5*lcgrand(0);
            // Start with a clean VRP object
            V.reset();
            CW.Construct(&V, lambda, false);
//.........这里部分代码省略.........

OsiSolverInterface* cpropagation_preprocess(CPropagation *cp, int nindexes[])
{
    if(cp->varsToFix == 0)
    {
        /* printf("There are no variables to remove from the problem!\n"); */
        return NULL; /* returns a pointer to original solver */
    }

    const double *colLb = problem_vars_lower_bound(cp->problem), *colUb = problem_vars_upper_bound(cp->problem);
    const double *objCoef = problem_vars_obj_coefs(cp->problem);
    const char *ctype = problem_vars_type(cp->problem);

    double sumFixedObj = 0.0; /* stores the sum of objective coefficients of all variables fixed to 1 */

    OsiSolverInterface *preProcSolver = new OsiClpSolverInterface();
    preProcSolver->setIntParam(OsiNameDiscipline, 2);
    preProcSolver->messageHandler()->setLogLevel(0);
    preProcSolver->setHintParam(OsiDoReducePrint,true,OsiHintTry);
    //preProcSolver->setObjName(cp->solver->getObjName());

    for(int i = 0, j = 0; i < problem_num_cols(cp->problem); i++)
    {
        nindexes[i] = -1;
        if(cp->isToFix[i] == UNFIXED)
        {
            preProcSolver->addCol(0, NULL, NULL, colLb[i], colUb[i], objCoef[i]);
            preProcSolver->setColName(j, problem_var_name(cp->problem, i));
            if(problem_var_type(cp->problem, i) == CONTINUOUS)
                preProcSolver->setContinuous(j);
            else 
                preProcSolver->setInteger(j);
            nindexes[i] = j++;
        }
        else if(cp->isToFix[i] == ACTIVATE)
            sumFixedObj += objCoef[i];
    }

    if(fabs(sumFixedObj) > EPS)
    {
        /* adding a variable with cost equals to the sum of all coefficients of variables fixed to 1 */
        preProcSolver->addCol(0, NULL, NULL, 1.0, 1.0, sumFixedObj);
        preProcSolver->setColName(preProcSolver->getNumCols()-1, "sumFixedObj");
        preProcSolver->setInteger(preProcSolver->getNumCols()-1);
    }

    for(int idxRow = 0; idxRow < problem_num_rows(cp->problem); idxRow++)
    {
        const int nElements = problem_row_size(cp->problem, idxRow);
        const int *idxs = problem_row_idxs(cp->problem, idxRow);
        const double *coefs = problem_row_coefs(cp->problem, idxRow);
        vector< int > vidx; vidx.reserve(problem_num_cols(cp->problem));
        vector< double > vcoef; vcoef.reserve(problem_num_cols(cp->problem));
        double activeCoefs = 0.0;

        for(int i = 0; i < nElements; i++)
        {
            if(cp->isToFix[idxs[i]] == UNFIXED)
            {
                assert(nindexes[idxs[i]] >= 0 && nindexes[idxs[i]] < problem_num_cols(cp->problem));
                vidx.push_back(nindexes[idxs[i]]);
                vcoef.push_back(coefs[i]);
            }
            else if(cp->isToFix[idxs[i]] == ACTIVATE)
            	activeCoefs += coefs[i];
        }

        if(!vidx.empty())
        {
        	double rlb, rub;
        	const char sense = problem_row_sense(cp->problem, idxRow);
        	
        	if(sense == 'E')
            {
                rlb = problem_row_rhs(cp->problem, idxRow) - activeCoefs;
                rub = problem_row_rhs(cp->problem, idxRow) - activeCoefs;
            }
        	else if(sense == 'L')
            {
                rlb = preProcSolver->getInfinity();
                rub = problem_row_rhs(cp->problem, idxRow) - activeCoefs;
            }
        	else if(sense == 'G')
            {
                rlb = problem_row_rhs(cp->problem, idxRow) - activeCoefs;
                rub = preProcSolver->getInfinity();
            }
        	else
        	{
        		fprintf(stderr, "Error: invalid type of constraint!\n");
        		exit(EXIT_FAILURE);
        	}

        	preProcSolver->addRow((int)vcoef.size(), &vidx[0], &vcoef[0], rlb, rub);
            preProcSolver->setRowName(idxRow, problem_row_name(cp->problem, idxRow));
        }
	}

    return preProcSolver;
}

int doBaCParam (CoinParam *param)

{
    assert (param != 0) ;
    CbcGenParam *genParam = dynamic_cast<CbcGenParam *>(param) ;
    assert (genParam != 0) ;
    CbcGenCtlBlk *ctlBlk = genParam->obj() ;
    assert (ctlBlk != 0) ;
    CbcModel *model = ctlBlk->model_ ;
    assert (model != 0) ;
    /*
      Setup to return nonfatal/fatal error (1/-1) by default.
    */
    int retval ;
    if (CoinParamUtils::isInteractive()) {
        retval = 1 ;
    } else {
        retval = -1 ;
    }
    ctlBlk->setBaBStatus(CbcGenCtlBlk::BACAbandon, CbcGenCtlBlk::BACmInvalid,
                         CbcGenCtlBlk::BACwNotStarted, false, 0) ;
    /*
      We ain't gonna do squat without a good model.
    */
    if (!ctlBlk->goodModel_) {
        std::cout << "** Current model not valid!" << std::endl ;
        return (retval) ;
    }
    /*
      Start the clock ticking.
    */
    double time1 = CoinCpuTime() ;
    double time2 ;
    /*
      Create a clone of the model which we can modify with impunity. Extract
      the underlying solver for convenient access.
    */
    CbcModel babModel(*model) ;
    OsiSolverInterface *babSolver = babModel.solver() ;
    assert (babSolver != 0) ;
# if CBC_TRACK_SOLVERS > 0
    std::cout
        << "doBaCParam: initial babSolver is "
        << std::hex << babSolver << std::dec
        << ", log level " << babSolver->messageHandler()->logLevel()
        << "." << std::endl ;
# endif
    /*
      Solve the root relaxation. Bail unless it solves to optimality.
    */
    if (!solveRelaxation(&babModel)) {
        ctlBlk->setBaBStatus(&babModel, CbcGenCtlBlk::BACwBareRoot) ;
        return (0) ;
    }
# if COIN_CBC_VERBOSITY > 0
    std::cout
        << "doBaCParam: initial relaxation z = "
        << babSolver->getObjValue() << "." << std::endl ;
# endif
    /*
      Are we up for fixing variables based on reduced cost alone?
    */
    if (ctlBlk->djFix_.action_ == true) {
        reducedCostHack(babSolver, ctlBlk->djFix_.threshold_) ;
    }
    /*
      Time to consider preprocessing. We'll do a bit of setup before getting to
      the meat of the issue.

      preIppSolver will hold a clone of the unpreprocessed constraint system.
      We'll need it when we postprocess. ippSolver holds the preprocessed
      constraint system.  Again, we clone it and give the clone to babModel for
      B&C. Presumably we need an unmodified copy of the preprocessed system to
      do postprocessing, but the copy itself is hidden inside the preprocess
      object.
    */
    OsiSolverInterface *preIppSolver = 0 ;
    CglPreProcess ippObj ;
    bool didIPP = false ;

    int numberChanged = 0 ;
    int numberOriginalColumns = babSolver->getNumCols() ;
    CbcGenCtlBlk::IPPControl ippAction = ctlBlk->getIPPAction() ;

    if (!(ippAction == CbcGenCtlBlk::IPPOff ||
            ippAction == CbcGenCtlBlk::IPPStrategy)) {
        double timeLeft = babModel.getMaximumSeconds() ;
        preIppSolver = babSolver->clone() ;
        OsiSolverInterface *ippSolver ;
#   if CBC_TRACK_SOLVERS > 0
        std::cout
            << "doBaCParam: clone made prior to IPP is "
            << std::hex << preIppSolver << std::dec
            << ", log level " << preIppSolver->messageHandler()->logLevel()
            << "." << std::endl ;
#   endif

        preIppSolver->setHintParam(OsiDoInBranchAndCut, true, OsiHintDo) ;
        ippObj.messageHandler()->setLogLevel(babModel.logLevel()) ;

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