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

Last change on this file since 2271 was 2070, checked in by forrest, 5 years ago

for some cbc ideas

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 44.1 KB
Line 
1/* $Id: ClpDualRowSteepest.cpp 2070 2014-11-18 11:12:54Z forrest $ */
2// Copyright (C) 2002, International Business Machines
3// Corporation and others.  All Rights Reserved.
4// This code is licensed under the terms of the Eclipse Public License (EPL).
5
6#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 CLP_DUAL_FIXED_COLUMN_MULTIPLIER 10.0
239     if (lastPivotRow >= 0 && lastPivotRow < model_->numberRows()) {
240#ifdef CLP_DUAL_COLUMN_MULTIPLIER
241          int numberColumns = model_->numberColumns();
242#endif
243          int iPivot = pivotVariable[lastPivotRow];
244          double value = solution[iPivot];
245          double lower = model_->lower(iPivot);
246          double upper = model_->upper(iPivot);
247          if (value > upper + tolerance) {
248               value -= upper;
249               value *= value;
250#ifdef CLP_DUAL_COLUMN_MULTIPLIER
251               if (iPivot < numberColumns)
252                    value *= CLP_DUAL_COLUMN_MULTIPLIER; // bias towards columns
253#endif
254               // store square in list
255               if (infeas[lastPivotRow])
256                    infeas[lastPivotRow] = value; // already there
257               else
258                    infeasible_->quickAdd(lastPivotRow, value);
259          } else if (value < lower - tolerance) {
260               value -= lower;
261               value *= value;
262#ifdef CLP_DUAL_COLUMN_MULTIPLIER
263               if (iPivot < numberColumns)
264                    value *= CLP_DUAL_COLUMN_MULTIPLIER; // bias towards columns
265#endif
266               // store square in list
267               if (infeas[lastPivotRow])
268                    infeas[lastPivotRow] = value; // already there
269               else
270                    infeasible_->add(lastPivotRow, value);
271          } else {
272               // feasible - was it infeasible - if so set tiny
273               if (infeas[lastPivotRow])
274                    infeas[lastPivotRow] = COIN_INDEXED_REALLY_TINY_ELEMENT;
275          }
276          number = infeasible_->getNumElements();
277     }
278     if(model_->numberIterations() < model_->lastBadIteration() + 200) {
279          // we can't really trust infeasibilities if there is dual error
280          if (model_->largestDualError() > model_->largestPrimalError()) {
281               tolerance *= CoinMin(model_->largestDualError() / model_->largestPrimalError(), 1000.0);
282               toleranceChanged = true;
283          }
284     }
285     int numberWanted;
286     if (mode_ < 2 ) {
287          numberWanted = number + 1;
288     } else if (mode_ == 2) {
289          numberWanted = CoinMax(2000, number / 8);
290     } else {
291          int numberElements = model_->factorization()->numberElements();
292          double ratio = static_cast<double> (numberElements) /
293                         static_cast<double> (model_->numberRows());
294          numberWanted = CoinMax(2000, number / 8);
295          if (ratio < 1.0) {
296               numberWanted = CoinMax(2000, number / 20);
297          } else if (ratio > 10.0) {
298               ratio = number * (ratio / 80.0);
299               if (ratio > number)
300                    numberWanted = number + 1;
301               else
302                    numberWanted = CoinMax(2000, static_cast<int> (ratio));
303          }
304     }
305     if (model_->largestPrimalError() > 1.0e-3)
306          numberWanted = number + 1; // be safe
307     int iPass;
308     // Setup two passes
309     int start[4];
310     start[1] = number;
311     start[2] = 0;
312     double dstart = static_cast<double> (number) * model_->randomNumberGenerator()->randomDouble();
313     start[0] = static_cast<int> (dstart);
314     start[3] = start[0];
315     //double largestWeight=0.0;
316     //double smallestWeight=1.0e100;
317     for (iPass = 0; iPass < 2; iPass++) {
318          int end = start[2*iPass+1];
319          for (i = start[2*iPass]; i < end; i++) {
320               iRow = index[i];
321               double value = infeas[iRow];
322               if (value > tolerance) {
323                    //#define OUT_EQ
324#ifdef OUT_EQ
325                    {
326                         int iSequence = pivotVariable[iRow];
327                         if (upper[iSequence] == lower[iSequence])
328                              value *= 2.0;
329                    }
330#endif
331                    double weight = CoinMin(weights_[iRow], 1.0e50);
332                    //largestWeight = CoinMax(largestWeight,weight);
333                    //smallestWeight = CoinMin(smallestWeight,weight);
334                    //double dubious = dubiousWeights_[iRow];
335                    //weight *= dubious;
336                    //if (value>2.0*largest*weight||(value>0.5*largest*weight&&value*largestWeight>dubious*largest*weight)) {
337                    if (value > largest * weight) {
338                         // make last pivot row last resort choice
339                         if (iRow == lastPivotRow) {
340                              if (value * 1.0e-10 < largest * weight)
341                                   continue;
342                              else
343                                   value *= 1.0e-10;
344                         }
345                         int iSequence = pivotVariable[iRow];
346                         if (!model_->flagged(iSequence)) {
347                              //#define CLP_DEBUG 3
348#ifdef CLP_DEBUG
349                              double value2 = 0.0;
350                              if (solution[iSequence] > upper[iSequence] + tolerance)
351                                   value2 = solution[iSequence] - upper[iSequence];
352                              else if (solution[iSequence] < lower[iSequence] - tolerance)
353                                   value2 = solution[iSequence] - lower[iSequence];
354                              assert(fabs(value2 * value2 - infeas[iRow]) < 1.0e-8 * CoinMin(value2 * value2, infeas[iRow]));
355#endif
356                              if (solution[iSequence] > upper[iSequence] + tolerance ||
357                                        solution[iSequence] < lower[iSequence] - tolerance) {
358                                   chosenRow = iRow;
359                                   largest = value / weight;
360                                   //largestWeight = dubious;
361                              }
362                         } else {
363                              // just to make sure we don't exit before got something
364                              numberWanted++;
365                         }
366                    }
367                    numberWanted--;
368                    if (!numberWanted)
369                         break;
370               }
371          }
372          if (!numberWanted)
373               break;
374     }
375     //printf("smallest %g largest %g\n",smallestWeight,largestWeight);
376     if (chosenRow < 0 && toleranceChanged) {
377          // won't line up with checkPrimalSolution - do again
378          double saveError = model_->largestDualError();
379          model_->setLargestDualError(0.0);
380          // can't loop
381          chosenRow = pivotRow();
382          model_->setLargestDualError(saveError);
383     }
384     if (chosenRow < 0 && lastPivotRow < 0) {
385          int nLeft = 0;
386          for (int i = 0; i < number; i++) {
387               int iRow = index[i];
388               if (fabs(infeas[iRow]) > 1.0e-50) {
389                    index[nLeft++] = iRow;
390               } else {
391                    infeas[iRow] = 0.0;
392               }
393          }
394          infeasible_->setNumElements(nLeft);
395          model_->setNumberPrimalInfeasibilities(nLeft);
396     }
397     return chosenRow;
398}
399#if 0
400static double ft_count = 0.0;
401static double up_count = 0.0;
402static double ft_count_in = 0.0;
403static double up_count_in = 0.0;
404static int xx_count = 0;
405#endif
406/* Updates weights and returns pivot alpha.
407   Also does FT update */
408double
409ClpDualRowSteepest::updateWeights(CoinIndexedVector * input,
410                                  CoinIndexedVector * spare,
411                                  CoinIndexedVector * spare2,
412                                  CoinIndexedVector * updatedColumn)
413{
414     //#define CLP_DEBUG 3
415#if CLP_DEBUG>2
416     // Very expensive debug
417     {
418          int numberRows = model_->numberRows();
419          CoinIndexedVector * temp = new CoinIndexedVector();
420          temp->reserve(numberRows +
421                        model_->factorization()->maximumPivots());
422          double * array = alternateWeights_->denseVector();
423          int * which = alternateWeights_->getIndices();
424          alternateWeights_->clear();
425          int i;
426          for (i = 0; i < numberRows; i++) {
427               double value = 0.0;
428               array[i] = 1.0;
429               which[0] = i;
430               alternateWeights_->setNumElements(1);
431               model_->factorization()->updateColumnTranspose(temp,
432                         alternateWeights_);
433               int number = alternateWeights_->getNumElements();
434               int j;
435               for (j = 0; j < number; j++) {
436                    int iRow = which[j];
437                    value += array[iRow] * array[iRow];
438                    array[iRow] = 0.0;
439               }
440               alternateWeights_->setNumElements(0);
441               double w = CoinMax(weights_[i], value) * .1;
442               if (fabs(weights_[i] - value) > w) {
443                    printf("%d old %g, true %g\n", i, weights_[i], value);
444                    weights_[i] = value; // to reduce printout
445               }
446               //else
447               //printf("%d matches %g\n",i,value);
448          }
449          delete temp;
450     }
451#endif
452     assert (input->packedMode());
453     if (!updatedColumn->packedMode()) {
454          // I think this means empty
455#ifdef COIN_DEVELOP
456          printf("updatedColumn not packed mode ClpDualRowSteepest::updateWeights\n");
457#endif
458          return 0.0;
459     }
460     double alpha = 0.0;
461     if (!model_->factorization()->networkBasis()) {
462          // clear other region
463          alternateWeights_->clear();
464          double norm = 0.0;
465          int i;
466          double * work = input->denseVector();
467          int numberNonZero = input->getNumElements();
468          int * which = input->getIndices();
469          double * work2 = spare->denseVector();
470          int * which2 = spare->getIndices();
471          // ftran
472          //permute and move indices into index array
473          //also compute norm
474          //int *regionIndex = alternateWeights_->getIndices (  );
475          const int *permute = model_->factorization()->permute();
476          //double * region = alternateWeights_->denseVector();
477          if (permute) {
478               for ( i = 0; i < numberNonZero; i ++ ) {
479                    int iRow = which[i];
480                    double value = work[i];
481                    norm += value * value;
482                    iRow = permute[iRow];
483                    work2[iRow] = value;
484                    which2[i] = iRow;
485               }
486          } else {
487               for ( i = 0; i < numberNonZero; i ++ ) {
488                    int iRow = which[i];
489                    double value = work[i];
490                    norm += value * value;
491                    //iRow = permute[iRow];
492                    work2[iRow] = value;
493                    which2[i] = iRow;
494               }
495          }
496          spare->setNumElements ( numberNonZero );
497          // Do FT update
498#if 0
499          ft_count_in += updatedColumn->getNumElements();
500          up_count_in += spare->getNumElements();
501#endif
502          if (permute || true) {
503#if CLP_DEBUG>2
504               printf("REGION before %d els\n", spare->getNumElements());
505               spare->print();
506#endif
507               model_->factorization()->updateTwoColumnsFT(spare2, updatedColumn,
508                         spare, permute != NULL);
509#if CLP_DEBUG>2
510               printf("REGION after %d els\n", spare->getNumElements());
511               spare->print();
512#endif
513          } else {
514               // Leave as old way
515               model_->factorization()->updateColumnFT(spare2, updatedColumn);
516               model_->factorization()->updateColumn(spare2, spare, false);
517          }
518#undef CLP_DEBUG
519#if 0
520          ft_count += updatedColumn->getNumElements();
521          up_count += spare->getNumElements();
522          xx_count++;
523          if ((xx_count % 1000) == 0)
524               printf("zz %d ft %g up %g (in %g %g)\n", xx_count, ft_count, up_count,
525                      ft_count_in, up_count_in);
526#endif
527          numberNonZero = spare->getNumElements();
528          // alternateWeights_ should still be empty
529          int pivotRow = model_->pivotRow();
530#ifdef CLP_DEBUG
531          if ( model_->logLevel (  ) > 4  &&
532                    fabs(norm - weights_[pivotRow]) > 1.0e-3 * (1.0 + norm))
533               printf("on row %d, true weight %g, old %g\n",
534                      pivotRow, sqrt(norm), sqrt(weights_[pivotRow]));
535#endif
536          // could re-initialize here (could be expensive)
537          norm /= model_->alpha() * model_->alpha();
538          assert(model_->alpha());
539          assert(norm);
540          // pivot element
541          alpha = 0.0;
542          double multiplier = 2.0 / model_->alpha();
543          // look at updated column
544          work = updatedColumn->denseVector();
545          numberNonZero = updatedColumn->getNumElements();
546          which = updatedColumn->getIndices();
547
548          int nSave = 0;
549          double * work3 = alternateWeights_->denseVector();
550          int * which3 = alternateWeights_->getIndices();
551          const int * pivotColumn = model_->factorization()->pivotColumn();
552          for (i = 0; i < numberNonZero; i++) {
553               int iRow = which[i];
554               double theta = work[i];
555               if (iRow == pivotRow)
556                    alpha = theta;
557               double devex = weights_[iRow];
558               work3[nSave] = devex; // save old
559               which3[nSave++] = iRow;
560               // transform to match spare
561               int jRow = permute ? pivotColumn[iRow] : iRow;
562               double value = work2[jRow];
563               devex +=  theta * (theta * norm + value * multiplier);
564               if (devex < DEVEX_TRY_NORM)
565                    devex = DEVEX_TRY_NORM;
566               weights_[iRow] = devex;
567          }
568          alternateWeights_->setPackedMode(true);
569          alternateWeights_->setNumElements(nSave);
570          if (norm < DEVEX_TRY_NORM)
571               norm = DEVEX_TRY_NORM;
572          // Try this to make less likely will happen again and stop cycling
573          //norm *= 1.02;
574          weights_[pivotRow] = norm;
575          spare->clear();
576#ifdef CLP_DEBUG
577          spare->checkClear();
578#endif
579     } else {
580          // Do FT update
581          model_->factorization()->updateColumnFT(spare, updatedColumn);
582          // clear other region
583          alternateWeights_->clear();
584          double norm = 0.0;
585          int i;
586          double * work = input->denseVector();
587          int number = input->getNumElements();
588          int * which = input->getIndices();
589          double * work2 = spare->denseVector();
590          int * which2 = spare->getIndices();
591          for (i = 0; i < number; i++) {
592               int iRow = which[i];
593               double value = work[i];
594               norm += value * value;
595               work2[iRow] = value;
596               which2[i] = iRow;
597          }
598          spare->setNumElements(number);
599          // ftran
600#ifndef NDEBUG
601          alternateWeights_->checkClear();
602#endif
603          model_->factorization()->updateColumn(alternateWeights_, spare);
604          // alternateWeights_ should still be empty
605#ifndef NDEBUG
606          alternateWeights_->checkClear();
607#endif
608          int pivotRow = model_->pivotRow();
609#ifdef CLP_DEBUG
610          if ( model_->logLevel (  ) > 4  &&
611                    fabs(norm - weights_[pivotRow]) > 1.0e-3 * (1.0 + norm))
612               printf("on row %d, true weight %g, old %g\n",
613                      pivotRow, sqrt(norm), sqrt(weights_[pivotRow]));
614#endif
615          // could re-initialize here (could be expensive)
616          norm /= model_->alpha() * model_->alpha();
617
618          assert(norm);
619          //if (norm < DEVEX_TRY_NORM)
620          //norm = DEVEX_TRY_NORM;
621          // pivot element
622          alpha = 0.0;
623          double multiplier = 2.0 / model_->alpha();
624          // look at updated column
625          work = updatedColumn->denseVector();
626          number = updatedColumn->getNumElements();
627          which = updatedColumn->getIndices();
628
629          int nSave = 0;
630          double * work3 = alternateWeights_->denseVector();
631          int * which3 = alternateWeights_->getIndices();
632          for (i = 0; i < number; i++) {
633               int iRow = which[i];
634               double theta = work[i];
635               if (iRow == pivotRow)
636                    alpha = theta;
637               double devex = weights_[iRow];
638               work3[nSave] = devex; // save old
639               which3[nSave++] = iRow;
640               double value = work2[iRow];
641               devex +=  theta * (theta * norm + value * multiplier);
642               if (devex < DEVEX_TRY_NORM)
643                    devex = DEVEX_TRY_NORM;
644               weights_[iRow] = devex;
645          }
646          if (!alpha) {
647               // error - but carry on
648               alpha = 1.0e-50;
649          }
650          alternateWeights_->setPackedMode(true);
651          alternateWeights_->setNumElements(nSave);
652          if (norm < DEVEX_TRY_NORM)
653               norm = DEVEX_TRY_NORM;
654          weights_[pivotRow] = norm;
655          spare->clear();
656     }
657#ifdef CLP_DEBUG
658     spare->checkClear();
659#endif
660     return alpha;
661}
662
663/* Updates primal solution (and maybe list of candidates)
664   Uses input vector which it deletes
665   Computes change in objective function
666*/
667void
668ClpDualRowSteepest::updatePrimalSolution(
669     CoinIndexedVector * primalUpdate,
670     double primalRatio,
671     double & objectiveChange)
672{
673     double * COIN_RESTRICT work = primalUpdate->denseVector();
674     int number = primalUpdate->getNumElements();
675     int * COIN_RESTRICT which = primalUpdate->getIndices();
676     int i;
677     double changeObj = 0.0;
678     double tolerance = model_->currentPrimalTolerance();
679     const int * COIN_RESTRICT pivotVariable = model_->pivotVariable();
680     double * COIN_RESTRICT infeas = infeasible_->denseVector();
681     double * COIN_RESTRICT solution = model_->solutionRegion();
682     const double * COIN_RESTRICT costModel = model_->costRegion();
683     const double * COIN_RESTRICT lowerModel = model_->lowerRegion();
684     const double * COIN_RESTRICT upperModel = model_->upperRegion();
685#ifdef CLP_DUAL_COLUMN_MULTIPLIER
686     int numberColumns = model_->numberColumns();
687#endif
688     if (primalUpdate->packedMode()) {
689          for (i = 0; i < number; i++) {
690               int iRow = which[i];
691               int iPivot = pivotVariable[iRow];
692               double value = solution[iPivot];
693               double cost = costModel[iPivot];
694               double change = primalRatio * work[i];
695               work[i] = 0.0;
696               value -= change;
697               changeObj -= change * cost;
698               double lower = lowerModel[iPivot];
699               double upper = upperModel[iPivot];
700               solution[iPivot] = value;
701               if (value < lower - tolerance) {
702                    value -= lower;
703                    value *= value;
704#ifdef CLP_DUAL_COLUMN_MULTIPLIER
705                    if (iPivot < numberColumns)
706                         value *= CLP_DUAL_COLUMN_MULTIPLIER; // bias towards columns
707#endif
708#ifdef CLP_DUAL_FIXED_COLUMN_MULTIPLIER
709                    if (lower == upper)
710                         value *= CLP_DUAL_FIXED_COLUMN_MULTIPLIER; // bias towards taking out fixed variables
711#endif
712                    // store square in list
713                    if (infeas[iRow])
714                         infeas[iRow] = value; // already there
715                    else
716                         infeasible_->quickAdd(iRow, value);
717               } else if (value > upper + tolerance) {
718                    value -= upper;
719                    value *= value;
720#ifdef CLP_DUAL_COLUMN_MULTIPLIER
721                    if (iPivot < numberColumns)
722                         value *= CLP_DUAL_COLUMN_MULTIPLIER; // bias towards columns
723#endif
724#ifdef CLP_DUAL_FIXED_COLUMN_MULTIPLIER
725                    if (lower == upper)
726                         value *= CLP_DUAL_FIXED_COLUMN_MULTIPLIER; // bias towards taking out fixed variables
727#endif
728                    // store square in list
729                    if (infeas[iRow])
730                         infeas[iRow] = value; // already there
731                    else
732                         infeasible_->quickAdd(iRow, value);
733               } else {
734                    // feasible - was it infeasible - if so set tiny
735                    if (infeas[iRow])
736                         infeas[iRow] = COIN_INDEXED_REALLY_TINY_ELEMENT;
737               }
738          }
739     } else {
740          for (i = 0; i < number; i++) {
741               int iRow = which[i];
742               int iPivot = pivotVariable[iRow];
743               double value = solution[iPivot];
744               double cost = costModel[iPivot];
745               double change = primalRatio * work[iRow];
746               value -= change;
747               changeObj -= change * cost;
748               double lower = lowerModel[iPivot];
749               double upper = upperModel[iPivot];
750               solution[iPivot] = value;
751               if (value < lower - tolerance) {
752                    value -= lower;
753                    value *= value;
754#ifdef CLP_DUAL_COLUMN_MULTIPLIER
755                    if (iPivot < numberColumns)
756                         value *= CLP_DUAL_COLUMN_MULTIPLIER; // bias towards columns
757#endif
758#ifdef CLP_DUAL_FIXED_COLUMN_MULTIPLIER
759                    if (lower == upper)
760                         value *= CLP_DUAL_FIXED_COLUMN_MULTIPLIER; // bias towards taking out fixed variables
761#endif
762                    // store square in list
763                    if (infeas[iRow])
764                         infeas[iRow] = value; // already there
765                    else
766                         infeasible_->quickAdd(iRow, value);
767               } else if (value > upper + tolerance) {
768                    value -= upper;
769                    value *= value;
770#ifdef CLP_DUAL_COLUMN_MULTIPLIER
771                    if (iPivot < numberColumns)
772                         value *= CLP_DUAL_COLUMN_MULTIPLIER; // bias towards columns
773#endif
774#ifdef CLP_DUAL_FIXED_COLUMN_MULTIPLIER
775                    if (lower == upper)
776                         value *= CLP_DUAL_FIXED_COLUMN_MULTIPLIER; // bias towards taking out fixed variables
777#endif
778                    // store square in list
779                    if (infeas[iRow])
780                         infeas[iRow] = value; // already there
781                    else
782                         infeasible_->quickAdd(iRow, value);
783               } else {
784                    // feasible - was it infeasible - if so set tiny
785                    if (infeas[iRow])
786                         infeas[iRow] = COIN_INDEXED_REALLY_TINY_ELEMENT;
787               }
788               work[iRow] = 0.0;
789          }
790     }
791     // Do pivot row
792     {
793          int iRow = model_->pivotRow();
794          // feasible - was it infeasible - if so set tiny
795          //assert (infeas[iRow]);
796          if (infeas[iRow])
797               infeas[iRow] = COIN_INDEXED_REALLY_TINY_ELEMENT;
798     }
799     primalUpdate->setNumElements(0);
800     objectiveChange += changeObj;
801}
802/* Saves any weights round factorization as pivot rows may change
803   1) before factorization
804   2) after factorization
805   3) just redo infeasibilities
806   4) as 2 but restore weights from previous snapshot
807   5) for strong branching - initialize to 1 , infeasibilities
808   6) scale back
809   7) for strong branching - initialize full weights , infeasibilities
810*/
811void
812ClpDualRowSteepest::saveWeights(ClpSimplex * model, int mode)
813{
814     // alternateWeights_ is defined as indexed but is treated oddly
815     model_ = model;
816     int numberRows = model_->numberRows();
817     int numberColumns = model_->numberColumns();
818     const int * pivotVariable = model_->pivotVariable();
819     int i;
820     if (mode == 1) {
821          if(weights_) {
822               // Check if size has changed
823               if (infeasible_->capacity() == numberRows) {
824                    alternateWeights_->clear();
825                    // change from row numbers to sequence numbers
826                    int * which = alternateWeights_->getIndices();
827                    for (i = 0; i < numberRows; i++) {
828                         int iPivot = pivotVariable[i];
829                         which[i] = iPivot;
830                    }
831                    state_ = 1;
832                    // clear marker on savedWeights_
833                    assert (savedWeights_);
834                    if (savedWeights_->packedMode()) {
835                      savedWeights_->setPackedMode(false);
836                      savedWeights_->setNumElements(0);
837                    }
838               } else {
839                    // size has changed - clear everything
840                    delete [] weights_;
841                    weights_ = NULL;
842                    delete [] dubiousWeights_;
843                    dubiousWeights_ = NULL;
844                    delete infeasible_;
845                    infeasible_ = NULL;
846                    delete alternateWeights_;
847                    alternateWeights_ = NULL;
848                    delete savedWeights_;
849                    savedWeights_ = NULL;
850                    state_ = -1;
851               }
852          }
853     } else if (mode == 2 || mode == 4 || mode >= 5) {
854          // restore
855          if (!weights_ || state_ == -1 || mode == 5 || mode == 7) {
856               // initialize weights
857               delete [] weights_;
858               delete alternateWeights_;
859               weights_ = new double[numberRows];
860               // initialize to 1.0 (can we do better?)
861               for (i = 0; i < numberRows; i++) {
862                 weights_[i] = 1.0;
863               }
864               alternateWeights_ = new CoinIndexedVector();
865               // enough space so can use it for factorization
866               alternateWeights_->reserve(numberRows +
867                                          model_->factorization()->maximumPivots());
868               if (mode_ != 1 || mode == 5) {
869               } else {
870                    CoinIndexedVector * temp = new CoinIndexedVector();
871                    temp->reserve(numberRows +
872                                  model_->factorization()->maximumPivots());
873                    double * array = alternateWeights_->denseVector();
874                    int * which = alternateWeights_->getIndices();
875                    int firstRow=0;
876                    int lastRow=numberRows;
877                    if (mode==7) {
878                      // use info passed in
879                      firstRow=model->spareIntArray_[0];
880                      lastRow=model->spareIntArray_[1];
881                    }
882                    for (i = firstRow; i < lastRow; i++) {
883                         double value = 0.0;
884                         array[0] = 1.0;
885                         which[0] = i;
886                         alternateWeights_->setNumElements(1);
887                         alternateWeights_->setPackedMode(true);
888                         model_->factorization()->updateColumnTranspose(temp,
889                                   alternateWeights_);
890                         int number = alternateWeights_->getNumElements();
891                         int j;
892                         for (j = 0; j < number; j++) {
893                              value += array[j] * array[j];
894                              array[j] = 0.0;
895                         }
896                         alternateWeights_->setNumElements(0);
897                         weights_[i] = value;
898                    }
899                    delete temp;
900               }
901               // create saved weights (not really indexedvector)
902               savedWeights_ = new CoinIndexedVector();
903               savedWeights_->reserve(numberRows);
904               for (int i=0;i<model_->numberRows();i++)
905                 savedWeights_->denseVector()[i]=1.0;
906
907               double * array = savedWeights_->denseVector();
908               int * which = savedWeights_->getIndices();
909               for (int i = 0; i < numberRows; i++) {
910                    array[i] = weights_[i];
911                    which[i] = pivotVariable[i];
912               }
913               if (mode==7) {
914                 savedWeights_->setNumElements(numberRows); //flag as special
915                 savedWeights_->setPackedMode(true);
916               }
917          } else if (mode != 6) {
918               int * which = alternateWeights_->getIndices();
919               CoinIndexedVector * rowArray3 = model_->rowArray(3);
920               rowArray3->clear();
921               int * back = rowArray3->getIndices();
922               // In case something went wrong
923               for (i = 0; i < numberRows + numberColumns; i++)
924                    back[i] = -1;
925               if (mode != 4) {
926                    // save
927                    CoinMemcpyN(which,  numberRows, savedWeights_->getIndices());
928                    CoinMemcpyN(weights_,       numberRows, savedWeights_->denseVector());
929               } else {
930                    // restore
931                    //memcpy(which,savedWeights_->getIndices(),
932                    //     numberRows*sizeof(int));
933                    //memcpy(weights_,savedWeights_->denseVector(),
934                    //     numberRows*sizeof(double));
935                    which = savedWeights_->getIndices();
936               }
937               // restore (a bit slow - but only every re-factorization)
938               double * array = savedWeights_->denseVector();
939               for (i = 0; i < numberRows; i++) {
940                    int iSeq = which[i];
941                    back[iSeq] = i;
942               }
943               int firstRow=0;
944               int lastRow=numberRows;
945               if (mode==7) {
946                 // use info passed in
947                 firstRow=model->spareIntArray_[0];
948                 lastRow=model->spareIntArray_[1];
949               }
950               for (i = firstRow; i < lastRow; i++) {
951                    int iPivot = pivotVariable[i];
952                    iPivot = back[iPivot];
953                    if (iPivot >= 0) {
954                         weights_[i] = array[iPivot];
955                         if (weights_[i] < DEVEX_TRY_NORM)
956                              weights_[i] = DEVEX_TRY_NORM; // may need to check more
957                    } else {
958                         // odd
959                         weights_[i] = 1.0;
960                         //printf("bad pivot row %d (%d->%d) iPivot %d\n",i,firstRow,lastRow,iPivot);
961                    }
962               }
963#if 0
964               printf("mode %d mode_ %d state_ %d\n",mode,mode_,state_);
965               if (!model_->numberIterations()) {
966                 for (int k=0;k<numberRows;k+=10)
967                   printf("%d %g ",k,weights_[k]);
968                 printf("\n");
969               }
970#endif
971          } else {
972               // mode 6 - scale back weights as primal errors
973               double primalError = model_->largestPrimalError();
974               double allowed ;
975               if (primalError > 1.0e3)
976                    allowed = 10.0;
977               else if (primalError > 1.0e2)
978                    allowed = 50.0;
979               else if (primalError > 1.0e1)
980                    allowed = 100.0;
981               else
982                    allowed = 1000.0;
983               double allowedInv = 1.0 / allowed;
984               for (i = 0; i < numberRows; i++) {
985                    double value = weights_[i];
986                    if (value < allowedInv)
987                         value = allowedInv;
988                    else if (value > allowed)
989                         value = allowed;
990                    weights_[i] = allowed;
991               }
992          }
993          state_ = 0;
994          // set up infeasibilities
995          if (!infeasible_) {
996               infeasible_ = new CoinIndexedVector();
997               infeasible_->reserve(numberRows);
998          }
999     }
1000     if (mode >= 2) {
1001          // Get dubious weights
1002          //if (!dubiousWeights_)
1003          //dubiousWeights_=new int[numberRows];
1004          //model_->factorization()->getWeights(dubiousWeights_);
1005          infeasible_->clear();
1006          int iRow;
1007          const int * pivotVariable = model_->pivotVariable();
1008          double tolerance = model_->currentPrimalTolerance();
1009          for (iRow = 0; iRow < numberRows; iRow++) {
1010               int iPivot = pivotVariable[iRow];
1011               double value = model_->solution(iPivot);
1012               double lower = model_->lower(iPivot);
1013               double upper = model_->upper(iPivot);
1014               if (value < lower - tolerance) {
1015                    value -= lower;
1016                    value *= value;
1017#ifdef CLP_DUAL_COLUMN_MULTIPLIER
1018                    if (iPivot < numberColumns)
1019                         value *= CLP_DUAL_COLUMN_MULTIPLIER; // bias towards columns
1020#endif
1021#ifdef CLP_DUAL_FIXED_COLUMN_MULTIPLIER
1022                    if (lower == upper)
1023                         value *= CLP_DUAL_FIXED_COLUMN_MULTIPLIER; // bias towards taking out fixed variables
1024#endif
1025                    // store square in list
1026                    infeasible_->quickAdd(iRow, value);
1027               } else if (value > upper + tolerance) {
1028                    value -= upper;
1029                    value *= value;
1030#ifdef CLP_DUAL_COLUMN_MULTIPLIER
1031                    if (iPivot < numberColumns)
1032                         value *= CLP_DUAL_COLUMN_MULTIPLIER; // bias towards columns
1033#endif
1034#ifdef CLP_DUAL_FIXED_COLUMN_MULTIPLIER
1035                    if (lower == upper)
1036                         value *= CLP_DUAL_FIXED_COLUMN_MULTIPLIER; // bias towards taking out fixed variables
1037#endif
1038                    // store square in list
1039                    infeasible_->quickAdd(iRow, value);
1040               }
1041          }
1042     }
1043     // see where coming from
1044     if (mode==2&&!model->numberIterations()) {
1045       int options=model->specialOptions();
1046       if ((options&16384)!=0&&true) {
1047         // fast of some sort - could be clever???
1048         // for now initialize
1049         if ((options&524288)!=0&&false) {
1050           // fathom
1051           for (int i = 0; i < numberRows; i++) 
1052             weights_[i] = 1.0;
1053         } else if (true) {
1054           // strong
1055           for (int i = 0; i < numberRows; i++) 
1056             weights_[i] = 1.0;
1057         }
1058       }
1059     }
1060}
1061// Pass in saved weights
1062void 
1063ClpDualRowSteepest::passInSavedWeights(const CoinIndexedVector * saved)
1064{
1065  delete savedWeights_;
1066  savedWeights_=new CoinIndexedVector(*saved);
1067}
1068// Gets rid of last update
1069void
1070ClpDualRowSteepest::unrollWeights()
1071{
1072     double * saved = alternateWeights_->denseVector();
1073     int number = alternateWeights_->getNumElements();
1074     int * which = alternateWeights_->getIndices();
1075     int i;
1076     if (alternateWeights_->packedMode()) {
1077          for (i = 0; i < number; i++) {
1078               int iRow = which[i];
1079               weights_[iRow] = saved[i];
1080               saved[i] = 0.0;
1081          }
1082     } else {
1083          for (i = 0; i < number; i++) {
1084               int iRow = which[i];
1085               weights_[iRow] = saved[iRow];
1086               saved[iRow] = 0.0;
1087          }
1088     }
1089     alternateWeights_->setNumElements(0);
1090}
1091//-------------------------------------------------------------------
1092// Clone
1093//-------------------------------------------------------------------
1094ClpDualRowPivot * ClpDualRowSteepest::clone(bool CopyData) const
1095{
1096     if (CopyData) {
1097          return new ClpDualRowSteepest(*this);
1098     } else {
1099          return new ClpDualRowSteepest();
1100     }
1101}
1102// Gets rid of all arrays
1103void
1104ClpDualRowSteepest::clearArrays()
1105{
1106     if (persistence_ == normal) {
1107          delete [] weights_;
1108          weights_ = NULL;
1109          delete [] dubiousWeights_;
1110          dubiousWeights_ = NULL;
1111          delete infeasible_;
1112          infeasible_ = NULL;
1113          delete alternateWeights_;
1114          alternateWeights_ = NULL;
1115          delete savedWeights_;
1116          savedWeights_ = NULL;
1117     }
1118     state_ = -1;
1119}
1120// Returns true if would not find any row
1121bool
1122ClpDualRowSteepest::looksOptimal() const
1123{
1124     int iRow;
1125     const int * pivotVariable = model_->pivotVariable();
1126     double tolerance = model_->currentPrimalTolerance();
1127     // we can't really trust infeasibilities if there is primal error
1128     // this coding has to mimic coding in checkPrimalSolution
1129     double error = CoinMin(1.0e-2, model_->largestPrimalError());
1130     // allow tolerance at least slightly bigger than standard
1131     tolerance = tolerance +  error;
1132     // But cap
1133     tolerance = CoinMin(1000.0, tolerance);
1134     int numberRows = model_->numberRows();
1135     int numberInfeasible = 0;
1136     for (iRow = 0; iRow < numberRows; iRow++) {
1137          int iPivot = pivotVariable[iRow];
1138          double value = model_->solution(iPivot);
1139          double lower = model_->lower(iPivot);
1140          double upper = model_->upper(iPivot);
1141          if (value < lower - tolerance) {
1142               numberInfeasible++;
1143          } else if (value > upper + tolerance) {
1144               numberInfeasible++;
1145          }
1146     }
1147     return (numberInfeasible == 0);
1148}
1149// Called when maximum pivots changes
1150void
1151ClpDualRowSteepest::maximumPivotsChanged()
1152{
1153     if (alternateWeights_ &&
1154               alternateWeights_->capacity() != model_->numberRows() +
1155               model_->factorization()->maximumPivots()) {
1156          delete alternateWeights_;
1157          alternateWeights_ = new CoinIndexedVector();
1158          // enough space so can use it for factorization
1159          alternateWeights_->reserve(model_->numberRows() +
1160                                     model_->factorization()->maximumPivots());
1161     }
1162}
1163
Note: See TracBrowser for help on using the repository browser.