source: trunk/Clp/examples/dualCuts.cpp

Last change on this file was 2278, checked in by forrest, 9 months ago

COIN_BIG_INDEX 2 changes

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