source: trunk/Cbc/examples/CbcSolver2.cpp @ 1173

Last change on this file since 1173 was 1173, checked in by forrest, 11 years ago

add $id

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