source: trunk/Clp/examples/iis.cpp

Last change on this file was 2317, checked in by forrest, 4 months ago

try at iis

File size: 9.4 KB
Line 
1/* $Id: iis.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/* Find an irreducible infeasible subsytem.
14   This is not trying to be fast - this implementation is
15   really trying to find the maximum feasible subsystem
16*/
17#include "ClpSimplex.hpp"
18#include "CoinSort.hpp"
19#include "CoinHelperFunctions.hpp"
20#include "CoinTime.hpp"
21#include <iomanip>
22
23int main(int argc, const char *argv[])
24{
25  ClpSimplex  model;
26  int status;
27  if (argc < 2) {
28    printf("please give a model\n");
29    status=1;
30  } else {
31    status = model.readMps(argv[1], true);
32  }
33  if (status)
34    exit(10);
35  // check infeasible
36  model.initialSolve();
37
38  if (model.problemStatus()==0) {
39    printf("Problem is feasible - nothing to do\n");
40    exit(11);
41  }
42  double * rowLower = model.rowLower();
43  double * rowUpper = model.rowUpper();
44  int numberColumns = model.numberColumns();
45  int numberRows = model.numberRows();
46  if (model.numberPrimalInfeasibilities()==1) {
47    const double * rowActivity = model.primalRowSolution();
48    for (int iRow = 0; iRow < numberRows; iRow++) {
49      if (rowActivity[iRow]<rowLower[iRow]-1.0e-5||
50          rowActivity[iRow]>rowUpper[iRow]+1.0e-5) {
51        printf("The following one row is the IIS\n");
52        printf("%d %s\n",iRow,model.getRowName(iRow).c_str());
53        return 0;
54      }
55    }
56  }
57  // This example only looks at constraints
58  // add infeasibility slacks
59  ClpSimplex model2(model);
60  double time1 = CoinCpuTime();
61
62  int nAdded=0;
63
64  CoinBigIndex * starts = new CoinBigIndex [2*numberRows+1];
65  int * rows = new int [2*numberRows+1];
66  double * elements = new double [2*numberRows];
67  double * newObjective = new double [2*numberRows];
68  memset(model2.objective(),0,numberColumns*sizeof(double));
69  for (int iRow = 0; iRow < numberRows; iRow++) {
70    if (rowLower[iRow]>-1.0e30) {
71      rows[nAdded]=iRow;
72      elements[nAdded++]=1.0;
73    }
74    if (rowUpper[iRow]<1.0e30) {
75      rows[nAdded]=iRow;
76      elements[nAdded++]=-1.0;
77    }
78  }
79  rows[nAdded]=-1; // so can test if two slacks for a row
80  starts[0]=0;
81#define OBJ 10000.0
82  for (int i=0;i<nAdded;i++) {
83    starts[i+1]=i+1;
84    newObjective[i]=OBJ;
85  }
86  model2.addColumns(nAdded,NULL,NULL,newObjective,starts,rows,elements);
87  delete [] newObjective;
88  delete [] starts;
89  delete [] elements;
90  int numberColumns2=numberColumns+nAdded;
91  rowLower = model2.rowLower();
92  rowUpper = model2.rowUpper();
93  // solve
94  model2.allSlackBasis();
95  model2.dual();
96  printf("Initial sum of infeasibilities is %g\n",model2.objectiveValue()/OBJ);
97  if (model2.numberPrimalInfeasibilities()) {
98    printf("ouch\n");
99    model2.writeMps("bad.mps");
100    abort();
101  }
102  model2.setLogLevel(0);
103  model2.setInfeasibilityCost(1.0e15);
104  //#define SWITCH_OFF_SCALING
105#ifdef SWITCH_OFF_SCALING
106  model2.scaling(0);
107#endif
108  int * coverRows = new int [numberRows];
109  int * candidateRows = new int [numberRows];
110  int * nextRows = new int [numberRows];
111  int nCover=0;
112  int nCandidate=0;
113  int outRow=9999;
114  const double * duals = model2.dualRowSolution();
115  for (int iRow=0;iRow<numberRows;iRow++) {
116    if (fabs(duals[iRow])>1.0e-5)
117      nextRows[nCandidate++]=iRow;
118  }
119  assert (nCandidate);
120  /* save basis info -
121     we could reduce problem size each time but
122     normally not many passes needed */
123  unsigned char * statusArray = new unsigned char [numberColumns2+numberRows];
124  double * solution = new double[numberColumns2];
125  memcpy(statusArray,model2.statusArray(),numberColumns2+numberRows);
126  memcpy(solution,model2.primalColumnSolution(),numberColumns2*sizeof(double));
127  double lastObjectiveValue = model2.objectiveValue();
128  double * objective = model2.objective();
129  int nSolves=0;
130  int nPass=0;
131  while (outRow>=0) {
132    memcpy(candidateRows,nextRows,nCandidate*sizeof(int));
133    double sumInf=COIN_DBL_MAX;
134    int nInf=numberColumns2;
135    bool badAccuracy=false;
136    outRow=-1;
137    nPass++;
138    for (int j=0;j<nCandidate;j++) {
139      memcpy(model2.statusArray(),statusArray,numberColumns2+numberRows);
140      memcpy(model2.primalColumnSolution(),solution,numberColumns2*sizeof(double));
141      int iRow=candidateRows[j];
142      /* One can free rows or zero out costs -
143         zeroing out costs is cleaner from basis */
144#define ZAP_COSTS
145#ifdef ZAP_COSTS
146      // zero cost and solve
147      int iColumn0=-1;
148      int iColumn1=-1;
149      for (int i=0;i<nAdded;i++) {
150        if (rows[i]==iRow) {
151          iColumn0=i+numberColumns;
152          objective[iColumn0]=0.0;
153          if (rows[i+1]==iRow) {
154            iColumn1=i+1+numberColumns;
155            objective[iColumn1]=0.0;
156          }
157          break;
158        }
159      }
160      model2.primal(1);
161      if (model2.objectiveValue()>lastObjectiveValue+1.0e-1) {
162        if (!badAccuracy)
163          printf("Sum infeasibilities increased from %g to %g! - problems with accuracy\n",
164                 lastObjectiveValue/OBJ,model2.objectiveValue()/OBJ);
165        badAccuracy=true;
166      }
167      nSolves++;
168      objective[iColumn0]=OBJ;
169      if (iColumn1>=0)
170        objective[iColumn1]=OBJ;
171#else
172      // delete and solve
173      double lower=rowLower[iRow];
174      double upper=rowUpper[iRow];
175      rowLower[iRow]=-COIN_DBL_MAX;
176      rowUpper[iRow]=COIN_DBL_MAX;
177      model2.primal(1);
178      if (model2.objectiveValue()>lastObjectiveValue+1.0e-1) {
179        if (!badAccuracy)
180          printf("Sum infeasibilities increased from %g to %g! - problems with accuracy\n",
181                 lastObjectiveValue/OBJ,model2.objectiveValue()/OBJ);
182        badAccuracy=true;
183      }
184      nSolves++;
185      rowLower[iRow]=lower;
186      rowUpper[iRow]=upper;
187#endif
188      double objectiveValue = model2.objectiveValue();
189      if (objectiveValue<sumInf+1.0e-3) {
190        int n=0;
191        for (int i=numberColumns;i<numberColumns2;i++) {
192          if (solution[i]>1.0e-5)
193            n++;
194        }
195        if (objectiveValue>1.0e-6) {
196          if (objectiveValue<sumInf-1.0e-3||n<nInf) {
197            sumInf=objectiveValue;
198            outRow=iRow;
199            nInf=n;
200          }
201        } else {
202          coverRows[nCover++]=iRow;
203          printf("Pass %d (%d solves, %.2f seconds) - candidate %d (%s) added - Final as feasible\n",nPass,nSolves,CoinCpuTime()-time1,iRow,
204                 model2.getRowName(iRow).c_str());
205          outRow=-1;
206          break; // finished
207        }
208      }
209    }
210    if (outRow>=0) {
211#ifdef ZAP_COSTS
212      int iColumn0=-1;
213      int iColumn1=-1;
214      for (int i=0;i<nAdded;i++) {
215        if (rows[i]==outRow) {
216          iColumn0=i+numberColumns;
217          objective[iColumn0]=0.0;
218          if (rows[i+1]==outRow) {
219            iColumn1=i+1+numberColumns;
220            objective[iColumn1]=0.0;
221          }
222          break;
223        }
224      }
225#else
226      rowLower[outRow]=-COIN_DBL_MAX;
227      rowUpper[outRow]=COIN_DBL_MAX;
228#endif
229      coverRows[nCover++]=outRow;
230      memcpy(model2.statusArray(),statusArray,numberColumns2+numberRows);
231      memcpy(model2.primalColumnSolution(),solution,numberColumns2*sizeof(double));
232      model2.primal(1);
233      if (model2.objectiveValue()>lastObjectiveValue+1.0e-1) {
234        if (!badAccuracy)
235          printf("Sum infeasibilities increased from %g to %g on cover solve! - problems with accuracy\n",
236                 lastObjectiveValue/OBJ,model2.objectiveValue()/OBJ);
237        badAccuracy=true;
238      }
239      lastObjectiveValue=model2.objectiveValue();
240      memcpy(statusArray,model2.statusArray(),numberColumns2+numberRows);
241      memcpy(solution,model2.primalColumnSolution(),numberColumns2*sizeof(double));
242      nCandidate=0;
243      for (int iRow=0;iRow<numberRows;iRow++) {
244        if (fabs(duals[iRow])>1.0e-6)
245          nextRows[nCandidate++]=iRow;
246      }
247      printf("Pass %d (%d solves, %.2f seconds) - candidate %d (%s) added - infeasibility %g (%d)\n",
248             nPass,nSolves,CoinCpuTime()-time1,outRow,
249             model2.getRowName(outRow).c_str(),
250             model2.objectiveValue()/OBJ,nInf);
251    }
252  }
253  model2=model;
254  model2.deleteRows(nCover,coverRows);
255  // make sure not unbounded
256  memset(model2.objective(),0,numberColumns*sizeof(double));
257  model2.dual();
258  printf("The following %d rows cover the IIS\n",nCover);
259  for (int i=0;i<nCover;i++) {
260    int iRow=coverRows[i];
261    printf("%d %s\n",iRow,model.getRowName(iRow).c_str());
262  }
263  if (model2.problemStatus()) {
264    printf("We seem to have an accuracy problem??\n");
265  } else {
266    // see if we can do better
267    CoinPackedMatrix * matrix = model.matrix();
268    // get row copy
269    CoinPackedMatrix rowCopy = *matrix;
270    rowCopy.reverseOrdering();
271    const int * column = rowCopy.getIndices();
272    const int * rowLength = rowCopy.getVectorLengths();
273    const CoinBigIndex * rowStart = rowCopy.getVectorStarts();
274    const double * element = rowCopy.getElements();
275    numberRows=model2.numberRows();
276    memcpy(statusArray,model2.statusArray(),numberColumns+numberRows);
277    memcpy(solution,model2.primalColumnSolution(),numberColumns*sizeof(double));
278    int lastRow=numberRows-1;
279    for (int i=0;i<nCover;i++) {
280      memcpy(model2.statusArray(),statusArray,numberColumns+numberRows);
281      memcpy(model2.primalColumnSolution(),solution,numberColumns*sizeof(double));
282      int iRow=coverRows[i];
283      model2.addRow(rowLength[iRow],
284                    column+rowStart[iRow],element+rowStart[iRow],
285                    model.rowLower()[iRow],model.rowUpper()[iRow]);
286      model2.dual();
287      if (!model2.problemStatus()) {
288        printf("%d %s should not be in cover\n",iRow,model.getRowName(iRow).c_str());
289      }
290      model2.deleteRows(1,&lastRow);
291    }
292  }
293  delete [] rows;
294  delete [] statusArray;
295  delete [] solution;
296  delete [] candidateRows;
297  delete [] nextRows;
298  delete [] coverRows;
299  return 0;
300}
Note: See TracBrowser for help on using the repository browser.