C++ OsiClpSolverInterface类代码示例

double
CbcGeneralBranchingObject::branch()
{
    double cutoff = model_->getCutoff();
    //printf("GenB %x whichNode %d numberLeft %d which %d\n",
    // this,whichNode_,numberBranchesLeft(),branchIndex());
    if (whichNode_ < 0) {
        assert (node_);
        bool applied = false;
        while (numberBranchesLeft()) {
            int which = branchIndex();
            decrementNumberBranchesLeft();
            CbcSubProblem * thisProb = subProblems_ + which;
            if (thisProb->objectiveValue_ < cutoff) {
                //printf("branch %x (sub %x) which now %d\n",this,
                //     subProblems_,which);
                OsiSolverInterface * solver = model_->solver();
                thisProb->apply(solver);
                OsiClpSolverInterface * clpSolver
                = dynamic_cast<OsiClpSolverInterface *> (solver);
                assert (clpSolver);
                // Move status to basis
                clpSolver->setWarmStart(NULL);
                //ClpSimplex * simplex = clpSolver->getModelPtr();
                node_->setObjectiveValue(thisProb->objectiveValue_);
                node_->setSumInfeasibilities(thisProb->sumInfeasibilities_);
                node_->setNumberUnsatisfied(thisProb->numberInfeasibilities_);
                applied = true;
                doingDoneBranch = true;
                break;
            } else if (numberBranchesLeft()) {
                node_->nodeInfo()->branchedOn() ;
            }
        }
        if (!applied) {
            // no good one
            node_->setObjectiveValue(cutoff + 1.0e20);
            node_->setSumInfeasibilities(1.0);
            node_->setNumberUnsatisfied(1);
            assert (whichNode_ < 0);
        }
    } else {
        decrementNumberBranchesLeft();
        CbcSubProblem * thisProb = subProblems_ + whichNode_;
        assert (thisProb->objectiveValue_ < cutoff);
        OsiSolverInterface * solver = model_->solver();
        thisProb->apply(solver);
        //OsiClpSolverInterface * clpSolver
        //= dynamic_cast<OsiClpSolverInterface *> (solver);
        //assert (clpSolver);
        // Move status to basis
        //clpSolver->setWarmStart(NULL);
    }
    return 0.0;
}

  void 
  LinearCutsGenerator::generateCuts(const OsiSolverInterface &solver, OsiCuts &cs,
                     const CglTreeInfo info) const {

    //const OsiTMINLPInterface * tmp = dynamic_cast<const OsiTMINLPInterface *>(&solver);
    OsiTMINLPInterface * nlp = dynamic_cast<OsiTMINLPInterface *>(solver.clone());//const_cast<OsiTMINLPInterface *>(tmp);
    assert(nlp);
    OuterApprox oa;
    //si.writeMps("toto");
    int numberRows = nlp->getNumRows();
    for(int i = 0 ; i < 5 ; i++){
      nlp->resolve();
      OsiClpSolverInterface si;
      oa(*nlp, &si, solver.getColSolution(), true); 
      si.resolve();
      OsiCuts cuts;
      for(std::list<Coin::SmartPtr<CuttingMethod> >::const_iterator i = methods_.begin() ;
          i != methods_.end() ; i++){
         (*i)->cgl->generateCuts(si, cuts, info);
      }
      std::vector<OsiRowCut *> mycuts(cuts.sizeRowCuts());
      for(int i = 0 ; i < cuts.sizeRowCuts() ; i++){
        mycuts[i] = cuts.rowCutPtr(i);
        cs.insert(*mycuts[i]);
      }
      nlp->applyRowCuts(mycuts.size(), const_cast<const OsiRowCut **> (&mycuts[0]));
    }

    // Take off slack cuts
    std::vector<int> kept;
    int numberRowsNow = nlp->getNumRows();
    int * del = new int [numberRowsNow-numberRows];
    nlp->resolve();
    
    const double * activity = nlp->getRowActivity();
    const double * lb = nlp->getRowLower();
    const double * ub = nlp->getRowUpper();
    CoinRelFltEq eq(1e-06);
    //int nDelete=0;
    for (int i=numberRowsNow -1;i>=numberRows;i--) {
      if ( !(eq(activity[i], lb[i]) || eq(activity[i], ub[i])) )
        cs.eraseRowCut(i - numberRows);
    }
    delete [] del;
    delete nlp;
  }

    // nothing to inherit from GC really
    ColoredGraph solve(const Graph& gr) override {

        CoinModel coinModel;
        for (auto i = 0; i < gr.nodeCount(); ++i) {
            coinModel.addCol(0, nullptr, nullptr, 0, gr.nodeCount()-1, 1, nullptr, true);
        }
        OsiClpSolverInterface solver;
        solver.loadFromCoinModel(coinModel);

        CbcModel model(solver);
        model.setLogLevel(0);
        model.passInEventHandler(make_unique<GC_LP_EventHandler>(recoveryPath).get());

        if (use_heuristic) {
            GC_LP_Heuristic heuristic;
            model.addHeuristic(&heuristic);
        }

        AddRules(model, gr);

        if (use_parallel) {
            model.setNumberThreads(std::thread::hardware_concurrency());
        }

        if (max_seconds != 0) model.setMaximumSeconds(max_seconds);

        model.initialSolve();
        model.branchAndBound();

        if (model.maximumSecondsReached()) Println(cout, "max seconds reached");
        if (model.isSecondsLimitReached()) Println(cout, "seconds limit reached");

        seconds_passed_ = model.getCurrentSeconds();
        iterations_passed_ = model.getIterationCount();

        const double *solution = model.bestSolution();

        if (solution == nullptr) {
            vector<Color> colors(gr.nodeCount());
            iota(colors.begin(), colors.end(), 0);
            return {gr, colors};
        }
        return ColoredGraph(gr, {solution, solution+gr.nodeCount()});
    }

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

	try{
		// Set up lp solver
		OsiClpSolverInterface lpSolver;
		lpSolver.getModelPtr()->setDualBound(1.0e10);
		lpSolver.messageHandler()->setLogLevel(0);
	
		// Create BLIS model 
		BlisModel model;
		model.setSolver(&lpSolver);
	
#ifdef  COIN_HAS_MPI
		AlpsKnowledgeBrokerMPI broker(argc, argv, model);
#else
		AlpsKnowledgeBrokerSerial broker(argc, argv, model); 
#endif

		// Search for best solution
		broker.search(&model);
	
		// Report the best solution found and its ojective value
		broker.printBestSolution();
 }	

	catch(CoinError& er) {
 	std::cerr << "\nBLIS ERROR: \"" << er.message() 
		  << "\""<< std::endl
		  << "             from function \"" << er.methodName()
		  << "\""<< std::endl
		  << "             from class \"" << er.className()
		  << "\"" << std::endl;
	}
  catch(...) {
		std::cerr << "Something went wrong!" << std::endl;
  }
    
    
  return 0;
}

// Infeasibility - large is 0.5
double
CbcGeneralDepth::infeasibility(const OsiBranchingInformation * /*info*/,
                               int &/*preferredWay*/) const
{
    whichSolution_ = -1;
    // should use genuine OsiBranchingInformation usefulInfo = model_->usefulInformation();
    // for now assume only called when correct
    //if (usefulInfo.depth_>=4&&!model_->parentModel()
    //     &&(usefulInfo.depth_%2)==0) {
    if (true) {
        OsiSolverInterface * solver = model_->solver();
        OsiClpSolverInterface * clpSolver
        = dynamic_cast<OsiClpSolverInterface *> (solver);
        if (clpSolver) {
	  if ((model_->moreSpecialOptions()&33554432)==0) {
            ClpNodeStuff * info = nodeInfo_;
            info->integerTolerance_ = model_->getIntegerTolerance();
            info->integerIncrement_ = model_->getCutoffIncrement();
            info->numberBeforeTrust_ = model_->numberBeforeTrust();
            info->stateOfSearch_ = model_->stateOfSearch();
            // Compute "small" change in branch
            int nBranches = model_->getIntParam(CbcModel::CbcNumberBranches);
            if (nBranches) {
                double average = model_->getDblParam(CbcModel::CbcSumChange) / static_cast<double>(nBranches);
                info->smallChange_ =
                    CoinMax(average * 1.0e-5, model_->getDblParam(CbcModel::CbcSmallestChange));
                info->smallChange_ = CoinMax(info->smallChange_, 1.0e-8);
            } else {
                info->smallChange_ = 1.0e-8;
            }
            int numberIntegers = model_->numberIntegers();
            double * down = new double[numberIntegers];
            double * up = new double[numberIntegers];
            int * priority = new int[numberIntegers];
            int * numberDown = new int[numberIntegers];
            int * numberUp = new int[numberIntegers];
            int * numberDownInfeasible = new int[numberIntegers];
            int * numberUpInfeasible = new int[numberIntegers];
            model_->fillPseudoCosts(down, up, priority, numberDown, numberUp,
                                    numberDownInfeasible, numberUpInfeasible);
            info->fillPseudoCosts(down, up, priority, numberDown, numberUp,
                                  numberDownInfeasible,
                                  numberUpInfeasible, numberIntegers);
            info->presolveType_ = 1;
            delete [] down;
            delete [] up;
            delete [] numberDown;
            delete [] numberUp;
            delete [] numberDownInfeasible;
            delete [] numberUpInfeasible;
            bool takeHint;
            OsiHintStrength strength;
            solver->getHintParam(OsiDoReducePrint, takeHint, strength);
            ClpSimplex * simplex = clpSolver->getModelPtr();
            int saveLevel = simplex->logLevel();
            if (strength != OsiHintIgnore && takeHint && saveLevel == 1)
                simplex->setLogLevel(0);
            clpSolver->setBasis();
            whichSolution_ = simplex->fathomMany(info);
            //printf("FAT %d nodes, %d iterations\n",
            //info->numberNodesExplored_,info->numberIterations_);
            //printf("CbcBranch %d rows, %d columns\n",clpSolver->getNumRows(),
            //     clpSolver->getNumCols());
            model_->incrementExtra(info->numberNodesExplored_,
                                   info->numberIterations_);
            // update pseudo costs
            double smallest = 1.0e50;
            double largest = -1.0;
            OsiObject ** objects = model_->objects();
#ifndef NDEBUG
            const int * integerVariable = model_->integerVariable();
#endif
            for (int i = 0; i < numberIntegers; i++) {
#ifndef NDEBUG
                CbcSimpleIntegerDynamicPseudoCost * obj =
                    dynamic_cast <CbcSimpleIntegerDynamicPseudoCost *>(objects[i]) ;
                assert (obj && obj->columnNumber() == integerVariable[i]);
#else
                CbcSimpleIntegerDynamicPseudoCost * obj =
                    static_cast <CbcSimpleIntegerDynamicPseudoCost *>(objects[i]) ;
#endif
                if (info->numberUp_[i] > 0) {
                    if (info->downPseudo_[i] > largest)
                        largest = info->downPseudo_[i];
                    if (info->downPseudo_[i] < smallest)
                        smallest = info->downPseudo_[i];
                    if (info->upPseudo_[i] > largest)
                        largest = info->upPseudo_[i];
                    if (info->upPseudo_[i] < smallest)
                        smallest = info->upPseudo_[i];
                    obj->updateAfterMini(info->numberDown_[i],
                                         info->numberDownInfeasible_[i],
                                         info->downPseudo_[i],
                                         info->numberUp_[i],
                                         info->numberUpInfeasible_[i],
                                         info->upPseudo_[i]);
                }
            }
            //printf("range of costs %g to %g\n",smallest,largest);
//.........这里部分代码省略.........

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

void ModelScenario(const char * const name )
{
	OsiClpSolverInterface *osiClp1 = new OsiClpSolverInterface();
	double INF=osiClp1->getInfinity();

	// example of direct interfaces for scenario generation

    /* Model dimensions */
    int nels=44; // ncol=27, nrow=9

	/* Sparse matrix data...organized by row */
    int mrow[]={ 0, 0, 0, 0, 0,
		1, 1, 1, 1,
		2, 2, 2,
		3, 3, 3, 3, 3,
		4, 4, 4, 4,
		5, 5, 5, 5, 5, 5,
		6, 6, 6, 6, 6,
		7, 7, 7, 7, 7, 7,
		8, 8, 8, 8, 8, 8 };
	  int mcol[]={ 0, 1, 2, 3, 4,
		5, 6, 7, 8,
		9,10, 11,
		12, 13, 14, 15, 16,
		0,        12, 17, 18,
		1, 5, 9,  13, 19, 20,
		2, 6,     14, 21, 22,
		3, 7, 10, 15, 23, 24,
		4, 8, 11, 16, 25, 26 };

    double dels[] = { 1.0, 1.0, 1.0, 1.0, 1.0,
		1.0, 1.0, 1.0, 1.0,
		1.0, 1.0, 1.0,
		1.0, 1.0, 1.0, 1.0, 1.0,
		16.0,              9.0, -1.0, 1.0,
		15.0, 10.0,  5.0, 11.0, -1.0, 1.0,
		28.0, 14.0,       22.0, -1.0, 1.0,
		23.0, 15.0,  7.0, 17.0, -1.0, 1.0,
		81.0, 57.0, 29.0, 55.0, -1.0, 1.0 };

    /* Objective */
    double dobj[]={ 18.0, 21.0, 18.0, 16.0, 10.0, 15.0, 16.0, 14.0, 9.0,
		10.0,  9.0,  6.0, 17.0, 16.0, 17.0, 15.0, 10.0, 0.0,
		13.0,  0.0, 13.0,  0.0,  7.0,  0.0,  7.0,  0.0, 1.0 };

    /* Column bounds */
    double dclo[]={ 0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,
		0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,
		0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0 };



    double dcup[]={ INF,  INF,  INF,  INF,  INF,  INF,  INF,  INF,  INF,
		INF,  INF,  INF,  INF,  INF,  INF,  INF,  INF,  INF,
		INF,  INF,  INF,  INF,  INF,  INF,  INF,  INF,  INF };

    /* Row bounds */
    double drlo[]={ -INF, -INF, -INF, -INF,  0.0, 4.0, 0.0, 8.0, 10.0 };
    double drup[]={ 10.0, 19.0, 25.0, 15.0,  0.0, 7.0, 0.0, 8.0, 90.0 };

    /* Stages */
    int nstg=2;
    int n_first_stg_rows=4;
    int rstg[]={ 0,0,0,0,1,1,1,1,1 };
    int cstg[]={ 0,0,0,0,0,0,0,0,0,
		0,0,0,0,0,0,0,0,1,
		1,1,1,1,1,1,1,1,1 };

    /* Stochastic data */
    int nindp=5;
    int nsamp[]={ 5, 2, 5, 5, 3 };
    double demand[]={ 200, 220, 250, 270, 300,
		50, 150,
		140, 160, 180, 200, 220,
		10, 50, 80, 100, 340,
		580, 600, 620 };
    double dprobs[]={ 0.2, 0.05, 0.35, 0.2, 0.2,
		0.3, 0.7,
		0.1, 0.2, 0.4, 0.2, 0.1,
		0.2, 0.2, 0.3, 0.2, 0.1,
		0.1, 0.8, 0.1 };

    /* local variables */
    int ns=1,ii,iii,jj,*indx,*incr;
    double dp=1.0;

    for (ii=0;ii<nindp;ii++) ns *= nsamp[ii];     /* Compute number of scenarios */


	// initialize SmiModel
	SmiScnModel *smiModel = new SmiScnModel();
	smiModel->setOsiSolverHandle(*osiClp1);

	// set core model using Osi interface
	OsiClpSolverInterface ocsi;
	ocsi.loadProblem(CoinPackedMatrix( 1,mrow,mcol,dels,nels),dclo,dcup,dobj,drlo,drup);
	SmiCoreData *osiCore = new SmiCoreData(&ocsi,nstg,cstg,rstg);

	cout << "ModelScenario: Created CoreData" << endl;

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

int CbcHeuristicDive::reducedCostFix (OsiSolverInterface* solver)

{
    //return 0; // temp
#ifndef JJF_ONE
    if (!model_->solverCharacteristics()->reducedCostsAccurate())
        return 0; //NLP
#endif
    double cutoff = model_->getCutoff() ;
    if (cutoff > 1.0e20)
        return 0;
#ifdef DIVE_DEBUG
    std::cout << "cutoff = " << cutoff << std::endl;
#endif
    double direction = solver->getObjSense() ;
    double gap = cutoff - solver->getObjValue() * direction ;
    gap *= 0.5; // Fix more
    double tolerance;
    solver->getDblParam(OsiDualTolerance, tolerance) ;
    if (gap <= 0.0)
        gap = tolerance; //return 0;
    gap += 100.0 * tolerance;
    double integerTolerance = model_->getDblParam(CbcModel::CbcIntegerTolerance);

    const double *lower = solver->getColLower() ;
    const double *upper = solver->getColUpper() ;
    const double *solution = solver->getColSolution() ;
    const double *reducedCost = solver->getReducedCost() ;

    int numberIntegers = model_->numberIntegers();
    const int * integerVariable = model_->integerVariable();

    int numberFixed = 0 ;

# ifdef COIN_HAS_CLP
    OsiClpSolverInterface * clpSolver
    = dynamic_cast<OsiClpSolverInterface *> (solver);
    ClpSimplex * clpSimplex = NULL;
    if (clpSolver)
        clpSimplex = clpSolver->getModelPtr();
# endif
    for (int i = 0 ; i < numberIntegers ; i++) {
        int iColumn = integerVariable[i] ;
        double djValue = direction * reducedCost[iColumn] ;
        if (upper[iColumn] - lower[iColumn] > integerTolerance) {
            if (solution[iColumn] < lower[iColumn] + integerTolerance && djValue > gap) {
#ifdef COIN_HAS_CLP
                // may just have been fixed before
                if (clpSimplex) {
                    if (clpSimplex->getColumnStatus(iColumn) == ClpSimplex::basic) {
#ifdef COIN_DEVELOP
                        printf("DJfix %d has status of %d, dj of %g gap %g, bounds %g %g\n",
                               iColumn, clpSimplex->getColumnStatus(iColumn),
                               djValue, gap, lower[iColumn], upper[iColumn]);
#endif
                    } else {
                        assert(clpSimplex->getColumnStatus(iColumn) == ClpSimplex::atLowerBound ||
                               clpSimplex->getColumnStatus(iColumn) == ClpSimplex::isFixed);
                    }
                }
#endif
                solver->setColUpper(iColumn, lower[iColumn]) ;
                numberFixed++ ;
            } else if (solution[iColumn] > upper[iColumn] - integerTolerance && -djValue > gap) {
#ifdef COIN_HAS_CLP
                // may just have been fixed before
                if (clpSimplex) {
                    if (clpSimplex->getColumnStatus(iColumn) == ClpSimplex::basic) {
#ifdef COIN_DEVELOP
                        printf("DJfix %d has status of %d, dj of %g gap %g, bounds %g %g\n",
                               iColumn, clpSimplex->getColumnStatus(iColumn),
                               djValue, gap, lower[iColumn], upper[iColumn]);
#endif
                    } else {
                        assert(clpSimplex->getColumnStatus(iColumn) == ClpSimplex::atUpperBound ||
                               clpSimplex->getColumnStatus(iColumn) == ClpSimplex::isFixed);
                    }
                }
#endif
                solver->setColLower(iColumn, upper[iColumn]) ;
                numberFixed++ ;
            }
        }
    }
    return numberFixed;
}

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

  // Define your favorite OsiSolver
  
  OsiClpSolverInterface solver1;

  // Read in model using argv[1]
  // and assert that it is a clean model
  std::string dirsep(1,CoinFindDirSeparator());
  std::string mpsFileName;
# if defined(SAMPLEDIR)
  mpsFileName = SAMPLEDIR ;
  mpsFileName += dirsep+"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(),"");
  if( numMpsReadErrors != 0 )
  {
     printf("%d errors reading MPS file\n", numMpsReadErrors);
     return numMpsReadErrors;
  }
  double time1 = CoinCpuTime();

  /* Options are:
     preprocess to do preprocessing
     time in minutes
     if 2 parameters and numeric taken as time
  */
  bool preProcess=false;
  double minutes=-1.0;
  int nGoodParam=0;
  for (int iParam=2; iParam<argc;iParam++) {
    if (!strcmp(argv[iParam],"preprocess")) {
      preProcess=true;
      nGoodParam++;
    } else if (!strcmp(argv[iParam],"time")) {
      if (iParam+1<argc&&isdigit(argv[iParam+1][0])) {
        minutes=atof(argv[iParam+1]);
        if (minutes>=0.0) {
          nGoodParam+=2;
          iParam++; // skip time
        }
      }
    }
  }
  if (nGoodParam==0&&argc==3&&isdigit(argv[2][0])) {
    // If time is given then stop after that number of minutes
    minutes = atof(argv[2]);
    if (minutes>=0.0) 
      nGoodParam=1;
  }
  if (nGoodParam!=argc-2&&argc>=2) {
    printf("Usage <file> [preprocess] [time <minutes>] or <file> <minutes>\n");
    exit(1);
  }
  solver1.initialSolve();
  // Reduce printout
  solver1.setHintParam(OsiDoReducePrint,true,OsiHintTry);
  // See if we want preprocessing
  OsiSolverInterface * solver2=&solver1;
#if PREPROCESS==1
  CglPreProcess process;
  if (preProcess) {
    /* Do not try and produce equality cliques and
       do up to 5 passes */
    solver2 = process.preProcess(solver1,false,5);
    if (!solver2) {
      printf("Pre-processing says infeasible\n");
      exit(2);
    }
    solver2->resolve();
  }
#endif
  CbcModel model(*solver2);
  model.solver()->setHintParam(OsiDoReducePrint,true,OsiHintTry);
  // Set up some cut generators and defaults
  // Probing first as gets tight bounds on continuous

  CglProbing generator1;
  generator1.setUsingObjective(true);
  generator1.setMaxPass(1);
  generator1.setMaxPassRoot(5);
  // Number of unsatisfied variables to look at
  generator1.setMaxProbe(10);
  generator1.setMaxProbeRoot(1000);
  // How far to follow the consequences
  generator1.setMaxLook(50);
  generator1.setMaxLookRoot(500);
  // Only look at rows with fewer than this number of elements
  generator1.setMaxElements(200);
  generator1.setRowCuts(3);

  CglGomory generator2;
  // try larger limit
//.........这里部分代码省略.........

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

  /* Define your favorite OsiSolver.

     CbcModel clones the solver so use solver1 up to the time you pass it
     to CbcModel then use a pointer to cloned solver (model.solver())
  */
  
  OsiClpSolverInterface solver1;
  /* From now on we can build model in a solver independent way.
     You can add rows one at a time but for large problems this is slow so
     this example uses CoinBuild or CoinModel
  */
  OsiSolverInterface * solver = &solver1;
  // Data (is exmip1.mps in Mps/Sample
  // Objective 
  double objValue[]={1.0,2.0,0.0,0.0,0.0,0.0,0.0,-1.0};
  // Lower bounds for columns
  double columnLower[]={2.5,0.0,0.0,0.0,0.5,0.0,0.0,0.0};
  // Upper bounds for columns
  double columnUpper[]={COIN_DBL_MAX,4.1,1.0,1.0,4.0,
                  COIN_DBL_MAX,COIN_DBL_MAX,4.3};
  // Lower bounds for row activities
  double rowLower[]={2.5,-COIN_DBL_MAX,-COIN_DBL_MAX,1.8,3.0};
  // Upper bounds for row activities
  double rowUpper[]={COIN_DBL_MAX,2.1,4.0,5.0,15.0};
  // Matrix stored packed
  int column[] = {0,1,3,4,7,
                  1,2,
                  2,5,
                  3,6,
                  4,7};
  double element[] = {3.0,1.0,-2.0,-1.0,-1.0,
                      2.0,1.1,
                      1.0,1.0,
                      2.8,-1.2,
                      1.0,1.9};
  int starts[]={0,5,7,9,11,13};
  // Integer variables (note upper bound already 1.0)
  int whichInt[]={2,3};
  int numberRows=(int) (sizeof(rowLower)/sizeof(double));
  int numberColumns=(int) (sizeof(columnLower)/sizeof(double));
#define BUILD 2
#if BUILD==1
  // Using CoinBuild 
  // First do columns (objective and bounds)
  int i;
  // We are not adding elements 
  for (i=0;i<numberColumns;i++) {
    solver->addCol(0,NULL,NULL,columnLower[i],columnUpper[i],
                    objValue[i]);
  }
  // mark as integer
  for (i=0;i<(int) (sizeof(whichInt)/sizeof(int));i++)
    solver->setInteger(whichInt[i]);
  // Now build rows
  CoinBuild build;
  for (i=0;i<numberRows;i++) {
    int startRow = starts[i];
    int numberInRow = starts[i+1]-starts[i];
    build.addRow(numberInRow,column+startRow,element+startRow,
                 rowLower[i],rowUpper[i]);
  }  
  // add rows into solver
  solver->addRows(build);
#else
  /* using CoinModel - more flexible but still beta.
     Can do exactly same way but can mix and match much more.
     Also all operations are on building object
  */
  CoinModel build;
  // First do columns (objective and bounds)
  int i;
  for (i=0;i<numberColumns;i++) {
    build.setColumnBounds(i,columnLower[i],columnUpper[i]);
    build.setObjective(i,objValue[i]);
  }
  // mark as integer
  for (i=0;i<(int) (sizeof(whichInt)/sizeof(int));i++)
    build.setInteger(whichInt[i]);
  // Now build rows
  for (i=0;i<numberRows;i++) {
    int startRow = starts[i];
    int numberInRow = starts[i+1]-starts[i];
    build.addRow(numberInRow,column+startRow,element+startRow,
                 rowLower[i],rowUpper[i]);
  }  
  // add rows into solver
  solver->loadFromCoinModel(build);
#endif

  // Pass to solver
  CbcModel model(*solver);
  model.solver()->setHintParam(OsiDoReducePrint,true,OsiHintTry);


  // Set up some cut generators and defaults
  // Probing first as gets tight bounds on continuous

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

//.........这里部分代码省略.........
    //   int numRows = m->getNumRows();
    CoinFillN(integers_ ,n, true);
    for (int i = 0 ;  i < numCols ; i++)
    {
        if (si.isContinuous(i))
            integers_[i] = false;
    }
    bool * integerSlacks = integers_ + numCols;
    for (int i = 0 ; i < nBasics ; i++)
    {
        if (rowLower[i] > -1e50 && INT_INFEAS(rowLower[i]) > 1e-15)
            integerSlacks[i] = false;
        if (rowUpper[i] < 1e50 && INT_INFEAS(rowUpper[i]) > 1e-15)
            integerSlacks[i] = false;
    }
    for (int i = 0 ;  i < numCols ; i++)
    {
        CoinBigIndex end = starts[i] + lengths[i];
        if (integers_[i])
        {
            for (CoinBigIndex k=starts[i] ; k < end; k++)
            {
                if (integerSlacks[inds[k]] && INT_INFEAS(elems[k])>1e-15 )
                    integerSlacks[inds[k]] = false;
            }
        }
        else
        {
            for (CoinBigIndex k=starts[i] ; k < end; k++)
            {
                if (integerSlacks[inds[k]])
                    integerSlacks[inds[k]] = false;
            }
        }
    }

    CoinCopyN(si.getColSolution(), si.getNumCols(), colsol_);
    CoinCopyN(si.getRowActivity(), si.getNumRows(), slacks_);
    for (int i = 0 ; i < si.getNumRows() ; i++)
    {
        slacks_[i]*=-1;
        if (rowLower[i]>-1e50)
        {
            slacks_[i] += rowLower[i];
        }
        else
        {
            slacks_[i] += rowUpper[i];
        }
    }
    //Now get the fill the arrays;
    nNonBasics = 0;
    nBasics = 0;



    //For having the index variables correctly ordered we need to access to OsiSimplexInterface
    {
        OsiSolverInterface * ncSi = (const_cast<OsiSolverInterface *>(&si));
        ncSi->enableSimplexInterface(0);
        ncSi->getBasics(basics_);
	// Save enabled solver
	solver_ = si.clone();
#ifdef COIN_HAS_OSICLP
	OsiClpSolverInterface * clpSi = dynamic_cast<OsiClpSolverInterface *>(solver_);
	const OsiClpSolverInterface * clpSiRhs = dynamic_cast<const OsiClpSolverInterface *>(&si);
	if (clpSi)
	  clpSi->getModelPtr()->copyEnabledStuff(clpSiRhs->getModelPtr());;
#endif
        ncSi->disableSimplexInterface();
    }

    int numStructural = basis_->getNumStructural();
    for (int i = 0 ; i < numStructural ; i++)
    {
        if (basis_->getStructStatus(i)== CoinWarmStartBasis::basic)
        {
            nBasics++;
            //Basically do nothing
        }
        else
        {
            nonBasics_[nNonBasics++] = i;
        }
    }

    int numArtificial = basis_->getNumArtificial();
    for (int i = 0 ; i < numArtificial ; i++)
    {
        if (basis_->getArtifStatus(i)== CoinWarmStartBasis::basic)
        {
            //Just check number of basics
            nBasics++;
        }
        else
        {
            nonBasics_[nNonBasics++] = i + basis_->getNumStructural();
        }
    }
}

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

  // Define your favorite OsiSolver
  
  OsiClpSolverInterface solver1;

  // 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);
  double time1 = CoinCpuTime();

  /* Options are:
     preprocess to do preprocessing
     time in minutes
     if 2 parameters and numeric taken as time
  */
  bool preProcess=false;
  double minutes=-1.0;
  int nGoodParam=0;
  for (int iParam=2; iParam<argc;iParam++) {
    if (!strcmp(argv[iParam],"preprocess")) {
      preProcess=true;
      nGoodParam++;
    } else if (!strcmp(argv[iParam],"time")) {
      if (iParam+1<argc&&isdigit(argv[iParam+1][0])) {
        minutes=atof(argv[iParam+1]);
        if (minutes>=0.0) {
          nGoodParam+=2;
          iParam++; // skip time
        }
      }
    }
  }
  if (nGoodParam==0&&argc==3&&isdigit(argv[2][0])) {
    // If time is given then stop after that number of minutes
    minutes = atof(argv[2]);
    if (minutes>=0.0) 
      nGoodParam=1;
  }
  if (nGoodParam!=argc-2&&argc>=2) {
    printf("Usage <file> [preprocess] [time <minutes>] or <file> <minutes>\n");
    exit(1);
  }
  //solver1.getModelPtr()->setLogLevel(0);
  solver1.messageHandler()->setLogLevel(0);
  solver1.initialSolve();
  // Reduce printout
  solver1.setHintParam(OsiDoReducePrint,true,OsiHintTry);
  CbcModel model(solver1);
  model.solver()->setHintParam(OsiDoReducePrint,true,OsiHintTry);
  // Set up some cut generators and defaults
  // Probing first as gets tight bounds on continuous

  CglProbing generator1;
  generator1.setUsingObjective(true);
  generator1.setMaxPass(1);
  generator1.setMaxPassRoot(5);
  // Number of unsatisfied variables to look at
  generator1.setMaxProbe(10);
  generator1.setMaxProbeRoot(1000);
  // How far to follow the consequences
  generator1.setMaxLook(50);
  generator1.setMaxLookRoot(500);
  // Only look at rows with fewer than this number of elements
  generator1.setMaxElements(200);
  generator1.setRowCuts(3);

  CglGomory generator2;
  // try larger limit
  generator2.setLimit(300);

  CglKnapsackCover generator3;

  CglRedSplit generator4;
  // try larger limit
  generator4.setLimit(200);

  CglClique generator5;
  generator5.setStarCliqueReport(false);
  generator5.setRowCliqueReport(false);

  CglMixedIntegerRounding2 mixedGen;
  CglFlowCover flowGen;
  
  // Add in generators
  // Experiment with -1 and -99 etc
  model.addCutGenerator(&generator1,-1,"Probing");
  model.addCutGenerator(&generator2,-1,"Gomory");
//.........这里部分代码省略.........

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

  // Define your favorite OsiSolver
  
  OsiClpSolverInterface solver1;

  // 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);
  double time1 = CoinCpuTime();
  OsiClpSolverInterface solverSave = solver1;

  /* Options are:
     preprocess to do preprocessing
     time in minutes
     if 2 parameters and numeric taken as time
  */
  bool preProcess=false;
  double minutes=-1.0;
  int nGoodParam=0;
  for (int iParam=2; iParam<argc;iParam++) {
    if (!strcmp(argv[iParam],"preprocess")) {
      preProcess=true;
      nGoodParam++;
    } else if (!strcmp(argv[iParam],"time")) {
      if (iParam+1<argc&&isdigit(argv[iParam+1][0])) {
        minutes=atof(argv[iParam+1]);
        if (minutes>=0.0) {
          nGoodParam+=2;
          iParam++; // skip time
        }
      }
    }
  }
  if (nGoodParam==0&&argc==3&&isdigit(argv[2][0])) {
    // If time is given then stop after that number of minutes
    minutes = atof(argv[2]);
    if (minutes>=0.0) 
      nGoodParam=1;
  }
  if (nGoodParam!=argc-2&&argc>=2) {
    printf("Usage <file> [preprocess] [time <minutes>] or <file> <minutes>\n");
    exit(1);
  }
  // Reduce printout
  solver1.setHintParam(OsiDoReducePrint,true,OsiHintTry);
  // See if we want preprocessing
  OsiSolverInterface * solver2=&solver1;
  CglPreProcess process;
  // Never do preprocessing until dual tests out as can fix incorrectly
  preProcess=false;
  if (preProcess) {
    /* Do not try and produce equality cliques and
       do up to 5 passes */
    solver2 = process.preProcess(solver1,false,5);
    if (!solver2) {
      printf("Pre-processing says infeasible\n");
      exit(2);
    }
    solver2->resolve();
  }
  // Turn L rows into cuts
  CglStoredUser stored;
 {
  int numberRows = solver2->getNumRows();

  int * whichRow = new int[numberRows];
  // get row copy
  const CoinPackedMatrix * rowCopy = solver2->getMatrixByRow();
  const int * column = rowCopy->getIndices();
  const int * rowLength = rowCopy->getVectorLengths();
  const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
  const double * rowLower = solver2->getRowLower();
  const double * rowUpper = solver2->getRowUpper();
  const double * element = rowCopy->getElements();
  int iRow,nDelete=0;
  for (iRow=0;iRow<numberRows;iRow++) {
    if (rowLower[iRow]<-1.0e20||rowUpper[iRow]>1.0e20) {
      // take out
      whichRow[nDelete++]=iRow;
    }
  }
  // leave some rows to avoid empty problem (Gomory does not like)
  nDelete = CoinMax(CoinMin(nDelete,numberRows-5),0);
  for (int jRow=0;jRow<nDelete;jRow++) {
    iRow=whichRow[jRow];
    int start = rowStart[iRow];
    stored.addCut(rowLower[iRow],rowUpper[iRow],rowLength[iRow],
//.........这里部分代码省略.........

int main(int argc, char **argv) 
{
  char *f_name_lp, *last_dot_pos, f_name[256], *f_name_pos;
  int i, ncol;

  if((argc < 2) || (argc > 2)) {
    printf("### ERROR: main(): Usage: One of the following\ncgl_data_test input_file_name.mps\ncgl_data_test input_file_name.lp\n");
    exit(1);
  }

  f_name_lp = strdup(argv[1]);
  f_name_pos = strrchr(f_name_lp, '/');
  if(f_name_pos != NULL) {
    strcpy(f_name, &(f_name_pos[1]));
  }
  else {
    strcpy(f_name, f_name_lp);
  }
  last_dot_pos = strrchr(f_name, '.');
  if(last_dot_pos != NULL) {
    last_dot_pos = '\0';
  }

  OsiClpSolverInterface *clp = new OsiClpSolverInterface;
  clp->messageHandler()->setLogLevel(0);
  if(strcmp(&(f_name_lp[strlen(f_name_lp)-3]), ".lp") == 0) {
    clp->readLp(f_name_lp);    
  }
  else {
    if(strcmp(&(f_name_lp[strlen(f_name_lp)-4]), ".mps") == 0) {
      clp->readMps(f_name_lp);    
    }
    else {
      printf("### ERROR: unrecognized file type\n");
      exit(1);
    }
  }
  ncol = clp->getNumCols();
  clp->initialSolve();

  printf("LP value: %12.2f\n", clp->getObjValue());

  OsiCuts cuts;

  // Define parameters for CglRedSplit generator
  CglParam cpar;
  cpar.setMAX_SUPPORT(ncol+1);
  CglRedSplitParam rspar(cpar);

  // Create a cut generator with the given parameters
  CglRedSplit cutGen(rspar);

  char *colType = new char[ncol];
  for(i=0; i<ncol; i++) {
    if(clp->isContinuous(i)) {
      colType[i] = 'C';
    }
    else {
      colType[i] = 'I';
    }
  }

  int round, max_rounds = 10;
  for(round=0; round<max_rounds; round++) {
    cutGen.generateCuts(*clp, cuts);

    int ncuts = cuts.sizeRowCuts();

    const OsiRowCut **newRowCuts = new const OsiRowCut * [ncuts];
    for(i=0; i<ncuts; i++) {
      newRowCuts[i] = &cuts.rowCut(i); 
    }
    clp->applyRowCuts(ncuts, newRowCuts);
    delete[] newRowCuts;

    printf("round %4d: %4d generated cuts  new objective value: %12.2f\n", 
	   round, ncuts, clp->getObjValue());

    clp->resolve();  

    if(clp->isAbandoned()) {
      printf("###ERROR: Numerical difficulties in Solver\n");
      exit(1);
    }
  
    if(clp->isProvenPrimalInfeasible()) {
      printf("### WARNING: Problem is infeasible\n");
      exit(1);
    }
  }

  delete clp;
  free(f_name_lp);
  delete[] colType;

  return(0);
}

/*
  Randomized Rounding Heuristic
  Returns 1 if solution, 0 if not
*/
int
CbcHeuristicRandRound::solution(double & solutionValue,
                                double * betterSolution)
{
    // rlh: Todo: Memory Cleanup

    //  std::cout << "Entering the Randomized Rounding Heuristic" << std::endl;

    setWhen(1);  // setWhen(1) didn't have the effect I expected (e.g., run once).

    // Run only once.
    //
    //    See if at root node
    bool atRoot = model_->getNodeCount() == 0;
    int passNumber = model_->getCurrentPassNumber();
    //    Just do once
    if (!atRoot || passNumber > 1) {
        // std::cout << "Leaving the Randomized Rounding Heuristic" << std::endl;
        return 0;
    }

    std::cout << "Entering the Randomized Rounding Heuristic" << std::endl;
    typedef struct {
        int numberSolutions;
        int maximumSolutions;
        int numberColumns;
        double ** solution;
        int * numberUnsatisfied;
    } clpSolution;

    double start = CoinCpuTime();
    numCouldRun_++; //
#ifdef HEURISTIC_INFORM
    printf("Entering heuristic %s - nRuns %d numCould %d when %d\n",
	   heuristicName(),numRuns_,numCouldRun_,when_);
#endif
    // Todo: Ask JJHF what "number of times
    // the heuristic could run" means.

    OsiSolverInterface * solver = model_->solver()->clone();
    double primalTolerance ;
    solver->getDblParam(OsiPrimalTolerance, primalTolerance) ;
    OsiClpSolverInterface * clpSolver = dynamic_cast<OsiClpSolverInterface *> (solver);
    assert (clpSolver);
    ClpSimplex * simplex = clpSolver->getModelPtr();

    // Initialize the structure holding the solutions for the Simplex iterations
    clpSolution solutions;
    // Set typeStruct field of ClpTrustedData struct to 1 to indicate
    // desired behavior for  RandRound heuristic (which is what?)
    ClpTrustedData trustedSolutions;
    trustedSolutions.typeStruct = 1;
    trustedSolutions.data = &solutions;
    solutions.numberSolutions = 0;
    solutions.maximumSolutions = 0;
    solutions.numberColumns = simplex->numberColumns();
    solutions.solution = NULL;
    solutions.numberUnsatisfied = NULL;
    simplex->setTrustedUserPointer(&trustedSolutions);

    // Solve from all slack to get some points
    simplex->allSlackBasis();

    // Calling primal() invalidates pointers to some rim vectors,
    // like...row sense (!)
    simplex->primal();

    // 1. Okay - so a workaround would be to copy the data I want BEFORE
    // calling primal.
    // 2. Another approach is to ask the simplex solvers NOT to mess up my
    // rims.
    // 3. See freeCachedResults() for what is getting
    // deleted. Everything else points into the structure.
    // ...or use collower and colupper rather than rowsense.
    // ..store address of where one of these

    // Store the basic problem information
    // -Get the number of columns, rows and rhs vector
    int numCols = clpSolver->getNumCols();
    int numRows = clpSolver->getNumRows();

    // Find the integer variables (use columnType(?))
    // One if not continuous, that is binary or general integer)
    // columnType() = 0 continuous
    //              = 1 binary
    //              = 2 general integer
    bool * varClassInt = new bool[numCols];
    const char* columnType = clpSolver->columnType();
    int numGenInt = 0;
    for (int i = 0; i < numCols; i++) {
        if (clpSolver->isContinuous(i))
            varClassInt[i] = 0;
        else
            varClassInt[i] = 1;
        if (columnType[i] == 2) numGenInt++;
    }
//.........这里部分代码省略.........