CBC User Guide

The COIN ^[1] Branch and Cut solver (CBC) is an open-source mixed-integer program (MIP) solver written in C++. CBC is intended to be used primarily as a callable library to create customized branch-and-cut solvers. A basic, stand-alone executable version is also available. CBC is an active open-source project led by John Forrest at www.coin-or.org.

The primary users of CBC are expected to be developers implementing customized branch-and-cut algorithms in C++ using CBC as a library. Consequently, this document assumes a working knowledge of C++, including basic object-oriented programming terminology, and familiarity with the fundamental concepts of linear programming (LP) and mixed integer programming (MIP).

CBC relies on other parts of the COIN repository. CBC needs a LP solver and relies on the COIN Open Solver Inteface (OSI) to communicate with the user's choice of solver. Any LP solver with an OSI interface can be used with CBC. The LP solver expected to be used most commonly is COIN's native linear program solver, CLP. For cut generators, CBC relies on the COIN Cut Generation Library (CGL). Any cut generator written to CGL standards can be used with CBC. Some of the cut generators in CGL rely on other parts of COIN, e.g., CGL's Gomory cut generator rely on the factorization functionality of CoinFactorization. This document assumes basic familiarity with OSI and CGL.

Technically speaking, CBC accesses the solver (and sometime the model and data it contains) through an OSISolverInterface. For the sake of simplicity, we will refer to the OsiSolverInterface as "the solver" in this document, rather than "the standard application programming interface to the solver." We hope any confusion caused by blurring this distinction will be mitigated by the shorter sentences.

In summary, readers should have the following prerequisites:

Before examining CBC in more detail, we tersely describe the basic branch-and-cut algorithm by way of example, (which should really be called branch-and-cut-and-bound) and show the major C++ class(es) in CBC related to each step. The major CBC classes, labeled (A) through (F), are described in Table 1.1.

Step 1. (Bound) Given a MIP model to minimize where some variables must take on integer values (e.g., 0, 1, or 2), relax the integrality requirements (e.g., consider each "integer" variable to be continuous with a lower bound of 0.0 and an upper bound of 2.0). Solve the resulting linear model with an LP solver to obtain a lower bound on the MIP's objective function value. If the optimal LP solution has integer values for the MIP's integer variables, we are finished. Any MIP-feasible solution provides an upper bound on the objective value. The upper bound equals the lower bound; the solution is optimal.

Step 2. (Branch) Otherwise, there exists an "integer" variable with a non-integral value. Choose one non-integral variable (e.g., with value 1.3) (A)(B) and branch. Create two ^[2] nodes, one with the branching variable having an upper bound of 1.0, and the other with the branching variable having a lower bound of 2.0. Add the two nodes to the search tree.

While (search tree is not empty) {

Step 3. (Choose Node) Pick a node off the tree (C)(D)

Step 4. (Re-optimize LP) Create an LP relaxation and solve.

Step 5. (Bound) Interrogate the optimal LP solution, and try to prune the node by one of the following.

Step 6. (Branch) If we were unable to prune the node, then branch. Choose one non-integral variable to branch on (A)(B). Create two nodes and add them to the search tree. }

This is the outline of a "branch-and-bound" algorithm. If in optimizing the linear programs, we use cuts to tighten the LP relaxations (E)(F), then we have a "branch-and-cut" algorithm. (Note, if cuts are only used in Step 1, the method is called a "cut-and-branch" algorithm.)

There are a number of resources available to help new CBC users get started. This document is designed to be used in conjunction with the files in the Samples subdirectory of the main CBC directory (COIN/Cbc/Samples). The Samples illustrate how to use CBC and may also serve as useful starting points for user projects. In the event that either this document or the available Doxygen content conflicts with the observed behavior of the source code, the comments in the header files, found in COIN/Cbc/include, are the ultimate reference.

The main class in CBC is CbcModel. The CbcModel class is where most of the parameter setting is done. The absolute minimum number of actions taken with CbcModel is two,

The first sample program shows how to perform simple branch-and-bound with CBC. This program is short enough to present in full. Most of the remaining examples will take the form of small code fragments. The complete code for all the examples in this Guide can be found in the CBC Samples directory, COIN/Cbc/Samples.

    
// Copyright (C) 2005, International Business Machines
// Corporation and others.  All Rights Reserved.

#include "CbcModel.hpp"

// Using CLP as the solver
#include "OsiClpSolverInterface.hpp"

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

  // Read in example model in MPS file format
  // and assert that it is a clean model
  int numMpsReadErrors = solver1.readMps("../../Mps/Sample/p0033.mps","");
  assert(numMpsReadErrors==0);

  // Pass the solver with the problem to be solved to CbcModel 
  CbcModel model(solver1);

  // Do complete search
  model.branchAndBound();

  /* Print the solution.  CbcModel clones the solver so we
     need to get current copy from the CbcModel */
  int numberColumns = model.solver()->getNumCols();
    
  const double * solution = model.bestSolution();
    
  for (int iColumn=0;iColumn<numberColumns;iColumn++) {
    double value=solution[iColumn];
    if (fabs(value)>1.0e-7&&model.solver()->isInteger(iColumn)) 
      printf("%d has value %g\n",iColumn,value);
   }
  return 0;
}    
     
  

The program in Example 2.1 creates a OsiClpSolverInterface solver interface (i.e., solver1), and reads an MPS file. If there are no errors, the program passes the problem to CbcModel which solves the problem using the branch-and-bound algorithm. The part of the program which solves the problem is very small (one line!) but before that one line, the LP solver (i.e., solver1) had to be created and populated with the problem. After that one line, the results were printed out.

The program in Example 2.1 illustrates the dependency of CBC on the OsiSolverInterface class. The constructor of CbcModel takes a pointer to an OsiSolverInterface (i.e., a solver). The CbcModel clones the solver, and uses its own instance of the solver. The CbcModel's solver and the original solver (e.g., solver1) are not in sync unless the user synchronizes them. The user can always access the CbcModel's solver through the model() class. To synchronize the two solvers, explicitly refreshing the original, e.g.,

  solver1 = model.solver();

CbcModel's method solver() returns a pointer to CBC's cloned solver.

For convenience, many of the OSI methods to access problem data have identical method names in CbcModel. (It's just more convenient to type model.getNumCols() rather than model.solver()->getNumCols()). The CbcModel refreshes its solver at certain logical points during the algorithm. At these points, the information from the CbcModel model will match the information from the model.solver(). Elsewhere, the information may vary. For instance, the method CbcModel::bestSolution() will contain the best solution so far, the OSI method getColSolution() may not. In this case, it is safer to use CbcModel::bestSolution().

While all the OSI methods used in minimum.cpp have equivalent methods in CbcModel, there are some OSI methods which do not. For example, if the program produced a lot of undesired output, one might add the line

  model.solver()->setHintParam(OsiDoReducePrint,true,OsiHintTry); 

to reduce the output. There is no setHintParam() method in CbcModel.

Optimality can be checked through a call to model.isProvenOptimal(). Also available are isProvenInfeasible(), isSolutionLimitReached(), isNodeLimitReached() or the feared isAbandoned(). There is also int status() which returns 0 if finished (which includes the case when the algorithm is finished because it has been proved infeasible), 1 if stopped by user, and 2 if difficulties arose.

In addition to these CbcModel methods, solution values can be accessed via OSI methods. The OSI methods pick up the current solution in the CBCModel. The current solution will match the best solution found so far if called after branchAndBound() and a solution was found.

Most of the parameter setting in CBC is done through CbcModel methods. The most commonly used set and get methods are listed in Table 2.2.

CbcModel is extremely flexible and customizable. The class structure of CBC is designed to make the most commonly desired customizations of branch-and-cut possible. These include:

To enable this flexibility, CbcModel uses other classes in CBC (some of which are virtual and may have multiple instances). Not all classes are created equal. The two tables below list in alphabetical order the classes used by CbcModel that are of most interest and of least interest.

There is not much about the classes listed in Table 2.4 that the average user needs to know about.

The order in which the nodes of the search tree are explored can strongly influence the performance of branch-and-cut algorithms. CBC give users complete control over the search order, including the ability to dynamically change the node selection logic as the search progresses. The search order is controlled via the CbcCompare... class, and its method test(). Dynamic changes can be made whenever

CBC provides an abstract base class, CbcCompareBase, and implementations of several commonly used node selection strategies as Compare Classes, see Table 3.1.

It is relatively simple for a user to create a customized node selection by creating a new compare class instances. The code in Example 3.1 describes how to build a new comparison class and the reasoning behind it. The complete source can be found in CbcCompareUser.hpp and CbcCompareUser.cpp, located in the CBC Samples directory. Besides the constructor, the only method the user -must- implement in CbcCompare is bool test(CbcNode* x, CbcNode* y)) which returns true if node y is preferred over node x. In the test() method, information from CbcNode can easily be used. Table 3.2 lists some commonly used methods to access information at a node.

The node desired in the tree is often a function of the how the search is progressing. In the design of CBC, there is no information on the state of the tree. The CBC is designed so that the method newSolution() is called whenever a solution is found and the method every1000Nodes() is called every 1000 nodes. When these methods are called, the user has the opportunity to modify the behavior of test() by adjusting their common variables (e.g., weight_). Because CbcNode has a pointer to the model, the user can also influence the search through actions such as changing the maximum time CBC is allowed, once a solution has been found (e.g., CbcModel::setMaximumSeconds(double value)). In CbcCompareUser.cpp of the COIN/Cbc/Samples directory, four items of data are used.

Initially, weight_ is -1.0 and the search is biased towards depth first. In fact, test() prefers y if y has fewer unsatisfied variables. In the case of a tie, test() prefers the node with the greater depth in tree. The full code for the CbcCompareUser::test() method is given in Example 3.1.

    
// Returns true if y better than x
bool 
CbcCompareUser::test (CbcNode * x, CbcNode * y)
{
  if (weight_==-1.0) {
    // before solution
    if (x->numberUnsatisfied() > y->numberUnsatisfied())
      return true;
    else if (x->numberUnsatisfied() < y->numberUnsatisfied())
      return false;
    else
      return x->depth() < y->depth();
  } else {
    // after solution. 
    // note: if weight_=0, comparison is based 
    //       solely on objective value
    double weight = CoinMax(weight_,0.0);
    return x->objectiveValue()+ weight*x->numberUnsatisfied() > 
      y->objectiveValue() + weight*y->numberUnsatisfied();
  }
}
     
  

CBC calls the method newSolution() after a new solution is found. The method newSolution() interacts with test() by means of the variable weight_. If the solution was achieved by branching, a calculation is made to determine the cost per unsatisfied integer variable to go from the continuous solution to an integer solution. The variable weight_ is then set to aim at a slightly better solution. From then on, test() returns true if it seems that y will lead to a better solution than x. This source for newSolution() in given in Example 3.2.

    
// This allows the test() method to change behavior by resetting weight_.
// It is called after each new solution is found.
void 
CbcCompareUser::newSolution(CbcModel * model,
			       double objectiveAtContinuous,
			       int numberInfeasibilitiesAtContinuous) 
{
  if (model->getSolutionCount()==model->getNumberHeuristicSolutions())
    return; // The number of solutions found by any means equals the 
            // number of solutions, so this solution was found by rounding. 
            // Ignore it.

  // set weight_ to get close to this solution
  double costPerInteger = 
    (model->getObjValue()-objectiveAtContinuous)/
    ((double) numberInfeasibilitiesAtContinuous);
  weight_ = 0.98*costPerInteger; // this aims for a solution 
                                 // slightly better than known.
                                 // why 0.98? why not?! Experiment yourself.                                 
  saveWeight_=weight_; // We're going to switching between depth-first and breadth-first
                       // branching strategies, depending on what we find in the tree.
                       // When doing depth first, we'll want to retrieve this weight. 
                       // So, let's save it. 
  numberSolutions_++;
  if (numberSolutions_>5)
    weight_ =0.0; // comparison in test() will be 
                  // based strictly on objective value.
}
     
  

As the search progresses, the comparison can be modified. If many nodes (or many solutions) have been generated, then weight_ is set to 0.0 leading to a breadth-first search. Breadth-first search can lead to an enormous tree. If the tree size is exceeds 10000, it may be desirable to return to a search biased towards depth first. Changing the behavior in this manner is done by the method every1000Nodes shown in Example 3.3.

    
// This allows the test() method to change behavior every 1000 nodes.
bool 
CbcCompareUser::every1000Nodes(CbcModel * model, int numberNodes)
{
  if (numberNodes>10000) 
    weight_ =0.0; // compare nodes based on objective value
    // get size of tree
  treeSize_ = model->tree()->size();
  if (treeSize_>10000) {
    // set weight to reduce size most of time
    if (treeSize_>20000)
      weight_=-1.0;
    else if ((numberNodes%4000)!=0)  // Flip-flop between the strategies.
                                     // Why 4000? Why not? Experiment yourself.
      weight_=-1.0;
    else
      weight_=saveWeight_;
  }
  return numberNodes==11000; // resort if first time
}
     
  

In practice, it is very useful to get a good solution reasonably fast. Any MIP-feasible solution produces an upper bound, and a good bound will greatly reduce the run time. Good solutions can satisfy the user on very large problems where a complete search is impossible. Obviously, heuristics are problem dependent, although some do have more general use. At present there is only one heuristic in CBC itself, CbcRounding. Hopefully, the number will grow. Other heuristics are in the COIN/Cbc/Samples directory. A heuristic tries to obtain a solution to the original problem so it only needs to consider the original rows and does not have to use the current bounds. CBC provides an abstract base class CbcHeuristic and a rounding heuristic in CBC.

This chapter describes how to build a greedy heuristic for a set covering problem, e.g., the miplib problem fast0507. A more general (and efficient) version of the heuristic is in CbcHeuristicGreedy.hpp and CbcHeuristicGreedy.cpp located in the COIN/Cbc/Samples directory, see Chapter 8, More Samples .

The greedy heuristic will leave all variables taking value one at this node of the tree at value one, and will initially set all other variables to value zero. All variables are then sorted in order of their cost divided by the number of entries in rows which are not yet covered. (We may randomize that value a bit so that ties will be broken in different ways on different runs of the heuristic.) The best one is choosen, and set to one. The process is repeated. Because this is a set covering problem (i.e., all constraints are ≥), the heuristic is guaranteed to find a solution (but not necessarily an improved solution). The speed of the heuristic could be improved by just redoing those affected, but for illustrative purposes we will keep it simple. (The speed could also be improved if all elements are 1.0).

The key CbcHeuristic method is int solution(double & solutionValue, double * betterSolution). The solution() method returns 0 if no solution found, and returns 1 if a solution is found, in which case it fills in the objective value and primal solution. The code in CbcHeuristicGreedy.cpp is a little more complicated than this following example. For instance, the code here assumes all variables are integer. The important bit of data is a copy of the matrix (stored by column) before any cuts have been made. The data used are bounds, objective and the matrix plus two work arrays.

    
  OsiSolverInterface * solver = model_->solver(); // Get solver from CbcModel
  const double * columnLower = solver->getColLower(); // Column Bounds
  const double * columnUpper = solver->getColUpper();
  const double * rowLower = solver->getRowLower(); // We know we only need lower bounds
  const double * solution = solver->getColSolution();
  const double * objective = solver->getObjCoefficients(); // In code we also use min/max
  double integerTolerance = model_->getDblParam(CbcModel::CbcIntegerTolerance);
  double primalTolerance;
  solver->getDblParam(OsiPrimalTolerance,primalTolerance);
  int numberRows = originalNumberRows_; // This is number of rows when matrix was passed in
  // Column copy of matrix (before cuts)
  const double * element = matrix_.getElements();
  const int * row = matrix_.getIndices();
  const CoinBigIndex * columnStart = matrix_.getVectorStarts();
  const int * columnLength = matrix_.getVectorLengths();

  // Get solution array for heuristic solution
  int numberColumns = solver->getNumCols();
  double * newSolution = new double [numberColumns];
  // And to sum row activities
  double * rowActivity = new double[numberRows];
     
  

The newSolution is then initialized to the rounded down solution.

    
  for (iColumn=0;iColumn<numberColumns;iColumn++) {
    CoinBigIndex j;
    double value = solution[iColumn];
    // Round down integer
    if (fabs(floor(value+0.5)-value)<integerTolerance) 
      value=floor(CoinMax(value+1.0e-3,columnLower[iColumn]));
    // make sure clean
    value = CoinMin(value,columnUpper[iColumn]);
    value = CoinMax(value,columnLower[iColumn]);
    newSolution[iColumn]=value;
    if (value) {
      double cost = objective[iColumn];
      newSolutionValue += value*cost;
      for (j=columnStart[iColumn];
           j<columnStart[iColumn]+columnLength[iColumn];j++) {
        int iRow=row[j];
        rowActivity[iRow] += value*element[j];
      }
    }
  }
     
  

At this point some row activities are below their lower bound. To correct the infeasibility, the variable which is cheapest in reducing the sum of infeasibilities is found and updated, and the process repeats. This is a finite process. (The implementation could be faster, but is kept simple for illustrative purposes.)

    
  while (true) {
    // Get column with best ratio
    int bestColumn=-1;
    double bestRatio=COIN_DBL_MAX;
    for (int iColumn=0;iColumn<numberColumns;iColumn++) {
      CoinBigIndex j;
      double value = newSolution[iColumn];
      double cost = direction * objective[iColumn];
      // we could use original upper rather than current
      if (value+0.99<columnUpper[iColumn]) {
        double sum=0.0; // Compute how much we will reduce infeasibility by
        for (j=columnStart[iColumn];
             j<columnStart[iColumn]+columnLength[iColumn];j++) {
          int iRow=row[j];
          double gap = rowLower[iRow]-rowActivity[iRow];
          if (gap>1.0e-7) {
            sum += CoinMin(element[j],gap);
          if (element[j]+rowActivity[iRow]<rowLower[iRow]+1.0e-7) {
            sum += element[j];
	  }
	}
        if (sum>0.0) {
	  double ratio = (cost/sum)*(1.0+0.1*CoinDrand48());
	  if (ratio<bestRatio) {
            bestRatio=ratio;
            bestColumn=iColumn;
          }
	}
      }
    }
    if (bestColumn<0)
      break; // we have finished
    // Increase chosen column
    newSolution[bestColumn] += 1.0;
    double cost = direction * objective[bestColumn];
    newSolutionValue += cost;
    for (CoinBigIndex j=columnStart[bestColumn];
         j<columnStart[bestColumn]+columnLength[bestColumn];j++) {
      int iRow = row[j];
      rowActivity[iRow] += element[j];
    }
  }
     
  

A solution value of newSolution is compared to the best solution value. If newSolution is an improvement, its feasibility is validated. We expect newSolution to be feasible, and are trapping for unexpected numerical errors.

    
  returnCode=0; // 0 means no good solution
  if (newSolutionValue<solutionValue) { // minimization
    // check feasible
    memset(rowActivity,0,numberRows*sizeof(double));
    for (iColumn=0;iColumn<numberColumns;iColumn++) {
      CoinBigIndex j;
      double value = newSolution[iColumn];
      if (value) {
	for (j=columnStart[iColumn];
	     j<columnStart[iColumn]+columnLength[iColumn];j++) {
	  int iRow=row[j];
	  rowActivity[iRow] += value*element[j];
	}
      }
    }
    // check was approximately feasible
    bool feasible=true;
    for (iRow=0;iRow<numberRows;iRow++) {
      if(rowActivity[iRow]<rowLower[iRow]) {
	if (rowActivity[iRow]<rowLower[iRow]-10.0*primalTolerance)
	  feasible = false;
      }
    }
    if (feasible) {
      // new solution
      memcpy(betterSolution,newSolution,numberColumns*sizeof(double));
      solutionValue = newSolutionValue;
      // We have good solution
      returnCode=1;
    }
  }
     
  

CBC's concept of branching is based on the idea of an "object". An object has (i) a feasible region, (ii) can be evaluated for infeasibility, (iii) can be branched on, e.g., a method of generating a branching object, which defines an up branch and a down branch, and (iv) allows comparsion of the effect of branching. Instances of objects include.

In Chapter 5, Branching , we give examples of how to use existing branching objects. (The next revision of this Guide should include an example of how to write your own branching object; Contributions of examples are welcome.)

If the user declares variables as integer but does no more, then CBC will treat them as simple integer variables. In many cases the user would like to do some more fine tuning. This section shows how to create integer variables with pseudo costs. When pseudo costs are given then it is assumed that if a variable is at 1.3 then the cost of branching that variable down will be 0.3 times the down pseudo cost and the cost of branching up would be 0.7 times the up pseudo cost. Pseudo costs can be used both for branching and for choosing a node. The full code is in longthin.cpp located in the CBC Samples directory, see Chapter 8, More Samples .

The idea is simple for set covering problems. Branching up gets us much closer to an integer solution so we will encourage that direction by branch up if variable value > 0.333333. The expected cost of going up obviously depends on the cost of the variable. The pseudo costs are choosen to reflect that fact.

    
  int iColumn;
  int numberColumns = solver3->getNumCols();
  // do pseudo costs
  CbcObject ** objects = new CbcObject * [numberColumns];
  // Point to objective
  const double * objective = model.getObjCoefficients();
  int numberIntegers=0;
  for (iColumn=0;iColumn<numberColumns;iColumn++) {
    if (solver3->isInteger(iColumn)) {
      double cost = objective[iColumn];
      CbcSimpleIntegerPseudoCost * newObject =
        new CbcSimpleIntegerPseudoCost(&model,numberIntegers,iColumn,
                                       2.0*cost,cost);
      newObject->setMethod(3);
      objects[numberIntegers++]= newObject;
    }
  }
  // Now add in objects (they will replace simple integers)
  model.addObjects(numberIntegers,objects);
  for (iColumn=0;iColumn<numberIntegers;iColumn++)
    delete objects[iColumn];
  delete [] objects;
     
  

The code in Example 5.1 also tries to give more importance to variables with more coefficients. Whether this sort of thing is worthwhile should be the subject of experimentation.

In crew scheduling, the problems are long and thin. A problem may have a few rows but many thousands of variables. Branching a variable to 1 is very powerful as it fixes many other variables to zero, but branching to zero is very weak as thousands of variables can increase from zero. In crew scheduling problems, each constraint is a flight leg, e.g., JFK airport to DFW airport. From DFW there may be several flights the crew could take next - suppose one flight is the 9:30 flight from DFW to LAX airport. A binary branch is that the crew arriving at DFW either take the 9:30 flight to LAX or they don't. This "follow-on" branching does not fix individual variables. Instead this branching divides all the variables with entries in the JFK-DFW constraint into two groups - those with entries in the DFW-LAX constraint and those without entries.

The full sample code for follow-on brancing is in crew.cpp located in the CBC Samples directory, see Chapter 8, More Samples ). In this case, the simple integer variables are left which may be necessary if other sorts of constraints exist. Follow-on branching rules are to be considered first, so the priorities are set to indicate the follow-on rules take precedence. Priority 1 is the highest priority.

    
  int iColumn;
  int numberColumns = solver3->getNumCols();
  /* We are going to add a single follow-on object but we
     want to give low priority to existing integers 
     As the default priority is 1000 we don't actually need to give
     integer priorities but it is here to show how.
  */
  // Normal integer priorities
  int * priority = new int [numberColumns];
  int numberIntegers=0;
  for (iColumn=0;iColumn<numberColumns;iColumn++) {
    if (solver3->isInteger(iColumn)) {
      priority[numberIntegers++]= 100; // low priority
    }
  }
  /* Second parameter is set to true for objects,
     and false for integers. This indicates integers */
  model.passInPriorities(priority,false);
  delete [] priority;
  /* Add in objects before we can give them a priority.
     In this case just one object
     - but it shows the general method
  */
  CbcObject ** objects = new CbcObject * [1];
  objects[0]=new CbcFollowOn(&model);
  model.addObjects(1,objects);
  delete objects[0];
  delete [] objects;
  // High priority
  int followPriority=1;
  model.passInPriorities(&followPriority,true);
     
  

In the next version of this Guide, we need to have an example illustrating how to use COIN's CGL with CBC. Contribtions are welcome.

CBC uses a generic OsiSolverInterface and its resolve capability. This does not give much flexibility so advanced users can inherit from their interface of choice. This section illustrates how to implement a specialized solver for a long thin problem, e.g., fast0507 again. As with the other examples in the Guide, the sample code is not guaranteed to be the fastest way to solve the problem. The main purpose of the example is to illustrate techniques. The full source is in CbcSolver2.hpp and CbcSolver2.cpp located in the CBC Samples directory, see Chapter 8, More Samples .

The method initialSolve is called a few times in CBC, and provides a convenient starting point. The modelPtr_ derives from OsiClpSolverInterface.

    
  // modelPtr_ is of type ClpSimplex *
  modelPtr_->setLogLevel(1); // switch on a bit of printout
  modelPtr_->scaling(0); // We don't want scaling for fast0507
  setBasis(basis_,modelPtr_); // Put basis into ClpSimplex
  // Do long thin by sprint
  ClpSolve options;
  options.setSolveType(ClpSolve::usePrimalorSprint);
  options.setPresolveType(ClpSolve::presolveOff);
  options.setSpecialOption(1,3,15); // Do 15 sprint iterations
  modelPtr_->initialSolve(options); // solve problem
  basis_ = getBasis(modelPtr_); // save basis
  modelPtr_->setLogLevel(0); // switch off printout
     
  

The resolve() method is more complicated than initialSolve(). The main pieces of data are a counter count_ which is incremented each solve and an integer array node_ which stores the last time a variable was active in a solution. For the first few solves, the normal Dual Simplex is called and node_ array is updated.

    
  if (count_<10) {
    OsiClpSolverInterface::resolve(); // Normal resolve
    if (modelPtr_->status()==0) {
      count_++; // feasible - save any nonzero or basic
      const double * solution = modelPtr_->primalColumnSolution();
      for (int i=0;i<numberColumns;i++) {
	if (solution[i]>1.0e-6||modelPtr_->getStatus(i)==ClpSimplex::basic) {
	  node_[i]=CoinMax(count_,node_[i]);
	  howMany_[i]++;
	}
      }
    } else {
      printf("infeasible early on\n");
    }
  }
     
  

After the first few solves, only those variables which took part in a solution in the last so many solves are used. As fast0507 is a set covering problem, any rows which are already covered can be taken out.

    
    int * whichRow = new int[numberRows]; // Array to say which rows used
    int * whichColumn = new int [numberColumns]; // Array to say which columns used
    int i;
    const double * lower = modelPtr_->columnLower();
    const double * upper = modelPtr_->columnUpper();
    setBasis(basis_,modelPtr_); // Set basis
    int nNewCol=0; // Number of columns in small model
    // Column copy of matrix
    const double * element = modelPtr_->matrix()->getElements();
    const int * row = modelPtr_->matrix()->getIndices();
    const CoinBigIndex * columnStart = modelPtr_->matrix()->getVectorStarts();
    const int * columnLength = modelPtr_->matrix()->getVectorLengths();
    
    int * rowActivity = new int[numberRows]; // Number of columns with entries in each row
    memset(rowActivity,0,numberRows*sizeof(int));
    int * rowActivity2 = new int[numberRows]; // Lower bound on row activity for each row
    memset(rowActivity2,0,numberRows*sizeof(int));
    char * mark = (char *) modelPtr_->dualColumnSolution(); // Get some space to mark columns
    memset(mark,0,numberColumns);
    for (i=0;i<numberColumns;i++) {
      bool choose = (node_[i]>count_-memory_&&node_[i]>0); // Choose if used recently
      // Take if used recently or active in some sense
      if ((choose&&upper[i])
	  ||(modelPtr_->getStatus(i)!=ClpSimplex::atLowerBound&&
             modelPtr_->getStatus(i)!=ClpSimplex::isFixed)
	  ||lower[i]>0.0) {
        mark[i]=1; // mark as used
	whichColumn[nNewCol++]=i; // add to list
        CoinBigIndex j;
        double value = upper[i];
        if (value) {
          for (j=columnStart[i];
               j<columnStart[i]+columnLength[i];j++) {
            int iRow=row[j];
            assert (element[j]==1.0);
            rowActivity[iRow] ++; // This variable can cover this row
          }
          if (lower[i]>0.0) {
            for (j=columnStart[i];
                 j<columnStart[i]+columnLength[i];j++) {
              int iRow=row[j];
              rowActivity2[iRow] ++; // This row redundant
            }
          }
        }
      }
    }
    int nOK=0; // Use to count rows which can be covered
    int nNewRow=0; // Use to make list of rows needed
    for (i=0;i<numberRows;i++) {
      if (rowActivity[i])
        nOK++;
      if (!rowActivity2[i])
        whichRow[nNewRow++]=i; // not satisfied
      else
        modelPtr_->setRowStatus(i,ClpSimplex::basic); // make slack basic
    }
    if (nOK<numberRows) {
      // The variables we have do not cover rows - see if we can find any that do
      for (i=0;i<numberColumns;i++) {
        if (!mark[i]&&upper[i]) {
          CoinBigIndex j;
          int good=0;
          for (j=columnStart[i];
               j<columnStart[i]+columnLength[i];j++) {
            int iRow=row[j];
            if (!rowActivity[iRow]) {
              rowActivity[iRow] ++;
              good++;
            }
          }
          if (good) {
            nOK+=good; // This covers - put in list
            whichColumn[nNewCol++]=i;
          }
        }
      }
    }
    delete [] rowActivity;
    delete [] rowActivity2;
    if (nOK<numberRows) {
      // By inspection the problem is infeasible - no need to solve
      modelPtr_->setProblemStatus(1);
      delete [] whichRow;
      delete [] whichColumn;
      printf("infeasible by inspection\n");
      return;
    }
    // Now make up a small model with the right rows and columns
    ClpSimplex *  temp = new ClpSimplex(modelPtr_,nNewRow,whichRow,nNewCol,whichColumn);
     
  

If the variables cover the rows, then the problem is feasible (no cuts are being used). (If the rows were equality constraints, then this might not be the case. More work would be needed.) After the solution to the subproblem, the reduced costs of the full problem are checked. If the reduced cost of any variable not in the subproblem is negative, the code goes back to the full problem and cleans up with Primal Simplex.

    
    temp->setDualObjectiveLimit(1.0e50); // Switch off dual cutoff as problem is restricted
    temp->dual(); // solve
    double * solution = modelPtr_->primalColumnSolution(); // put back solution
    const double * solution2 = temp->primalColumnSolution();
    memset(solution,0,numberColumns*sizeof(double));
    for (i=0;i<nNewCol;i++) {
      int iColumn = whichColumn[i];
      solution[iColumn]=solution2[i];
      modelPtr_->setStatus(iColumn,temp->getStatus(i));
    }
    double * rowSolution = modelPtr_->primalRowSolution();
    const double * rowSolution2 = temp->primalRowSolution();
    double * dual = modelPtr_->dualRowSolution();
    const double * dual2 = temp->dualRowSolution();
    memset(dual,0,numberRows*sizeof(double));
    for (i=0;i<nNewRow;i++) {
      int iRow=whichRow[i];
      modelPtr_->setRowStatus(iRow,temp->getRowStatus(i));
      rowSolution[iRow]=rowSolution2[i];
      dual[iRow]=dual2[i];
    }
    // See if optimal
    double * dj = modelPtr_->dualColumnSolution();
    // get reduced cost for large problem
    // this assumes minimization
    memcpy(dj,modelPtr_->objective(),numberColumns*sizeof(double));
    modelPtr_->transposeTimes(-1.0,dual,dj);
    modelPtr_->setObjectiveValue(temp->objectiveValue());
    modelPtr_->setProblemStatus(0);
    int nBad=0;
      
    for (i=0;i<numberColumns;i++) {
      if (modelPtr_->getStatus(i)==ClpSimplex::atLowerBound
          &&upper[i]>lower[i]&&dj[i]<-1.0e-5)
        nBad++;
    }
    // If necessary clean up with primal (and save some statistics)
    if (nBad) {
      timesBad_++;
      modelPtr_->primal(1);
      iterationsBad_ += modelPtr_->numberIterations();
    }
     
  

The array node_ is updated, as for the first few solves. To give some idea of the effect of this tactic, the problem fast0507 has 63,009 variables but the small problem never has more than 4,000 variables. In only about ten percent of solves was it necessary to resolve, and then the average number of iterations on full problem was less than 20.

To give another example - again only for illustrative purposes -- it is possible to do quadratic MIP with CBC. In this case, we make resolve the same as initialSolve. The full code is in ClpQuadInterface.hpp and ClpQuadInterface.cpp located in the CBC Samples directory, see Chapter 8, More Samples ).

    
  // save cutoff
  double cutoff = modelPtr_->dualObjectiveLimit();
  modelPtr_->setDualObjectiveLimit(1.0e50);
  modelPtr_->scaling(0);
  modelPtr_->setLogLevel(0);
  // solve with no objective to get feasible solution
  setBasis(basis_,modelPtr_);
  modelPtr_->dual();
  basis_ = getBasis(modelPtr_);
  modelPtr_->setDualObjectiveLimit(cutoff);
  if (modelPtr_->problemStatus()) 
    return; // problem was infeasible 
  // Now pass in quadratic objective
  ClpObjective * saveObjective  = modelPtr_->objectiveAsObject();
  modelPtr_->setObjectivePointer(quadraticObjective_);
  modelPtr_->primal(); // Th model has a quadratic objective, 
                       // so this invokes quadratic primal.
  modelPtr_->setDualObjectiveLimit(cutoff);
  if (modelPtr_->objectiveValue()>cutoff)
    modelPtr_->setProblemStatus(1);
  modelPtr_->setObjectivePointer(saveObjective);
     
  

Rather than implementing all the method from scratch, we based the quadratic solver ClpQuadInteface on the linear programming solver OsiClpSolverInterface. This is a convenient approach to take when prototyping ideas. After the merit of an idea is proven, the user can decide is a more serious implementation is warranted.

The CBC distribution includes a number of .cpp sample files. Users are encouraged to use them as starting points for their own CBC projects. The files can be found in the COIN/Cbc/Samples/ directory. For the latest information on compiling and running these samples, please see the file COIN/Cbc/Samples/INSTALL. Most of them can be built by

make DRIVER=name

which produces an executable testit. Below is a list of some of the most useful sample files with a short description for each file.