source: trunk/Clp/examples/sprint2.cpp @ 1597

Last change on this file since 1597 was 1597, checked in by forrest, 11 years ago

fix stupidity in example

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 6.2 KB
Line 
1/* $Id: sprint2.cpp 1597 2010-08-31 11:16:24Z forrest $ */
2// Copyright (C) 2003, International Business Machines
3// Corporation and others.  All Rights Reserved.
4
5#include "ClpSimplex.hpp"
6#include "ClpPresolve.hpp"
7#include "CoinSort.hpp"
8#include <iomanip>
9
10int main(int argc, const char *argv[])
11{
12     ClpSimplex  model;
13     int status;
14     // Keep names
15     if (argc < 2) {
16          status = model.readMps("small.mps", true);
17     } else {
18          status = model.readMps(argv[1], true);
19     }
20     if (status)
21          exit(10);
22     /*
23       This driver implements the presolve variation of Sprint.
24       This assumes we can get feasible easily
25     */
26
27     int numberRows = model.numberRows();
28     int numberColumns = model.numberColumns();
29
30     // We will need arrays to choose variables.  These are too big but ..
31     double * weight = new double [numberRows+numberColumns];
32     int * sort = new int [numberRows+numberColumns];
33
34     double * columnLower = model.columnLower();
35     double * columnUpper = model.columnUpper();
36     double * saveLower = new double [numberColumns];
37     memcpy(saveLower, columnLower, numberColumns * sizeof(double));
38     double * saveUpper = new double [numberColumns];
39     memcpy(saveUpper, columnUpper, numberColumns * sizeof(double));
40     // Fix in some magical way so remaining problem is easy
41#if 0
42     // This is from a real-world problem
43     for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
44          char firstCharacter = model.columnName(iColumn)[0];
45          if (firstCharacter == 'F' || firstCharacter == 'P'
46                    || firstCharacter == 'L' || firstCharacter == 'T') {
47               columnUpper[iColumn] = columnLower[iColumn];
48          }
49     }
50#else
51     double * obj = model->objective();
52     double * saveObj = new double [numberColumns];
53     memcpy(saveObj, obj, numberColumns * sizeof(double));
54     memset(obj, 0, numberColumns * sizeof(double));
55     model->dual();
56     memcpy(obj, saveObj, numberColumns * sizeof(double));
57     delete [] saveObj;
58     for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
59         if (solution[iColumn]<columnLower[iColumn]+1.0e-8) {
60             columnUpper[iColumn] = columnLower[iColumn];
61       }
62     }
63#endif
64     double * solution = model.primalColumnSolution();
65
66     // Just do this number of passes
67     int maxPass = 100;
68     int iPass;
69     double lastObjective = 1.0e31;
70
71     // Just take this number of columns in small problem
72     int smallNumberColumns = 3 * numberRows;
73     // And we want number of rows to be this
74     int smallNumberRows = numberRows / 4;
75
76     for (iPass = 0; iPass < maxPass; iPass++) {
77          printf("Start of pass %d\n", iPass);
78          ClpSimplex * model2;
79          ClpPresolve pinfo;
80          int numberPasses = 1; // can change this
81          model2 = pinfo.presolvedModel(model, 1.0e-8, false, numberPasses, false);
82          if (!model2) {
83               fprintf(stdout, "ClpPresolve says %s is infeasible with tolerance of %g\n",
84                       argv[1], 1.0e-8);
85               // model was infeasible - maybe try again with looser tolerances
86               model2 = pinfo.presolvedModel(model, 1.0e-7, false, numberPasses, false);
87               if (!model2) {
88                    fprintf(stdout, "ClpPresolve says %s is infeasible with tolerance of %g\n",
89                            argv[1], 1.0e-7);
90                    exit(2);
91               }
92          }
93          // change factorization frequency from 200
94          model2->setFactorizationFrequency(100 + model2->numberRows() / 50);
95          model2->primal();
96          pinfo.postsolve(true);
97
98          // adjust smallNumberColumns if necessary
99          if (iPass) {
100               double ratio = ((double) smallNumberRows) / ((double) model2->numberRows());
101               smallNumberColumns = (int)(smallNumberColumns * ratio);
102          }
103          delete model2;
104          /* After this postsolve model should be optimal.
105             We can use checkSolution and test feasibility */
106          model.checkSolution();
107          if (model.numberDualInfeasibilities() ||
108                    model.numberPrimalInfeasibilities())
109               printf("%g dual %g(%d) Primal %g(%d)\n",
110                      model.objectiveValue(),
111                      model.sumDualInfeasibilities(),
112                      model.numberDualInfeasibilities(),
113                      model.sumPrimalInfeasibilities(),
114                      model.numberPrimalInfeasibilities());
115          // Put back true bounds
116          memcpy(columnLower, saveLower, numberColumns * sizeof(double));
117          memcpy(columnUpper, saveUpper, numberColumns * sizeof(double));
118          if ((model.objectiveValue() > lastObjective - 1.0e-7 && iPass > 5) ||
119                    iPass == maxPass - 1) {
120               break; // finished
121          } else {
122               lastObjective = model.objectiveValue();
123               // now massage weight so all basic in plus good djs
124               const double * djs = model->primalDualSolution();
125               for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
126                    double dj = djs[iColumn];
127                    double value = solution[iColumn];
128                    if (model.getStatus(iColumn) == ClpSimplex::basic)
129                         dj = -1.0e50;
130                    else if (dj < 0.0 && value < columnUpper[iColumn])
131                         dj = dj;
132                    else if (dj > 0.0 && value > columnLower[iColumn])
133                         dj = -dj;
134                    else if (columnUpper[iColumn] > columnLower[iColumn])
135                         dj = fabs(dj);
136                    else
137                         dj = 1.0e50;
138                    weight[iColumn] = dj;
139                    sort[iColumn] = iColumn;
140               }
141               // sort
142               CoinSort_2(weight, weight + numberColumns, sort);
143               // and fix others
144               for (int iColumn = smallNumberColumns; iColumn < numberColumns; iColumn++) {
145                    int kColumn = sort[iColumn];
146                    double value = solution[kColumn];
147                    columnLower[kColumn] = value;
148                    columnUpper[kColumn] = value;
149               }
150          }
151     }
152     delete [] weight;
153     delete [] sort;
154     delete [] saveLower;
155     delete [] saveUpper;
156     model.primal(1);
157     return 0;
158}
Note: See TracBrowser for help on using the repository browser.