source: trunk/Cbc/examples/CbcSolver2.cpp

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

Change to EPL license notice.

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