source: trunk/Clp/examples/dualCuts.cpp @ 778

Last change on this file since 778 was 673, checked in by forrest, 14 years ago

mods

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