source: branches/pre/ClpPrimalColumnSteepest.cpp @ 180

Last change on this file since 180 was 180, checked in by forrest, 18 years ago

new stuff

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 36.7 KB
Line 
1// Copyright (C) 2002, International Business Machines
2// Corporation and others.  All Rights Reserved.
3
4#include "CoinPragma.hpp"
5
6#include "ClpSimplex.hpp"
7#include "ClpPrimalColumnSteepest.hpp"
8#include "CoinIndexedVector.hpp"
9#include "ClpFactorization.hpp"
10#include "ClpNonLinearCost.hpp"
11#include "ClpMessage.hpp"
12#include "CoinHelperFunctions.hpp"
13#include <stdio.h>
14
15
16// bias for free variables
17#define FREE_BIAS 1.0e1
18// Acceptance criteria for free variables
19#define FREE_ACCEPT 1.0e2
20//#############################################################################
21// Constructors / Destructor / Assignment
22//#############################################################################
23
24//-------------------------------------------------------------------
25// Default Constructor
26//-------------------------------------------------------------------
27ClpPrimalColumnSteepest::ClpPrimalColumnSteepest (int mode) 
28  : ClpPrimalColumnPivot(),
29    devex_(0.0),
30    weights_(NULL),
31    infeasible_(NULL),
32    alternateWeights_(NULL),
33    savedWeights_(NULL),
34    reference_(NULL),
35    state_(-1),
36    mode_(mode),
37    pivotSequence_(-1),
38    savedPivotSequence_(-1),
39    savedSequenceOut_(-1)
40{
41  type_=2+64*mode;
42}
43
44//-------------------------------------------------------------------
45// Copy constructor
46//-------------------------------------------------------------------
47ClpPrimalColumnSteepest::ClpPrimalColumnSteepest (const ClpPrimalColumnSteepest & rhs) 
48: ClpPrimalColumnPivot(rhs)
49{ 
50  state_=rhs.state_;
51  mode_ = rhs.mode_;
52  model_ = rhs.model_;
53  pivotSequence_ = rhs.pivotSequence_;
54  savedPivotSequence_ = rhs.savedPivotSequence_;
55  savedSequenceOut_ = rhs.savedSequenceOut_;
56  devex_ = rhs.devex_;
57  if (rhs.infeasible_) {
58    infeasible_= new CoinIndexedVector(rhs.infeasible_);
59  } else {
60    infeasible_=NULL;
61  }
62  reference_=NULL;
63  if (rhs.weights_) {
64    assert(model_);
65    int number = model_->numberRows()+model_->numberColumns();
66    weights_= new double[number];
67    ClpDisjointCopyN(rhs.weights_,number,weights_);
68    savedWeights_= new double[number];
69    ClpDisjointCopyN(rhs.savedWeights_,number,savedWeights_);
70    if (!mode_) {
71      reference_ = new unsigned int[(number+31)>>5];
72      memcpy(reference_,rhs.reference_,((number+31)>>5)*sizeof(unsigned int));
73    }
74  } else {
75    weights_=NULL;
76    savedWeights_=NULL;
77  }
78  if (rhs.alternateWeights_) {
79    alternateWeights_= new CoinIndexedVector(rhs.alternateWeights_);
80  } else {
81    alternateWeights_=NULL;
82  }
83}
84
85//-------------------------------------------------------------------
86// Destructor
87//-------------------------------------------------------------------
88ClpPrimalColumnSteepest::~ClpPrimalColumnSteepest ()
89{
90  delete [] weights_;
91  delete infeasible_;
92  delete alternateWeights_;
93  delete [] savedWeights_;
94  delete [] reference_;
95}
96
97//----------------------------------------------------------------
98// Assignment operator
99//-------------------------------------------------------------------
100ClpPrimalColumnSteepest &
101ClpPrimalColumnSteepest::operator=(const ClpPrimalColumnSteepest& rhs)
102{
103  if (this != &rhs) {
104    ClpPrimalColumnPivot::operator=(rhs);
105    state_=rhs.state_;
106    mode_ = rhs.mode_;
107    model_ = rhs.model_;
108    pivotSequence_ = rhs.pivotSequence_;
109    savedPivotSequence_ = rhs.savedPivotSequence_;
110    savedSequenceOut_ = rhs.savedSequenceOut_;
111    devex_ = rhs.devex_;
112    delete [] weights_;
113    delete [] reference_;
114    reference_=NULL;
115    delete infeasible_;
116    delete alternateWeights_;
117    delete [] savedWeights_;
118    savedWeights_ = NULL;
119    if (rhs.infeasible_!=NULL) {
120      infeasible_= new CoinIndexedVector(rhs.infeasible_);
121    } else {
122      infeasible_=NULL;
123    }
124    if (rhs.weights_!=NULL) {
125      assert(model_);
126      int number = model_->numberRows()+model_->numberColumns();
127      weights_= new double[number];
128      ClpDisjointCopyN(rhs.weights_,number,weights_);
129      savedWeights_= new double[number];
130      ClpDisjointCopyN(rhs.savedWeights_,number,savedWeights_);
131      if (!mode_) {
132        reference_ = new unsigned int[(number+31)>>5];
133        memcpy(reference_,rhs.reference_,
134               ((number+31)>>5)*sizeof(unsigned int));
135      }
136    } else {
137      weights_=NULL;
138    }
139    if (rhs.alternateWeights_!=NULL) {
140      alternateWeights_= new CoinIndexedVector(rhs.alternateWeights_);
141    } else {
142      alternateWeights_=NULL;
143    }
144  }
145  return *this;
146}
147
148#define TRY_NORM 1.0e-4
149#define ADD_ONE 1.0
150// Returns pivot column, -1 if none
151int 
152ClpPrimalColumnSteepest::pivotColumn(CoinIndexedVector * updates,
153                                    CoinIndexedVector * spareRow1,
154                                    CoinIndexedVector * spareRow2,
155                                    CoinIndexedVector * spareColumn1,
156                                    CoinIndexedVector * spareColumn2)
157{
158  assert(model_);
159  int iSection,j;
160  int number;
161  int * index;
162  double * updateBy;
163  double * reducedCost;
164  // dj could be very small (or even zero - take care)
165  double dj = model_->dualIn();
166  double * infeas = infeasible_->denseVector();
167  double tolerance=model_->currentDualTolerance();
168  // we can't really trust infeasibilities if there is dual error
169  // this coding has to mimic coding in checkDualSolution
170  double error = min(1.0e-3,model_->largestDualError());
171  // allow tolerance at least slightly bigger than standard
172  tolerance = tolerance +  error;
173  int pivotRow = model_->pivotRow();
174
175  int anyUpdates;
176
177  if (updates->getNumElements()) {
178    // would have to have two goes for devex, three for steepest
179    anyUpdates=2;
180    // add in pivot contribution
181    if (pivotRow>=0) 
182      updates->add(pivotRow,-dj);
183  } else if (pivotRow>=0) {
184    if (fabs(dj)>1.0e-15) {
185      // some dj
186      updates->insert(pivotRow,-dj);
187      if (fabs(dj)>1.0e-6) {
188        // reasonable size
189        anyUpdates=1;
190      } else {
191        // too small
192        anyUpdates=2;
193      }
194    } else {
195      // zero dj
196      anyUpdates=-1;
197    }
198  } else if (pivotSequence_>=0){
199    // just after re-factorization
200    anyUpdates=-1;
201  } else {
202    // sub flip - nothing to do
203    anyUpdates=0;
204  }
205
206  if (anyUpdates>0) {
207    model_->factorization()->updateColumnTranspose(spareRow2,updates);
208   
209    // put row of tableau in rowArray and columnArray
210    model_->clpMatrix()->transposeTimes(model_,-1.0,
211                                        updates,spareColumn2,spareColumn1);
212    // normal
213    for (iSection=0;iSection<2;iSection++) {
214     
215      reducedCost=model_->djRegion(iSection);
216      int addSequence;
217     
218      if (!iSection) {
219        number = updates->getNumElements();
220        index = updates->getIndices();
221        updateBy = updates->denseVector();
222        addSequence = model_->numberColumns();
223      } else {
224        number = spareColumn1->getNumElements();
225        index = spareColumn1->getIndices();
226        updateBy = spareColumn1->denseVector();
227        addSequence = 0;
228      }
229      if (!model_->nonLinearCost()->lookBothWays()) {
230
231        for (j=0;j<number;j++) {
232          int iSequence = index[j];
233          double value = reducedCost[iSequence];
234          value -= updateBy[iSequence];
235          reducedCost[iSequence] = value;
236          ClpSimplex::Status status = model_->getStatus(iSequence+addSequence);
237         
238          switch(status) {
239           
240          case ClpSimplex::basic:
241            infeasible_->zero(iSequence+addSequence);
242          case ClpSimplex::isFixed:
243            break;
244          case ClpSimplex::isFree:
245          case ClpSimplex::superBasic:
246            if (fabs(value)>FREE_ACCEPT*tolerance) {
247              // we are going to bias towards free (but only if reasonable)
248              value *= FREE_BIAS;
249              // store square in list
250              if (infeas[iSequence+addSequence])
251                infeas[iSequence+addSequence] = value*value; // already there
252              else
253                infeasible_->quickAdd(iSequence+addSequence,value*value);
254            } else {
255              infeasible_->zero(iSequence+addSequence);
256            }
257            break;
258          case ClpSimplex::atUpperBound:
259            if (value>tolerance) {
260              // store square in list
261              if (infeas[iSequence+addSequence])
262                infeas[iSequence+addSequence] = value*value; // already there
263              else
264                infeasible_->quickAdd(iSequence+addSequence,value*value);
265            } else {
266              infeasible_->zero(iSequence+addSequence);
267            }
268            break;
269          case ClpSimplex::atLowerBound:
270            if (value<-tolerance) {
271              // store square in list
272              if (infeas[iSequence+addSequence])
273                infeas[iSequence+addSequence] = value*value; // already there
274              else
275                infeasible_->quickAdd(iSequence+addSequence,value*value);
276            } else {
277              infeasible_->zero(iSequence+addSequence);
278            }
279          }
280        }
281      } else {
282        ClpNonLinearCost * nonLinear = model_->nonLinearCost(); 
283        // We can go up OR down
284        for (j=0;j<number;j++) {
285          int iSequence = index[j];
286          double value = reducedCost[iSequence];
287          value -= updateBy[iSequence];
288          reducedCost[iSequence] = value;
289          ClpSimplex::Status status = model_->getStatus(iSequence+addSequence);
290         
291          switch(status) {
292           
293          case ClpSimplex::basic:
294            infeasible_->zero(iSequence+addSequence);
295          case ClpSimplex::isFixed:
296            break;
297          case ClpSimplex::isFree:
298          case ClpSimplex::superBasic:
299            if (fabs(value)>FREE_ACCEPT*tolerance) {
300              // we are going to bias towards free (but only if reasonable)
301              value *= FREE_BIAS;
302              // store square in list
303              if (infeas[iSequence+addSequence])
304                infeas[iSequence+addSequence] = value*value; // already there
305              else
306                infeasible_->quickAdd(iSequence+addSequence,value*value);
307            } else {
308              infeasible_->zero(iSequence+addSequence);
309            }
310            break;
311          case ClpSimplex::atUpperBound:
312            if (value>tolerance) {
313              // store square in list
314              if (infeas[iSequence+addSequence])
315                infeas[iSequence+addSequence] = value*value; // already there
316              else
317                infeasible_->quickAdd(iSequence+addSequence,value*value);
318            } else {
319              // look other way - change up should be negative
320              value -= nonLinear->changeUpInCost(iSequence+addSequence);
321              if (value>-tolerance) {
322                infeasible_->zero(iSequence+addSequence);
323              } else {
324                // store square in list
325                if (infeas[iSequence+addSequence])
326                  infeas[iSequence+addSequence] = value*value; // already there
327                else
328                  infeasible_->quickAdd(iSequence+addSequence,value*value);
329              }
330            }
331            break;
332          case ClpSimplex::atLowerBound:
333            if (value<-tolerance) {
334              // store square in list
335              if (infeas[iSequence+addSequence])
336                infeas[iSequence+addSequence] = value*value; // already there
337              else
338                infeasible_->quickAdd(iSequence+addSequence,value*value);
339            } else {
340              // look other way - change down should be positive
341              value -= nonLinear->changeDownInCost(iSequence+addSequence);
342              if (value<tolerance) {
343                infeasible_->zero(iSequence+addSequence);
344              } else {
345                // store square in list
346                if (infeas[iSequence+addSequence])
347                  infeas[iSequence+addSequence] = value*value; // already there
348                else
349                  infeasible_->quickAdd(iSequence+addSequence,value*value);
350              }
351            }
352          }
353        }
354      }
355    }
356    if (anyUpdates==2) {
357      // we can zero out as will have to get pivot row
358      updates->clear();
359      spareColumn1->clear();
360    }
361    if (pivotRow>=0) {
362      // make sure infeasibility on incoming is 0.0
363      int sequenceIn = model_->sequenceIn();
364      infeasible_->zero(sequenceIn);
365    }
366  }
367  // make sure outgoing from last iteration okay
368  int sequenceOut = model_->sequenceOut();
369  if (sequenceOut>=0) {
370    ClpSimplex::Status status = model_->getStatus(sequenceOut);
371    double value = model_->reducedCost(sequenceOut);
372   
373    switch(status) {
374     
375    case ClpSimplex::basic:
376    case ClpSimplex::isFixed:
377      break;
378    case ClpSimplex::isFree:
379    case ClpSimplex::superBasic:
380      if (fabs(value)>FREE_ACCEPT*tolerance) { 
381        // we are going to bias towards free (but only if reasonable)
382        value *= FREE_BIAS;
383        // store square in list
384        if (infeas[sequenceOut])
385          infeas[sequenceOut] = value*value; // already there
386        else
387          infeasible_->quickAdd(sequenceOut,value*value);
388      } else {
389        infeasible_->zero(sequenceOut);
390      }
391      break;
392    case ClpSimplex::atUpperBound:
393      if (value>tolerance) {
394        // store square in list
395        if (infeas[sequenceOut])
396          infeas[sequenceOut] = value*value; // already there
397        else
398          infeasible_->quickAdd(sequenceOut,value*value);
399      } else {
400        infeasible_->zero(sequenceOut);
401      }
402      break;
403    case ClpSimplex::atLowerBound:
404      if (value<-tolerance) {
405        // store square in list
406        if (infeas[sequenceOut])
407          infeas[sequenceOut] = value*value; // already there
408        else
409          infeasible_->quickAdd(sequenceOut,value*value);
410      } else {
411        infeasible_->zero(sequenceOut);
412      }
413    }
414  }
415
416#ifdef RECOMPUTE_ALL_DJS
417  for (iSection=0;iSection<2;iSection++) {
418     
419    reducedCost=model_->djRegion(iSection);
420    int addSequence;
421    int iSequence;
422     
423    if (!iSection) {
424      number = model_->numberRows();
425      addSequence = model_->numberColumns();
426    } else {
427      number = model_->numberColumns();
428      addSequence = 0;
429    }
430
431    if (!model_->nonLinearCost()->lookBothWays()) {
432      for (iSequence=0;iSequence<number;iSequence++) {
433        double value = reducedCost[iSequence];
434        ClpSimplex::Status status = model_->getStatus(iSequence+addSequence);
435
436        switch(status) {
437         
438        case ClpSimplex::basic:
439          infeasible_->zero(iSequence+addSequence);
440        case ClpSimplex::isFixed:
441          break;
442        case ClpSimplex::isFree:
443          if (fabs(value)>tolerance) {
444            // we are going to bias towards free (but only if reasonable)
445            value *= FREE_BIAS;
446            // store square in list
447            if (infeas[iSequence+addSequence])
448              infeas[iSequence+addSequence] = value*value; // already there
449            else
450              infeasible_->quickAdd(iSequence+addSequence,value*value);
451          } else {
452            infeasible_->zero(iSequence+addSequence);
453          }
454          break;
455        case ClpSimplex::atUpperBound:
456          if (value>tolerance) {
457            // store square in list
458            if (infeas[iSequence+addSequence])
459              infeas[iSequence+addSequence] = value*value; // already there
460            else
461              infeasible_->quickAdd(iSequence+addSequence,value*value);
462          } else {
463            infeasible_->zero(iSequence+addSequence);
464          }
465          break;
466        case ClpSimplex::atLowerBound:
467          if (value<-tolerance) {
468            // store square in list
469            if (infeas[iSequence+addSequence])
470              infeas[iSequence+addSequence] = value*value; // already there
471            else
472              infeasible_->quickAdd(iSequence+addSequence,value*value);
473          } else {
474            infeasible_->zero(iSequence+addSequence);
475          }
476        }
477      }
478    } else {
479      // we can go both ways
480      ClpNonLinearCost * nonLinear = model_->nonLinearCost(); 
481      for (iSequence=0;iSequence<number;iSequence++) {
482        double value = reducedCost[iSequence];
483        ClpSimplex::Status status = model_->getStatus(iSequence+addSequence);
484
485        switch(status) {
486         
487        case ClpSimplex::basic:
488          infeasible_->zero(iSequence+addSequence);
489        case ClpSimplex::isFixed:
490          break;
491        case ClpSimplex::isFree:
492          if (fabs(value)>tolerance) {
493            // we are going to bias towards free (but only if reasonable)
494            value *= FREE_BIAS;
495            // store square in list
496            if (infeas[iSequence+addSequence])
497              infeas[iSequence+addSequence] = value*value; // already there
498            else
499              infeasible_->quickAdd(iSequence+addSequence,value*value);
500          } else {
501            infeasible_->zero(iSequence+addSequence);
502          }
503          break;
504        case ClpSimplex::atUpperBound:
505          if (value>tolerance) {
506            // store square in list
507            if (infeas[iSequence+addSequence])
508              infeas[iSequence+addSequence] = value*value; // already there
509            else
510              infeasible_->quickAdd(iSequence+addSequence,value*value);
511          } else {
512            // look other way - change up should be negative
513            value -= nonLinear->changeUpInCost(iSequence+addSequence);
514            if (value>-tolerance) {
515              infeasible_->zero(iSequence+addSequence);
516            } else {
517              // store square in list
518              if (infeas[iSequence+addSequence])
519                infeas[iSequence+addSequence] = value*value; // already there
520              else
521                infeasible_->quickAdd(iSequence+addSequence,value*value);
522            }
523          }
524          break;
525        case ClpSimplex::atLowerBound:
526          if (value<-tolerance) {
527            // store square in list
528            if (infeas[iSequence+addSequence])
529              infeas[iSequence+addSequence] = value*value; // already there
530            else
531              infeasible_->quickAdd(iSequence+addSequence,value*value);
532          } else {
533            // look other way - change down should be positive
534            value -= nonLinear->changeDownInCost(iSequence+addSequence);
535            if (value<tolerance) {
536              infeasible_->zero(iSequence+addSequence);
537            } else {
538              // store square in list
539              if (infeas[iSequence+addSequence])
540                infeas[iSequence+addSequence] = value*value; // already there
541              else
542                infeasible_->quickAdd(iSequence+addSequence,value*value);
543            }
544          }
545        }
546      }
547    }
548  }
549#endif
550  // for weights update we use pivotSequence
551  pivotRow = pivotSequence_;
552  // unset in case sub flip
553  pivotSequence_=-1;
554  if (pivotRow>=0) {
555    // make sure infeasibility on incoming is 0.0
556    const int * pivotVariable = model_->pivotVariable();
557    int sequenceIn = pivotVariable[pivotRow];
558    infeasible_->zero(sequenceIn);
559    // and we can see if reference
560    double referenceIn=0.0;
561    if (!mode_&&reference(sequenceIn))
562      referenceIn=1.0;
563    // save outgoing weight round update
564    double outgoingWeight=weights_[sequenceOut];
565    // update weights
566    if (anyUpdates!=1) {
567      updates->setNumElements(0);
568      spareColumn1->setNumElements(0);
569      // might as well set dj to 1
570      dj = 1.0;
571      updates->insert(pivotRow,-dj);
572      model_->factorization()->updateColumnTranspose(spareRow2,updates);
573      // put row of tableau in rowArray and columnArray
574      model_->clpMatrix()->transposeTimes(model_,-1.0,
575                                          updates,spareColumn2,spareColumn1);
576    }
577    // now update weight update array
578    model_->factorization()->updateColumnTranspose(spareRow2,
579                                                   alternateWeights_);
580    double * other = alternateWeights_->denseVector();
581    double * weight;
582    int numberColumns = model_->numberColumns();
583    double scaleFactor = -1.0/dj; // as formula is with 1.0
584    // rows
585    number = updates->getNumElements();
586    index = updates->getIndices();
587    updateBy = updates->denseVector();
588    weight = weights_+numberColumns;
589   
590    for (j=0;j<number;j++) {
591      int iSequence = index[j];
592      double thisWeight = weight[iSequence];
593      // row has -1
594      double pivot = updateBy[iSequence]*scaleFactor;
595      updateBy[iSequence]=0.0;
596      double modification = other[iSequence];
597      double pivotSquared = pivot * pivot;
598     
599      thisWeight += pivotSquared * devex_ + pivot * modification;
600      if (thisWeight<TRY_NORM) {
601        if (mode_) {
602          // steepest
603          thisWeight = max(TRY_NORM,ADD_ONE+pivotSquared);
604        } else {
605          // exact
606          thisWeight = referenceIn*pivotSquared;
607          if (reference(iSequence+numberColumns))
608            thisWeight += 1.0;
609          thisWeight = max(thisWeight,TRY_NORM);
610        }
611      }
612      weight[iSequence] = thisWeight;
613    }
614   
615    // columns
616    weight = weights_;
617   
618    scaleFactor = -scaleFactor;
619   
620    number = spareColumn1->getNumElements();
621    index = spareColumn1->getIndices();
622    updateBy = spareColumn1->denseVector();
623    // get subset which have nonzero tableau elements
624    model_->clpMatrix()->subsetTransposeTimes(model_,alternateWeights_,
625                                              spareColumn1,
626                                              spareColumn2);
627    double * updateBy2 = spareColumn2->denseVector();
628    for (j=0;j<number;j++) {
629      int iSequence = index[j];
630      double thisWeight = weight[iSequence];
631      double pivot = updateBy[iSequence]*scaleFactor;
632      updateBy[iSequence]=0.0;
633      double modification = updateBy2[iSequence];
634      updateBy2[iSequence]=0.0;
635      double pivotSquared = pivot * pivot;
636     
637      thisWeight += pivotSquared * devex_ + pivot * modification;
638      if (thisWeight<TRY_NORM) {
639        if (mode_) {
640          // steepest
641          thisWeight = max(TRY_NORM,ADD_ONE+pivotSquared);
642        } else {
643          // exact
644          thisWeight = referenceIn*pivotSquared;
645          if (reference(iSequence))
646            thisWeight += 1.0;
647          thisWeight = max(thisWeight,TRY_NORM);
648        }
649      }
650      weight[iSequence] = thisWeight;
651    }
652    // restore outgoing weight
653    weights_[sequenceOut]=outgoingWeight;
654    alternateWeights_->clear();
655    spareColumn2->setNumElements(0);
656#ifdef SOME_DEBUG_1
657    // check for accuracy
658    int iCheck=229;
659    printf("weight for iCheck is %g\n",weights_[iCheck]);
660    if (model_->getStatus(iCheck)!=ClpSimplex::basic&&
661        !model_->getStatus(iCheck)!=isFixed)
662      checkAccuracy(iCheck,1.0e-1,updates,spareRow2);
663#endif
664    updates->setNumElements(0);
665    spareColumn1->setNumElements(0);
666  }
667
668  // update of duals finished - now do pricing
669
670
671  double bestDj = 1.0e-30;
672  int bestSequence=-1;
673
674  int i,iSequence;
675  index = infeasible_->getIndices();
676  number = infeasible_->getNumElements();
677  if(model_->numberIterations()<model_->lastBadIteration()+200) {
678    // we can't really trust infeasibilities if there is dual error
679    double checkTolerance = 1.0e-8;
680    if (!model_->factorization()->pivots())
681      checkTolerance = 1.0e-6;
682    if (model_->largestDualError()>checkTolerance)
683      tolerance *= model_->largestDualError()/checkTolerance;
684    // But cap
685    tolerance = min(1000.0,tolerance);
686  }
687#ifdef CLP_DEBUG
688  if (model_->numberDualInfeasibilities()==1) 
689    printf("** %g %g %g %x %x %d\n",tolerance,model_->dualTolerance(),
690           model_->largestDualError(),model_,model_->messageHandler(),
691           number);
692#endif
693  // stop last one coming immediately
694  double saveOutInfeasibility=0.0;
695  if (sequenceOut>=0) {
696    saveOutInfeasibility = infeas[sequenceOut];
697    infeas[sequenceOut]=0.0;
698  }
699  tolerance *= tolerance; // as we are using squares
700  for (i=0;i<number;i++) {
701    iSequence = index[i];
702    double value = infeas[iSequence];
703    double weight = weights_[iSequence];
704    //weight=1.0;
705    if (value>bestDj*weight&&value>tolerance) {
706      // check flagged variable
707      if (!model_->flagged(iSequence)) {
708        /*if (model_->numberIterations()%100==0)
709          printf("chosen %g => %g\n",bestDj,value/weight);*/
710        bestDj=value/weight;
711        bestSequence = iSequence;
712      }
713    }
714  }
715  if (sequenceOut>=0) {
716    infeas[sequenceOut]=saveOutInfeasibility;
717  }
718  /*if (model_->numberIterations()%100==0)
719    printf("%d best %g\n",bestSequence,bestDj);*/
720 
721#ifdef CLP_DEBUG
722  if (bestSequence>=0) {
723    if (model_->getStatus(bestSequence)==ClpSimplex::atLowerBound)
724      assert(model_->reducedCost(bestSequence)<0.0);
725    if (model_->getStatus(bestSequence)==ClpSimplex::atUpperBound)
726      assert(model_->reducedCost(bestSequence)>0.0);
727  }
728#endif
729  return bestSequence;
730}
731/*
732   1) before factorization
733   2) after factorization
734   3) just redo infeasibilities
735   4) restore weights
736*/
737void 
738ClpPrimalColumnSteepest::saveWeights(ClpSimplex * model,int mode)
739{
740  // alternateWeights_ is defined as indexed but is treated oddly
741  // at times
742  model_ = model;
743  int numberRows = model_->numberRows();
744  int numberColumns = model_->numberColumns();
745  const int * pivotVariable = model_->pivotVariable();
746  if (mode==1&&weights_) {
747    if (pivotSequence_>=0) {
748      // save pivot order
749      memcpy(alternateWeights_->getIndices(),pivotVariable,
750             numberRows*sizeof(int));
751      // change from pivot row number to sequence number
752      pivotSequence_=pivotVariable[pivotSequence_];
753    } else {
754      alternateWeights_->clear();
755    }
756    state_=1;
757  } else if (mode==2||mode==4||mode==5) {
758    // restore
759    if (!weights_||state_==-1||mode==5) {
760      // initialize weights
761      delete [] weights_;
762      delete alternateWeights_;
763      weights_ = new double[numberRows+numberColumns];
764      alternateWeights_ = new CoinIndexedVector();
765      // enough space so can use it for factorization
766      alternateWeights_->reserve(numberRows+
767                                 model_->factorization()->maximumPivots());
768      initializeWeights();
769      // create saved weights
770      delete [] savedWeights_;
771      savedWeights_ = new double[numberRows+numberColumns];
772      memcpy(savedWeights_,weights_,(numberRows+numberColumns)*
773             sizeof(double));
774      savedPivotSequence_=-2;
775      savedSequenceOut_ = -2;
776     
777    } else {
778      if (mode!=4) {
779        // save
780        memcpy(savedWeights_,weights_,(numberRows+numberColumns)*
781               sizeof(double));
782        savedPivotSequence_= pivotSequence_;
783        savedSequenceOut_ = model_->sequenceOut();
784      } else {
785        // restore
786        memcpy(weights_,savedWeights_,(numberRows+numberColumns)*
787               sizeof(double));
788        // was - but I think should not be
789        //pivotSequence_= savedPivotSequence_;
790        //model_->setSequenceOut(savedSequenceOut_);
791        pivotSequence_= -1;
792        model_->setSequenceOut(-1); 
793      }
794    }
795    state_=0;
796    // set up infeasibilities
797    if (!infeasible_) {
798      infeasible_ = new CoinIndexedVector();
799      infeasible_->reserve(numberColumns+numberRows);
800    }
801  }
802  if (mode>=2&&mode!=5) {
803    if (mode!=3&&pivotSequence_>=0) {
804      // restore pivot row
805      int iRow;
806      // permute alternateWeights
807      double * temp = new double[numberRows+numberColumns];
808      double * work = alternateWeights_->denseVector();
809      int * oldPivotOrder = alternateWeights_->getIndices();
810      for (iRow=0;iRow<numberRows;iRow++) {
811        int iPivot=oldPivotOrder[iRow];
812        temp[iPivot]=work[iRow];
813      }
814      int number=0;
815      int found=-1;
816      int * which = oldPivotOrder;
817      // find pivot row and re-create alternateWeights
818      for (iRow=0;iRow<numberRows;iRow++) {
819        int iPivot=pivotVariable[iRow];
820        if (iPivot==pivotSequence_) 
821          found=iRow;
822        work[iRow]=temp[iPivot];
823        if (work[iRow])
824          which[number++]=iRow;
825      }
826      alternateWeights_->setNumElements(number);
827#ifdef CLP_DEBUG
828      // Can happen but I should clean up
829      assert(found>=0);
830#endif
831      pivotSequence_ = found;
832      delete [] temp;
833    }
834    infeasible_->clear();
835    double tolerance=model_->currentDualTolerance();
836    int number = model_->numberRows() + model_->numberColumns();
837    int iSequence;
838   
839    double * reducedCost = model_->djRegion();
840     
841    if (!model_->nonLinearCost()->lookBothWays()) {
842      for (iSequence=0;iSequence<number;iSequence++) {
843        double value = reducedCost[iSequence];
844        ClpSimplex::Status status = model_->getStatus(iSequence);
845
846        switch(status) {
847         
848        case ClpSimplex::basic:
849        case ClpSimplex::isFixed:
850          break;
851        case ClpSimplex::isFree:
852        case ClpSimplex::superBasic:
853          if (fabs(value)>FREE_ACCEPT*tolerance) { 
854            // we are going to bias towards free (but only if reasonable)
855            value *= FREE_BIAS;
856            // store square in list
857            infeasible_->quickAdd(iSequence,value*value);
858          }
859          break;
860        case ClpSimplex::atUpperBound:
861          if (value>tolerance) {
862            infeasible_->quickAdd(iSequence,value*value);
863          }
864          break;
865        case ClpSimplex::atLowerBound:
866          if (value<-tolerance) {
867            infeasible_->quickAdd(iSequence,value*value);
868          }
869        }
870      }
871    } else {
872      ClpNonLinearCost * nonLinear = model_->nonLinearCost(); 
873      // can go both ways
874      for (iSequence=0;iSequence<number;iSequence++) {
875        double value = reducedCost[iSequence];
876        ClpSimplex::Status status = model_->getStatus(iSequence);
877
878        switch(status) {
879         
880        case ClpSimplex::basic:
881        case ClpSimplex::isFixed:
882          break;
883        case ClpSimplex::isFree:
884        case ClpSimplex::superBasic:
885          if (fabs(value)>FREE_ACCEPT*tolerance) { 
886            // we are going to bias towards free (but only if reasonable)
887            value *= FREE_BIAS;
888            // store square in list
889            infeasible_->quickAdd(iSequence,value*value);
890          }
891          break;
892        case ClpSimplex::atUpperBound:
893          if (value>tolerance) {
894            infeasible_->quickAdd(iSequence,value*value);
895          } else {
896            // look other way - change up should be negative
897            value -= nonLinear->changeUpInCost(iSequence);
898            if (value<-tolerance) {
899              // store square in list
900              infeasible_->quickAdd(iSequence,value*value);
901            }
902          }
903          break;
904        case ClpSimplex::atLowerBound:
905          if (value<-tolerance) {
906            infeasible_->quickAdd(iSequence,value*value);
907          } else {
908            // look other way - change down should be positive
909            value -= nonLinear->changeDownInCost(iSequence);
910            if (value>tolerance) {
911              // store square in list
912              infeasible_->quickAdd(iSequence,value*value);
913            }
914          }
915        }
916      }
917    }
918  }
919}
920// Gets rid of last update
921void 
922ClpPrimalColumnSteepest::unrollWeights()
923{
924  double * saved = alternateWeights_->denseVector();
925  int number = alternateWeights_->getNumElements();
926  int * which = alternateWeights_->getIndices();
927  int i;
928  for (i=0;i<number;i++) {
929    int iRow = which[i];
930    weights_[iRow]=saved[iRow];
931    saved[iRow]=0.0;
932  }
933  alternateWeights_->setNumElements(0);
934}
935 
936//-------------------------------------------------------------------
937// Clone
938//-------------------------------------------------------------------
939ClpPrimalColumnPivot * ClpPrimalColumnSteepest::clone(bool CopyData) const
940{
941  if (CopyData) {
942    return new ClpPrimalColumnSteepest(*this);
943  } else {
944    return new ClpPrimalColumnSteepest();
945  }
946}
947void
948ClpPrimalColumnSteepest::updateWeights(CoinIndexedVector * input)
949{
950  int number=input->getNumElements();
951  int * which = input->getIndices();
952  double * work = input->denseVector();
953  int newNumber = 0;
954  int * newWhich = alternateWeights_->getIndices();
955  double * newWork = alternateWeights_->denseVector();
956  int i;
957  int sequenceIn = model_->sequenceIn();
958  int sequenceOut = model_->sequenceOut();
959  const int * pivotVariable = model_->pivotVariable();
960
961  int pivotRow = model_->pivotRow();
962  pivotSequence_ = pivotRow;
963
964  devex_ =0.0;
965
966  if (pivotRow>=0) {
967    if (mode_) {
968      for (i=0;i<number;i++) {
969        int iRow = which[i];
970        devex_ += work[iRow]*work[iRow];
971        newWork[iRow] = -2.0*work[iRow];
972      }
973      newWork[pivotRow] = -2.0*max(devex_,0.0);
974      devex_+=ADD_ONE;
975      weights_[sequenceOut]=1.0+ADD_ONE;
976      memcpy(newWhich,which,number*sizeof(int));
977      alternateWeights_->setNumElements(number);
978    } else {
979      for (i=0;i<number;i++) {
980        int iRow = which[i];
981        int iPivot = pivotVariable[iRow];
982        if (reference(iPivot)) {
983          devex_ += work[iRow]*work[iRow];
984          newWork[iRow] = -2.0*work[iRow];
985          newWhich[newNumber++]=iRow;
986        }
987      }
988      if (!newWork[pivotRow])
989        newWhich[newNumber++]=pivotRow; // add if not already in
990      newWork[pivotRow] = -2.0*max(devex_,0.0);
991      if (reference(sequenceIn)) {
992        devex_+=1.0;
993      } else {
994      }
995      if (reference(sequenceOut)) {
996        weights_[sequenceOut]=1.0+1.0;
997      } else {
998        weights_[sequenceOut]=1.0;
999      }
1000      alternateWeights_->setNumElements(newNumber);
1001    }
1002  } else {
1003    if (mode_) {
1004      for (i=0;i<number;i++) {
1005        int iRow = which[i];
1006        devex_ += work[iRow]*work[iRow];
1007      }
1008      devex_ += ADD_ONE;
1009    } else {
1010      for (i=0;i<number;i++) {
1011        int iRow = which[i];
1012        int iPivot = pivotVariable[iRow];
1013        if (reference(iPivot)) {
1014          devex_ += work[iRow]*work[iRow];
1015        }
1016      }
1017      if (reference(sequenceIn)) 
1018        devex_+=1.0;
1019    }
1020  }
1021  double oldDevex = weights_[sequenceIn];
1022#ifdef CLP_DEBUG
1023  if ((model_->messageHandler()->logLevel()&32))
1024    printf("old weight %g, new %g\n",oldDevex,devex_);
1025#endif
1026  double check = max(devex_,oldDevex);;
1027  weights_[sequenceIn] = devex_;
1028  if ( fabs ( devex_ - oldDevex ) > 1.0e-1 * check ) {
1029#ifdef CLP_DEBUG
1030    if ((model_->messageHandler()->logLevel()&48)==16)
1031      printf("old weight %g, new %g\n",oldDevex,devex_);
1032#endif
1033    if ( fabs ( devex_ - oldDevex ) > 1.0e5 * check ) {
1034      // need to redo
1035      model_->messageHandler()->message(CLP_INITIALIZE_STEEP,
1036                                        *model_->messagesPointer())
1037                                          <<oldDevex<<devex_
1038                                          <<CoinMessageEol;
1039      initializeWeights();
1040    }
1041  }
1042  if (pivotRow>=0) {
1043    // set outgoing weight here
1044    weights_[model_->sequenceOut()]=devex_/(model_->alpha()*model_->alpha());
1045  }
1046}
1047// Checks accuracy - just for debug
1048void 
1049ClpPrimalColumnSteepest::checkAccuracy(int sequence,
1050                                       double relativeTolerance,
1051                                       CoinIndexedVector * rowArray1,
1052                                       CoinIndexedVector * rowArray2)
1053{
1054  model_->unpack(rowArray1,sequence);
1055  model_->factorization()->updateColumn(rowArray2,rowArray1);
1056  int number=rowArray1->getNumElements();
1057  int * which = rowArray1->getIndices();
1058  double * work = rowArray1->denseVector();
1059  const int * pivotVariable = model_->pivotVariable();
1060 
1061  double devex =0.0;
1062  int i;
1063
1064  if (mode_) {
1065    for (i=0;i<number;i++) {
1066      int iRow = which[i];
1067      devex += work[iRow]*work[iRow];
1068      work[iRow]=0.0;
1069    }
1070    devex += ADD_ONE;
1071  } else {
1072    for (i=0;i<number;i++) {
1073      int iRow = which[i];
1074      int iPivot = pivotVariable[iRow];
1075      if (reference(iPivot)) {
1076        devex += work[iRow]*work[iRow];
1077      }
1078      work[iRow]=0.0;
1079    }
1080    if (reference(sequence)) 
1081      devex+=1.0;
1082  }
1083
1084  double oldDevex = weights_[sequence];
1085  double check = max(devex,oldDevex);;
1086  if ( fabs ( devex - oldDevex ) > relativeTolerance * check ) {
1087    printf("check %d old weight %g, new %g\n",sequence,oldDevex,devex);
1088    // update so won't print again
1089    weights_[sequence]=devex;
1090  }
1091  rowArray1->setNumElements(0);
1092}
1093
1094// Initialize weights
1095void 
1096ClpPrimalColumnSteepest::initializeWeights()
1097{
1098  int numberRows = model_->numberRows();
1099  int numberColumns = model_->numberColumns();
1100  int number = numberRows + numberColumns;
1101  int iSequence;
1102  if (!mode_) {
1103    // initialize to 1.0
1104    // and set reference framework
1105    if (!reference_) {
1106      int nWords = (number+31)>>5;
1107      reference_ = new unsigned int[nWords];
1108      assert (sizeof(unsigned int)==4);
1109      // tiny overhead to zero out (stops valgrind complaining)
1110      memset(reference_,0,nWords*sizeof(int));
1111    }
1112   
1113    for (iSequence=0;iSequence<number;iSequence++) {
1114      weights_[iSequence]=1.0;
1115      if (model_->getStatus(iSequence)==ClpSimplex::basic) {
1116        setReference(iSequence,false);
1117      } else {
1118        setReference(iSequence,true);
1119      }
1120    }
1121  } else {
1122    CoinIndexedVector * temp = new CoinIndexedVector();
1123    temp->reserve(numberRows+
1124                  model_->factorization()->maximumPivots());
1125    double * array = alternateWeights_->denseVector();
1126    int * which = alternateWeights_->getIndices();
1127     
1128    for (iSequence=0;iSequence<number;iSequence++) {
1129      weights_[iSequence]=1.0+ADD_ONE;
1130      if (model_->getStatus(iSequence)!=ClpSimplex::basic&&
1131          model_->getStatus(iSequence)!=ClpSimplex::isFixed) {
1132        model_->unpack(alternateWeights_,iSequence);
1133        double value=ADD_ONE;
1134        model_->factorization()->updateColumn(temp,alternateWeights_);
1135        int number = alternateWeights_->getNumElements();       int j;
1136        for (j=0;j<number;j++) {
1137          int iRow=which[j];
1138          value+=array[iRow]*array[iRow];
1139          array[iRow]=0.0;
1140        }
1141        alternateWeights_->setNumElements(0);
1142        weights_[iSequence] = value;
1143      }
1144    }
1145    delete temp;
1146  }
1147}
1148// Gets rid of all arrays
1149void 
1150ClpPrimalColumnSteepest::clearArrays()
1151{
1152  delete [] weights_;
1153  weights_=NULL;
1154  delete infeasible_;
1155  infeasible_ = NULL;
1156  delete alternateWeights_;
1157  alternateWeights_ = NULL;
1158  delete [] savedWeights_;
1159  savedWeights_ = NULL;
1160  delete [] reference_;
1161  reference_ = NULL;
1162  pivotSequence_=-1;
1163  state_ = -1;
1164  savedPivotSequence_ = -1;
1165  savedSequenceOut_ = -1; 
1166  devex_ = 0.0;
1167}
1168// Returns true if would not find any column
1169bool 
1170ClpPrimalColumnSteepest::looksOptimal() const
1171{
1172  //**** THIS MUST MATCH the action coding above
1173  double tolerance=model_->currentDualTolerance();
1174  // we can't really trust infeasibilities if there is dual error
1175  // this coding has to mimic coding in checkDualSolution
1176  double error = min(1.0e-3,model_->largestDualError());
1177  // allow tolerance at least slightly bigger than standard
1178  tolerance = tolerance +  error;
1179  if(model_->numberIterations()<model_->lastBadIteration()+200) {
1180    // we can't really trust infeasibilities if there is dual error
1181    double checkTolerance = 1.0e-8;
1182    if (!model_->factorization()->pivots())
1183      checkTolerance = 1.0e-6;
1184    if (model_->largestDualError()>checkTolerance)
1185      tolerance *= model_->largestDualError()/checkTolerance;
1186    // But cap
1187    tolerance = min(1000.0,tolerance);
1188  }
1189  int number = model_->numberRows() + model_->numberColumns();
1190  int iSequence;
1191 
1192  double * reducedCost = model_->djRegion();
1193  int numberInfeasible=0;
1194  if (!model_->nonLinearCost()->lookBothWays()) {
1195    for (iSequence=0;iSequence<number;iSequence++) {
1196      double value = reducedCost[iSequence];
1197      ClpSimplex::Status status = model_->getStatus(iSequence);
1198     
1199      switch(status) {
1200
1201      case ClpSimplex::basic:
1202      case ClpSimplex::isFixed:
1203        break;
1204      case ClpSimplex::isFree:
1205      case ClpSimplex::superBasic:
1206        if (fabs(value)>FREE_ACCEPT*tolerance) 
1207          numberInfeasible++;
1208        break;
1209      case ClpSimplex::atUpperBound:
1210        if (value>tolerance) 
1211          numberInfeasible++;
1212        break;
1213      case ClpSimplex::atLowerBound:
1214        if (value<-tolerance) 
1215          numberInfeasible++;
1216      }
1217    }
1218  } else {
1219    ClpNonLinearCost * nonLinear = model_->nonLinearCost(); 
1220    // can go both ways
1221    for (iSequence=0;iSequence<number;iSequence++) {
1222      double value = reducedCost[iSequence];
1223      ClpSimplex::Status status = model_->getStatus(iSequence);
1224     
1225      switch(status) {
1226
1227      case ClpSimplex::basic:
1228      case ClpSimplex::isFixed:
1229        break;
1230      case ClpSimplex::isFree:
1231      case ClpSimplex::superBasic:
1232        if (fabs(value)>FREE_ACCEPT*tolerance) 
1233          numberInfeasible++;
1234        break;
1235      case ClpSimplex::atUpperBound:
1236        if (value>tolerance) {
1237          numberInfeasible++;
1238        } else {
1239          // look other way - change up should be negative
1240          value -= nonLinear->changeUpInCost(iSequence);
1241          if (value<-tolerance) 
1242            numberInfeasible++;
1243        }
1244        break;
1245      case ClpSimplex::atLowerBound:
1246        if (value<-tolerance) {
1247          numberInfeasible++;
1248        } else {
1249          // look other way - change down should be positive
1250          value -= nonLinear->changeDownInCost(iSequence);
1251          if (value>tolerance) 
1252            numberInfeasible++;
1253        }
1254      }
1255    }
1256  }
1257  return numberInfeasible==0;
1258}
Note: See TracBrowser for help on using the repository browser.