source: trunk/Clp/src/ClpPEPrimalColumnSteepest.cpp

Last change on this file was 2385, checked in by unxusr, 9 months ago

formatting

File size: 20.6 KB
Line 
1/* $Id: ClpPrimalColumnSteepest.cpp 1955 2013-05-14 10:10:07Z forrest $ */
2// Copyright (C) 2002, 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   Authors
7
8   Jeremy Omer, Mehdi Towhidi
9
10   Last update: april 10, 2015
11
12 */
13
14#include "CoinPragma.hpp"
15
16#include "ClpPEPrimalColumnSteepest.hpp"
17
18#include "ClpSimplex.hpp"
19#include "ClpMessage.hpp"
20#include "ClpNonLinearCost.hpp"
21
22#include <stdio.h>
23//#define CLP_DEBUG
24//#############################################################################
25// Constructors / Destructor / Assignment
26//#############################################################################
27
28//-------------------------------------------------------------------
29// Default Constructor
30//-------------------------------------------------------------------
31ClpPEPrimalColumnSteepest::ClpPEPrimalColumnSteepest(double psi, int mode)
32  : ClpPrimalColumnSteepest(mode)
33  , modelPE_(NULL)
34  , psi_(psi)
35  , iCurrent_(0)
36  , iInterval_(100)
37  , coDegenCompatibles_(0)
38  , coConsecutiveCompatibles_(0)
39  , updateCompatibles_(true)
40{
41}
42
43//-------------------------------------------------------------------
44// Copy constructor
45//-------------------------------------------------------------------
46ClpPEPrimalColumnSteepest::ClpPEPrimalColumnSteepest(const ClpPEPrimalColumnSteepest &source)
47  : ClpPrimalColumnSteepest(source)
48{
49  modelPE_ = NULL;
50  psi_ = source.psi_;
51  iCurrent_ = source.iCurrent_;
52  iInterval_ = source.iInterval_;
53  updateCompatibles_ = source.updateCompatibles_;
54  coDegenCompatibles_ = source.coDegenCompatibles_;
55  coConsecutiveCompatibles_ = source.coConsecutiveCompatibles_;
56}
57
58//-------------------------------------------------------------------
59// Destructor
60//-------------------------------------------------------------------
61ClpPEPrimalColumnSteepest::~ClpPEPrimalColumnSteepest()
62{
63  delete modelPE_;
64}
65
66//----------------------------------------------------------------
67// Assignment operator
68//-------------------------------------------------------------------
69ClpPEPrimalColumnSteepest &
70ClpPEPrimalColumnSteepest::operator=(const ClpPEPrimalColumnSteepest &rhs)
71{
72  if (this != &rhs) {
73    ClpPrimalColumnSteepest::operator=(rhs);
74    delete modelPE_;
75    modelPE_ = NULL;
76  }
77  return *this;
78}
79
80//-------------------------------------------------------------------
81// Clone
82//-------------------------------------------------------------------
83ClpPrimalColumnPivot *ClpPEPrimalColumnSteepest::clone(bool CopyData) const
84{
85  if (CopyData) {
86    return new ClpPEPrimalColumnSteepest(*this);
87  } else {
88    return new ClpPEPrimalColumnSteepest(psi_);
89  }
90}
91
92// These have to match ClpPackedMatrix version
93#define TRY_NORM 1.0e-4
94#define ADD_ONE 1.0
95
96// Returns pivot column, -1 if none
97/*   The Packed CoinIndexedVector updates has cost updates - for normal LP
98that is just +-weight where a feasibility changed.  It also has
99reduced cost from last iteration in pivot row*/
100int ClpPEPrimalColumnSteepest::pivotColumn(CoinIndexedVector *updates,
101  CoinIndexedVector *spareRow1,
102  CoinIndexedVector *spareRow2,
103  CoinIndexedVector *spareColumn1,
104  CoinIndexedVector *spareColumn2)
105{
106  // check that the model exists and has the proper form
107  assert(model_);
108  assert(!model_->nonLinearCost()->lookBothWays() && model_->algorithm() != 2);
109
110#if PE_DEBUG >= 1
111  std::cout << "@@@@@@@@@@ClpPEPrimalColumnSteepest::pivotColumn\n";
112#endif
113
114  int number = 0;
115  int *index;
116
117  double tolerance = model_->currentDualTolerance();
118  // we can't really trust infeasibilities if there is dual error
119  // this coding has to mimic coding in checkDualSolution
120  double error = CoinMin(1.0e-2, model_->largestDualError());
121  // allow tolerance at least slightly bigger than standard
122  tolerance = tolerance + error;
123  int pivotRow = model_->pivotRow();
124  int anyUpdates;
125  double *infeas = infeasible_->denseVector();
126
127  // Local copy of mode so can decide what to do
128  int switchType;
129  if (mode_ == 4)
130    switchType = 5 - numberSwitched_;
131  else if (mode_ >= 10)
132    switchType = 3;
133  else
134    switchType = mode_;
135
136  /* switchType -
137        0 - all exact devex
138        1 - all steepest
139        2 - some exact devex
140        3 - auto some exact devex
141        4 - devex
142        Default value is switchType = 3
143        I removed the mode >= 5 out of simplicity
144        */
145
146  // Definition in ClpGubMatrix, mode is set to 4
147  model_->clpMatrix()->dualExpanded(model_, updates, NULL, 4);
148
149  if (updates->getNumElements() > 1) {
150    // would have to have two goes for devex, three for steepest
151    anyUpdates = 2;
152  } else if (updates->getNumElements()) {
153    if (updates->getIndices()[0] == pivotRow && fabs(updates->denseVector()[0]) > 1.0e-6) {
154      // reasonable size
155      anyUpdates = 1;
156    } else {
157      // too small
158      anyUpdates = 2;
159    }
160  } else if (pivotSequence_ >= 0) {
161    // just after re-factorization
162    anyUpdates = -1;
163  } else {
164    // sub flip - nothing to do
165    anyUpdates = 0;
166  }
167  int sequenceOut = model_->sequenceOut();
168
169  if (anyUpdates == 1) {
170    if (switchType < 4) {
171      // exact etc when can use dj
172      djsAndSteepest(updates, spareRow2,
173        spareColumn1, spareColumn2);
174    } else {
175      // devex etc when can use dj
176      djsAndDevex(updates, spareRow2,
177        spareColumn1, spareColumn2);
178    }
179  } else if (anyUpdates == -1) {
180    if (switchType < 4) {
181      // exact etc when djs okay
182      justSteepest(updates, spareRow2,
183        spareColumn1, spareColumn2);
184    } else {
185      // devex etc when djs okay
186      justDevex(updates, spareRow2,
187        spareColumn1, spareColumn2);
188    }
189  } else if (anyUpdates == 2) {
190    if (switchType < 4) {
191      // exact etc when have to use pivot
192      djsAndSteepest2(updates, spareRow2,
193        spareColumn1, spareColumn2);
194    } else {
195      // devex etc when have to use pivot
196      djsAndDevex2(updates, spareRow2,
197        spareColumn1, spareColumn2);
198    }
199  }
200
201  // make sure outgoing from last iteration okay
202  if (sequenceOut >= 0) {
203    ClpSimplex::Status status = model_->getStatus(sequenceOut);
204    double value = model_->reducedCost(sequenceOut);
205
206    switch (status) {
207
208    case ClpSimplex::basic:
209    case ClpSimplex::isFixed:
210      break;
211    case ClpSimplex::isFree:
212    case ClpSimplex::superBasic:
213      if (fabs(value) > FREE_ACCEPT * tolerance) {
214        // we are going to bias towards free (but only if reasonable)
215        value *= FREE_BIAS;
216        // store square in list
217        if (infeas[sequenceOut])
218          infeas[sequenceOut] = value * value; // already there
219        else
220          infeasible_->quickAdd(sequenceOut, value * value);
221      } else {
222        infeasible_->zero(sequenceOut);
223      }
224      break;
225    case ClpSimplex::atUpperBound:
226      if (value > tolerance) {
227        // store square in list
228        if (infeas[sequenceOut])
229          infeas[sequenceOut] = value * value; // already there
230        else
231          infeasible_->quickAdd(sequenceOut, value * value);
232      } else {
233        infeasible_->zero(sequenceOut);
234      }
235      break;
236    case ClpSimplex::atLowerBound:
237      if (value < -tolerance) {
238        // store square in list
239        if (infeas[sequenceOut])
240          infeas[sequenceOut] = value * value; // already there
241        else
242          infeasible_->quickAdd(sequenceOut, value * value);
243      } else {
244        infeasible_->zero(sequenceOut);
245      }
246    }
247  }
248
249  /* Determine whether the set of compatible variables should be updated
250        */
251  // store the number of degenerate pivots on compatible variables and the
252  // overall number of degenerate pivots
253  //double progress = fabs(modelPE_->lastObjectiveValue() - model_->objectiveValue());
254  //bool isLastDegenerate = progress <= 1.0e-12*fabs(model_->objectiveValue()) ? true:false;
255  bool isLastDegenerate = fabs(model_->theta()) < 1.01e-7;
256  bool isLastDegenerate2;
257  // could fine tune a bit e.g. use 0.0 rather than tolerance
258  if (model_->directionOut() < 0) {
259    // going out to lower bound
260    isLastDegenerate2 = (model_->valueOut() - model_->lowerOut() < model_->primalTolerance());
261  } else {
262    // going out to upper bound
263    isLastDegenerate2 = (model_->upperOut() - model_->valueOut() < model_->primalTolerance());
264  }
265#if 0
266        if ( isLastDegenerate2 &&!isLastDegenerate)
267          printf("PE2 thinks degenerate - theta %g value %g - to %s bound\n",
268                 model_->theta(),model_->valueOut(),
269                 (model_->directionOut()<0) ? "lower" : "upper");
270        else if ( !isLastDegenerate2 &&isLastDegenerate)
271          printf("PE2 thinks not degenerate - theta %g value %g - to %s bound\n",
272                 model_->theta(),model_->valueOut(),
273                 (model_->directionOut()<0) ? "lower" : "upper");
274#else
275  isLastDegenerate = isLastDegenerate2;
276#endif
277  if (isLastDegenerate) {
278    modelPE_->addDegeneratePivot();
279    modelPE_->addDegeneratePivotConsecutive();
280
281    if (modelPE_->isLastPivotCompatible()) {
282      modelPE_->addDegenerateCompatiblePivot();
283    }
284  } else {
285    modelPE_->resetDegeneratePivotsConsecutive();
286  }
287
288  // the compatible variables are updated when the number of degenerate pivots
289  // on compatible variables is more than 20% the overall number of degenerate pivots
290  if (modelPE_->isLastPivotCompatible()) {
291    coConsecutiveCompatibles_++;
292    if (isLastDegenerate) {
293      coDegenCompatibles_++;
294      if (coConsecutiveCompatibles_ >= 10 && 5 * coDegenCompatibles_ * model_->numberIterations() > modelPE_->coDegeneratePivots() * coConsecutiveCompatibles_) {
295        updateCompatibles_ = true;
296      }
297    }
298  }
299
300  if (modelPE_->doStatistics()) {
301    /* For results comparison.
302        count the number of degenerate variables if psi==1
303        add the time spent in counting to the time limit
304        */
305    modelPE_->startTimer();
306    if (psi_ >= 1 && iCurrent_ >= 100) {
307      modelPE_->updateDualDegenerates();
308      modelPE_->updateDualDegeneratesAvg(100);
309      //model_->setMaximumSeconds(36000.0+modelPE_->timeCompatibility()-CoinCpuTime());
310      iCurrent_ = 0;
311    }
312    modelPE_->stopTimer();
313  }
314
315  /* Update the set of compatible variables if necessary
316        */
317  if (modelPE_->doStatistics())
318    modelPE_->startTimer();
319  double psiTmp = psi_;
320  if ((psi_ < 1.0) && (iCurrent_ >= iInterval_) && (updateCompatibles_ || iCurrent_ >= 1000)) {
321    // the compatible variables are never updated if the last pivot is non degenerate
322    // this could be counterproductive
323    if (isLastDegenerate) {
324      modelPE_->updatePrimalDegenerates();
325      for (int i = 0; i < 4; i++)
326        assert(!model_->rowArray(i)->getNumElements());
327      modelPE_->identifyCompatibleCols(model_->numberRows() + model_->numberColumns(), NULL, spareRow2, spareRow1);
328      for (int i = 0; i < 4; i++)
329        assert(!model_->rowArray(i)->getNumElements());
330
331      if (modelPE_->doStatistics()) {
332        // update the average number of degenerates and compatibles for statistics
333        modelPE_->updatePrimalDegeneratesAvg(iCurrent_);
334        modelPE_->updateCompatibleColsAvg(iCurrent_);
335        if (modelPE_->doStatistics() > 3) {
336          char generalPrint[100];
337
338          sprintf(generalPrint, "coDegen = %d; coComp = %d; iCurrent_ = %d; compatibleColumns = %d",
339            coDegenCompatibles_, coConsecutiveCompatibles_,
340            iCurrent_, modelPE_->coCompatibleCols());
341          model_->messageHandler()->message(CLP_GENERAL,
342            *model_->messagesPointer())
343            << generalPrint << CoinMessageEol;
344        }
345      }
346      // switch off ? - maybe just until drops back
347      if (modelPE_->coCompatibleCols() * 1.0 > model_->numberColumns()) {
348        printf("switching off pe\n");
349        psi_ = 1.0;
350      }
351
352      // the checking frequency is modified to reflect what appears to be needed
353      if (iCurrent_ == iInterval_)
354        iInterval_ = std::max(50, iInterval_ - 50);
355      else
356        iInterval_ = std::min(300, iInterval_ + 50);
357
358      // reset all the indicators that are used to decide whether the compatible
359      // variables must be updated
360      iCurrent_ = 0;
361      updateCompatibles_ = false;
362      coConsecutiveCompatibles_ = 0;
363      coDegenCompatibles_ = 0;
364    } else
365      iInterval_++;
366  }
367  /* otherwise, update the value of the priority factor depending on the number of
368         * consecutive degenerate pivots
369         */
370  else {
371    // the idea is that when a lot of consecutive pivots are degenerate, there is
372    // a potentially high added value in putting a very high priority on compatible
373    // variables
374    if (modelPE_->coDegeneratePivotsConsecutive() >= 10) {
375      psiTmp = 0.01;
376      psiTmp = 0.25 * psi_;
377
378#if PE_DEBUG >= 1
379      std::cout << "Lower the priority factor to " << psiTmp << std::endl;
380      std::cout << "Consecutive degenerate pivots=" << modelPE_->coDegeneratePivotsConsecutive() << std::endl;
381#endif
382    }
383  }
384  iCurrent_++;
385  if (modelPE_->doStatistics())
386    modelPE_->stopTimer();
387
388  /* Update of the dual solution and compatible variables is finished,
389        ** now do the pricing */
390
391  // See what sort of pricing
392  int numberWanted = 10;
393  number = infeasible_->getNumElements();
394  int numberColumns = model_->numberColumns();
395  int numberRows = model_->numberRows();
396  // ratio is done on number of rows here
397  double ratio = static_cast< double >(sizeFactorization_) / static_cast< double >(numberRows);
398  if (switchType == 4) {
399    // Still in devex mode
400    // Go to steepest if lot of iterations?
401    if (ratio < 5.0) {
402      numberWanted = CoinMax(2000, number / 10);
403      numberWanted = CoinMax(numberWanted, numberColumns / 20);
404    } else if (ratio < 7.0) {
405      numberWanted = CoinMax(2000, number / 5);
406      numberWanted = CoinMax(numberWanted, numberColumns / 10);
407    } else {
408      // we can zero out
409      updates->clear();
410      spareColumn1->clear();
411      switchType = 3;
412      // initialize
413      pivotSequence_ = -1;
414      pivotRow = -1;
415      numberSwitched_++;
416      // Make sure will re-do
417      delete[] weights_;
418      weights_ = NULL;
419      saveWeights(model_, 4);
420      COIN_DETAIL_PRINT(printf("switching to exact %d nel ratio %g\n", sizeFactorization_, ratio));
421      updates->clear();
422    }
423    if (model_->numberIterations() % 1000 == 0)
424      COIN_DETAIL_PRINT(printf("numels %d ratio %g wanted %d type x\n", sizeFactorization_, ratio, numberWanted));
425  }
426  if (switchType < 4) {
427    if (switchType < 2) {
428      numberWanted = number + 1;
429    } else if (switchType == 2) {
430      numberWanted = CoinMax(2000, number / 8);
431    } else {
432      if (ratio < 1.0) {
433        numberWanted = CoinMax(2000, number / 20);
434      } else if (ratio < 5.0) {
435        numberWanted = CoinMax(2000, number / 10);
436        numberWanted = CoinMax(numberWanted, numberColumns / 40);
437      } else if (ratio < 10.0) {
438        numberWanted = CoinMax(2000, number / 8);
439        numberWanted = CoinMax(numberWanted, numberColumns / 20);
440      } else {
441        ratio = number * (ratio / 80.0);
442        if (ratio > number) {
443          numberWanted = number + 1;
444        } else {
445          numberWanted = CoinMax(2000, static_cast< int >(ratio));
446          numberWanted = CoinMax(numberWanted, numberColumns / 10);
447        }
448      }
449    }
450  }
451
452  // initialize the best reduced cost values
453  double bestDj = 1.0e-30;
454  int bestSequence = -1;
455  double bestDjComp = 1.0e-30;
456  int bestSequenceComp = -1;
457
458  int i, iSequence;
459  index = infeasible_->getIndices();
460  number = infeasible_->getNumElements();
461  // Re-sort infeasible every 100 pivots
462  // Not a good idea
463  if (0 && model_->factorization()->pivots() > 0 && (model_->factorization()->pivots() % 100) == 0) {
464    int nLook = model_->numberRows() + numberColumns;
465    number = 0;
466    for (i = 0; i < nLook; i++) {
467      if (infeas[i]) {
468        if (fabs(infeas[i]) > COIN_INDEXED_TINY_ELEMENT)
469          index[number++] = i;
470        else
471          infeas[i] = 0.0;
472      }
473    }
474    infeasible_->setNumElements(number);
475  }
476  if (model_->numberIterations() < model_->lastBadIteration() + 200 && model_->factorization()->pivots() > 10) {
477    // we can't really trust infeasibilities if there is dual error
478    double checkTolerance = 1.0e-8;
479    if (model_->largestDualError() > checkTolerance)
480      tolerance *= model_->largestDualError() / checkTolerance;
481    // But cap
482    tolerance = CoinMin(1000.0, tolerance);
483  }
484
485  // stop last one coming immediately
486  double saveOutInfeasibility = 0.0;
487  if (sequenceOut >= 0) {
488    saveOutInfeasibility = infeas[sequenceOut];
489    infeas[sequenceOut] = 0.0;
490  }
491  if (model_->factorization()->pivots() && model_->numberPrimalInfeasibilities())
492    tolerance = CoinMax(tolerance, 1.0e-10 * model_->infeasibilityCost());
493  tolerance *= tolerance; // as we are using squares
494
495  // only check the compatible variables when the bidimensional factor is less than 1
496  // and the ratio of compatible variables is larger than 0.01
497  bool checkCompatibles = true;
498  double ratioCompatibles = static_cast< double >(modelPE_->coCompatibleCols()) / static_cast< double >((model_->numberRows() + model_->numberColumns()));
499  double ratioCompatibles2 = static_cast< double >(modelPE_->coCompatibleCols()) / static_cast< double >(model_->numberColumns());
500
501  if (psi_ >= 1.0 || ratioCompatibles < 0.01)
502    checkCompatibles = false;
503  if (ratioCompatibles2 > 0.5)
504    checkCompatibles = false;
505  // proceed to a partial pricing in two phases checking the weighted reduced
506  // costs of the variables until numberWanted variables have been checked,
507  // the number of variables explored in the first phase is randomly generated
508  int iPass;
509  int start[4];
510  start[1] = number;
511  start[2] = 0;
512  double dstart = static_cast< double >(number) * model_->randomNumberGenerator()->randomDouble();
513  start[0] = static_cast< int >(dstart);
514  start[3] = start[0];
515  for (iPass = 0; iPass < 2; iPass++) {
516    int end = start[2 * iPass + 1];
517
518    for (i = start[2 * iPass]; i < end; i++) {
519      iSequence = index[i];
520      double value = infeas[iSequence];
521      double weight = weights_[iSequence];
522      double weightedDj = weight * bestDj;
523      double largestWeightedDj = std::max(psi_ * weightedDj, weight * bestDjComp);
524      if (value > tolerance) {
525        if (value > largestWeightedDj) {
526          if (model_->flagged(iSequence)) {
527            // just to make sure we don't exit before got something
528            numberWanted++;
529          } else if (checkCompatibles && modelPE_->isCompatibleCol(iSequence)) {
530            // the condition on value is sufficient if the variable is compatible
531            bestDjComp = value / weight;
532            bestSequenceComp = iSequence;
533          } else if (value > weightedDj) {
534            // the usual condition is checked if the variable is not compatible
535            bestDj = value / weight;
536            bestSequence = iSequence;
537          }
538        }
539        numberWanted--;
540      }
541      if (!numberWanted)
542        break;
543    }
544    if (!numberWanted)
545      break;
546  }
547  if (bestSequence >= 0 && model_->getStatus(bestSequence) == ClpSimplex::isFree) {
548    printf("Free in %d compat %c dj %g\n",
549      bestSequence, modelPE_->isCompatibleCol(bestSequence) ? 'y' : 'n', bestDj);
550    bestDjComp = 0.0;
551  }
552  // choose a compatible or an incompatible row depending on the
553  // largest values and on the factor of preference psi_
554  if (bestSequenceComp >= 0 && bestDjComp >= psiTmp * bestDj) {
555
556#if 0
557        if (bestDjComp < bestDj)
558          printf("pivoting on column %d bestDjComp %g bestDj %g\n",
559                 bestSequenceComp,bestDjComp,bestDj);
560        else if (bestSequenceComp==bestSequence)
561          printf("pivoting on column %d bestDjComp %g == bestDj %g\n",
562                 bestSequenceComp,bestDjComp,bestDj);
563#endif
564    // record the number of pivots done on compatible variables
565    // that would not have been done without positive edge
566    if (modelPE_->doStatistics()) {
567      if (bestDjComp < bestDj)
568        modelPE_->addPriorityPivot();
569    }
570    bestSequence = bestSequenceComp;
571  }
572  if (bestSequence >= 0 && psi_ < 1.0 && modelPE_->isCompatibleCol(bestSequence)) {
573    modelPE_->isLastPivotCompatible(true);
574    modelPE_->addCompatiblePivot();
575  } else
576    modelPE_->isLastPivotCompatible(false);
577
578  model_->clpMatrix()->setSavedBestSequence(bestSequence);
579  if (bestSequence >= 0)
580    model_->clpMatrix()->setSavedBestDj(model_->djRegion()[bestSequence]);
581  if (sequenceOut >= 0) {
582    infeas[sequenceOut] = saveOutInfeasibility;
583  }
584
585  modelPE_->updateLastObjectiveValue();
586  return bestSequence;
587}
588/* Save weights - this may initialize weights as well
589   This is as parent but may initialize ClpPESimplex
590*/
591void ClpPEPrimalColumnSteepest::saveWeights(ClpSimplex *model, int mode)
592{
593  // See if we need to initialize ClpPESimplex
594  if (!modelPE_ || model != modelPE_->clpModel() || !modelPE_->checkSize()) {
595    delete modelPE_;
596    modelPE_ = new ClpPESimplex(model);
597  }
598  ClpPrimalColumnSteepest::saveWeights(model, mode);
599}
600// Updates weights - as ordinary but checks for zero moves
601void ClpPEPrimalColumnSteepest::updateWeights(CoinIndexedVector *input)
602{
603  //if (modelPE_->isLastPivotCompatible()) {
604  //double theta = modelPE_->model()->theta();
605  //if (theta
606  //}
607  ClpPrimalColumnSteepest::updateWeights(input);
608}
609
610/* vi: softtabstop=2 shiftwidth=2 expandtab tabstop=2
611*/
Note: See TracBrowser for help on using the repository browser.