source: trunk/Clp/src/ClpDualRowSteepest.cpp @ 1665

Last change on this file since 1665 was 1665, checked in by lou, 9 years ago

Add EPL license notice in src.

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 41.7 KB
Line 
1/* $Id: ClpDualRowSteepest.cpp 1665 2011-01-04 17:55:54Z lou $ */
2// Copyright (C) 2002, International Business Machines
3// Corporation and others.  All Rights Reserved.
4// This code is licensed under the terms of the Eclipse Public License (EPL).
5
6#include "CoinPragma.hpp"
7#include "ClpSimplex.hpp"
8#include "ClpDualRowSteepest.hpp"
9#include "CoinIndexedVector.hpp"
10#include "ClpFactorization.hpp"
11#include "CoinHelperFunctions.hpp"
12#include <cstdio>
13//#############################################################################
14// Constructors / Destructor / Assignment
15//#############################################################################
16//#define CLP_DEBUG 4
17//-------------------------------------------------------------------
18// Default Constructor
19//-------------------------------------------------------------------
20ClpDualRowSteepest::ClpDualRowSteepest (int mode)
21     : ClpDualRowPivot(),
22       state_(-1),
23       mode_(mode),
24       persistence_(normal),
25       weights_(NULL),
26       infeasible_(NULL),
27       alternateWeights_(NULL),
28       savedWeights_(NULL),
29       dubiousWeights_(NULL)
30{
31     type_ = 2 + 64 * mode;
32}
33
34//-------------------------------------------------------------------
35// Copy constructor
36//-------------------------------------------------------------------
37ClpDualRowSteepest::ClpDualRowSteepest (const ClpDualRowSteepest & rhs)
38     : ClpDualRowPivot(rhs)
39{
40     state_ = rhs.state_;
41     mode_ = rhs.mode_;
42     persistence_ = rhs.persistence_;
43     model_ = rhs.model_;
44     if ((model_ && model_->whatsChanged() & 1) != 0) {
45          int number = model_->numberRows();
46          if (rhs.savedWeights_)
47               number = CoinMin(number, rhs.savedWeights_->capacity());
48          if (rhs.infeasible_) {
49               infeasible_ = new CoinIndexedVector(rhs.infeasible_);
50          } else {
51               infeasible_ = NULL;
52          }
53          if (rhs.weights_) {
54               weights_ = new double[number];
55               ClpDisjointCopyN(rhs.weights_, number, weights_);
56          } else {
57               weights_ = NULL;
58          }
59          if (rhs.alternateWeights_) {
60               alternateWeights_ = new CoinIndexedVector(rhs.alternateWeights_);
61          } else {
62               alternateWeights_ = NULL;
63          }
64          if (rhs.savedWeights_) {
65               savedWeights_ = new CoinIndexedVector(rhs.savedWeights_);
66          } else {
67               savedWeights_ = NULL;
68          }
69          if (rhs.dubiousWeights_) {
70               assert(model_);
71               int number = model_->numberRows();
72               dubiousWeights_ = new int[number];
73               ClpDisjointCopyN(rhs.dubiousWeights_, number, dubiousWeights_);
74          } else {
75               dubiousWeights_ = NULL;
76          }
77     } else {
78          infeasible_ = NULL;
79          weights_ = NULL;
80          alternateWeights_ = NULL;
81          savedWeights_ = NULL;
82          dubiousWeights_ = NULL;
83     }
84}
85
86//-------------------------------------------------------------------
87// Destructor
88//-------------------------------------------------------------------
89ClpDualRowSteepest::~ClpDualRowSteepest ()
90{
91     delete [] weights_;
92     delete [] dubiousWeights_;
93     delete infeasible_;
94     delete alternateWeights_;
95     delete savedWeights_;
96}
97
98//----------------------------------------------------------------
99// Assignment operator
100//-------------------------------------------------------------------
101ClpDualRowSteepest &
102ClpDualRowSteepest::operator=(const ClpDualRowSteepest& rhs)
103{
104     if (this != &rhs) {
105          ClpDualRowPivot::operator=(rhs);
106          state_ = rhs.state_;
107          mode_ = rhs.mode_;
108          persistence_ = rhs.persistence_;
109          model_ = rhs.model_;
110          delete [] weights_;
111          delete [] dubiousWeights_;
112          delete infeasible_;
113          delete alternateWeights_;
114          delete savedWeights_;
115          assert(model_);
116          int number = model_->numberRows();
117          if (rhs.savedWeights_)
118               number = CoinMin(number, rhs.savedWeights_->capacity());
119          if (rhs.infeasible_ != NULL) {
120               infeasible_ = new CoinIndexedVector(rhs.infeasible_);
121          } else {
122               infeasible_ = NULL;
123          }
124          if (rhs.weights_ != NULL) {
125               weights_ = new double[number];
126               ClpDisjointCopyN(rhs.weights_, number, weights_);
127          } else {
128               weights_ = NULL;
129          }
130          if (rhs.alternateWeights_ != NULL) {
131               alternateWeights_ = new CoinIndexedVector(rhs.alternateWeights_);
132          } else {
133               alternateWeights_ = NULL;
134          }
135          if (rhs.savedWeights_ != NULL) {
136               savedWeights_ = new CoinIndexedVector(rhs.savedWeights_);
137          } else {
138               savedWeights_ = NULL;
139          }
140          if (rhs.dubiousWeights_) {
141               assert(model_);
142               int number = model_->numberRows();
143               dubiousWeights_ = new int[number];
144               ClpDisjointCopyN(rhs.dubiousWeights_, number, dubiousWeights_);
145          } else {
146               dubiousWeights_ = NULL;
147          }
148     }
149     return *this;
150}
151// Fill most values
152void
153ClpDualRowSteepest::fill(const ClpDualRowSteepest& rhs)
154{
155     state_ = rhs.state_;
156     mode_ = rhs.mode_;
157     persistence_ = rhs.persistence_;
158     assert (model_->numberRows() == rhs.model_->numberRows());
159     model_ = rhs.model_;
160     assert(model_);
161     int number = model_->numberRows();
162     if (rhs.savedWeights_)
163          number = CoinMin(number, rhs.savedWeights_->capacity());
164     if (rhs.infeasible_ != NULL) {
165          if (!infeasible_)
166               infeasible_ = new CoinIndexedVector(rhs.infeasible_);
167          else
168               *infeasible_ = *rhs.infeasible_;
169     } else {
170          delete infeasible_;
171          infeasible_ = NULL;
172     }
173     if (rhs.weights_ != NULL) {
174          if (!weights_)
175               weights_ = new double[number];
176          ClpDisjointCopyN(rhs.weights_, number, weights_);
177     } else {
178          delete [] weights_;
179          weights_ = NULL;
180     }
181     if (rhs.alternateWeights_ != NULL) {
182          if (!alternateWeights_)
183               alternateWeights_ = new CoinIndexedVector(rhs.alternateWeights_);
184          else
185               *alternateWeights_ = *rhs.alternateWeights_;
186     } else {
187          delete alternateWeights_;
188          alternateWeights_ = NULL;
189     }
190     if (rhs.savedWeights_ != NULL) {
191          if (!savedWeights_)
192               savedWeights_ = new CoinIndexedVector(rhs.savedWeights_);
193          else
194               *savedWeights_ = *rhs.savedWeights_;
195     } else {
196          delete savedWeights_;
197          savedWeights_ = NULL;
198     }
199     if (rhs.dubiousWeights_) {
200          assert(model_);
201          int number = model_->numberRows();
202          if (!dubiousWeights_)
203               dubiousWeights_ = new int[number];
204          ClpDisjointCopyN(rhs.dubiousWeights_, number, dubiousWeights_);
205     } else {
206          delete [] dubiousWeights_;
207          dubiousWeights_ = NULL;
208     }
209}
210// Returns pivot row, -1 if none
211int
212ClpDualRowSteepest::pivotRow()
213{
214     assert(model_);
215     int i, iRow;
216     double * infeas = infeasible_->denseVector();
217     double largest = 0.0;
218     int * index = infeasible_->getIndices();
219     int number = infeasible_->getNumElements();
220     const int * pivotVariable = model_->pivotVariable();
221     int chosenRow = -1;
222     int lastPivotRow = model_->pivotRow();
223     assert (lastPivotRow < model_->numberRows());
224     double tolerance = model_->currentPrimalTolerance();
225     // we can't really trust infeasibilities if there is primal error
226     // this coding has to mimic coding in checkPrimalSolution
227     double error = CoinMin(1.0e-2, model_->largestPrimalError());
228     // allow tolerance at least slightly bigger than standard
229     tolerance = tolerance +  error;
230     // But cap
231     tolerance = CoinMin(1000.0, tolerance);
232     tolerance *= tolerance; // as we are using squares
233     bool toleranceChanged = false;
234     double * solution = model_->solutionRegion();
235     double * lower = model_->lowerRegion();
236     double * upper = model_->upperRegion();
237     // do last pivot row one here
238     //#define COLUMN_BIAS 4.0
239     //#define FIXED_BIAS 10.0
240     if (lastPivotRow >= 0 && lastPivotRow < model_->numberRows()) {
241#ifdef COLUMN_BIAS
242          int numberColumns = model_->numberColumns();
243#endif
244          int iPivot = pivotVariable[lastPivotRow];
245          double value = solution[iPivot];
246          double lower = model_->lower(iPivot);
247          double upper = model_->upper(iPivot);
248          if (value > upper + tolerance) {
249               value -= upper;
250               value *= value;
251#ifdef COLUMN_BIAS
252               if (iPivot < numberColumns)
253                    value *= COLUMN_BIAS; // bias towards columns
254               k
255#endif
256               // store square in list
257               if (infeas[lastPivotRow])
258                    infeas[lastPivotRow] = value; // already there
259               else
260                    infeasible_->quickAdd(lastPivotRow, value);
261          } else if (value < lower - tolerance) {
262               value -= lower;
263               value *= value;
264#ifdef COLUMN_BIAS
265               if (iPivot < numberColumns)
266                    value *= COLUMN_BIAS; // bias towards columns
267#endif
268               // store square in list
269               if (infeas[lastPivotRow])
270                    infeas[lastPivotRow] = value; // already there
271               else
272                    infeasible_->add(lastPivotRow, value);
273          } else {
274               // feasible - was it infeasible - if so set tiny
275               if (infeas[lastPivotRow])
276                    infeas[lastPivotRow] = COIN_INDEXED_REALLY_TINY_ELEMENT;
277          }
278          number = infeasible_->getNumElements();
279     }
280     if(model_->numberIterations() < model_->lastBadIteration() + 200) {
281          // we can't really trust infeasibilities if there is dual error
282          if (model_->largestDualError() > model_->largestPrimalError()) {
283               tolerance *= CoinMin(model_->largestDualError() / model_->largestPrimalError(), 1000.0);
284               toleranceChanged = true;
285          }
286     }
287     int numberWanted;
288     if (mode_ < 2 ) {
289          numberWanted = number + 1;
290     } else if (mode_ == 2) {
291          numberWanted = CoinMax(2000, number / 8);
292     } else {
293          int numberElements = model_->factorization()->numberElements();
294          double ratio = static_cast<double> (numberElements) /
295                         static_cast<double> (model_->numberRows());
296          numberWanted = CoinMax(2000, number / 8);
297          if (ratio < 1.0) {
298               numberWanted = CoinMax(2000, number / 20);
299          } else if (ratio > 10.0) {
300               ratio = number * (ratio / 80.0);
301               if (ratio > number)
302                    numberWanted = number + 1;
303               else
304                    numberWanted = CoinMax(2000, static_cast<int> (ratio));
305          }
306     }
307     if (model_->largestPrimalError() > 1.0e-3)
308          numberWanted = number + 1; // be safe
309     int iPass;
310     // Setup two passes
311     int start[4];
312     start[1] = number;
313     start[2] = 0;
314     double dstart = static_cast<double> (number) * model_->randomNumberGenerator()->randomDouble();
315     start[0] = static_cast<int> (dstart);
316     start[3] = start[0];
317     //double largestWeight=0.0;
318     //double smallestWeight=1.0e100;
319     for (iPass = 0; iPass < 2; iPass++) {
320          int end = start[2*iPass+1];
321          for (i = start[2*iPass]; i < end; i++) {
322               iRow = index[i];
323               double value = infeas[iRow];
324               if (value > tolerance) {
325                    //#define OUT_EQ
326#ifdef OUT_EQ
327                    {
328                         int iSequence = pivotVariable[iRow];
329                         if (upper[iSequence] == lower[iSequence])
330                              value *= 2.0;
331                    }
332#endif
333                    double weight = CoinMin(weights_[iRow], 1.0e50);
334                    //largestWeight = CoinMax(largestWeight,weight);
335                    //smallestWeight = CoinMin(smallestWeight,weight);
336                    //double dubious = dubiousWeights_[iRow];
337                    //weight *= dubious;
338                    //if (value>2.0*largest*weight||(value>0.5*largest*weight&&value*largestWeight>dubious*largest*weight)) {
339                    if (value > largest * weight) {
340                         // make last pivot row last resort choice
341                         if (iRow == lastPivotRow) {
342                              if (value * 1.0e-10 < largest * weight)
343                                   continue;
344                              else
345                                   value *= 1.0e-10;
346                         }
347                         int iSequence = pivotVariable[iRow];
348                         if (!model_->flagged(iSequence)) {
349                              //#define CLP_DEBUG 3
350#ifdef CLP_DEBUG
351                              double value2 = 0.0;
352                              if (solution[iSequence] > upper[iSequence] + tolerance)
353                                   value2 = solution[iSequence] - upper[iSequence];
354                              else if (solution[iSequence] < lower[iSequence] - tolerance)
355                                   value2 = solution[iSequence] - lower[iSequence];
356                              assert(fabs(value2 * value2 - infeas[iRow]) < 1.0e-8 * CoinMin(value2 * value2, infeas[iRow]));
357#endif
358                              if (solution[iSequence] > upper[iSequence] + tolerance ||
359                                        solution[iSequence] < lower[iSequence] - tolerance) {
360                                   chosenRow = iRow;
361                                   largest = value / weight;
362                                   //largestWeight = dubious;
363                              }
364                         } else {
365                              // just to make sure we don't exit before got something
366                              numberWanted++;
367                         }
368                    }
369                    numberWanted--;
370                    if (!numberWanted)
371                         break;
372               }
373          }
374          if (!numberWanted)
375               break;
376     }
377     //printf("smallest %g largest %g\n",smallestWeight,largestWeight);
378     if (chosenRow < 0 && toleranceChanged) {
379          // won't line up with checkPrimalSolution - do again
380          double saveError = model_->largestDualError();
381          model_->setLargestDualError(0.0);
382          // can't loop
383          chosenRow = pivotRow();
384          model_->setLargestDualError(saveError);
385     }
386     if (chosenRow < 0 && lastPivotRow < 0) {
387          int nLeft = 0;
388          for (int i = 0; i < number; i++) {
389               int iRow = index[i];
390               if (fabs(infeas[iRow]) > 1.0e-50) {
391                    index[nLeft++] = iRow;
392               } else {
393                    infeas[iRow] = 0.0;
394               }
395          }
396          infeasible_->setNumElements(nLeft);
397          model_->setNumberPrimalInfeasibilities(nLeft);
398     }
399     return chosenRow;
400}
401#if 0
402static double ft_count = 0.0;
403static double up_count = 0.0;
404static double ft_count_in = 0.0;
405static double up_count_in = 0.0;
406static int xx_count = 0;
407#endif
408/* Updates weights and returns pivot alpha.
409   Also does FT update */
410double
411ClpDualRowSteepest::updateWeights(CoinIndexedVector * input,
412                                  CoinIndexedVector * spare,
413                                  CoinIndexedVector * spare2,
414                                  CoinIndexedVector * updatedColumn)
415{
416     //#define CLP_DEBUG 3
417#if CLP_DEBUG>2
418     // Very expensive debug
419     {
420          int numberRows = model_->numberRows();
421          CoinIndexedVector * temp = new CoinIndexedVector();
422          temp->reserve(numberRows +
423                        model_->factorization()->maximumPivots());
424          double * array = alternateWeights_->denseVector();
425          int * which = alternateWeights_->getIndices();
426          alternateWeights_->clear();
427          int i;
428          for (i = 0; i < numberRows; i++) {
429               double value = 0.0;
430               array[i] = 1.0;
431               which[0] = i;
432               alternateWeights_->setNumElements(1);
433               model_->factorization()->updateColumnTranspose(temp,
434                         alternateWeights_);
435               int number = alternateWeights_->getNumElements();
436               int j;
437               for (j = 0; j < number; j++) {
438                    int iRow = which[j];
439                    value += array[iRow] * array[iRow];
440                    array[iRow] = 0.0;
441               }
442               alternateWeights_->setNumElements(0);
443               double w = CoinMax(weights_[i], value) * .1;
444               if (fabs(weights_[i] - value) > w) {
445                    printf("%d old %g, true %g\n", i, weights_[i], value);
446                    weights_[i] = value; // to reduce printout
447               }
448               //else
449               //printf("%d matches %g\n",i,value);
450          }
451          delete temp;
452     }
453#endif
454     assert (input->packedMode());
455     if (!updatedColumn->packedMode()) {
456          // I think this means empty
457#ifdef COIN_DEVELOP
458          printf("updatedColumn not packed mode ClpDualRowSteepest::updateWeights\n");
459#endif
460          return 0.0;
461     }
462     double alpha = 0.0;
463     if (!model_->factorization()->networkBasis()) {
464          // clear other region
465          alternateWeights_->clear();
466          double norm = 0.0;
467          int i;
468          double * work = input->denseVector();
469          int numberNonZero = input->getNumElements();
470          int * which = input->getIndices();
471          double * work2 = spare->denseVector();
472          int * which2 = spare->getIndices();
473          // ftran
474          //permute and move indices into index array
475          //also compute norm
476          //int *regionIndex = alternateWeights_->getIndices (  );
477          const int *permute = model_->factorization()->permute();
478          //double * region = alternateWeights_->denseVector();
479          if (permute) {
480               for ( i = 0; i < numberNonZero; i ++ ) {
481                    int iRow = which[i];
482                    double value = work[i];
483                    norm += value * value;
484                    iRow = permute[iRow];
485                    work2[iRow] = value;
486                    which2[i] = iRow;
487               }
488          } else {
489               for ( i = 0; i < numberNonZero; i ++ ) {
490                    int iRow = which[i];
491                    double value = work[i];
492                    norm += value * value;
493                    //iRow = permute[iRow];
494                    work2[iRow] = value;
495                    which2[i] = iRow;
496               }
497          }
498          spare->setNumElements ( numberNonZero );
499          // Do FT update
500#if 0
501          ft_count_in += updatedColumn->getNumElements();
502          up_count_in += spare->getNumElements();
503#endif
504          if (permute || true) {
505#if CLP_DEBUG>2
506               printf("REGION before %d els\n", spare->getNumElements());
507               spare->print();
508#endif
509               model_->factorization()->updateTwoColumnsFT(spare2, updatedColumn,
510                         spare, permute != NULL);
511#if CLP_DEBUG>2
512               printf("REGION after %d els\n", spare->getNumElements());
513               spare->print();
514#endif
515          } else {
516               // Leave as old way
517               model_->factorization()->updateColumnFT(spare2, updatedColumn);
518               model_->factorization()->updateColumn(spare2, spare, false);
519          }
520#undef CLP_DEBUG
521#if 0
522          ft_count += updatedColumn->getNumElements();
523          up_count += spare->getNumElements();
524          xx_count++;
525          if ((xx_count % 1000) == 0)
526               printf("zz %d ft %g up %g (in %g %g)\n", xx_count, ft_count, up_count,
527                      ft_count_in, up_count_in);
528#endif
529          numberNonZero = spare->getNumElements();
530          // alternateWeights_ should still be empty
531          int pivotRow = model_->pivotRow();
532#ifdef CLP_DEBUG
533          if ( model_->logLevel (  ) > 4  &&
534                    fabs(norm - weights_[pivotRow]) > 1.0e-3 * (1.0 + norm))
535               printf("on row %d, true weight %g, old %g\n",
536                      pivotRow, sqrt(norm), sqrt(weights_[pivotRow]));
537#endif
538          // could re-initialize here (could be expensive)
539          norm /= model_->alpha() * model_->alpha();
540          assert(model_->alpha());
541          assert(norm);
542          // pivot element
543          alpha = 0.0;
544          double multiplier = 2.0 / model_->alpha();
545          // look at updated column
546          work = updatedColumn->denseVector();
547          numberNonZero = updatedColumn->getNumElements();
548          which = updatedColumn->getIndices();
549
550          int nSave = 0;
551          double * work3 = alternateWeights_->denseVector();
552          int * which3 = alternateWeights_->getIndices();
553          const int * pivotColumn = model_->factorization()->pivotColumn();
554          for (i = 0; i < numberNonZero; i++) {
555               int iRow = which[i];
556               double theta = work[i];
557               if (iRow == pivotRow)
558                    alpha = theta;
559               double devex = weights_[iRow];
560               work3[nSave] = devex; // save old
561               which3[nSave++] = iRow;
562               // transform to match spare
563               int jRow = permute ? pivotColumn[iRow] : iRow;
564               double value = work2[jRow];
565               devex +=  theta * (theta * norm + value * multiplier);
566               if (devex < DEVEX_TRY_NORM)
567                    devex = DEVEX_TRY_NORM;
568               weights_[iRow] = devex;
569          }
570          alternateWeights_->setPackedMode(true);
571          alternateWeights_->setNumElements(nSave);
572          if (norm < DEVEX_TRY_NORM)
573               norm = DEVEX_TRY_NORM;
574          // Try this to make less likely will happen again and stop cycling
575          //norm *= 1.02;
576          weights_[pivotRow] = norm;
577          spare->clear();
578#ifdef CLP_DEBUG
579          spare->checkClear();
580#endif
581     } else {
582          // Do FT update
583          model_->factorization()->updateColumnFT(spare, updatedColumn);
584          // clear other region
585          alternateWeights_->clear();
586          double norm = 0.0;
587          int i;
588          double * work = input->denseVector();
589          int number = input->getNumElements();
590          int * which = input->getIndices();
591          double * work2 = spare->denseVector();
592          int * which2 = spare->getIndices();
593          for (i = 0; i < number; i++) {
594               int iRow = which[i];
595               double value = work[i];
596               norm += value * value;
597               work2[iRow] = value;
598               which2[i] = iRow;
599          }
600          spare->setNumElements(number);
601          // ftran
602#ifndef NDEBUG
603          alternateWeights_->checkClear();
604#endif
605          model_->factorization()->updateColumn(alternateWeights_, spare);
606          // alternateWeights_ should still be empty
607#ifndef NDEBUG
608          alternateWeights_->checkClear();
609#endif
610          int pivotRow = model_->pivotRow();
611#ifdef CLP_DEBUG
612          if ( model_->logLevel (  ) > 4  &&
613                    fabs(norm - weights_[pivotRow]) > 1.0e-3 * (1.0 + norm))
614               printf("on row %d, true weight %g, old %g\n",
615                      pivotRow, sqrt(norm), sqrt(weights_[pivotRow]));
616#endif
617          // could re-initialize here (could be expensive)
618          norm /= model_->alpha() * model_->alpha();
619
620          assert(norm);
621          //if (norm < DEVEX_TRY_NORM)
622          //norm = DEVEX_TRY_NORM;
623          // pivot element
624          alpha = 0.0;
625          double multiplier = 2.0 / model_->alpha();
626          // look at updated column
627          work = updatedColumn->denseVector();
628          number = updatedColumn->getNumElements();
629          which = updatedColumn->getIndices();
630
631          int nSave = 0;
632          double * work3 = alternateWeights_->denseVector();
633          int * which3 = alternateWeights_->getIndices();
634          for (i = 0; i < number; i++) {
635               int iRow = which[i];
636               double theta = work[i];
637               if (iRow == pivotRow)
638                    alpha = theta;
639               double devex = weights_[iRow];
640               work3[nSave] = devex; // save old
641               which3[nSave++] = iRow;
642               double value = work2[iRow];
643               devex +=  theta * (theta * norm + value * multiplier);
644               if (devex < DEVEX_TRY_NORM)
645                    devex = DEVEX_TRY_NORM;
646               weights_[iRow] = devex;
647          }
648          if (!alpha) {
649               // error - but carry on
650               alpha = 1.0e-50;
651          }
652          alternateWeights_->setPackedMode(true);
653          alternateWeights_->setNumElements(nSave);
654          if (norm < DEVEX_TRY_NORM)
655               norm = DEVEX_TRY_NORM;
656          weights_[pivotRow] = norm;
657          spare->clear();
658     }
659#ifdef CLP_DEBUG
660     spare->checkClear();
661#endif
662     return alpha;
663}
664
665/* Updates primal solution (and maybe list of candidates)
666   Uses input vector which it deletes
667   Computes change in objective function
668*/
669void
670ClpDualRowSteepest::updatePrimalSolution(
671     CoinIndexedVector * primalUpdate,
672     double primalRatio,
673     double & objectiveChange)
674{
675     double * COIN_RESTRICT work = primalUpdate->denseVector();
676     int number = primalUpdate->getNumElements();
677     int * COIN_RESTRICT which = primalUpdate->getIndices();
678     int i;
679     double changeObj = 0.0;
680     double tolerance = model_->currentPrimalTolerance();
681     const int * COIN_RESTRICT pivotVariable = model_->pivotVariable();
682     double * COIN_RESTRICT infeas = infeasible_->denseVector();
683     double * COIN_RESTRICT solution = model_->solutionRegion();
684     const double * COIN_RESTRICT costModel = model_->costRegion();
685     const double * COIN_RESTRICT lowerModel = model_->lowerRegion();
686     const double * COIN_RESTRICT upperModel = model_->upperRegion();
687#ifdef COLUMN_BIAS
688     int numberColumns = model_->numberColumns();
689#endif
690     if (primalUpdate->packedMode()) {
691          for (i = 0; i < number; i++) {
692               int iRow = which[i];
693               int iPivot = pivotVariable[iRow];
694               double value = solution[iPivot];
695               double cost = costModel[iPivot];
696               double change = primalRatio * work[i];
697               work[i] = 0.0;
698               value -= change;
699               changeObj -= change * cost;
700               double lower = lowerModel[iPivot];
701               double upper = upperModel[iPivot];
702               solution[iPivot] = value;
703               if (value < lower - tolerance) {
704                    value -= lower;
705                    value *= value;
706#ifdef COLUMN_BIAS
707                    if (iPivot < numberColumns)
708                         value *= COLUMN_BIAS; // bias towards columns
709#endif
710#ifdef FIXED_BIAS
711                    if (lower == upper)
712                         value *= FIXED_BIAS; // bias towards taking out fixed variables
713#endif
714                    // store square in list
715                    if (infeas[iRow])
716                         infeas[iRow] = value; // already there
717                    else
718                         infeasible_->quickAdd(iRow, value);
719               } else if (value > upper + tolerance) {
720                    value -= upper;
721                    value *= value;
722#ifdef COLUMN_BIAS
723                    if (iPivot < numberColumns)
724                         value *= COLUMN_BIAS; // bias towards columns
725#endif
726#ifdef FIXED_BIAS
727                    if (lower == upper)
728                         value *= FIXED_BIAS; // bias towards taking out fixed variables
729#endif
730                    // store square in list
731                    if (infeas[iRow])
732                         infeas[iRow] = value; // already there
733                    else
734                         infeasible_->quickAdd(iRow, value);
735               } else {
736                    // feasible - was it infeasible - if so set tiny
737                    if (infeas[iRow])
738                         infeas[iRow] = COIN_INDEXED_REALLY_TINY_ELEMENT;
739               }
740          }
741     } else {
742          for (i = 0; i < number; i++) {
743               int iRow = which[i];
744               int iPivot = pivotVariable[iRow];
745               double value = solution[iPivot];
746               double cost = costModel[iPivot];
747               double change = primalRatio * work[iRow];
748               value -= change;
749               changeObj -= change * cost;
750               double lower = lowerModel[iPivot];
751               double upper = upperModel[iPivot];
752               solution[iPivot] = value;
753               if (value < lower - tolerance) {
754                    value -= lower;
755                    value *= value;
756#ifdef COLUMN_BIAS
757                    if (iPivot < numberColumns)
758                         value *= COLUMN_BIAS; // bias towards columns
759#endif
760#ifdef FIXED_BIAS
761                    if (lower == upper)
762                         value *= FIXED_BIAS; // bias towards taking out fixed variables
763#endif
764                    // store square in list
765                    if (infeas[iRow])
766                         infeas[iRow] = value; // already there
767                    else
768                         infeasible_->quickAdd(iRow, value);
769               } else if (value > upper + tolerance) {
770                    value -= upper;
771                    value *= value;
772#ifdef COLUMN_BIAS
773                    if (iPivot < numberColumns)
774                         value *= COLUMN_BIAS; // bias towards columns
775#endif
776#ifdef FIXED_BIAS
777                    if (lower == upper)
778                         value *= FIXED_BIAS; // bias towards taking out fixed variables
779#endif
780                    // store square in list
781                    if (infeas[iRow])
782                         infeas[iRow] = value; // already there
783                    else
784                         infeasible_->quickAdd(iRow, value);
785               } else {
786                    // feasible - was it infeasible - if so set tiny
787                    if (infeas[iRow])
788                         infeas[iRow] = COIN_INDEXED_REALLY_TINY_ELEMENT;
789               }
790               work[iRow] = 0.0;
791          }
792     }
793     // Do pivot row
794     {
795          int iRow = model_->pivotRow();
796          // feasible - was it infeasible - if so set tiny
797          //assert (infeas[iRow]);
798          if (infeas[iRow])
799               infeas[iRow] = COIN_INDEXED_REALLY_TINY_ELEMENT;
800     }
801     primalUpdate->setNumElements(0);
802     objectiveChange += changeObj;
803}
804/* Saves any weights round factorization as pivot rows may change
805   1) before factorization
806   2) after factorization
807   3) just redo infeasibilities
808   4) restore weights
809*/
810void
811ClpDualRowSteepest::saveWeights(ClpSimplex * model, int mode)
812{
813     // alternateWeights_ is defined as indexed but is treated oddly
814     model_ = model;
815     int numberRows = model_->numberRows();
816     int numberColumns = model_->numberColumns();
817     const int * pivotVariable = model_->pivotVariable();
818     int i;
819     if (mode == 1) {
820          if(weights_) {
821               // Check if size has changed
822               if (infeasible_->capacity() == numberRows) {
823                    alternateWeights_->clear();
824                    // change from row numbers to sequence numbers
825                    int * which = alternateWeights_->getIndices();
826                    for (i = 0; i < numberRows; i++) {
827                         int iPivot = pivotVariable[i];
828                         which[i] = iPivot;
829                    }
830                    state_ = 1;
831               } else {
832                    // size has changed - clear everything
833                    delete [] weights_;
834                    weights_ = NULL;
835                    delete [] dubiousWeights_;
836                    dubiousWeights_ = NULL;
837                    delete infeasible_;
838                    infeasible_ = NULL;
839                    delete alternateWeights_;
840                    alternateWeights_ = NULL;
841                    delete savedWeights_;
842                    savedWeights_ = NULL;
843                    state_ = -1;
844               }
845          }
846     } else if (mode == 2 || mode == 4 || mode >= 5) {
847          // restore
848          if (!weights_ || state_ == -1 || mode == 5) {
849               // initialize weights
850               delete [] weights_;
851               delete alternateWeights_;
852               weights_ = new double[numberRows];
853               alternateWeights_ = new CoinIndexedVector();
854               // enough space so can use it for factorization
855               alternateWeights_->reserve(numberRows +
856                                          model_->factorization()->maximumPivots());
857               if (mode_ != 1 || mode == 5) {
858                    // initialize to 1.0 (can we do better?)
859                    for (i = 0; i < numberRows; i++) {
860                         weights_[i] = 1.0;
861                    }
862               } else {
863                    CoinIndexedVector * temp = new CoinIndexedVector();
864                    temp->reserve(numberRows +
865                                  model_->factorization()->maximumPivots());
866                    double * array = alternateWeights_->denseVector();
867                    int * which = alternateWeights_->getIndices();
868                    for (i = 0; i < numberRows; i++) {
869                         double value = 0.0;
870                         array[0] = 1.0;
871                         which[0] = i;
872                         alternateWeights_->setNumElements(1);
873                         alternateWeights_->setPackedMode(true);
874                         model_->factorization()->updateColumnTranspose(temp,
875                                   alternateWeights_);
876                         int number = alternateWeights_->getNumElements();
877                         int j;
878                         for (j = 0; j < number; j++) {
879                              value += array[j] * array[j];
880                              array[j] = 0.0;
881                         }
882                         alternateWeights_->setNumElements(0);
883                         weights_[i] = value;
884                    }
885                    delete temp;
886               }
887               // create saved weights (not really indexedvector)
888               savedWeights_ = new CoinIndexedVector();
889               savedWeights_->reserve(numberRows);
890
891               double * array = savedWeights_->denseVector();
892               int * which = savedWeights_->getIndices();
893               for (i = 0; i < numberRows; i++) {
894                    array[i] = weights_[i];
895                    which[i] = pivotVariable[i];
896               }
897          } else if (mode != 6) {
898               int * which = alternateWeights_->getIndices();
899               CoinIndexedVector * rowArray3 = model_->rowArray(3);
900               rowArray3->clear();
901               int * back = rowArray3->getIndices();
902               // In case something went wrong
903               for (i = 0; i < numberRows + numberColumns; i++)
904                    back[i] = -1;
905               if (mode != 4) {
906                    // save
907                    CoinMemcpyN(which,  numberRows, savedWeights_->getIndices());
908                    CoinMemcpyN(weights_,       numberRows, savedWeights_->denseVector());
909               } else {
910                    // restore
911                    //memcpy(which,savedWeights_->getIndices(),
912                    //     numberRows*sizeof(int));
913                    //memcpy(weights_,savedWeights_->denseVector(),
914                    //     numberRows*sizeof(double));
915                    which = savedWeights_->getIndices();
916               }
917               // restore (a bit slow - but only every re-factorization)
918               double * array = savedWeights_->denseVector();
919               for (i = 0; i < numberRows; i++) {
920                    int iSeq = which[i];
921                    back[iSeq] = i;
922               }
923               for (i = 0; i < numberRows; i++) {
924                    int iPivot = pivotVariable[i];
925                    iPivot = back[iPivot];
926                    if (iPivot >= 0) {
927                         weights_[i] = array[iPivot];
928                         if (weights_[i] < DEVEX_TRY_NORM)
929                              weights_[i] = DEVEX_TRY_NORM; // may need to check more
930                    } else {
931                         // odd
932                         weights_[i] = 1.0;
933                    }
934               }
935          } else {
936               // mode 6 - scale back weights as primal errors
937               double primalError = model_->largestPrimalError();
938               double allowed ;
939               if (primalError > 1.0e3)
940                    allowed = 10.0;
941               else if (primalError > 1.0e2)
942                    allowed = 50.0;
943               else if (primalError > 1.0e1)
944                    allowed = 100.0;
945               else
946                    allowed = 1000.0;
947               double allowedInv = 1.0 / allowed;
948               for (i = 0; i < numberRows; i++) {
949                    double value = weights_[i];
950                    if (value < allowedInv)
951                         value = allowedInv;
952                    else if (value > allowed)
953                         value = allowed;
954                    weights_[i] = allowed;
955               }
956          }
957          state_ = 0;
958          // set up infeasibilities
959          if (!infeasible_) {
960               infeasible_ = new CoinIndexedVector();
961               infeasible_->reserve(numberRows);
962          }
963     }
964     if (mode >= 2) {
965          // Get dubious weights
966          //if (!dubiousWeights_)
967          //dubiousWeights_=new int[numberRows];
968          //model_->factorization()->getWeights(dubiousWeights_);
969          infeasible_->clear();
970          int iRow;
971          const int * pivotVariable = model_->pivotVariable();
972          double tolerance = model_->currentPrimalTolerance();
973          for (iRow = 0; iRow < numberRows; iRow++) {
974               int iPivot = pivotVariable[iRow];
975               double value = model_->solution(iPivot);
976               double lower = model_->lower(iPivot);
977               double upper = model_->upper(iPivot);
978               if (value < lower - tolerance) {
979                    value -= lower;
980                    value *= value;
981#ifdef COLUMN_BIAS
982                    if (iPivot < numberColumns)
983                         value *= COLUMN_BIAS; // bias towards columns
984#endif
985#ifdef FIXED_BIAS
986                    if (lower == upper)
987                         value *= FIXED_BIAS; // bias towards taking out fixed variables
988#endif
989                    // store square in list
990                    infeasible_->quickAdd(iRow, value);
991               } else if (value > upper + tolerance) {
992                    value -= upper;
993                    value *= value;
994#ifdef COLUMN_BIAS
995                    if (iPivot < numberColumns)
996                         value *= COLUMN_BIAS; // bias towards columns
997#endif
998#ifdef FIXED_BIAS
999                    if (lower == upper)
1000                         value *= FIXED_BIAS; // bias towards taking out fixed variables
1001#endif
1002                    // store square in list
1003                    infeasible_->quickAdd(iRow, value);
1004               }
1005          }
1006     }
1007}
1008// Gets rid of last update
1009void
1010ClpDualRowSteepest::unrollWeights()
1011{
1012     double * saved = alternateWeights_->denseVector();
1013     int number = alternateWeights_->getNumElements();
1014     int * which = alternateWeights_->getIndices();
1015     int i;
1016     if (alternateWeights_->packedMode()) {
1017          for (i = 0; i < number; i++) {
1018               int iRow = which[i];
1019               weights_[iRow] = saved[i];
1020               saved[i] = 0.0;
1021          }
1022     } else {
1023          for (i = 0; i < number; i++) {
1024               int iRow = which[i];
1025               weights_[iRow] = saved[iRow];
1026               saved[iRow] = 0.0;
1027          }
1028     }
1029     alternateWeights_->setNumElements(0);
1030}
1031//-------------------------------------------------------------------
1032// Clone
1033//-------------------------------------------------------------------
1034ClpDualRowPivot * ClpDualRowSteepest::clone(bool CopyData) const
1035{
1036     if (CopyData) {
1037          return new ClpDualRowSteepest(*this);
1038     } else {
1039          return new ClpDualRowSteepest();
1040     }
1041}
1042// Gets rid of all arrays
1043void
1044ClpDualRowSteepest::clearArrays()
1045{
1046     if (persistence_ == normal) {
1047          delete [] weights_;
1048          weights_ = NULL;
1049          delete [] dubiousWeights_;
1050          dubiousWeights_ = NULL;
1051          delete infeasible_;
1052          infeasible_ = NULL;
1053          delete alternateWeights_;
1054          alternateWeights_ = NULL;
1055          delete savedWeights_;
1056          savedWeights_ = NULL;
1057     }
1058     state_ = -1;
1059}
1060// Returns true if would not find any row
1061bool
1062ClpDualRowSteepest::looksOptimal() const
1063{
1064     int iRow;
1065     const int * pivotVariable = model_->pivotVariable();
1066     double tolerance = model_->currentPrimalTolerance();
1067     // we can't really trust infeasibilities if there is primal error
1068     // this coding has to mimic coding in checkPrimalSolution
1069     double error = CoinMin(1.0e-2, model_->largestPrimalError());
1070     // allow tolerance at least slightly bigger than standard
1071     tolerance = tolerance +  error;
1072     // But cap
1073     tolerance = CoinMin(1000.0, tolerance);
1074     int numberRows = model_->numberRows();
1075     int numberInfeasible = 0;
1076     for (iRow = 0; iRow < numberRows; iRow++) {
1077          int iPivot = pivotVariable[iRow];
1078          double value = model_->solution(iPivot);
1079          double lower = model_->lower(iPivot);
1080          double upper = model_->upper(iPivot);
1081          if (value < lower - tolerance) {
1082               numberInfeasible++;
1083          } else if (value > upper + tolerance) {
1084               numberInfeasible++;
1085          }
1086     }
1087     return (numberInfeasible == 0);
1088}
1089// Called when maximum pivots changes
1090void
1091ClpDualRowSteepest::maximumPivotsChanged()
1092{
1093     if (alternateWeights_ &&
1094               alternateWeights_->capacity() != model_->numberRows() +
1095               model_->factorization()->maximumPivots()) {
1096          delete alternateWeights_;
1097          alternateWeights_ = new CoinIndexedVector();
1098          // enough space so can use it for factorization
1099          alternateWeights_->reserve(model_->numberRows() +
1100                                     model_->factorization()->maximumPivots());
1101     }
1102}
1103
Note: See TracBrowser for help on using the repository browser.