source: trunk/Clp/src/ClpPEPrimalColumnDantzig.cpp @ 2385

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

formatting

File size: 12.6 KB
Line 
1// Copyright (C) 2002, International Business Machines
2// Corporation and others.  All Rights Reserved.
3
4/*
5   Authors
6
7   Jeremy Omer, Mehdi Towhidi
8
9   Last update: april 10, 2015
10
11 */
12
13#include "CoinPragma.hpp"
14
15#include <cstdio>
16
17#include "ClpPEPrimalColumnDantzig.hpp"
18
19//-------------------------------------------------------------------
20// Default Constructor
21//-------------------------------------------------------------------
22ClpPEPrimalColumnDantzig::ClpPEPrimalColumnDantzig(double psi)
23  : ClpPrimalColumnDantzig()
24  , modelPE_(NULL)
25  , psi_(psi)
26  , iCurrent_(0)
27  , iInterval_(100)
28  , coDegenCompatibles_(0)
29  , coConsecutiveCompatibles_(0)
30  , updateCompatibles_(true)
31{
32}
33
34//-------------------------------------------------------------------
35// Copy constructor
36//-------------------------------------------------------------------
37ClpPEPrimalColumnDantzig::ClpPEPrimalColumnDantzig(const ClpPEPrimalColumnDantzig &source)
38  : ClpPrimalColumnDantzig(source)
39{
40  modelPE_ = NULL;
41  psi_ = source.psi_;
42  iCurrent_ = source.iCurrent_;
43  iInterval_ = source.iInterval_;
44  updateCompatibles_ = source.updateCompatibles_;
45  coDegenCompatibles_ = source.coDegenCompatibles_;
46  coConsecutiveCompatibles_ = source.coConsecutiveCompatibles_;
47}
48
49//-------------------------------------------------------------------
50// Destructor
51//-------------------------------------------------------------------
52ClpPEPrimalColumnDantzig::~ClpPEPrimalColumnDantzig()
53{
54  delete modelPE_;
55}
56
57//----------------------------------------------------------------
58// Assignment operator
59//-------------------------------------------------------------------
60ClpPEPrimalColumnDantzig &
61ClpPEPrimalColumnDantzig::operator=(const ClpPEPrimalColumnDantzig &rhs)
62{
63  if (this != &rhs) {
64    ClpPrimalColumnDantzig::operator=(rhs);
65    delete modelPE_;
66    modelPE_ = NULL;
67  }
68  return *this;
69}
70
71//-------------------------------------------------------------------
72//// Clone
73////-------------------------------------------------------------------
74ClpPrimalColumnPivot *ClpPEPrimalColumnDantzig::clone(bool CopyData) const
75{
76  if (CopyData) {
77    return new ClpPEPrimalColumnDantzig(*this);
78  } else {
79    return new ClpPEPrimalColumnDantzig(psi_);
80  }
81}
82
83//-------------------------------------------------------------------
84// pivotColumn
85// Returns pivot column, -1 if none
86// With this class, a two-dimensional pricing is applied to choose
87// between the most negative reduced cost and the most negative
88// recuded cost of the compatible variables.
89////-------------------------------------------------------------------
90int ClpPEPrimalColumnDantzig::pivotColumn(CoinIndexedVector *updates,
91  CoinIndexedVector *spareRow1,
92  CoinIndexedVector *spareRow2,
93  CoinIndexedVector *spareColumn1,
94  CoinIndexedVector *spareColumn2)
95{
96  assert(model_);
97
98#if PE_DEBUG >= 1
99  std::cout << "@@@@@@@@@@PEPrimalColumnDantzig::pivotColumn\n";
100#endif
101
102  int iSection, j;
103  int number;
104  int *index;
105  double *updateBy;
106  double *reducedCost;
107
108  bool anyUpdates;
109
110  // If the input vector updates is not empty, the dual solution and
111  // the reduced cost needs to be updated
112  //
113  if (updates->getNumElements()) {
114    anyUpdates = true;
115  } else {
116    // sub flip - nothing to do
117    anyUpdates = false;
118  }
119  if (anyUpdates) {
120// compute y:  y^T * A = updates^T and updates <- spareRow2 + updates
121#ifdef PE_STATISTICS
122    model_->factorization()->doStatistics(false);
123#endif
124    model_->factorization()->updateColumnTranspose(spareRow2, updates);
125#ifdef PE_STATISTICS
126    model_->factorization()->doStatistics(true);
127#endif
128    // put row of tableau in rowArray and columnArray
129    model_->clpMatrix()->transposeTimes(model_, -1.0,
130      updates, spareColumn2, spareColumn1);
131    for (iSection = 0; iSection < 2; iSection++) {
132
133      reducedCost = model_->djRegion(iSection);
134
135      if (!iSection) {
136        number = updates->getNumElements();
137        index = updates->getIndices();
138        updateBy = updates->denseVector();
139      } else {
140        number = spareColumn1->getNumElements();
141        index = spareColumn1->getIndices();
142        updateBy = spareColumn1->denseVector();
143      }
144
145      for (j = 0; j < number; j++) {
146        int iSequence = index[j];
147        double value = reducedCost[iSequence];
148        value -= updateBy[j];
149        updateBy[j] = 0.0;
150        reducedCost[iSequence] = value;
151      }
152    }
153    updates->setNumElements(0);
154    spareColumn1->setNumElements(0);
155  }
156
157  // Determine whether the set of compatible variables should be updated
158  //
159  // store the number of degenerate pivots on compatible variables and the
160  // overal number of degenerate pivots
161  double progress = fabs(modelPE_->lastObjectiveValue() - model_->objectiveValue());
162  bool isLastDegenerate = progress <= 1.0e-12 * fabs(model_->objectiveValue()) ? true : false;
163  if (isLastDegenerate) {
164    modelPE_->addDegeneratePivot();
165    modelPE_->addDegeneratePivotConsecutive();
166
167    if (modelPE_->isLastPivotCompatible()) {
168      modelPE_->addDegenerateCompatiblePivot();
169    }
170  } else {
171    modelPE_->resetDegeneratePivotsConsecutive();
172  }
173
174  // the compatible variables are updated when the number of degenerate pivots
175  // on compatible variables is more than 20% the overall number of degenerate pivots
176  if (modelPE_->isLastPivotCompatible()) {
177    coConsecutiveCompatibles_++;
178    if (isLastDegenerate) {
179      coDegenCompatibles_++;
180      if (coConsecutiveCompatibles_ >= 10 && 5 * coDegenCompatibles_ * model_->numberIterations() > modelPE_->coDegeneratePivots() * coConsecutiveCompatibles_) {
181        updateCompatibles_ = true;
182      }
183    }
184  }
185
186  if (modelPE_->doStatistics()) {
187    // For results comparison.
188    // count the number of degenerate variables if psi==1
189    // add the time spent in counting to the time limit
190    //
191    modelPE_->startTimer();
192    if (psi_ >= 1 && iCurrent_ >= 100) {
193      modelPE_->updateDualDegenerates();
194      modelPE_->updateDualDegeneratesAvg(100);
195      model_->setMaximumSeconds(36000.0 + modelPE_->timeCompatibility() - CoinCpuTime());
196      iCurrent_ = 0;
197    }
198    modelPE_->stopTimer();
199  }
200
201  // Update the set of compatible variables if necessary
202  //
203  if (modelPE_->doStatistics())
204    modelPE_->startTimer();
205  double psiTmp = psi_;
206  if ((psi_ < 1.0) && (iCurrent_ >= iInterval_) && (updateCompatibles_ || iCurrent_ >= 1000)) {
207    // the compatible variables are never updated if the last pivot is non degenerate
208    // this could be counterproductive
209    if (isLastDegenerate) {
210      modelPE_->updatePrimalDegenerates();
211      modelPE_->identifyCompatibleCols(model_->numberRows() + model_->numberColumns(), NULL, spareRow2, spareRow1);
212
213      if (modelPE_->doStatistics()) {
214        // update the average number of degenerates and compatibles for statistics
215        modelPE_->updatePrimalDegeneratesAvg(iCurrent_);
216        modelPE_->updateCompatibleColsAvg(iCurrent_);
217      }
218
219#ifdef PE_TEST
220      std::cout << "coDegen=" << coDegenCompatibles_ << " ; coComp = " << coConsecutiveCompatibles_
221                << " ; iCurrent_ = " << iCurrent_ << " ;compatibleColumns = "
222                << modelPE_->coCompatibleCols() << std::endl;
223#endif
224
225      // the checking frequency is modified to reflect what appears to be needed
226      if (iCurrent_ == iInterval_)
227        iInterval_ = std::max(50, iInterval_ - 50);
228      else
229        iInterval_ = std::min(300, iInterval_ + 50);
230
231      // reset all the indicators that are used to decide whether the compatible
232      // variables must be updated
233      iCurrent_ = 0;
234      updateCompatibles_ = false;
235      coConsecutiveCompatibles_ = 0;
236      coDegenCompatibles_ = 0;
237    } else
238      iInterval_++;
239  }
240  // otherwise, update the value of the priority factor depending on the number of
241  // consecutive degenerate pivots
242  //
243  else {
244    // the idea is that when a lot of consecutive pivots are degenerate, there is
245    // a potentially hign added value in putting a very high priroity on compatible
246    // variables
247    if (modelPE_->coDegeneratePivotsConsecutive() >= 10) {
248      psiTmp = 0.01;
249
250#if PE_DEBUG >= 1
251      std::cout << "Lower the priority factor to " << psiTmp << std::endl;
252      std::cout << "Consecutive degenerate pivots=" << modelPE_->coDegeneratePivotsConsecutive() << std::endl;
253#endif
254    }
255  }
256  iCurrent_++;
257  if (modelPE_->doStatistics())
258    modelPE_->stopTimer();
259
260  // Updates of dual solutions and compatible variables finished,
261  // now do the pricing
262  //
263
264  // we can't really trust infeasibilities if there is primal error
265  double largest = model_->currentPrimalTolerance();
266  if (model_->largestDualError() > 1.0e-8)
267    largest *= model_->largestDualError() / 1.0e-8;
268
269  // initialize the best reduced cost values
270  double dualTolerance = model_->dualTolerance();
271  double bestDj = 1.0e-30;
272  int bestSequence = -1;
273  double bestDjComp = 1.0e-30;
274  int bestSequenceComp = -1;
275
276  number = model_->numberRows() + model_->numberColumns();
277  reducedCost = model_->djRegion();
278
279  // only check the compatible variables when the bidimensional factor is less than 1
280  // and the ratio of compatible variables is larger than 0.01
281  bool checkCompatibles = true;
282  double ratioCompatibles = static_cast< double >(modelPE_->coCompatibleCols()) / static_cast< double >((model_->numberRows() + model_->numberColumns()));
283
284  if (psi_ >= 1.0 || ratioCompatibles < 0.01)
285    checkCompatibles = false;
286
287  int indicesToConsiderCount = 0;
288  for (int iSequence = 0; iSequence < number; iSequence++) {
289    // check flagged variable
290    if (!model_->flagged(iSequence)) {
291      double value = reducedCost[iSequence];
292      double largestDj = std::max(psi_ * bestDj, bestDjComp);
293      ClpSimplex::Status status = model_->getStatus(iSequence);
294
295      // we choose the nonbasic column whose reduced cost is either
296      // the most negative or the most positive one, depending on the
297      // status of the variable
298      // if a variable is compatible, a factor  0 < psi_ < 1 gives
299      // a greater priority to it
300      switch (status) {
301
302        // basic and fixed variable cannot be chosen to enter the basis
303      case ClpSimplex::basic:
304      case ClpSimplex::isFixed:
305        break;
306        // free and superbasic are given priority using the 0.1 factor
307        // since these variables never leave basis once they're in
308      case ClpSimplex::isFree:
309      case ClpSimplex::superBasic:
310        value = fabs(value);
311        if (checkCompatibles && modelPE_->isCompatibleCol(iSequence) && value > 0.1 * bestDjComp) {
312          bestDjComp = 10.0 * value;
313          bestSequenceComp = iSequence;
314        } else if (value > 0.1 * bestDj) {
315          bestDj = 10.0 * value;
316          bestSequence = iSequence;
317        }
318        break;
319        // just be careful with the sign of the reduced cost for the other
320        // variables and give priority to the compatible ones
321      case ClpSimplex::atUpperBound:
322        if (value > largestDj) {
323          if (checkCompatibles && modelPE_->isCompatibleCol(iSequence)) {
324            bestDjComp = value;
325            bestSequenceComp = iSequence;
326          } else if (value > bestDj) {
327            bestDj = value;
328            bestSequence = iSequence;
329          }
330        }
331        break;
332      case ClpSimplex::atLowerBound:
333        if (value < -largestDj) {
334          if (checkCompatibles && modelPE_->isCompatibleCol(iSequence)) {
335            bestDjComp = -value;
336            bestSequenceComp = iSequence;
337          } else if (value < -bestDj) {
338            bestDj = -value;
339            bestSequence = iSequence;
340          }
341        }
342        break;
343      }
344    }
345  }
346  // choose a compatible or an incompatible row depending on the
347  // largest values and on the factor of preference psi_
348  if (modelPE_->doStatistics())
349    modelPE_->startTimer();
350  if (bestSequenceComp >= 0 && bestDjComp >= psiTmp * bestDj) {
351    bestSequence = bestSequenceComp;
352
353    // record the number of pivots done on compatible variables
354    // that would not have been done without positive edge
355    if (modelPE_->doStatistics())
356      if (bestDjComp < bestDj)
357        modelPE_->addPriorityPivot();
358  }
359  if (psi_ < 1 && modelPE_->isCompatibleCol(bestSequence)) {
360    modelPE_->isLastPivotCompatible(true);
361    modelPE_->addCompatiblePivot();
362  } else
363    modelPE_->isLastPivotCompatible(false);
364  if (modelPE_->doStatistics())
365    modelPE_->stopTimer();
366
367  // save the current objective value
368  modelPE_->updateLastObjectiveValue();
369
370  return bestSequence;
371}
372/* Save weights - this may initialize weights as well
373   This is as parent but may initialize ClpPESimplex
374*/
375void ClpPEPrimalColumnDantzig::saveWeights(ClpSimplex *model, int mode)
376{
377  // See if we need to initialize ClpPESimplex
378  if (!modelPE_ || model != modelPE_->clpModel()) {
379    delete modelPE_;
380    modelPE_ = new ClpPESimplex(model);
381  }
382  ClpPrimalColumnDantzig::saveWeights(model, mode);
383}
384
385/* vi: softtabstop=2 shiftwidth=2 expandtab tabstop=2
386*/
Note: See TracBrowser for help on using the repository browser.