source: html/ch07.html @ 1591

Last change on this file since 1591 was 561, checked in by ladanyi, 14 years ago

Cbc static pages moved
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 11.3 KB
Line 
1<?xml version="1.0" encoding="ISO-8859-1" standalone="no"?>
2<html xmlns="http://www.w3.org/1999/xhtml"><head><title>Chapter 7. 
3  Advanced Solver Uses
4</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="index.html" title="CBC User Guide"/><link rel="previous" href="ch06.html" title="Chapter 6. Cutting planes"/><link rel="next" href="ch07s02.html" title="Quadratic MIP"/></head><body><div class="navheader"><table width="100%" summary="Navigation header"><tr><th colspan="3" align="center">Chapter 7. 
5  Advanced Solver Uses
6</th></tr><tr><td width="20%" align="left"><a accesskey="p" href="ch06.html">Prev</a> </td><th width="60%" align="center"> </th><td width="20%" align="right"> <a accesskey="n" href="ch07s02.html">Next</a></td></tr></table><hr/></div><div class="chapter" lang="en"><div class="titlepage"><div><div><h2 class="title"><a id="SolverChap"/>Chapter 7. 
7  Advanced Solver Uses
8</h2></div></div><div/></div><div class="toc"><p><b>Table of Contents</b></p><dl><dt><a href="ch07.html#solver">Creating a Solver via Inheritance</a></dt><dt><a href="ch07s02.html">Quadratic MIP</a></dt></dl></div><div class="section" lang="en"><div class="titlepage"><div><div><h2 class="title" style="clear: both"><a id="solver"/>Creating a Solver via Inheritance</h2></div></div><div/></div><p>
9  CBC uses a generic <tt class="classname">OsiSolverInterface</tt> and its <tt class="function">resolve</tt> capability.
10  This does not give much flexibility so advanced users can inherit from their interface
11  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 <tt class="filename">CbcSolver2.hpp</tt> and <tt class="filename">CbcSolver2.cpp</tt> located in the CBC Samples directory, see
12  <a href="ch08.html" title="Chapter 8. &#10;More Samples&#10;">Chapter 8, <i>
13More Samples
14</i></a>.
15</p><p>
16The method <tt class="function">initialSolve</tt> is called a few times in CBC, and provides a convenient starting point. The <tt class="varname">modelPtr_</tt> derives from <tt class="classname">OsiClpSolverInterface</tt>.
17  </p><div class="example"><a id="initialSolve"/><p class="title"><b>Example 7.1. <tt class="function">initialSolve()</tt></b></p><pre class="programlisting">
18   
19  // modelPtr_ is of type ClpSimplex *
20  modelPtr_-&gt;setLogLevel(1); // switch on a bit of printout
21  modelPtr_-&gt;scaling(0); // We don't want scaling for fast0507
22  setBasis(basis_,modelPtr_); // Put basis into ClpSimplex
23  // Do long thin by sprint
24  ClpSolve options;
25  options.setSolveType(ClpSolve::usePrimalorSprint);
26  options.setPresolveType(ClpSolve::presolveOff);
27  options.setSpecialOption(1,3,15); // Do 15 sprint iterations
28  modelPtr_-&gt;initialSolve(options); // solve problem
29  basis_ = getBasis(modelPtr_); // save basis
30  modelPtr_-&gt;setLogLevel(0); // switch off printout
31     
32  </pre></div><p>
33The <tt class="function">resolve()</tt> method is more complicated than <tt class="function">initialSolve()</tt>.  The main pieces of data are a counter <tt class="varname">count_</tt> which is incremented each solve and an integer array <tt class="varname">node_</tt> which stores the last time
34a variable was active in a solution.  For the first few solves, the normal Dual Simplex is called and
35<tt class="varname">node_</tt> array is updated.
36</p><div class="example"><a id="id3000834"/><p class="title"><b>Example 7.2. First Few Solves</b></p><pre class="programlisting">
37   
38  if (count_&lt;10) {
39    OsiClpSolverInterface::resolve(); // Normal resolve
40    if (modelPtr_-&gt;status()==0) {
41      count_++; // feasible - save any nonzero or basic
42      const double * solution = modelPtr_-&gt;primalColumnSolution();
43      for (int i=0;i&lt;numberColumns;i++) {
44        if (solution[i]&gt;1.0e-6||modelPtr_-&gt;getStatus(i)==ClpSimplex::basic) {
45          node_[i]=CoinMax(count_,node_[i]);
46          howMany_[i]++;
47        }
48      }
49    } else {
50      printf("infeasible early on\n");
51    }
52  }
53     
54  </pre></div><p>
55After the first few solves, only those variables which took part in a solution in the last so many
56solves are used.  As fast0507 is a set covering problem, any rows which are already covered can be taken out.
57  </p><div class="example"><a id="id3000862"/><p class="title"><b>Example 7.3. Create Small Sub-Problem</b></p><pre class="programlisting">
58   
59    int * whichRow = new int[numberRows]; // Array to say which rows used
60    int * whichColumn = new int [numberColumns]; // Array to say which columns used
61    int i;
62    const double * lower = modelPtr_-&gt;columnLower();
63    const double * upper = modelPtr_-&gt;columnUpper();
64    setBasis(basis_,modelPtr_); // Set basis
65    int nNewCol=0; // Number of columns in small model
66    // Column copy of matrix
67    const double * element = modelPtr_-&gt;matrix()-&gt;getElements();
68    const int * row = modelPtr_-&gt;matrix()-&gt;getIndices();
69    const CoinBigIndex * columnStart = modelPtr_-&gt;matrix()-&gt;getVectorStarts();
70    const int * columnLength = modelPtr_-&gt;matrix()-&gt;getVectorLengths();
71   
72    int * rowActivity = new int[numberRows]; // Number of columns with entries in each row
73    memset(rowActivity,0,numberRows*sizeof(int));
74    int * rowActivity2 = new int[numberRows]; // Lower bound on row activity for each row
75    memset(rowActivity2,0,numberRows*sizeof(int));
76    char * mark = (char *) modelPtr_-&gt;dualColumnSolution(); // Get some space to mark columns
77    memset(mark,0,numberColumns);
78    for (i=0;i&lt;numberColumns;i++) {
79      bool choose = (node_[i]&gt;count_-memory_&amp;&amp;node_[i]&gt;0); // Choose if used recently
80      // Take if used recently or active in some sense
81      if ((choose&amp;&amp;upper[i])
82          ||(modelPtr_-&gt;getStatus(i)!=ClpSimplex::atLowerBound&amp;&amp;
83             modelPtr_-&gt;getStatus(i)!=ClpSimplex::isFixed)
84          ||lower[i]&gt;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&lt;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]&gt;0.0) {
97            for (j=columnStart[i];
98                 j&lt;columnStart[i]+columnLength[i];j++) {
99              int iRow=row[j];
100              rowActivity2[iRow] ++; // This row redundant
101            }
102          }
103        }
104      }
105    }
106    int nOK=0; // Use to count rows which can be covered
107    int nNewRow=0; // Use to make list of rows needed
108    for (i=0;i&lt;numberRows;i++) {
109      if (rowActivity[i])
110        nOK++;
111      if (!rowActivity2[i])
112        whichRow[nNewRow++]=i; // not satisfied
113      else
114        modelPtr_-&gt;setRowStatus(i,ClpSimplex::basic); // make slack basic
115    }
116    if (nOK&lt;numberRows) {
117      // The variables we have do not cover rows - see if we can find any that do
118      for (i=0;i&lt;numberColumns;i++) {
119        if (!mark[i]&amp;&amp;upper[i]) {
120          CoinBigIndex j;
121          int good=0;
122          for (j=columnStart[i];
123               j&lt;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    }
137    delete [] rowActivity;
138    delete [] rowActivity2;
139    if (nOK&lt;numberRows) {
140      // By inspection the problem is infeasible - no need to solve
141      modelPtr_-&gt;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
148    ClpSimplex *  temp = new ClpSimplex(modelPtr_,nNewRow,whichRow,nNewCol,whichColumn);
149     
150  </pre></div><p>
151If the variables cover the rows, then the problem is feasible (no cuts are being used). (If the rows
152were 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.
153  </p><div class="example"><a id="id3000908"/><p class="title"><b>Example 7.4. Check Optimal Solution</b></p><pre class="programlisting">
154   
155    temp-&gt;setDualObjectiveLimit(1.0e50); // Switch off dual cutoff as problem is restricted
156    temp-&gt;dual(); // solve
157    double * solution = modelPtr_-&gt;primalColumnSolution(); // put back solution
158    const double * solution2 = temp-&gt;primalColumnSolution();
159    memset(solution,0,numberColumns*sizeof(double));
160    for (i=0;i&lt;nNewCol;i++) {
161      int iColumn = whichColumn[i];
162      solution[iColumn]=solution2[i];
163      modelPtr_-&gt;setStatus(iColumn,temp-&gt;getStatus(i));
164    }
165    double * rowSolution = modelPtr_-&gt;primalRowSolution();
166    const double * rowSolution2 = temp-&gt;primalRowSolution();
167    double * dual = modelPtr_-&gt;dualRowSolution();
168    const double * dual2 = temp-&gt;dualRowSolution();
169    memset(dual,0,numberRows*sizeof(double));
170    for (i=0;i&lt;nNewRow;i++) {
171      int iRow=whichRow[i];
172      modelPtr_-&gt;setRowStatus(iRow,temp-&gt;getRowStatus(i));
173      rowSolution[iRow]=rowSolution2[i];
174      dual[iRow]=dual2[i];
175    }
176    // See if optimal
177    double * dj = modelPtr_-&gt;dualColumnSolution();
178    // get reduced cost for large problem
179    // this assumes minimization
180    memcpy(dj,modelPtr_-&gt;objective(),numberColumns*sizeof(double));
181    modelPtr_-&gt;transposeTimes(-1.0,dual,dj);
182    modelPtr_-&gt;setObjectiveValue(temp-&gt;objectiveValue());
183    modelPtr_-&gt;setProblemStatus(0);
184    int nBad=0;
185     
186    for (i=0;i&lt;numberColumns;i++) {
187      if (modelPtr_-&gt;getStatus(i)==ClpSimplex::atLowerBound
188          &amp;&amp;upper[i]&gt;lower[i]&amp;&amp;dj[i]&lt;-1.0e-5)
189        nBad++;
190    }
191    // If necessary clean up with primal (and save some statistics)
192    if (nBad) {
193      timesBad_++;
194      modelPtr_-&gt;primal(1);
195      iterationsBad_ += modelPtr_-&gt;numberIterations();
196    }
197     
198  </pre></div><p>
199The array <tt class="varname">node_</tt> 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
200on full problem was less than 20.
201</p></div></div><div class="navfooter"><hr/><table width="100%" summary="Navigation footer"><tr><td width="40%" align="left"><a accesskey="p" href="ch06.html">Prev</a> </td><td width="20%" align="center"><a accesskey="u" href="index.html">Up</a></td><td width="40%" align="right"> <a accesskey="n" href="ch07s02.html">Next</a></td></tr><tr><td width="40%" align="left" valign="top">Chapter 6. Cutting planes </td><td width="20%" align="center"><a accesskey="h" href="index.html">Home</a></td><td width="40%" align="right" valign="top"> Quadratic MIP</td></tr></table></div></body></html>
Note: See TracBrowser for help on using the repository browser.