source: stable/1.13/Clp/examples/dualCuts.cpp @ 1898

Last change on this file since 1898 was 1559, checked in by stefan, 10 years ago

merge split branch into trunk

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 12.3 KB
Line 
1/* $Id: dualCuts.cpp 1559 2010-06-05 19:42:36Z stefan $ */
2/* Copyright (C) 2003, International Business Machines Corporation
3   and others.  All Rights Reserved.
4
5   This sample program is designed to illustrate programming
6   techniques using CoinLP, has not been thoroughly tested
7   and comes without any warranty whatsoever.
8
9   You may copy, modify and distribute this sample program without
10   any restrictions whatsoever and without any payment to anyone.
11*/
12
13#include "ClpSimplex.hpp"
14#include "ClpPresolve.hpp"
15#include "ClpFactorization.hpp"
16#include "CoinSort.hpp"
17#include "CoinHelperFunctions.hpp"
18#include "CoinTime.hpp"
19#include "CoinMpsIO.hpp"
20#include <iomanip>
21
22int main(int argc, const char *argv[])
23{
24     ClpSimplex  model;
25     int status;
26     // Keep names
27     if (argc < 2) {
28          status = model.readMps("small.mps", true);
29     } else {
30          status = model.readMps(argv[1], false);
31     }
32     if (status)
33          exit(10);
34     /*
35       This driver implements a method of treating a problem as all cuts.
36       So it adds in all E rows, solves and then adds in violated rows.
37     */
38
39     double time1 = CoinCpuTime();
40     ClpSimplex * model2;
41     ClpPresolve pinfo;
42     int numberPasses = 5; // can change this
43     /* Use a tolerance of 1.0e-8 for feasibility, treat problem as
44        not being integer, do "numberpasses" passes and throw away names
45        in presolved model */
46     model2 = pinfo.presolvedModel(model, 1.0e-8, false, numberPasses, false);
47     if (!model2) {
48          fprintf(stderr, "ClpPresolve says %s is infeasible with tolerance of %g\n",
49                  argv[1], 1.0e-8);
50          fprintf(stdout, "ClpPresolve says %s is infeasible with tolerance of %g\n",
51                  argv[1], 1.0e-8);
52          // model was infeasible - maybe try again with looser tolerances
53          model2 = pinfo.presolvedModel(model, 1.0e-7, false, numberPasses, false);
54          if (!model2) {
55               fprintf(stderr, "ClpPresolve says %s is infeasible with tolerance of %g\n",
56                       argv[1], 1.0e-7);
57               fprintf(stdout, "ClpPresolve says %s is infeasible with tolerance of %g\n",
58                       argv[1], 1.0e-7);
59               exit(2);
60          }
61     }
62     // change factorization frequency from 200
63     model2->setFactorizationFrequency(100 + model2->numberRows() / 50);
64
65     int numberColumns = model2->numberColumns();
66     int originalNumberRows = model2->numberRows();
67
68     // We will need arrays to choose rows to add
69     double * weight = new double [originalNumberRows];
70     int * sort = new int [originalNumberRows];
71     int numberSort = 0;
72     char * take = new char [originalNumberRows];
73
74     const double * rowLower = model2->rowLower();
75     const double * rowUpper = model2->rowUpper();
76     int iRow, iColumn;
77     // Set up initial list
78     numberSort = 0;
79     for (iRow = 0; iRow < originalNumberRows; iRow++) {
80          weight[iRow] = 1.123e50;
81          if (rowLower[iRow] == rowUpper[iRow]) {
82               sort[numberSort++] = iRow;
83               weight[iRow] = 0.0;
84          }
85     }
86     numberSort /= 2;
87     // Just add this number of rows each time in small problem
88     int smallNumberRows = 2 * numberColumns;
89     smallNumberRows = CoinMin(smallNumberRows, originalNumberRows / 20);
90     // and pad out with random rows
91     double ratio = ((double)(smallNumberRows - numberSort)) / ((double) originalNumberRows);
92     for (iRow = 0; iRow < originalNumberRows; iRow++) {
93          if (weight[iRow] == 1.123e50 && CoinDrand48() < ratio)
94               sort[numberSort++] = iRow;
95     }
96     /* This is optional.
97        The best thing to do is to miss out random rows and do a set which makes dual feasible.
98        If that is not possible then make sure variables have bounds.
99
100        One way that normally works is to automatically tighten bounds.
101     */
102#if 0
103     // However for some we need to do anyway
104     double * columnLower = model2->columnLower();
105     double * columnUpper = model2->columnUpper();
106     for (iColumn = 0; iColumn < numberColumns; iColumn++) {
107          columnLower[iColumn] = CoinMax(-1.0e6, columnLower[iColumn]);
108          columnUpper[iColumn] = CoinMin(1.0e6, columnUpper[iColumn]);
109     }
110#endif
111     model2->tightenPrimalBounds(-1.0e4, true);
112     printf("%d rows in initial problem\n", numberSort);
113     double * fullSolution = model2->primalRowSolution();
114
115     // Just do this number of passes
116     int maxPass = 50;
117     // And take out slack rows until this pass
118     int takeOutPass = 30;
119     int iPass;
120
121     const int * start = model2->clpMatrix()->getVectorStarts();
122     const int * length = model2->clpMatrix()->getVectorLengths();
123     const int * row = model2->clpMatrix()->getIndices();
124     int * whichColumns = new int [numberColumns];
125     for (iColumn = 0; iColumn < numberColumns; iColumn++)
126          whichColumns[iColumn] = iColumn;
127     int numberSmallColumns = numberColumns;
128     for (iPass = 0; iPass < maxPass; iPass++) {
129          printf("Start of pass %d\n", iPass);
130          // Cleaner this way
131          std::sort(sort, sort + numberSort);
132          // Create small problem
133          ClpSimplex small(model2, numberSort, sort, numberSmallColumns, whichColumns);
134          small.setFactorizationFrequency(100 + numberSort / 200);
135          //small.setPerturbation(50);
136          //small.setLogLevel(63);
137          // A variation is to just do N iterations
138          //if (iPass)
139          //small.setMaximumIterations(100);
140          // Solve
141          small.factorization()->messageLevel(8);
142          if (iPass) {
143               small.dual();
144          } else {
145               small.writeMps("continf.mps");
146               ClpSolve solveOptions;
147               solveOptions.setSolveType(ClpSolve::useDual);
148               //solveOptions.setSolveType(ClpSolve::usePrimalorSprint);
149               //solveOptions.setSpecialOption(1,2,200); // idiot
150               small.initialSolve(solveOptions);
151          }
152          bool dualInfeasible = (small.status() == 2);
153          // move solution back
154          double * solution = model2->primalColumnSolution();
155          const double * smallSolution = small.primalColumnSolution();
156          for (int j = 0; j < numberSmallColumns; j++) {
157               iColumn = whichColumns[j];
158               solution[iColumn] = smallSolution[j];
159               model2->setColumnStatus(iColumn, small.getColumnStatus(j));
160          }
161          for (iRow = 0; iRow < numberSort; iRow++) {
162               int kRow = sort[iRow];
163               model2->setRowStatus(kRow, small.getRowStatus(iRow));
164          }
165          // compute full solution
166          memset(fullSolution, 0, originalNumberRows * sizeof(double));
167          model2->times(1.0, model2->primalColumnSolution(), fullSolution);
168          if (iPass != maxPass - 1) {
169               // Mark row as not looked at
170               for (iRow = 0; iRow < originalNumberRows; iRow++)
171                    weight[iRow] = 1.123e50;
172               // Look at rows already in small problem
173               int iSort;
174               int numberDropped = 0;
175               int numberKept = 0;
176               int numberBinding = 0;
177               int numberInfeasibilities = 0;
178               double sumInfeasibilities = 0.0;
179               for (iSort = 0; iSort < numberSort; iSort++) {
180                    iRow = sort[iSort];
181                    //printf("%d %g %g\n",iRow,fullSolution[iRow],small.primalRowSolution()[iSort]);
182                    if (model2->getRowStatus(iRow) == ClpSimplex::basic) {
183                         // Basic - we can get rid of if early on
184                         if (iPass < takeOutPass && !dualInfeasible) {
185                              // may have hit max iterations so check
186                              double infeasibility = CoinMax(fullSolution[iRow] - rowUpper[iRow],
187                                                             rowLower[iRow] - fullSolution[iRow]);
188                              weight[iRow] = -infeasibility;
189                              if (infeasibility > 1.0e-8) {
190                                   numberInfeasibilities++;
191                                   sumInfeasibilities += infeasibility;
192                              } else {
193                                   weight[iRow] = 1.0;
194                                   numberDropped++;
195                              }
196                         } else {
197                              // keep
198                              weight[iRow] = -1.0e40;
199                              numberKept++;
200                         }
201                    } else {
202                         // keep
203                         weight[iRow] = -1.0e50;
204                         numberKept++;
205                         numberBinding++;
206                    }
207               }
208               // Now rest
209               for (iRow = 0; iRow < originalNumberRows; iRow++) {
210                    sort[iRow] = iRow;
211                    if (weight[iRow] == 1.123e50) {
212                         // not looked at yet
213                         double infeasibility = CoinMax(fullSolution[iRow] - rowUpper[iRow],
214                                                        rowLower[iRow] - fullSolution[iRow]);
215                         weight[iRow] = -infeasibility;
216                         if (infeasibility > 1.0e-8) {
217                              numberInfeasibilities++;
218                              sumInfeasibilities += infeasibility;
219                         }
220                    }
221               }
222               // sort
223               CoinSort_2(weight, weight + originalNumberRows, sort);
224               numberSort = CoinMin(originalNumberRows, smallNumberRows + numberKept);
225               memset(take, 0, originalNumberRows);
226               for (iRow = 0; iRow < numberSort; iRow++)
227                    take[sort[iRow]] = 1;
228               numberSmallColumns = 0;
229               for (iColumn = 0; iColumn < numberColumns; iColumn++) {
230                    int n = 0;
231                    for (int j = start[iColumn]; j < start[iColumn] + length[iColumn]; j++) {
232                         int iRow = row[j];
233                         if (take[iRow])
234                              n++;
235                    }
236                    if (n)
237                         whichColumns[numberSmallColumns++] = iColumn;
238               }
239               printf("%d rows binding, %d rows kept, %d rows dropped - new size %d rows, %d columns\n",
240                      numberBinding, numberKept, numberDropped, numberSort, numberSmallColumns);
241               printf("%d rows are infeasible - sum is %g\n", numberInfeasibilities,
242                      sumInfeasibilities);
243               if (!numberInfeasibilities) {
244                    printf("Exiting as looks optimal\n");
245                    break;
246               }
247               numberInfeasibilities = 0;
248               sumInfeasibilities = 0.0;
249               for (iSort = 0; iSort < numberSort; iSort++) {
250                    if (weight[iSort] > -1.0e30 && weight[iSort] < -1.0e-8) {
251                         numberInfeasibilities++;
252                         sumInfeasibilities += -weight[iSort];
253                    }
254               }
255               printf("in small model %d rows are infeasible - sum is %g\n", numberInfeasibilities,
256                      sumInfeasibilities);
257          }
258     }
259     delete [] weight;
260     delete [] sort;
261     delete [] whichColumns;
262     delete [] take;
263     // If problem is big you may wish to skip this
264     model2->dual();
265     int numberBinding = 0;
266     for (iRow = 0; iRow < originalNumberRows; iRow++) {
267          if (model2->getRowStatus(iRow) != ClpSimplex::basic)
268               numberBinding++;
269     }
270     printf("%d binding rows at end\n", numberBinding);
271     pinfo.postsolve(true);
272
273     int numberIterations = model2->numberIterations();;
274     delete model2;
275     /* After this postsolve model should be optimal.
276        We can use checkSolution and test feasibility */
277     model.checkSolution();
278     if (model.numberDualInfeasibilities() ||
279               model.numberPrimalInfeasibilities())
280          printf("%g dual %g(%d) Primal %g(%d)\n",
281                 model.objectiveValue(),
282                 model.sumDualInfeasibilities(),
283                 model.numberDualInfeasibilities(),
284                 model.sumPrimalInfeasibilities(),
285                 model.numberPrimalInfeasibilities());
286     // But resolve for safety
287     model.primal(1);
288
289     numberIterations += model.numberIterations();;
290     printf("Solve took %g seconds\n", CoinCpuTime() - time1);
291     return 0;
292}
Note: See TracBrowser for help on using the repository browser.