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

Last change on this file since 2496 was 2469, checked in by unxusr, 7 months ago

formatting

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 13.9 KB
Line 
1// $Id: CbcSolver2.cpp 2469 2019-01-06 23:17:46Z 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// 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 CbcSolver2::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.