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

Last change on this file since 1898 was 1599, checked in by forrest, 10 years ago

example is rubbish - but at least it should compile!

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 6.2 KB
Line 
1/* $Id: sprint2.cpp 1599 2010-08-31 13:36:05Z stefan $ */
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     double * solution = model.primalColumnSolution();
41     // Fix in some magical way so remaining problem is easy
42#if 0
43     // This is from a real-world problem
44     for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
45          char firstCharacter = model.columnName(iColumn)[0];
46          if (firstCharacter == 'F' || firstCharacter == 'P'
47                    || firstCharacter == 'L' || firstCharacter == 'T') {
48               columnUpper[iColumn] = columnLower[iColumn];
49          }
50     }
51#else
52     double * obj = model.objective();
53     double * saveObj = new double [numberColumns];
54     memcpy(saveObj, obj, numberColumns * sizeof(double));
55     memset(obj, 0, numberColumns * sizeof(double));
56     model.dual();
57     memcpy(obj, saveObj, numberColumns * sizeof(double));
58     delete [] saveObj;
59     for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
60         if (solution[iColumn]<columnLower[iColumn]+1.0e-8) {
61             columnUpper[iColumn] = columnLower[iColumn];
62       }
63     }
64#endif
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.dualColumnSolution();
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.