Advanced use of solver

Coin Branch and Cut uses a generic OsiSolverInterface and its resolve capability. This does not give much flexibility so advanced users can inherit from the interface of choice. This describes such a solver for a long thin problem e.g. fast0507 again. As with all these examples it is not guaranteed that this is the fastest way to solve any of these problems - they are to illustrate techniques. The full code is in CbcSolver2.hpp and CbcSolver2.cpp (this code can be found in the CBC Samples directory, see Chapter 5, More Samples ). initialSolve is called a few times so although we will not gain much this is a simpler place to start. The example derives from OsiClpSolverInterface and the code is:

The resolve method is more complicated. The main pieces of data are a counter count_ which is incremented each solve and an int array node_ which stores the last time a variable was active in a solution. For the first few times normal dual is called and node_ array is updated.

After the first few solves we only use those which took part in a solution in the last so many solves. As fast0507 is a set covering problem we can also take out any rows which are already covered.

Example 3.12. Create small problem

    
    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 we know that the problem is feasible (We are not using cuts). If the rows were E rows then this might not be the case and we would have to do more work. When we have solved then we see if there are any negative reduced costs and if there are then we have to go to the full problem and use primal to clean up.

Example 3.13. Check optimal solution

    
    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 claen up with primal (and save some statistics)
    if (nBad) {
      timesBad_++;
      modelPtr_->primal(1);
      iterationsBad_ += modelPtr_->numberIterations();
    }
     
  

We then update node_ array as for the first few solves. To give some idea of the effect of this tactic fast0507 has 63,009 variables and the small problem never has more than 4,000 variables. In just over ten percent of solves did we have 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. In this case we make resolve the same as initialSolve. The full code is in ClpQuadInterface.hpp and ClpQuadInterface.cpp (this code can be found in the CBC Samples directory, see Chapter 5, More Samples ).