source: stable/1.17/Clp/src/ClpPEDualRowSteepest.cpp @ 2439

Last change on this file since 2439 was 2385, checked in by unxusr, 11 months ago

formatting

File size: 16.9 KB
Line 
1// Corporation and others.  All Rights Reserved.
2/*
3   Authors
4
5   Jeremy Omer
6
7   Last update: april 10, 2015
8
9 */
10
11#include "CoinPragma.hpp"
12#include "ClpSimplex.hpp"
13#include "ClpPEDualRowSteepest.hpp"
14#include "CoinIndexedVector.hpp"
15#include "ClpFactorization.hpp"
16#include "CoinHelperFunctions.hpp"
17#include "ClpMessage.hpp"
18#include <cstdio>
19
20//#############################################################################
21// Constructors / Destructor / Assignment
22//#############################################################################
23
24//-------------------------------------------------------------------
25// Default Constructor
26//-------------------------------------------------------------------
27ClpPEDualRowSteepest::ClpPEDualRowSteepest(double psi, int mode)
28  : ClpDualRowSteepest(mode)
29  , modelPE_(NULL)
30  , psi_(psi)
31  , iCurrent_(0)
32  , iInterval_(100)
33  , updateCompatibles_(true)
34  , coDegenCompatibles_(0)
35  , coConsecutiveCompatibles_(0)
36{
37  iInterval_ = 100;
38}
39
40//-------------------------------------------------------------------
41// Copy constructor
42//-------------------------------------------------------------------
43ClpPEDualRowSteepest::ClpPEDualRowSteepest(const ClpPEDualRowSteepest &source)
44  : ClpDualRowSteepest(source)
45{
46  modelPE_ = NULL;
47  psi_ = source.psi_;
48  iCurrent_ = source.iCurrent_;
49  iInterval_ = source.iInterval_;
50  updateCompatibles_ = source.updateCompatibles_;
51  coDegenCompatibles_ = source.coDegenCompatibles_;
52  coConsecutiveCompatibles_ = source.coConsecutiveCompatibles_;
53}
54
55//-------------------------------------------------------------------
56// Destructor
57//-------------------------------------------------------------------
58ClpPEDualRowSteepest::~ClpPEDualRowSteepest()
59{
60  delete modelPE_;
61}
62
63//----------------------------------------------------------------
64// Assignment operator
65//-------------------------------------------------------------------
66ClpPEDualRowSteepest &
67ClpPEDualRowSteepest::operator=(const ClpPEDualRowSteepest &rhs)
68{
69  if (this != &rhs) {
70    ClpDualRowSteepest::operator=(rhs);
71    delete modelPE_;
72    modelPE_ = NULL;
73  }
74  return *this;
75}
76
77//-------------------------------------------------------------------
78// Clone
79//-------------------------------------------------------------------
80ClpDualRowPivot *ClpPEDualRowSteepest::clone(bool CopyData) const
81{
82  if (CopyData) {
83    return new ClpPEDualRowSteepest(*this);
84  } else {
85    return new ClpPEDualRowSteepest(psi_);
86  }
87}
88//extern int thinkDegenerate;
89// Returns pivot row, -1 if none
90int ClpPEDualRowSteepest::pivotRow()
91{
92  assert(model_);
93
94  // Determine whether the set of compatible variables should be updated
95  //
96  // store the number of degenerate pivots on compatible variables and the
97  // overal number of degenerate pivots
98  //double progress = fabs(modelPE_->lastObjectiveValue() - model_->objectiveValue());
99  //bool isLastDegenerate = progress <= 1.0e-12*fabs(model_->objectiveValue()) ? true:false;
100  bool isLastDegenerate = fabs(model_->theta()) < 1.0e-7;
101  bool isLastDegenerate2;
102  // could fine tune a bit e.g. use 0.0 rather than tolerance
103  if (model_->directionIn() > 0) {
104    // coming in from lower bound
105    isLastDegenerate2 = (model_->dualIn() < model_->dualTolerance());
106  } else {
107    // coming in from upper bound
108    isLastDegenerate2 = (model_->dualIn() > -model_->dualTolerance());
109  }
110#if 0
111        if ( isLastDegenerate2 &&!isLastDegenerate)
112          printf("PE2 thinks degenerate - theta %g dual %g - at %s bound\n",
113                 model_->theta(),model_->dualIn(),
114                 (fabs(model_->valueIn()-model_->lowerIn())<
115                  fabs(model_->valueIn()-model_->upperIn())) ? "lower" : "upper");
116        else if ( !isLastDegenerate2 &&isLastDegenerate)
117          printf("PE2 thinks not degenerate - theta %g dual %g - at %s bound\n",
118                 model_->theta(),model_->dualIn(),
119                 (fabs(model_->valueIn()-model_->lowerIn())<
120                  fabs(model_->valueIn()-model_->upperIn())) ? "lower" : "upper");
121        //else if ( !isLastDegenerate &&thinkDegenerate)
122        //printf("PE thinks non degenerate - dualColumn says yes\n");
123#else
124  isLastDegenerate = isLastDegenerate2;
125#endif
126  if (isLastDegenerate) {
127    modelPE_->addDegeneratePivot();
128    modelPE_->addDegeneratePivotConsecutive();
129
130    if (modelPE_->isLastPivotCompatible()) {
131      modelPE_->addDegenerateCompatiblePivot();
132    }
133  } else {
134    modelPE_->resetDegeneratePivotsConsecutive();
135  }
136
137  // the compatible variables are updated when the number of degenerate pivots
138  // on compatible variables is more than 20% the overall number of degenerate pivots
139  if (modelPE_->isLastPivotCompatible()) {
140    coConsecutiveCompatibles_++;
141    if (isLastDegenerate) {
142      coDegenCompatibles_++;
143      if (coConsecutiveCompatibles_ >= 10 && 5 * coDegenCompatibles_ * model_->numberIterations() > modelPE_->coDegeneratePivots() * coConsecutiveCompatibles_) {
144        updateCompatibles_ = true;
145      }
146    }
147  }
148
149  // For results comparison.
150  // count the number of degenerate variables if psi==1
151  // add the time spent in counting to the time limit
152  //
153  if (modelPE_->doStatistics()) {
154    modelPE_->startTimer();
155    if (psi_ >= 1 && iCurrent_ >= 100) {
156      modelPE_->updateDualDegenerates();
157      modelPE_->updateDualDegeneratesAvg(100);
158      model_->setMaximumSeconds(36000.0 + modelPE_->timeCompatibility() - CoinCpuTime());
159      iCurrent_ = 0;
160    }
161    modelPE_->stopTimer();
162  }
163
164  // Update the set of compatible variables if necessary
165  //
166  if (modelPE_->doStatistics())
167    modelPE_->startTimer();
168  double psiTmp = psi_;
169  if ((psi_ < 1.0) && (iCurrent_ >= iInterval_) && (updateCompatibles_ || iCurrent_ >= 1000 /*|| !model_->factorization()->pivots()*/)) {
170    // the compatible variables are never updated if the last pivot is non degenerate
171    if (isLastDegenerate) {
172      modelPE_->updateDualDegenerates();
173      if (modelPE_->doStatistics()) {
174        for (int i = 0; i < 4; i++)
175          assert(!model_->rowArray(i)->getNumElements());
176      }
177      modelPE_->identifyCompatibleRows(model_->rowArray(2),
178        model_->rowArray(1));
179      if (modelPE_->doStatistics() > 3) {
180        for (int i = 0; i < 4; i++)
181          assert(!model_->rowArray(i)->getNumElements());
182        char generalPrint[100];
183        sprintf(generalPrint, "updating - iCurrent,iInterval %d,%d degenerate pivots %d ? %d codegen since last %d",
184          iCurrent_, iInterval_,
185          modelPE_->coDegeneratePivots(),
186          modelPE_->coDualDegenerates(),
187          coDegenCompatibles_);
188        model_->messageHandler()->message(CLP_GENERAL,
189          *model_->messagesPointer())
190          << generalPrint << CoinMessageEol;
191        // update the average number of degenerates and compatibles for statistics
192        modelPE_->updateDualDegeneratesAvg(iCurrent_);
193        modelPE_->updateCompatibleRowsAvg(iCurrent_);
194      }
195#if PE_DEBUG >= 1
196      std::cout << "Number of compatible rows = " << modelPE_->coCompatibleRows() << " ; average  = " << modelPE_->coCompatibleRowsAvg() << std::endl;
197      std::cout << "Number of degenerates = " << modelPE_->coDualDegenerates() << " ; average  = " << modelPE_->coDualDegeneratesAvg() << std::endl;
198#endif
199
200      // the checking frequency is modified to reflect what appears to be needed
201      if (iCurrent_ == iInterval_)
202        iInterval_ = std::max(50, iInterval_ - 50);
203      else
204        iInterval_ = std::min(300, iInterval_ + 50);
205
206      // reset all the indicators that are used to decide whether the compatible
207      // variables must be updated
208      iCurrent_ = 0;
209      updateCompatibles_ = false;
210      coConsecutiveCompatibles_ = 0;
211      coDegenCompatibles_ = 0;
212    } else {
213#if 0
214                  if ((iCurrent_%500)==0)
215                        printf("not updating - ic,ii %d,%d degenerate pivots %d ? %d codegen since last %d\n",
216                               iCurrent_,iInterval_,
217                               modelPE_->coDegeneratePivots(),
218                               modelPE_->coDualDegenerates(),
219                               coDegenCompatibles_);
220#endif
221      iInterval_++;
222    }
223  }
224  // otherwise, update the value of the priority factor depending on the number of
225  // consecutive degenerate pivots
226  //
227  else {
228    // the idea is that when a lot of consecutive pivots are degenerate, there is
229    // a potentially high added value in putting a very high priroity on compatible
230    // variables
231    if (modelPE_->coDegeneratePivotsConsecutive() >= 10) {
232      //psiTmp = 0.01;
233      psiTmp = 0.25 * psi_;
234
235#if PE_DEBUG >= 1
236      std::cout << "Lower the priority factor to " << psiTmp << std::endl;
237      std::cout << "Consecutive degenerate pivots=" << modelPE_->coDegeneratePivotsConsecutive() << std::endl;
238#endif
239    }
240  }
241  iCurrent_++;
242  if (modelPE_->doStatistics())
243    modelPE_->stopTimer();
244
245  // Do the pricing and give priority to dual compatible variables
246  //
247  int i, iRow;
248  double *infeas = infeasible_->denseVector();
249  double largest = 0.0;
250  int *index = infeasible_->getIndices();
251  int number = infeasible_->getNumElements();
252  const int *pivotVariable = model_->pivotVariable();
253  int chosenRow = -1;
254  int lastPivotRow = model_->pivotRow();
255  assert(lastPivotRow < model_->numberRows());
256  double tolerance = model_->currentPrimalTolerance();
257  // we can't really trust infeasibilities if there is primal error
258  // this coding has to mimic coding in checkPrimalSolution
259  double error = CoinMin(1.0e-2, model_->largestPrimalError());
260  // allow tolerance at least slightly bigger than standard
261  tolerance = tolerance + error;
262  // But cap
263  tolerance = CoinMin(1000.0, tolerance);
264  tolerance *= tolerance; // as we are using squares
265  bool toleranceChanged = false;
266  double *solution = model_->solutionRegion();
267  double *lower = model_->lowerRegion();
268  double *upper = model_->upperRegion();
269  // do last pivot row one here
270  //#define CLP_DUAL_FIXED_COLUMN_MULTIPLIER 10.0
271  if (lastPivotRow >= 0 && lastPivotRow < model_->numberRows()) {
272#ifdef CLP_DUAL_COLUMN_MULTIPLIER
273    int numberColumns = model_->numberColumns();
274#endif
275    int iPivot = pivotVariable[lastPivotRow];
276    double value = solution[iPivot];
277    double lower = model_->lower(iPivot);
278    double upper = model_->upper(iPivot);
279    if (value > upper + tolerance) {
280      value -= upper;
281      value *= value;
282#ifdef CLP_DUAL_COLUMN_MULTIPLIER
283      if (iPivot < numberColumns)
284        value *= CLP_DUAL_COLUMN_MULTIPLIER; // bias towards columns
285#endif
286      // store square in list
287      if (infeas[lastPivotRow])
288        infeas[lastPivotRow] = value; // already there
289      else
290        infeasible_->quickAdd(lastPivotRow, value);
291    } else if (value < lower - tolerance) {
292      value -= lower;
293      value *= value;
294#ifdef CLP_DUAL_COLUMN_MULTIPLIER
295      if (iPivot < numberColumns)
296        value *= CLP_DUAL_COLUMN_MULTIPLIER; // bias towards columns
297#endif
298      // store square in list
299      if (infeas[lastPivotRow])
300        infeas[lastPivotRow] = value; // already there
301      else
302        infeasible_->add(lastPivotRow, value);
303    } else {
304      // feasible - was it infeasible - if so set tiny
305      if (infeas[lastPivotRow])
306        infeas[lastPivotRow] = COIN_INDEXED_REALLY_TINY_ELEMENT;
307    }
308    number = infeasible_->getNumElements();
309  }
310  if (model_->numberIterations() < model_->lastBadIteration() + 200) {
311    // we can't really trust infeasibilities if there is dual error
312    if (model_->largestDualError() > model_->largestPrimalError()) {
313      tolerance *= CoinMin(model_->largestDualError() / model_->largestPrimalError(), 1000.0);
314      toleranceChanged = true;
315    }
316  }
317  int numberWanted;
318  if (mode_ < 2) {
319    numberWanted = number + 1;
320  } else if (mode_ == 2) {
321    numberWanted = CoinMax(2000, number / 8);
322  } else {
323    int numberElements = model_->factorization()->numberElements();
324    double ratio = static_cast< double >(numberElements) / static_cast< double >(model_->numberRows());
325    numberWanted = CoinMax(2000, number / 8);
326    if (ratio < 1.0) {
327      numberWanted = CoinMax(2000, number / 20);
328    } else if (ratio > 10.0) {
329      ratio = number * (ratio / 80.0);
330      if (ratio > number)
331        numberWanted = number + 1;
332      else
333        numberWanted = CoinMax(2000, static_cast< int >(ratio));
334    }
335  }
336  if (model_->largestPrimalError() > 1.0e-3)
337    numberWanted = number + 1; // be safe
338  int iPass;
339  // Setup two passes
340  int start[4];
341  start[1] = number;
342  start[2] = 0;
343  double dstart = static_cast< double >(number) * model_->randomNumberGenerator()->randomDouble();
344  start[0] = static_cast< int >(dstart);
345  start[3] = start[0];
346
347  int chosenRowComp = -1;
348  double largestComp = 0.0;
349
350  // only check the compatible variables when the bidimensional factor is less than 1
351  // and the ratio of compatible variables is larger than 0.01
352  // the percentage of compatible variables is computed as the ratio to the
353  // smallest number among columns and rows
354  bool checkCompatibles = true;
355  double ratioCompatibles = static_cast< double >(modelPE_->coCompatibleRows()) / static_cast< double >(std::min(model_->numberRows(), model_->numberColumns()));
356
357  if (psi_ >= 1.0 || ratioCompatibles < 0.01)
358    checkCompatibles = false;
359
360  // potentially, a partial pricing may be done
361  for (iPass = 0; iPass < 2; iPass++) {
362    int end = start[2 * iPass + 1];
363    for (i = start[2 * iPass]; i < end; i++) {
364      iRow = index[i];
365      double value = infeas[iRow];
366      if (value > tolerance) {
367//#define OUT_EQ
368#ifdef OUT_EQ
369        {
370          int iSequence = pivotVariable[iRow];
371          if (upper[iSequence] == lower[iSequence])
372            value *= 2.0;
373        }
374#endif
375        double weight = CoinMin(weights_[iRow], 1.0e50);
376        double largestMax = std::max(psiTmp * largest, largestComp);
377        if (value > weight * largestMax) {
378          // make last pivot row last resort choice
379          if (iRow == lastPivotRow) {
380            if (value * 1.0e-10 < largestMax * weight)
381              continue;
382            else
383              value *= 1.0e-10;
384          }
385          int iSequence = pivotVariable[iRow];
386          if (model_->flagged(iSequence)) {
387            // just to make sure we don't exit before got something
388            numberWanted++;
389          } else if (checkCompatibles && modelPE_->isCompatibleRow(iRow) && value > weight * largestComp) {
390            chosenRowComp = iRow;
391            largestComp = value / weight;
392          } else if (value > largest * weight) {
393            chosenRow = iRow;
394            largest = value / weight;
395          }
396        }
397        numberWanted--;
398        if (!numberWanted)
399          break;
400      }
401    }
402    if (!numberWanted)
403      break;
404  }
405  // choose the a compatible or an incompatible row depending on the
406  // largest values and on the factor of preference psiTmp
407  if (modelPE_->doStatistics())
408    modelPE_->startTimer();
409  if (chosenRowComp >= 0 && largestComp >= psiTmp * largest) {
410
411    if (modelPE_->doStatistics()) {
412      // record the number of pivots done on compatible variables
413      // that would not have been done without positive edge
414      if (largestComp < largest)
415        modelPE_->addPriorityPivot();
416#if 0
417        if (largestComp < largest)
418          printf("pivoting on row %d largestComp %g largest %g\n",
419                 chosenRowComp,largestComp,largest);
420        else if (chosenRowComp==chosenRow)
421          printf("pivoting on row %d largestComp %g == largest %g\n",
422                 chosenRowComp,largestComp,largest);
423#endif
424    }
425    chosenRow = chosenRowComp;
426
427#if PE_DEBUG >= 1
428    modelPE_->checkCompatibilityRow(chosenRow);
429#endif
430  }
431  if (chosenRow >= 0 && psiTmp < 1.0 && modelPE_->isCompatibleRow(chosenRow)) {
432    modelPE_->isLastPivotCompatible(true);
433    modelPE_->addCompatiblePivot();
434  } else
435    modelPE_->isLastPivotCompatible(false);
436
437  if (modelPE_->doStatistics())
438    modelPE_->stopTimer();
439
440  if (chosenRow < 0 && toleranceChanged) {
441    // won't line up with checkPrimalSolution - do again
442    double saveError = model_->largestDualError();
443    model_->setLargestDualError(0.0);
444    // can't loop
445    chosenRow = pivotRow();
446    model_->setLargestDualError(saveError);
447  }
448  if (chosenRow < 0 && lastPivotRow < 0) {
449    int nLeft = 0;
450    for (int i = 0; i < number; i++) {
451      int iRow = index[i];
452      if (fabs(infeas[iRow]) > 1.0e-50) {
453        index[nLeft++] = iRow;
454      } else {
455        infeas[iRow] = 0.0;
456      }
457    }
458    infeasible_->setNumElements(nLeft);
459    model_->setNumberPrimalInfeasibilities(nLeft);
460  }
461
462  modelPE_->updateLastObjectiveValue();
463  return chosenRow;
464}
465/* Save weights - this may initialize weights as well
466   This is as parent but may initialize ClpPESimplex
467*/
468void ClpPEDualRowSteepest::saveWeights(ClpSimplex *model, int mode)
469{
470  // See if we need to initialize ClpPESimplex
471  if (!modelPE_ || model != modelPE_->clpModel() || !modelPE_->checkSize()) {
472    delete modelPE_;
473    modelPE_ = new ClpPESimplex(model);
474  }
475  ClpDualRowSteepest::saveWeights(model, mode);
476}
477/* Updates primal solution (and maybe list of candidates)
478   Uses input vector which it deletes
479   Computes change in objective function
480   As ordinary steepest but checks for zero moves
481*/
482void ClpPEDualRowSteepest::updatePrimalSolution(CoinIndexedVector *input,
483  double theta,
484  double &changeInObjective)
485{
486  int iColumn = model_->sequenceIn();
487  if (iColumn >= 0)
488    modelPE_->updateCompatibleRows(iColumn);
489  ClpDualRowSteepest::updatePrimalSolution(input, theta, changeInObjective);
490}
491
492/* vi: softtabstop=2 shiftwidth=2 expandtab tabstop=2
493*/
Note: See TracBrowser for help on using the repository browser.