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

Last change on this file since 1662 was 1662, checked in by lou, 8 years ago

Add EPL license notice in examples.

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