source: trunk/ClpDualRowSteepest.cpp @ 97

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

Take out references to stopQuickadd (as CoinIndexedVector? changed)

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 15.9 KB
Line 
1// Copyright (C) 2002, International Business Machines
2// Corporation and others.  All Rights Reserved.
3
4#include "CoinPragma.hpp"
5#include "ClpSimplex.hpp"
6#include "ClpDualRowSteepest.hpp"
7#include "CoinIndexedVector.hpp"
8#include "ClpFactorization.hpp"
9#include "CoinHelperFunctions.hpp"
10#include <cstdio>
11
12//#############################################################################
13// Constructors / Destructor / Assignment
14//#############################################################################
15
16//-------------------------------------------------------------------
17// Default Constructor
18//-------------------------------------------------------------------
19ClpDualRowSteepest::ClpDualRowSteepest (int mode) 
20  : ClpDualRowPivot(),
21    state_(-1),
22    mode_(mode),
23    weights_(NULL),
24    infeasible_(NULL),
25    alternateWeights_(NULL),
26    savedWeights_(NULL)
27{
28  type_=2+64*mode;
29}
30
31//-------------------------------------------------------------------
32// Copy constructor
33//-------------------------------------------------------------------
34ClpDualRowSteepest::ClpDualRowSteepest (const ClpDualRowSteepest & rhs) 
35: ClpDualRowPivot(rhs)
36{ 
37  state_=rhs.state_;
38  mode_ = rhs.mode_;
39  model_ = rhs.model_;
40  if (rhs.infeasible_) {
41    infeasible_= new CoinIndexedVector(rhs.infeasible_);
42  } else {
43    infeasible_=NULL;
44  }
45  if (rhs.weights_) {
46    assert(model_);
47    int number = model_->numberRows();
48    weights_= new double[number];
49    ClpDisjointCopyN(rhs.weights_,number,weights_);
50  } else {
51    weights_=NULL;
52  }
53  if (rhs.alternateWeights_) {
54    alternateWeights_= new CoinIndexedVector(rhs.alternateWeights_);
55  } else {
56    alternateWeights_=NULL;
57  }
58  if (rhs.savedWeights_) {
59    savedWeights_= new CoinIndexedVector(rhs.savedWeights_);
60  } else {
61    savedWeights_=NULL;
62  }
63}
64
65//-------------------------------------------------------------------
66// Destructor
67//-------------------------------------------------------------------
68ClpDualRowSteepest::~ClpDualRowSteepest ()
69{
70  delete [] weights_;
71  delete infeasible_;
72  delete alternateWeights_;
73  delete savedWeights_;
74}
75
76//----------------------------------------------------------------
77// Assignment operator
78//-------------------------------------------------------------------
79ClpDualRowSteepest &
80ClpDualRowSteepest::operator=(const ClpDualRowSteepest& rhs)
81{
82  if (this != &rhs) {
83    ClpDualRowPivot::operator=(rhs);
84    state_=rhs.state_;
85    mode_ = rhs.mode_;
86    model_ = rhs.model_;
87    delete [] weights_;
88    delete infeasible_;
89    delete alternateWeights_;
90    delete savedWeights_;
91    if (rhs.infeasible_!=NULL) {
92      infeasible_= new CoinIndexedVector(rhs.infeasible_);
93    } else {
94      infeasible_=NULL;
95    }
96    if (rhs.weights_!=NULL) {
97      assert(model_);
98      int number = model_->numberRows();
99      weights_= new double[number];
100      ClpDisjointCopyN(rhs.weights_,number,weights_);
101    } else {
102      weights_=NULL;
103    }
104    if (rhs.alternateWeights_!=NULL) {
105      alternateWeights_= new CoinIndexedVector(rhs.alternateWeights_);
106    } else {
107      alternateWeights_=NULL;
108    }
109    if (rhs.savedWeights_!=NULL) {
110      savedWeights_= new CoinIndexedVector(rhs.savedWeights_);
111    } else {
112      savedWeights_=NULL;
113    }
114  }
115  return *this;
116}
117
118// Returns pivot row, -1 if none
119int 
120ClpDualRowSteepest::pivotRow()
121{
122  assert(model_);
123  int i,iRow;
124  double * infeas = infeasible_->denseVector();
125  double largest=1.0e-50;
126  int * index = infeasible_->getIndices();
127  int number = infeasible_->getNumElements();
128  const int * pivotVariable =model_->pivotVariable();
129  int chosenRow=-1;
130  int lastPivotRow = model_->pivotRow();
131  double tolerance=model_->currentPrimalTolerance();
132  // we can't really trust infeasibilities if there is primal error
133  // this coding has to mimic coding in checkPrimalSolution
134  double error = min(1.0e-3,model_->largestPrimalError());
135  // allow tolerance at least slightly bigger than standard
136  tolerance = tolerance +  error;
137  tolerance *= tolerance; // as we are using squares
138  double * solution = model_->solutionRegion();
139  double * lower = model_->lowerRegion();
140  double * upper = model_->upperRegion();
141  // do last pivot row one here
142  if (lastPivotRow>=0) {
143    int iPivot=pivotVariable[lastPivotRow];
144    double value = solution[iPivot];
145    double lower = model_->lower(iPivot);
146    double upper = model_->upper(iPivot);
147    if (value>upper+tolerance) {
148      // store square in list
149      if (infeas[lastPivotRow])
150        infeas[lastPivotRow] = (value-upper)*(value-upper); // already there
151      else
152        infeasible_->quickAdd(lastPivotRow,(value-upper)*(value-upper));
153    } else if (value<lower-tolerance) {
154      // store square in list
155      if (infeas[lastPivotRow])
156        infeas[lastPivotRow] = (value-lower)*(value-lower); // already there
157      else
158        infeasible_->add(lastPivotRow,(value-lower)*(value-lower));
159    } else {
160      // feasible - was it infeasible - if so set tiny
161      if (infeas[lastPivotRow])
162        infeas[lastPivotRow] = 1.0e-100;
163    }
164    number = infeasible_->getNumElements();
165  }
166  for (i=0;i<number;i++) {
167    iRow = index[i];
168    double value = infeas[iRow];
169    if (value>largest*weights_[iRow]&&value>tolerance) {
170      // make last pivot row last resort choice
171      if (iRow==lastPivotRow) {
172        if (value*1.0e-10<largest*weights_[iRow]) 
173          continue;
174        else 
175          value *= 1.0e-10;
176      }
177      int iSequence = pivotVariable[iRow];
178      if (!model_->flagged(iSequence)) {
179        //#define CLP_DEBUG 1
180#ifdef CLP_DEBUG
181        double value2=0.0;
182        if (solution[iSequence]>upper[iSequence]+tolerance) 
183          value2=solution[iSequence]-upper[iSequence];
184        else if (solution[iSequence]<lower[iSequence]-tolerance) 
185          value2=solution[iSequence]-lower[iSequence];
186        assert(fabs(value2*value2-infeas[iRow])<1.0e-8*min(value2*value2,infeas[iRow]));
187#endif
188        if (solution[iSequence]>upper[iSequence]+tolerance||
189            solution[iSequence]<lower[iSequence]-tolerance) {
190          chosenRow=iRow;
191          largest=value/weights_[iRow];
192        }
193      }
194    }
195  }
196  return chosenRow;
197}
198#define TRY_NORM 1.0e-4
199// Updates weights
200void 
201ClpDualRowSteepest::updateWeights(CoinIndexedVector * input,
202                                  CoinIndexedVector * spare,
203                                  CoinIndexedVector * updatedColumn)
204{
205  // clear other region
206  alternateWeights_->clear();
207  double norm = 0.0;
208  double * work = input->denseVector();
209  int number = input->getNumElements();
210  int * which = input->getIndices();
211  double * work2 = spare->denseVector();
212  int * which2 = spare->getIndices();
213  double * work3 = alternateWeights_->denseVector();
214  int * which3 = alternateWeights_->getIndices();
215  int i;
216#if CLP_DEBUG>2
217  // Very expensive debug
218  {
219    int numberRows = model_->numberRows();
220    CoinIndexedVector * temp = new CoinIndexedVector();
221    temp->reserve(numberRows+
222                  model_->factorization()->maximumPivots());
223    double * array = alternateWeights_->denseVector();
224    int * which = alternateWeights_->getIndices();
225    for (i=0;i<numberRows;i++) {
226      double value=0.0;
227      array[i] = 1.0;
228      which[0] = i;
229      alternateWeights_->setNumElements(1);
230      model_->factorization()->updateColumnTranspose(temp,
231                                                     alternateWeights_);
232      int number = alternateWeights_->getNumElements();
233      int j;
234      for (j=0;j<number;j++) {
235        int iRow=which[j];
236        value+=array[iRow]*array[iRow];
237        array[iRow]=0.0;
238      }
239      alternateWeights_->setNumElements(0);
240      if (fabs(weights_[i]-value)>1.0e-4)
241        printf("%d old %g, true %g\n",i,weights_[i],value);
242    }
243    delete temp;
244  }
245#endif
246  for (i=0;i<number;i++) {
247    int iRow = which[i];
248    double value = work[iRow];
249    norm += value*value;
250    work2[iRow]=value;
251    which2[i]=iRow;
252  }
253  spare->setNumElements(number);
254  // ftran
255  model_->factorization()->updateColumn(alternateWeights_,spare);
256  // alternateWeights_ should still be empty
257  int pivotRow = model_->pivotRow();
258#ifdef CLP_DEBUG
259  if ( model_->logLevel (  ) >4  && 
260       fabs(norm-weights_[pivotRow])>1.0e-3*(1.0+norm)) 
261    printf("on row %d, true weight %g, old %g\n",
262           pivotRow,sqrt(norm),sqrt(weights_[pivotRow]));
263#endif
264  // could re-initialize here (could be expensive)
265  norm /= model_->alpha() * model_->alpha();
266
267  assert(norm);
268  if (norm < TRY_NORM) 
269    norm = TRY_NORM;
270  if (norm != 0.) {
271    double multiplier = 2.0 / model_->alpha();
272    // look at updated column
273    work = updatedColumn->denseVector();
274    number = updatedColumn->getNumElements();
275    which = updatedColumn->getIndices();
276
277    int nSave=0;
278
279    for (i =0; i < number; i++) {
280      int iRow = which[i];
281      double theta = work[iRow];
282      if (theta) {
283        double devex = weights_[iRow];
284        work3[iRow]=devex; // save old
285        which3[nSave++]=iRow;
286        double value = work2[iRow];
287        devex +=  theta * (theta*norm+value * multiplier);
288        if (devex < TRY_NORM) 
289          devex = TRY_NORM;
290        weights_[iRow]=devex;
291      }
292    }
293#ifdef CLP_DEBUG
294    assert(work3[pivotRow]&&work[pivotRow]);
295#endif
296    alternateWeights_->setNumElements(nSave);
297    if (norm < TRY_NORM) 
298      norm = TRY_NORM;
299    weights_[pivotRow] = norm;
300  }
301  spare->clear();
302#ifdef CLP_DEBUG
303  spare->checkClear();
304#endif
305}
306 
307/* Updates primal solution (and maybe list of candidates)
308   Uses input vector which it deletes
309   Computes change in objective function
310*/
311void 
312ClpDualRowSteepest::updatePrimalSolution(
313                                        CoinIndexedVector * primalUpdate,
314                                        double primalRatio,
315                                        double & objectiveChange)
316{
317  double * work = primalUpdate->denseVector();
318  int number = primalUpdate->getNumElements();
319  int * which = primalUpdate->getIndices();
320  int i;
321  double changeObj=0.0;
322  double tolerance=model_->currentPrimalTolerance();
323  const int * pivotVariable = model_->pivotVariable();
324  double * infeas = infeasible_->denseVector();
325  int pivotRow = model_->pivotRow();
326  double * solution = model_->solutionRegion();
327  for (i=0;i<number;i++) {
328    int iRow=which[i];
329    int iPivot=pivotVariable[iRow];
330    double value = solution[iPivot];
331    double cost = model_->cost(iPivot);
332    double change = primalRatio*work[iRow];
333    value -= change;
334    changeObj -= change*cost;
335    solution[iPivot] = value;
336    double lower = model_->lower(iPivot);
337    double upper = model_->upper(iPivot);
338    // But if pivot row then use value of incoming
339    // Although it is safer to recompute before next selection
340    // in case something odd happens
341    if (iRow==pivotRow) {
342      iPivot = model_->sequenceIn();
343      lower = model_->lower(iPivot);
344      upper = model_->upper(iPivot);
345      value = model_->valueIncomingDual();
346    }
347    if (value>upper+tolerance) {
348      // store square in list
349      if (infeas[iRow])
350        infeas[iRow] = (value-upper)*(value-upper); // already there
351      else
352        infeasible_->quickAdd(iRow,(value-upper)*(value-upper));
353    } else if (value<lower-tolerance) {
354      // store square in list
355      if (infeas[iRow])
356        infeas[iRow] = (value-lower)*(value-lower); // already there
357      else
358        infeasible_->quickAdd(iRow,(value-lower)*(value-lower));
359    } else {
360      // feasible - was it infeasible - if so set tiny
361      if (infeas[iRow])
362        infeas[iRow] = 1.0e-100;
363    }
364    work[iRow]=0.0;
365  }
366  primalUpdate->setNumElements(0);
367  objectiveChange += changeObj;
368}
369/* Saves any weights round factorization as pivot rows may change
370   1) before factorization
371   2) after factorization
372   3) just redo infeasibilities
373   4) restore weights
374*/
375void 
376ClpDualRowSteepest::saveWeights(ClpSimplex * model,int mode)
377{
378  // alternateWeights_ is defined as indexed but is treated oddly
379  model_ = model;
380  int numberRows = model_->numberRows();
381  int numberColumns = model_->numberColumns();
382  const int * pivotVariable = model_->pivotVariable();
383  int i;
384  if (mode==1&&weights_) {
385    alternateWeights_->clear();
386    // change from row numbers to sequence numbers
387    int * which = alternateWeights_->getIndices();
388    for (i=0;i<numberRows;i++) {
389      int iPivot=pivotVariable[i];
390      which[i]=iPivot;
391    }
392    state_=1;
393  } else if (mode==2||mode==4||mode==5) {
394    // restore
395    if (!weights_||state_==-1||mode==5) {
396      // initialize weights
397      delete [] weights_;
398      delete alternateWeights_;
399      weights_ = new double[numberRows];
400      alternateWeights_ = new CoinIndexedVector();
401      // enough space so can use it for factorization
402      alternateWeights_->reserve(numberRows+
403                                 model_->factorization()->maximumPivots());
404      if (!mode_||mode==5) {
405        // initialize to 1.0 (can we do better?)
406        for (i=0;i<numberRows;i++) {
407          weights_[i]=1.0;
408        }
409      } else {
410        CoinIndexedVector * temp = new CoinIndexedVector();
411        temp->reserve(numberRows+
412                      model_->factorization()->maximumPivots());
413        double * array = alternateWeights_->denseVector();
414        int * which = alternateWeights_->getIndices();
415        for (i=0;i<numberRows;i++) {
416          double value=0.0;
417          array[i] = 1.0;
418          which[0] = i;
419          alternateWeights_->setNumElements(1);
420          model_->factorization()->updateColumnTranspose(temp,
421                                                         alternateWeights_);
422          int number = alternateWeights_->getNumElements();
423          int j;
424          for (j=0;j<number;j++) {
425            int iRow=which[j];
426            value+=array[iRow]*array[iRow];
427            array[iRow]=0.0;
428          }
429          alternateWeights_->setNumElements(0);
430          weights_[i] = value;
431        }
432        delete temp;
433      }
434      // create saved weights (not really indexedvector)
435      savedWeights_ = new CoinIndexedVector();
436      savedWeights_->reserve(numberRows);
437     
438      double * array = savedWeights_->denseVector();
439      int * which = savedWeights_->getIndices();
440      for (i=0;i<numberRows;i++) {
441        array[i]=weights_[i];
442        which[i]=pivotVariable[i];
443      }
444    } else {
445      int * which = alternateWeights_->getIndices();
446      if (mode!=4) {
447        // save
448        memcpy(savedWeights_->getIndices(),which,
449               numberRows*sizeof(int));
450        memcpy(savedWeights_->denseVector(),weights_,
451               numberRows*sizeof(double));
452      } else {
453        // restore
454        memcpy(which,savedWeights_->getIndices(),
455               numberRows*sizeof(int));
456        memcpy(weights_,savedWeights_->denseVector(),
457               numberRows*sizeof(double));
458      }
459      // restore (a bit slow - but only every re-factorization)
460      double * array = new double[numberRows+numberColumns];
461      for (i=0;i<numberRows;i++) {
462        int iSeq=which[i];
463        array[iSeq]=weights_[i];
464      }
465      for (i=0;i<numberRows;i++) {
466        int iPivot=pivotVariable[i];
467        weights_[i]=array[iPivot];
468        if (weights_[i]<TRY_NORM)
469          weights_[i] = TRY_NORM; // may need to check more
470      }
471      delete [] array;
472    }
473    state_=0;
474    // set up infeasibilities
475    if (!infeasible_) {
476      infeasible_ = new CoinIndexedVector();
477      infeasible_->reserve(numberRows);
478    }
479  }
480  if (mode>=2) {
481    infeasible_->clear();
482    int iRow;
483    const int * pivotVariable = model_->pivotVariable();
484    double tolerance=model_->currentPrimalTolerance();
485    for (iRow=0;iRow<model_->numberRows();iRow++) {
486      int iPivot=pivotVariable[iRow];
487      double value = model_->solution(iPivot);
488      double lower = model_->lower(iPivot);
489      double upper = model_->upper(iPivot);
490      if (value>upper+tolerance) {
491        value -= upper;
492        // store square in list
493        infeasible_->quickAdd(iRow,value*value);
494      } else if (value<lower-tolerance) {
495        value -= lower;
496        // store square in list
497        infeasible_->quickAdd(iRow,value*value);
498      }
499    }
500  }
501}
502// Gets rid of last update
503void 
504ClpDualRowSteepest::unrollWeights()
505{
506  double * saved = alternateWeights_->denseVector();
507  int number = alternateWeights_->getNumElements();
508  int * which = alternateWeights_->getIndices();
509  int i;
510  for (i=0;i<number;i++) {
511    int iRow = which[i];
512    weights_[iRow]=saved[iRow];
513    saved[iRow]=0.0;
514  }
515  alternateWeights_->setNumElements(0);
516}
517//-------------------------------------------------------------------
518// Clone
519//-------------------------------------------------------------------
520ClpDualRowPivot * ClpDualRowSteepest::clone(bool CopyData) const
521{
522  if (CopyData) {
523    return new ClpDualRowSteepest(*this);
524  } else {
525    return new ClpDualRowSteepest();
526  }
527}
528// Gets rid of all arrays
529void 
530ClpDualRowSteepest::clearArrays()
531{
532  delete [] weights_;
533  weights_=NULL;
534  delete infeasible_;
535  infeasible_ = NULL;
536  delete alternateWeights_;
537  alternateWeights_ = NULL;
538  delete savedWeights_;
539  savedWeights_ = NULL;
540  state_ =-1;
541}
542
Note: See TracBrowser for help on using the repository browser.