source: trunk/Clp/src/ClpPEPrimalColumnSteepest.cpp @ 2149

Last change on this file since 2149 was 2149, checked in by forrest, 4 years ago

add PE code

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