source: html/trunk/Cbc/ch03s04.html @ 554

Last change on this file since 554 was 554, checked in by rlh, 16 years ago

initial import of Cbc documentation

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