source: branches/gh-pages/advanced.md @ 2521

Last change on this file since 2521 was 2521, checked in by stefan, 2 months ago

migrate Cbc documentation to gh-pages

  • converted docbook to markdown
  • source now in gh-pages branch, so remove Cbc/doc
  • some cleanup
  • update link in README.md
  • Property svn:eol-style set to native
File size: 8.6 KB
Line 
1# Advanced Solver Uses
2
3## Creating a Solver via Inheritance
4
5CBC uses a generic `OsiSolverInterface` and its `resolve` capability.
6This does not give much flexibility so advanced users can inherit from
7their interface of choice. This section illustrates how to implement
8such a solver for a long thin problem, e.g., fast0507 again. As with the
9other examples in the Guide, the sample code is not guaranteed to be the
10fastest way to solve the problem. The main purpose of the example is to
11illustrate techniques. The full source is in `CbcSolver2.hpp` and
12`CbcSolver2.cpp` located in the CBC Samples directory.
13
14The method `initialSolve` is called a few times in CBC, and provides a
15convenient starting point. The `modelPtr_` derives from
16`OsiClpSolverInterface`.
17
18```
19// modelPtr_ is of type ClpSimplex *
20modelPtr_->setLogLevel(1); // switch on a bit of printout
21modelPtr_->scaling(0); // We don't want scaling for fast0507
22setBasis(basis_,modelPtr_); // Put basis into ClpSimplex
23// Do long thin by sprint
24ClpSolve options;
25options.setSolveType(ClpSolve::usePrimalorSprint);
26options.setPresolveType(ClpSolve::presolveOff);
27options.setSpecialOption(1,3,15); // Do 15 sprint iterations
28modelPtr_->initialSolve(options); // solve problem
29basis_ = getBasis(modelPtr_); // save basis
30modelPtr_->setLogLevel(0); // switch off printout
31```
32
33The `resolve()` method is more complicated than `initialSolve()`. The
34main pieces of data are a counter `count_` which is incremented each
35solve and an integer array `node_` which stores the last time a variable
36was active in a solution. For the first few times, the normal Dual
37Simplex is called and `node_` array is updated.
38
39```
40if (count_<10) {
41  OsiClpSolverInterface::resolve(); // Normal resolve
42  if (modelPtr_->status()==0) {
43    count_++; // feasible - save any nonzero or basic
44    const double * solution = modelPtr_->primalColumnSolution();
45    for (int i=0;i<numberColumns;i++) {
46      if (solution[i]>1.0e-6||modelPtr_->getStatus(i)==ClpSimplex::basic) {
47        node_[i]=CoinMax(count_,node_[i]);
48        howMany_[i]++;
49      }
50    }
51  } else {
52    printf("infeasible early on\n");
53  }
54}
55```
56
57After the first few solves, only those variables which took part in a
58solution in the last so many solves are used. As fast0507 is a set
59covering problem, any rows which are already covered can be taken out.
60```
61int * whichRow = new int[numberRows]; // Array to say which rows used
62int * whichColumn = new int [numberColumns]; // Array to say which columns used
63int i;
64const double * lower = modelPtr_->columnLower();
65const double * upper = modelPtr_->columnUpper();
66setBasis(basis_,modelPtr_); // Set basis
67int nNewCol=0; // Number of columns in small model
68// Column copy of matrix
69const double * element = modelPtr_->matrix()->getElements();
70const int * row = modelPtr_->matrix()->getIndices();
71const CoinBigIndex * columnStart = modelPtr_->matrix()->getVectorStarts();
72const int * columnLength = modelPtr_->matrix()->getVectorLengths();
73
74int * rowActivity = new int[numberRows]; // Number of columns with entries in each row
75memset(rowActivity,0,numberRows*sizeof(int));
76int * rowActivity2 = new int[numberRows]; // Lower bound on row activity for each row
77memset(rowActivity2,0,numberRows*sizeof(int));
78char * mark = (char *) modelPtr_->dualColumnSolution(); // Get some space to mark columns
79memset(mark,0,numberColumns);
80for (i=0;i<numberColumns;i++) {
81  bool choose = (node_[i]>count_-memory_&&node_[i]>0); // Choose if used recently
82  // Take if used recently or active in some sense
83  if ((choose&&upper[i])||(modelPtr_->getStatus(i)!=ClpSimplex::atLowerBound&&
84       modelPtr_->getStatus(i)!=ClpSimplex::isFixed)||lower[i]>0.0) {
85    mark[i]=1; // mark as used
86    whichColumn[nNewCol++]=i; // add to list
87    CoinBigIndex j;
88    double value = upper[i];
89    if (value) {
90      for (j=columnStart[i];
91           j<columnStart[i]+columnLength[i];j++) {
92        int iRow=row[j];
93        assert (element[j]==1.0);
94        rowActivity[iRow] ++; // This variable can cover this row
95      }
96      if (lower[i]>0.0) {
97        for (j=columnStart[i];
98             j<columnStart[i]+columnLength[i];j++) {
99          int iRow=row[j];
100          rowActivity2[iRow] ++; // This row redundant
101        }
102      }
103    }
104  }
105}
106int nOK=0; // Use to count rows which can be covered
107int nNewRow=0; // Use to make list of rows needed
108for (i=0;i<numberRows;i++) {
109  if (rowActivity[i])
110    nOK++;
111  if (!rowActivity2[i])
112    whichRow[nNewRow++]=i; // not satisfied
113  else
114    modelPtr_->setRowStatus(i,ClpSimplex::basic); // make slack basic
115}
116if (nOK<numberRows) {
117  // The variables we have do not cover rows - see if we can find any that do
118  for (i=0;i<numberColumns;i++) {
119    if (!mark[i]&&upper[i]) {
120      CoinBigIndex j;
121      int good=0;
122      for (j=columnStart[i];
123           j<columnStart[i]+columnLength[i];j++) {
124        int iRow=row[j];
125        if (!rowActivity[iRow]) {
126          rowActivity[iRow] ++;
127          good++;
128        }
129      }
130      if (good) {
131        nOK+=good; // This covers - put in list
132        whichColumn[nNewCol++]=i;
133      }
134    }
135  }
136}
137delete [] rowActivity;
138delete [] rowActivity2;
139if (nOK<numberRows) {
140  // By inspection the problem is infeasible - no need to solve
141  modelPtr_->setProblemStatus(1);
142  delete [] whichRow;
143  delete [] whichColumn;
144  printf("infeasible by inspection\n");
145  return;
146}
147// Now make up a small model with the right rows and columns
148ClpSimplex *  temp = new ClpSimplex(modelPtr_,nNewRow,whichRow,nNewCol,whichColumn);
149```
150
151If the variables cover the rows, then the problem is feasible (no cuts
152are being used). If the rows were equality constraints, then this might
153not be the case. More work would be needed. After the solution, the
154reduct costs are checked. If any reduced costs are negative, the code
155goes back to the full problem and cleans up with Primal Simplex.
156
157```
158temp->setDualObjectiveLimit(1.0e50); // Switch off dual cutoff as problem is restricted
159temp->dual(); // solve
160double * solution = modelPtr_->primalColumnSolution(); // put back solution
161const double * solution2 = temp->primalColumnSolution();
162memset(solution,0,numberColumns*sizeof(double));
163for (i=0;i<nNewCol;i++) {
164  int iColumn = whichColumn[i];
165  solution[iColumn]=solution2[i];
166  modelPtr_->setStatus(iColumn,temp->getStatus(i));
167}
168double * rowSolution = modelPtr_->primalRowSolution();
169const double * rowSolution2 = temp->primalRowSolution();
170double * dual = modelPtr_->dualRowSolution();
171const double * dual2 = temp->dualRowSolution();
172memset(dual,0,numberRows*sizeof(double));
173for (i=0;i<nNewRow;i++) {
174  int iRow=whichRow[i];
175  modelPtr_->setRowStatus(iRow,temp->getRowStatus(i));
176  rowSolution[iRow]=rowSolution2[i];
177  dual[iRow]=dual2[i];
178}
179// See if optimal
180double * dj = modelPtr_->dualColumnSolution();
181// get reduced cost for large problem
182// this assumes minimization
183memcpy(dj,modelPtr_->objective(),numberColumns*sizeof(double));
184modelPtr_->transposeTimes(-1.0,dual,dj);
185modelPtr_->setObjectiveValue(temp->objectiveValue());
186modelPtr_->setProblemStatus(0);
187int nBad=0;
188
189for (i=0;i<numberColumns;i++) {
190  if (modelPtr_->getStatus(i)==ClpSimplex::atLowerBound
191      &&upper[i]>lower[i]&&dj[i]<-1.0e-5)
192    nBad++;
193}
194// If necessary clean up with primal (and save some statistics)
195if (nBad) {
196  timesBad_++;
197  modelPtr_->primal(1);
198  iterationsBad_ += modelPtr_->numberIterations();
199}
200```
201
202The array `node_` is updated, as for the first few solves. To give some
203idea of the effect of this tactic, the problem fast0507 has 63,009
204variables but the small problem never has more than 4,000 variables. In
205only about ten percent of solves was it necessary to resolve, and then
206the average number of iterations on full problem was less than 20.
207
208## Quadratic MIP
209
210To give another example - again only for illustrative purposes -- it is
211possible to do quadratic MIP with CBC. In this case, we make `resolve`
212the same as `initialSolve`. The full code is in `ClpQuadInterface.hpp`
213and `ClpQuadInterface.cpp` located in the CBC Samples directory.
214
215```
216// save cutoff
217double cutoff = modelPtr_->dualObjectiveLimit();
218modelPtr_->setDualObjectiveLimit(1.0e50);
219modelPtr_->scaling(0);
220modelPtr_->setLogLevel(0);
221// solve with no objective to get feasible solution
222setBasis(basis_,modelPtr_);
223modelPtr_->dual();
224basis_ = getBasis(modelPtr_);
225modelPtr_->setDualObjectiveLimit(cutoff);
226if (modelPtr_->problemStatus())
227  return; // problem was infeasible
228// Now pass in quadratic objective
229ClpObjective * saveObjective  = modelPtr_->objectiveAsObject();
230modelPtr_->setObjectivePointer(quadraticObjective_);
231modelPtr_->primal();
232modelPtr_->setDualObjectiveLimit(cutoff);
233if (modelPtr_->objectiveValue()>cutoff)
234  modelPtr_->setProblemStatus(1);
235modelPtr_->setObjectivePointer(saveObjective);
236```
Note: See TracBrowser for help on using the repository browser.