source: branches/sandbox/Cbc/examples/CbcSolver2.cpp @ 1273

Last change on this file since 1273 was 705, checked in by forrest, 12 years ago

INT_MAX to COIN_INT_MAX in Cbc/examples

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 13.4 KB
Line 
1// Copyright (C) 2005, International Business Machines
2// Corporation and others.  All Rights Reserved.
3
4#include <cassert>
5
6#include "CoinTime.hpp"
7
8#include "CoinHelperFunctions.hpp"
9#include "CoinIndexedVector.hpp"
10#include "ClpFactorization.hpp"
11#include "ClpObjective.hpp"
12#include "ClpSimplex.hpp"
13#include "ClpSolve.hpp"
14#include "CbcSolver2.hpp"
15#include "CbcModel.hpp"
16#include "CoinModel.hpp"
17
18static int timesBad_=0;
19static int iterationsBad_=0;
20//#############################################################################
21// Solve methods
22//#############################################################################
23void CbcSolver2::initialSolve()
24{
25  modelPtr_->scaling(0);
26  setBasis(basis_,modelPtr_);
27  // Do long thin by sprint
28  ClpSolve options;
29  options.setSolveType(ClpSolve::usePrimalorSprint);
30  options.setPresolveType(ClpSolve::presolveOff);
31  options.setSpecialOption(1,3,30);
32  modelPtr_->initialSolve(options);
33  basis_ = getBasis(modelPtr_);
34  modelPtr_->setLogLevel(0);
35}
36
37//-----------------------------------------------------------------------------
38void CbcSolver2::resolve()
39{
40  int numberColumns = modelPtr_->numberColumns();
41  if ((count_<10&&algorithm_==2)||!algorithm_) {
42    OsiClpSolverInterface::resolve();
43    if (modelPtr_->status()==0) {
44      count_++;
45      double * solution = modelPtr_->primalColumnSolution();
46      int i;
47      for (i=0;i<numberColumns;i++) {
48        if (solution[i]>1.0e-6||modelPtr_->getStatus(i)==ClpSimplex::basic) {
49          node_[i]=CoinMax(count_,node_[i]);
50          howMany_[i]++;
51        }
52      }
53    } else {
54      if (!algorithm_==2)
55        printf("infeasible early on\n");
56    }
57  } else {
58    // use counts
59    int numberRows=modelPtr_->numberRows();
60    int * whichRow = new int[numberRows];
61    int * whichColumn = new int [numberColumns];
62    int i;
63    const double * lower = modelPtr_->columnLower();
64    const double * upper = modelPtr_->columnUpper();
65    const double * rowUpper = modelPtr_->rowUpper();
66    bool equality=false;
67    for (i=0;i<numberRows;i++) {
68      if (rowUpper[i]==1.0) {
69        equality=true;
70        break;
71      }
72    }
73    setBasis(basis_,modelPtr_);
74    int nNewCol=0;
75    // Column copy
76    //const double * element = modelPtr_->matrix()->getElements();
77    const int * row = modelPtr_->matrix()->getIndices();
78    const CoinBigIndex * columnStart = modelPtr_->matrix()->getVectorStarts();
79    const int * columnLength = modelPtr_->matrix()->getVectorLengths();
80   
81    int * rowActivity = new int[numberRows];
82    memset(rowActivity,0,numberRows*sizeof(int));
83    int * rowActivity2 = new int[numberRows];
84    memset(rowActivity2,0,numberRows*sizeof(int));
85    char * mark = new char[numberColumns];
86    memset(mark,0,numberColumns);
87    // Get rows which are satisfied
88    for (i=0;i<numberColumns;i++) {
89      if (lower[i]>0.0) {
90        CoinBigIndex j;
91        for (j=columnStart[i];
92             j<columnStart[i]+columnLength[i];j++) {
93          int iRow=row[j];
94          rowActivity2[iRow] ++;
95        }
96      } else if (!upper[i]) {
97        mark[i]=2; // no good
98      }
99    }
100    // If equality - check not infeasible
101    if (equality) {
102      bool feasible=true;
103      for (i=0;i<numberRows;i++) {
104        if (rowActivity2[i]>1) {
105          feasible=false;
106          break;
107        }
108      }
109      if (!feasible) {
110        delete [] rowActivity;
111        delete [] rowActivity2;
112        modelPtr_->setProblemStatus(1);
113        delete [] whichRow;
114        delete [] whichColumn;
115        delete [] mark;
116        printf("infeasible by inspection (over)\n");
117        return;
118      }
119    }
120    int nNoGood=0;
121    for (i=0;i<numberColumns;i++) {
122      if (mark[i]==2) {
123        nNoGood++;
124        continue;
125      }
126      bool choose;
127      if (algorithm_==1)
128        choose = true;
129      else
130        choose = (node_[i]>count_-memory_&&node_[i]>0);
131      bool any;
132      if (equality) {
133        // See if forced to be zero
134        CoinBigIndex j;
135        any=true;
136        for (j=columnStart[i];
137             j<columnStart[i]+columnLength[i];j++) {
138          int iRow=row[j];
139          if (rowActivity2[iRow])
140            any=false; // can't be in
141        }
142      } else {
143        // See if not useful
144        CoinBigIndex j;
145        any=false;
146        for (j=columnStart[i];
147             j<columnStart[i]+columnLength[i];j++) {
148          int iRow=row[j];
149          if (!rowActivity2[iRow])
150            any=true; // useful
151        }
152      }
153      if (!any&&!lower[i]) {
154        choose=false;
155        // and say can't be useful
156        mark[i]=2;
157        nNoGood++;
158      }
159      if (strategy_&&modelPtr_->getColumnStatus(i)==ClpSimplex::basic)
160        choose=true;
161      if (choose||lower[i]>0.0) {
162        mark[i]=1;
163        whichColumn[nNewCol++]=i;
164        CoinBigIndex j;
165        double value = upper[i];
166        if (value) {
167          for (j=columnStart[i];
168               j<columnStart[i]+columnLength[i];j++) {
169            int iRow=row[j];
170            rowActivity[iRow] ++;
171          }
172        }
173      }
174    }
175    // If equality add in slacks
176    CoinModel build;
177    if (equality) {
178      int row=0;
179      for (i=0;i<numberRows;i++) {
180        // put in all rows if wanted
181        if(strategy_)
182          rowActivity2[i]=0;
183        if (!rowActivity2[i]) {
184          double element=1.0;
185          build.addColumn(1,&row,&element,0.0,1.0,1.0e8); // large cost
186          row++;
187        }
188      }
189    }
190    int nOK=0;
191    int nNewRow=0;
192    for (i=0;i<numberRows;i++) {
193      if (rowActivity[i])
194        nOK++;
195      if (!rowActivity2[i])
196        whichRow[nNewRow++]=i; // not satisfied
197      else
198        modelPtr_->setRowStatus(i,ClpSimplex::basic); // make slack basic
199    }
200    if (nOK<numberRows) {
201      for (i=0;i<numberColumns;i++) {
202        if (!mark[i]) {
203          CoinBigIndex j;
204          int good=0;
205          for (j=columnStart[i];
206               j<columnStart[i]+columnLength[i];j++) {
207            int iRow=row[j];
208            if (!rowActivity[iRow]) {
209              rowActivity[iRow] ++;
210              good++;
211            }
212          }
213          if (good) {
214            nOK+=good;
215            whichColumn[nNewCol++]=i;
216          }
217        }
218      }
219    }
220    delete [] rowActivity;
221    delete [] rowActivity2;
222    if (nOK<numberRows) {
223      modelPtr_->setProblemStatus(1);
224      delete [] whichRow;
225      delete [] whichColumn;
226      delete [] mark;
227      printf("infeasible by inspection\n");
228      return;
229    }
230    bool allIn=false;
231    if (nNewCol+nNoGood+numberRows>numberColumns) {
232      // add in all
233      allIn=true;
234      for (i=0;i<numberColumns;i++) {
235        if (!mark[i]) {
236          whichColumn[nNewCol++]=i;
237        }
238      }
239    }
240    ClpSimplex *  temp = new ClpSimplex(modelPtr_,nNewRow,whichRow,nNewCol,whichColumn);
241    if (equality)
242      temp->addColumns(build);
243    temp->setLogLevel(1);
244    printf("small has %d rows and %d columns (%d impossible to help) %s\n",
245           nNewRow,nNewCol,nNoGood,allIn ? "all in" : "");
246    temp->setSpecialOptions(128+512);
247    temp->setDualObjectiveLimit(1.0e50);
248    temp->dual();
249    assert (!temp->status());
250    double * solution = modelPtr_->primalColumnSolution();
251    const double * solution2 = temp->primalColumnSolution();
252    memset(solution,0,numberColumns*sizeof(double));
253    for (i=0;i<nNewCol;i++) {
254      int iColumn = whichColumn[i];
255      solution[iColumn]=solution2[i];
256      modelPtr_->setStatus(iColumn,temp->getStatus(i));
257    }
258    double * rowSolution = modelPtr_->primalRowSolution();
259    const double * rowSolution2 = temp->primalRowSolution();
260    double * dual = modelPtr_->dualRowSolution();
261    const double * dual2 = temp->dualRowSolution();
262    memset(dual,0,numberRows*sizeof(double));
263    for (i=0;i<nNewRow;i++) {
264      int iRow=whichRow[i];
265      modelPtr_->setRowStatus(iRow,temp->getRowStatus(i));
266      rowSolution[iRow]=rowSolution2[i];
267      dual[iRow]=dual2[i];
268    }
269    // See if optimal
270    double * dj = modelPtr_->dualColumnSolution();
271    // get reduced cost for large problem
272    // this assumes minimization
273    memcpy(dj,modelPtr_->objective(),numberColumns*sizeof(double));
274    modelPtr_->transposeTimes(-1.0,dual,dj);
275    modelPtr_->setObjectiveValue(temp->objectiveValue());
276    modelPtr_->setProblemStatus(0);
277    int nBad=0;
278    for (i=0;i<numberColumns;i++) {
279      if (modelPtr_->getStatus(i)==ClpSimplex::atLowerBound
280          &&upper[i]>lower[i]&&dj[i]<-1.0e-5)
281        nBad++;
282    }
283    //modelPtr_->writeMps("bada.mps");
284    //temp->writeMps("badb.mps");
285    delete temp;
286    if (nBad&&!allIn) {
287      assert (algorithm_==2);
288      //printf("%d bad\n",nBad);
289      timesBad_++;
290      // just non mark==2
291      int nAdded=0;
292      for (i=0;i<numberColumns;i++) {
293        if (!mark[i]) {
294          whichColumn[nNewCol++]=i;
295          nAdded++;
296        }
297      }
298      assert (nAdded);
299      {
300        temp = new ClpSimplex(modelPtr_,nNewRow,whichRow,nNewCol,whichColumn);
301        if (equality)
302          temp->addColumns(build);
303        temp->setLogLevel(2);
304        temp->setSpecialOptions(128+512);
305        temp->setDualObjectiveLimit(1.0e50);
306        temp->primal(1);
307        assert (!temp->status());
308        double * solution = modelPtr_->primalColumnSolution();
309        const double * solution2 = temp->primalColumnSolution();
310        memset(solution,0,numberColumns*sizeof(double));
311        for (i=0;i<nNewCol;i++) {
312          int iColumn = whichColumn[i];
313          solution[iColumn]=solution2[i];
314          modelPtr_->setStatus(iColumn,temp->getStatus(i));
315        }
316        double * rowSolution = modelPtr_->primalRowSolution();
317        const double * rowSolution2 = temp->primalRowSolution();
318        double * dual = modelPtr_->dualRowSolution();
319        const double * dual2 = temp->dualRowSolution();
320        memset(dual,0,numberRows*sizeof(double));
321        for (i=0;i<nNewRow;i++) {
322          int iRow=whichRow[i];
323          modelPtr_->setRowStatus(iRow,temp->getRowStatus(i));
324          rowSolution[iRow]=rowSolution2[i];
325          dual[iRow]=dual2[i];
326        }
327        modelPtr_->setObjectiveValue(temp->objectiveValue());
328        modelPtr_->setProblemStatus(0);
329        iterationsBad_ += temp->numberIterations();
330        printf("clean %d\n",temp->numberIterations());
331        delete temp;
332      }
333    }
334    delete [] mark;
335    delete [] whichRow;
336    delete [] whichColumn;
337    basis_ = getBasis(modelPtr_);
338    modelPtr_->setSpecialOptions(0);
339    count_++;
340    if ((count_%100)==0&&algorithm_==2)
341      printf("count %d, bad %d - iterations %d\n",count_,timesBad_,iterationsBad_);
342    for (i=0;i<numberColumns;i++) {
343      if (solution[i]>1.0e-6||modelPtr_->getStatus(i)==ClpSimplex::basic) {
344        node_[i]=CoinMax(count_,node_[i]);
345        howMany_[i]++;
346      }
347    }
348    if (modelPtr_->objectiveValue()>=modelPtr_->dualObjectiveLimit())
349      modelPtr_->setProblemStatus(1);
350  }
351}
352
353//#############################################################################
354// Constructors, destructors clone and assignment
355//#############################################################################
356
357//-------------------------------------------------------------------
358// Default Constructor
359//-------------------------------------------------------------------
360CbcSolver2::CbcSolver2 ()
361  : OsiClpSolverInterface()
362{
363  node_=NULL;
364  howMany_=NULL;
365  count_=0;
366  model_ = NULL;
367  memory_=300;
368  algorithm_=0;
369  strategy_=0;
370}
371
372//-------------------------------------------------------------------
373// Clone
374//-------------------------------------------------------------------
375OsiSolverInterface * 
376CbcSolver2::clone(bool CopyData) const
377{
378  if (CopyData) {
379    return new CbcSolver2(*this);
380  } else {
381    printf("warning CbcSolveUser clone with copyData false\n");
382    return new CbcSolver2();
383  }
384}
385
386
387//-------------------------------------------------------------------
388// Copy constructor
389//-------------------------------------------------------------------
390CbcSolver2::CbcSolver2 (
391                  const CbcSolver2 & rhs)
392  : OsiClpSolverInterface(rhs)
393{
394  model_ = rhs.model_;
395  int numberColumns = modelPtr_->numberColumns();
396  node_=CoinCopyOfArray(rhs.node_,numberColumns);
397  howMany_=CoinCopyOfArray(rhs.howMany_,numberColumns);
398  count_=rhs.count_;
399  memory_=rhs.memory_;
400  algorithm_=rhs.algorithm_;
401  strategy_=rhs.strategy_;
402}
403
404//-------------------------------------------------------------------
405// Destructor
406//-------------------------------------------------------------------
407CbcSolver2::~CbcSolver2 ()
408{
409  delete [] node_;
410  delete [] howMany_;
411}
412
413//-------------------------------------------------------------------
414// Assignment operator
415//-------------------------------------------------------------------
416CbcSolver2 &
417CbcSolver2::operator=(const CbcSolver2& rhs)
418{
419  if (this != &rhs) { 
420    OsiClpSolverInterface::operator=(rhs);
421    delete [] node_;
422    delete [] howMany_;
423    model_ = rhs.model_;
424    int numberColumns = modelPtr_->numberColumns();
425    node_=CoinCopyOfArray(rhs.node_,numberColumns);
426    howMany_=CoinCopyOfArray(rhs.howMany_,numberColumns);
427    count_=rhs.count_;
428    memory_=rhs.memory_;
429    algorithm_=rhs.algorithm_;
430    strategy_=rhs.strategy_;
431  }
432  return *this;
433}
434//-------------------------------------------------------------------
435// Real initializer
436//-------------------------------------------------------------------
437void
438CbcSolver2::initialize (CbcModel * model, const char * keep)
439{
440  model_=model;
441  int numberColumns = modelPtr_->numberColumns();
442  if (numberColumns) {
443    node_ = new int[numberColumns];
444    howMany_ = new int[numberColumns];
445    for (int i=0;i<numberColumns;i++) {
446      if (keep&&keep[i])
447        node_[i]=COIN_INT_MAX;
448      else
449        node_[i]=0;
450      howMany_[i]=0;
451    }
452  } else {
453    node_=NULL;
454    howMany_=NULL;
455  }
456}
Note: See TracBrowser for help on using the repository browser.