| 1 | <?xml version="1.0" encoding="ISO-8859-1" standalone="no"?> |
|---|
| 2 | <html xmlns="http://www.w3.org/1999/xhtml"><head><title>Advanced use of solver</title><meta name="generator" content="DocBook XSL Stylesheets V1.61.2"/><link rel="home" href="index.html" title="CBC User Guide"/><link rel="up" href="ch04.html" title="Chapter 4. Branching "/><link rel="previous" href="ch04.html" title="Chapter 4. Branching "/><link rel="next" href="ch05.html" title="Chapter 5. Building and Modifying a Model "/></head><body><div class="navheader"><table width="100%" summary="Navigation header"><tr><th colspan="3" align="center">Advanced use of solver</th></tr><tr><td width="20%" align="left"><a accesskey="p" href="ch04.html">Prev</a> </td><th width="60%" align="center">Chapter 4. |
|---|
| 3 | Branching |
|---|
| 4 | </th><td width="20%" align="right"> <a accesskey="n" href="ch05.html">Next</a></td></tr></table><hr/></div><div class="section" lang="en"><div class="titlepage"><div><div><h2 class="title" style="clear: both"><a id="solver"/>Advanced use of solver</h2></div></div><div/></div><p> |
|---|
| 5 | Coin Branch and Cut uses a generic OsiSolverInterface and its <tt class="function">resolve</tt> capability. |
|---|
| 6 | This does not give much flexibility so advanced users can inherit from the interface |
|---|
| 7 | of choice. This describes such a solver for a long thin problem e.g. fast0507 again. |
|---|
| 8 | As with all these examples it is not guaranteed that this is the fastest way to solve |
|---|
| 9 | any of these problems - they are to illustrate techniques. |
|---|
| 10 | The full code is in <tt class="filename">CbcSolver2.hpp</tt> and |
|---|
| 11 | <tt class="filename">CbcSolver2.cpp</tt> |
|---|
| 12 | (this code can be found in the CBC Samples directory, see |
|---|
| 13 | <a href="ch06.html" title="Chapter 6. More Samples ">Chapter 6, <i> |
|---|
| 14 | More Samples |
|---|
| 15 | </i></a>). |
|---|
| 16 | <tt class="function">initialSolve</tt> is called a few times so although we will not gain much |
|---|
| 17 | this is a simpler place to start. The example derives from OsiClpSolverInterface and the code |
|---|
| 18 | is: |
|---|
| 19 | </p><div class="example"><a id="id2985361"/><p class="title"><b>Example 4.3. initialSolve</b></p><pre class="programlisting"> |
|---|
| 20 | |
|---|
| 21 | // modelPtr_ is of type ClpSimplex * |
|---|
| 22 | modelPtr_->setLogLevel(1); // switch on a bit of printout |
|---|
| 23 | modelPtr_->scaling(0); // We don't want scaling for fast0507 |
|---|
| 24 | setBasis(basis_,modelPtr_); // Put basis into ClpSimplex |
|---|
| 25 | // Do long thin by sprint |
|---|
| 26 | ClpSolve options; |
|---|
| 27 | options.setSolveType(ClpSolve::usePrimalorSprint); |
|---|
| 28 | options.setPresolveType(ClpSolve::presolveOff); |
|---|
| 29 | options.setSpecialOption(1,3,15); // Do 15 sprint iterations |
|---|
| 30 | modelPtr_->initialSolve(options); // solve problem |
|---|
| 31 | basis_ = getBasis(modelPtr_); // save basis |
|---|
| 32 | modelPtr_->setLogLevel(0); // switch off printout |
|---|
| 33 | |
|---|
| 34 | </pre></div><p> |
|---|
| 35 | The <tt class="function">resolve</tt> method is more complicated. The main pieces of data are |
|---|
| 36 | a counter count_ which is incremented each solve and an int array node_ which stores the last time |
|---|
| 37 | a variable was active in a solution. For the first few times normal dual is called and |
|---|
| 38 | node_ array is updated. |
|---|
| 39 | </p><div class="example"><a id="id2985410"/><p class="title"><b>Example 4.4. First few solves</b></p><pre class="programlisting"> |
|---|
| 40 | |
|---|
| 41 | if (count_<10) { |
|---|
| 42 | OsiClpSolverInterface::resolve(); // Normal resolve |
|---|
| 43 | if (modelPtr_->status()==0) { |
|---|
| 44 | count_++; // feasible - save any nonzero or basic |
|---|
| 45 | const double * solution = modelPtr_->primalColumnSolution(); |
|---|
| 46 | for (int i=0;i<numberColumns;i++) { |
|---|
| 47 | if (solution[i]>1.0e-6||modelPtr_->getStatus(i)==ClpSimplex::basic) { |
|---|
| 48 | node_[i]=CoinMax(count_,node_[i]); |
|---|
| 49 | howMany_[i]++; |
|---|
| 50 | } |
|---|
| 51 | } |
|---|
| 52 | } else { |
|---|
| 53 | printf("infeasible early on\n"); |
|---|
| 54 | } |
|---|
| 55 | } |
|---|
| 56 | |
|---|
| 57 | </pre></div><p> |
|---|
| 58 | After the first few solves we only use those which took part in a solution in the last so many |
|---|
| 59 | solves. As fast0507 is a set covering problem we can also take out any rows which are |
|---|
| 60 | already covered. |
|---|
| 61 | </p><div class="example"><a id="id2985437"/><p class="title"><b>Example 4.5. Create small problem</b></p><pre class="programlisting"> |
|---|
| 62 | |
|---|
| 63 | int * whichRow = new int[numberRows]; // Array to say which rows used |
|---|
| 64 | int * whichColumn = new int [numberColumns]; // Array to say which columns used |
|---|
| 65 | int i; |
|---|
| 66 | const double * lower = modelPtr_->columnLower(); |
|---|
| 67 | const double * upper = modelPtr_->columnUpper(); |
|---|
| 68 | setBasis(basis_,modelPtr_); // Set basis |
|---|
| 69 | int nNewCol=0; // Number of columns in small model |
|---|
| 70 | // Column copy of matrix |
|---|
| 71 | const double * element = modelPtr_->matrix()->getElements(); |
|---|
| 72 | const int * row = modelPtr_->matrix()->getIndices(); |
|---|
| 73 | const CoinBigIndex * columnStart = modelPtr_->matrix()->getVectorStarts(); |
|---|
| 74 | const int * columnLength = modelPtr_->matrix()->getVectorLengths(); |
|---|
| 75 | |
|---|
| 76 | int * rowActivity = new int[numberRows]; // Number of columns with entries in each row |
|---|
| 77 | memset(rowActivity,0,numberRows*sizeof(int)); |
|---|
| 78 | int * rowActivity2 = new int[numberRows]; // Lower bound on row activity for each row |
|---|
| 79 | memset(rowActivity2,0,numberRows*sizeof(int)); |
|---|
| 80 | char * mark = (char *) modelPtr_->dualColumnSolution(); // Get some space to mark columns |
|---|
| 81 | memset(mark,0,numberColumns); |
|---|
| 82 | for (i=0;i<numberColumns;i++) { |
|---|
| 83 | bool choose = (node_[i]>count_-memory_&&node_[i]>0); // Choose if used recently |
|---|
| 84 | // Take if used recently or active in some sense |
|---|
| 85 | if ((choose&&upper[i]) |
|---|
| 86 | ||(modelPtr_->getStatus(i)!=ClpSimplex::atLowerBound&& |
|---|
| 87 | modelPtr_->getStatus(i)!=ClpSimplex::isFixed) |
|---|
| 88 | ||lower[i]>0.0) { |
|---|
| 89 | mark[i]=1; // mark as used |
|---|
| 90 | whichColumn[nNewCol++]=i; // add to list |
|---|
| 91 | CoinBigIndex j; |
|---|
| 92 | double value = upper[i]; |
|---|
| 93 | if (value) { |
|---|
| 94 | for (j=columnStart[i]; |
|---|
| 95 | j<columnStart[i]+columnLength[i];j++) { |
|---|
| 96 | int iRow=row[j]; |
|---|
| 97 | assert (element[j]==1.0); |
|---|
| 98 | rowActivity[iRow] ++; // This variable can cover this row |
|---|
| 99 | } |
|---|
| 100 | if (lower[i]>0.0) { |
|---|
| 101 | for (j=columnStart[i]; |
|---|
| 102 | j<columnStart[i]+columnLength[i];j++) { |
|---|
| 103 | int iRow=row[j]; |
|---|
| 104 | rowActivity2[iRow] ++; // This row redundant |
|---|
| 105 | } |
|---|
| 106 | } |
|---|
| 107 | } |
|---|
| 108 | } |
|---|
| 109 | } |
|---|
| 110 | int nOK=0; // Use to count rows which can be covered |
|---|
| 111 | int nNewRow=0; // Use to make list of rows needed |
|---|
| 112 | for (i=0;i<numberRows;i++) { |
|---|
| 113 | if (rowActivity[i]) |
|---|
| 114 | nOK++; |
|---|
| 115 | if (!rowActivity2[i]) |
|---|
| 116 | whichRow[nNewRow++]=i; // not satisfied |
|---|
| 117 | else |
|---|
| 118 | modelPtr_->setRowStatus(i,ClpSimplex::basic); // make slack basic |
|---|
| 119 | } |
|---|
| 120 | if (nOK<numberRows) { |
|---|
| 121 | // The variables we have do not cover rows - see if we can find any that do |
|---|
| 122 | for (i=0;i<numberColumns;i++) { |
|---|
| 123 | if (!mark[i]&&upper[i]) { |
|---|
| 124 | CoinBigIndex j; |
|---|
| 125 | int good=0; |
|---|
| 126 | for (j=columnStart[i]; |
|---|
| 127 | j<columnStart[i]+columnLength[i];j++) { |
|---|
| 128 | int iRow=row[j]; |
|---|
| 129 | if (!rowActivity[iRow]) { |
|---|
| 130 | rowActivity[iRow] ++; |
|---|
| 131 | good++; |
|---|
| 132 | } |
|---|
| 133 | } |
|---|
| 134 | if (good) { |
|---|
| 135 | nOK+=good; // This covers - put in list |
|---|
| 136 | whichColumn[nNewCol++]=i; |
|---|
| 137 | } |
|---|
| 138 | } |
|---|
| 139 | } |
|---|
| 140 | } |
|---|
| 141 | delete [] rowActivity; |
|---|
| 142 | delete [] rowActivity2; |
|---|
| 143 | if (nOK<numberRows) { |
|---|
| 144 | // By inspection the problem is infeasible - no need to solve |
|---|
| 145 | modelPtr_->setProblemStatus(1); |
|---|
| 146 | delete [] whichRow; |
|---|
| 147 | delete [] whichColumn; |
|---|
| 148 | printf("infeasible by inspection\n"); |
|---|
| 149 | return; |
|---|
| 150 | } |
|---|
| 151 | // Now make up a small model with the right rows and columns |
|---|
| 152 | ClpSimplex * temp = new ClpSimplex(modelPtr_,nNewRow,whichRow,nNewCol,whichColumn); |
|---|
| 153 | |
|---|
| 154 | </pre></div><p> |
|---|
| 155 | If the variables cover the rows then we know that the problem is feasible (We are not using |
|---|
| 156 | cuts). If the rows |
|---|
| 157 | were E rows then this might not be the case and we would have to do more work. When we have solved |
|---|
| 158 | then we see if there are any negative reduced costs and if there are then we have to go to the |
|---|
| 159 | full problem and use primal to clean up. |
|---|
| 160 | </p><div class="example"><a id="id2985478"/><p class="title"><b>Example 4.6. Check optimal solution</b></p><pre class="programlisting"> |
|---|
| 161 | |
|---|
| 162 | temp->setDualObjectiveLimit(1.0e50); // Switch off dual cutoff as problem is restricted |
|---|
| 163 | temp->dual(); // solve |
|---|
| 164 | double * solution = modelPtr_->primalColumnSolution(); // put back solution |
|---|
| 165 | const double * solution2 = temp->primalColumnSolution(); |
|---|
| 166 | memset(solution,0,numberColumns*sizeof(double)); |
|---|
| 167 | for (i=0;i<nNewCol;i++) { |
|---|
| 168 | int iColumn = whichColumn[i]; |
|---|
| 169 | solution[iColumn]=solution2[i]; |
|---|
| 170 | modelPtr_->setStatus(iColumn,temp->getStatus(i)); |
|---|
| 171 | } |
|---|
| 172 | double * rowSolution = modelPtr_->primalRowSolution(); |
|---|
| 173 | const double * rowSolution2 = temp->primalRowSolution(); |
|---|
| 174 | double * dual = modelPtr_->dualRowSolution(); |
|---|
| 175 | const double * dual2 = temp->dualRowSolution(); |
|---|
| 176 | memset(dual,0,numberRows*sizeof(double)); |
|---|
| 177 | for (i=0;i<nNewRow;i++) { |
|---|
| 178 | int iRow=whichRow[i]; |
|---|
| 179 | modelPtr_->setRowStatus(iRow,temp->getRowStatus(i)); |
|---|
| 180 | rowSolution[iRow]=rowSolution2[i]; |
|---|
| 181 | dual[iRow]=dual2[i]; |
|---|
| 182 | } |
|---|
| 183 | // See if optimal |
|---|
| 184 | double * dj = modelPtr_->dualColumnSolution(); |
|---|
| 185 | // get reduced cost for large problem |
|---|
| 186 | // this assumes minimization |
|---|
| 187 | memcpy(dj,modelPtr_->objective(),numberColumns*sizeof(double)); |
|---|
| 188 | modelPtr_->transposeTimes(-1.0,dual,dj); |
|---|
| 189 | modelPtr_->setObjectiveValue(temp->objectiveValue()); |
|---|
| 190 | modelPtr_->setProblemStatus(0); |
|---|
| 191 | int nBad=0; |
|---|
| 192 | |
|---|
| 193 | for (i=0;i<numberColumns;i++) { |
|---|
| 194 | if (modelPtr_->getStatus(i)==ClpSimplex::atLowerBound |
|---|
| 195 | &&upper[i]>lower[i]&&dj[i]<-1.0e-5) |
|---|
| 196 | nBad++; |
|---|
| 197 | } |
|---|
| 198 | // If necessary claen up with primal (and save some statistics) |
|---|
| 199 | if (nBad) { |
|---|
| 200 | timesBad_++; |
|---|
| 201 | modelPtr_->primal(1); |
|---|
| 202 | iterationsBad_ += modelPtr_->numberIterations(); |
|---|
| 203 | } |
|---|
| 204 | |
|---|
| 205 | </pre></div><p> |
|---|
| 206 | We then update node_ array as for the first few solves. To give some idea of the effect of this |
|---|
| 207 | tactic fast0507 has 63,009 variables and the small problem never has more than 4,000 variables. |
|---|
| 208 | In just over ten percent of solves did we have to resolve and then the average number of iterations |
|---|
| 209 | on full problem was less than 20. |
|---|
| 210 | To give another example - again only for illustrative purposes it is possible to do quadratic |
|---|
| 211 | mip. In this case we make <tt class="function">resolve</tt> the same as |
|---|
| 212 | <tt class="function">initialSolve</tt>. |
|---|
| 213 | The full code is in <tt class="filename">ClpQuadInterface.hpp</tt> and |
|---|
| 214 | <tt class="filename">ClpQuadInterface.cpp</tt> |
|---|
| 215 | (this code can be found in the CBC Samples directory, see |
|---|
| 216 | <a href="ch06.html" title="Chapter 6. More Samples ">Chapter 6, <i> |
|---|
| 217 | More Samples |
|---|
| 218 | </i></a>). |
|---|
| 219 | </p><div class="example"><a id="id2985520"/><p class="title"><b>Example 4.7. Solve a quadratic mip</b></p><pre class="programlisting"> |
|---|
| 220 | |
|---|
| 221 | // save cutoff |
|---|
| 222 | double cutoff = modelPtr_->dualObjectiveLimit(); |
|---|
| 223 | modelPtr_->setDualObjectiveLimit(1.0e50); |
|---|
| 224 | modelPtr_->scaling(0); |
|---|
| 225 | modelPtr_->setLogLevel(0); |
|---|
| 226 | // solve with no objective to get feasible solution |
|---|
| 227 | setBasis(basis_,modelPtr_); |
|---|
| 228 | modelPtr_->dual(); |
|---|
| 229 | basis_ = getBasis(modelPtr_); |
|---|
| 230 | modelPtr_->setDualObjectiveLimit(cutoff); |
|---|
| 231 | if (modelPtr_->problemStatus()) |
|---|
| 232 | return; // problem was infeasible |
|---|
| 233 | // Now pass in quadratic objective |
|---|
| 234 | ClpObjective * saveObjective = modelPtr_->objectiveAsObject(); |
|---|
| 235 | modelPtr_->setObjectivePointer(quadraticObjective_); |
|---|
| 236 | modelPtr_->primal(); |
|---|
| 237 | modelPtr_->setDualObjectiveLimit(cutoff); |
|---|
| 238 | if (modelPtr_->objectiveValue()>cutoff) |
|---|
| 239 | modelPtr_->setProblemStatus(1); |
|---|
| 240 | modelPtr_->setObjectivePointer(saveObjective); |
|---|
| 241 | |
|---|
| 242 | </pre></div></div><div class="navfooter"><hr/><table width="100%" summary="Navigation footer"><tr><td width="40%" align="left"><a accesskey="p" href="ch04.html">Prev</a> </td><td width="20%" align="center"><a accesskey="u" href="ch04.html">Up</a></td><td width="40%" align="right"> <a accesskey="n" href="ch05.html">Next</a></td></tr><tr><td width="40%" align="left" valign="top">Chapter 4. |
|---|
| 243 | Branching |
|---|
| 244 | </td><td width="20%" align="center"><a accesskey="h" href="index.html">Home</a></td><td width="40%" align="right" valign="top"> Chapter 5. |
|---|
| 245 | Building and Modifying a Model |
|---|
| 246 | </td></tr></table></div></body></html> |
|---|
