source: trunk/ClpPrimalColumnSteepest.cpp @ 162

Last change on this file since 162 was 162, checked in by forrest, 17 years ago

Piecewise linear

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