source: trunk/Clp/src/ClpNode.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:executable set to *
  • Property svn:keywords set to Id
File size: 46.2 KB
Line 
1/* $Id: ClpNode.cpp 1665 2011-01-04 17:55:54Z lou $ */
2// Copyright (C) 2008, 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 "ClpNode.hpp"
9#include "ClpFactorization.hpp"
10#include "ClpDualRowSteepest.hpp"
11
12//#############################################################################
13// Constructors / Destructor / Assignment
14//#############################################################################
15
16//-------------------------------------------------------------------
17// Default Constructor
18//-------------------------------------------------------------------
19ClpNode::ClpNode () :
20     branchingValue_(0.5),
21     objectiveValue_(0.0),
22     sumInfeasibilities_(0.0),
23     estimatedSolution_(0.0),
24     factorization_(NULL),
25     weights_(NULL),
26     status_(NULL),
27     primalSolution_(NULL),
28     dualSolution_(NULL),
29     lower_(NULL),
30     upper_(NULL),
31     pivotVariables_(NULL),
32     fixed_(NULL),
33     sequence_(1),
34     numberInfeasibilities_(0),
35     depth_(0),
36     numberFixed_(0),
37     flags_(0),
38     maximumFixed_(0),
39     maximumRows_(0),
40     maximumColumns_(0),
41     maximumIntegers_(0)
42{
43     branchState_.firstBranch = 0;
44     branchState_.branch = 0;
45}
46//-------------------------------------------------------------------
47// Useful Constructor from model
48//-------------------------------------------------------------------
49ClpNode::ClpNode (ClpSimplex * model, const ClpNodeStuff * stuff, int depth) :
50     branchingValue_(0.5),
51     objectiveValue_(0.0),
52     sumInfeasibilities_(0.0),
53     estimatedSolution_(0.0),
54     factorization_(NULL),
55     weights_(NULL),
56     status_(NULL),
57     primalSolution_(NULL),
58     dualSolution_(NULL),
59     lower_(NULL),
60     upper_(NULL),
61     pivotVariables_(NULL),
62     fixed_(NULL),
63     sequence_(1),
64     numberInfeasibilities_(0),
65     depth_(0),
66     numberFixed_(0),
67     flags_(0),
68     maximumFixed_(0),
69     maximumRows_(0),
70     maximumColumns_(0),
71     maximumIntegers_(0)
72{
73     branchState_.firstBranch = 0;
74     branchState_.branch = 0;
75     gutsOfConstructor(model, stuff, 0, depth);
76}
77
78//-------------------------------------------------------------------
79// Most of work of constructor from model
80//-------------------------------------------------------------------
81void
82ClpNode::gutsOfConstructor (ClpSimplex * model, const ClpNodeStuff * stuff,
83                            int arraysExist, int depth)
84{
85     int numberRows = model->numberRows();
86     int numberColumns = model->numberColumns();
87     int numberTotal = numberRows + numberColumns;
88     int maximumTotal = maximumRows_ + maximumColumns_;
89     depth_ = depth;
90     // save stuff
91     objectiveValue_ = model->objectiveValue() * model->optimizationDirection();
92     estimatedSolution_ = objectiveValue_;
93     flags_ = 1; //say scaled
94     if (!arraysExist) {
95          maximumRows_ = CoinMax(maximumRows_, numberRows);
96          maximumColumns_ = CoinMax(maximumColumns_, numberColumns);
97          maximumTotal = maximumRows_ + maximumColumns_;
98          assert (!factorization_);
99          factorization_ = new ClpFactorization(*model->factorization(), numberRows);
100          status_ = CoinCopyOfArrayPartial(model->statusArray(), maximumTotal, numberTotal);
101          primalSolution_ = CoinCopyOfArrayPartial(model->solutionRegion(), maximumTotal, numberTotal);
102          dualSolution_ = CoinCopyOfArrayPartial(model->djRegion(), maximumTotal, numberTotal); //? has duals as well?
103          pivotVariables_ = CoinCopyOfArrayPartial(model->pivotVariable(), maximumRows_, numberRows);
104          ClpDualRowSteepest* pivot =
105               dynamic_cast< ClpDualRowSteepest*>(model->dualRowPivot());
106          if (pivot) {
107               assert (!weights_);
108               weights_ = new ClpDualRowSteepest(*pivot);
109          }
110     } else {
111          if (arraysExist == 2)
112               assert(lower_);
113          if (numberRows <= maximumRows_ && numberColumns <= maximumColumns_) {
114               CoinMemcpyN(model->statusArray(), numberTotal, status_);
115               if (arraysExist == 1) {
116                    *factorization_ = *model->factorization();
117                    CoinMemcpyN(model->solutionRegion(), numberTotal, primalSolution_);
118                    CoinMemcpyN(model->djRegion(), numberTotal, dualSolution_); //? has duals as well?
119                    ClpDualRowSteepest* pivot =
120                         dynamic_cast< ClpDualRowSteepest*>(model->dualRowPivot());
121                    if (pivot) {
122                         if (weights_) {
123                              //if (weights_->numberRows()==pivot->numberRows()) {
124                              weights_->fill(*pivot);
125                              //} else {
126                              //delete weights_;
127                              //weights_ = new ClpDualRowSteepest(*pivot);
128                              //}
129                         } else {
130                              weights_ = new ClpDualRowSteepest(*pivot);
131                         }
132                    }
133                    CoinMemcpyN(model->pivotVariable(), numberRows, pivotVariables_);
134               } else {
135                    CoinMemcpyN(model->primalColumnSolution(), numberColumns, primalSolution_);
136                    CoinMemcpyN(model->dualColumnSolution(), numberColumns, dualSolution_);
137                    flags_ = 0;
138                    CoinMemcpyN(model->dualRowSolution(), numberRows, dualSolution_ + numberColumns);
139               }
140          } else {
141               // size has changed
142               maximumRows_ = CoinMax(maximumRows_, numberRows);
143               maximumColumns_ = CoinMax(maximumColumns_, numberColumns);
144               maximumTotal = maximumRows_ + maximumColumns_;
145               delete weights_;
146               weights_ = NULL;
147               delete [] status_;
148               delete [] primalSolution_;
149               delete [] dualSolution_;
150               delete [] pivotVariables_;
151               status_ = CoinCopyOfArrayPartial(model->statusArray(), maximumTotal, numberTotal);
152               primalSolution_ = new double [maximumTotal*sizeof(double)];
153               dualSolution_ = new double [maximumTotal*sizeof(double)];
154               if (arraysExist == 1) {
155                    *factorization_ = *model->factorization(); // I think this is OK
156                    CoinMemcpyN(model->solutionRegion(), numberTotal, primalSolution_);
157                    CoinMemcpyN(model->djRegion(), numberTotal, dualSolution_); //? has duals as well?
158                    ClpDualRowSteepest* pivot =
159                         dynamic_cast< ClpDualRowSteepest*>(model->dualRowPivot());
160                    if (pivot) {
161                         assert (!weights_);
162                         weights_ = new ClpDualRowSteepest(*pivot);
163                    }
164               } else {
165                    CoinMemcpyN(model->primalColumnSolution(), numberColumns, primalSolution_);
166                    CoinMemcpyN(model->dualColumnSolution(), numberColumns, dualSolution_);
167                    flags_ = 0;
168                    CoinMemcpyN(model->dualRowSolution(), numberRows, dualSolution_ + numberColumns);
169               }
170               pivotVariables_ = new int [maximumRows_];
171               if (model->pivotVariable() && model->numberRows() == numberRows)
172                    CoinMemcpyN(model->pivotVariable(), numberRows, pivotVariables_);
173               else
174                    CoinFillN(pivotVariables_, numberRows, -1);
175          }
176     }
177     numberFixed_ = 0;
178     const double * lower = model->columnLower();
179     const double * upper = model->columnUpper();
180     const double * solution = model->primalColumnSolution();
181     const char * integerType = model->integerInformation();
182     const double * columnScale = model->columnScale();
183     if (!flags_)
184          columnScale = NULL; // as duals correct
185     int iColumn;
186     sequence_ = -1;
187     double integerTolerance = stuff->integerTolerance_;
188     double mostAway = 0.0;
189     int bestPriority = COIN_INT_MAX;
190     sumInfeasibilities_ = 0.0;
191     numberInfeasibilities_ = 0;
192     int nFix = 0;
193     double gap = CoinMax(model->dualObjectiveLimit() - objectiveValue_, 1.0e-4);
194#define PSEUDO 3
195#if PSEUDO==1||PSEUDO==2
196     // Column copy of matrix
197     ClpPackedMatrix * matrix = model->clpScaledMatrix();
198     const double *objective = model->costRegion() ;
199     if (!objective) {
200          objective = model->objective();
201          //if (!matrix)
202          matrix = dynamic_cast< ClpPackedMatrix*>(model->clpMatrix());
203     } else if (!matrix) {
204          matrix = dynamic_cast< ClpPackedMatrix*>(model->clpMatrix());
205     }
206     const double * element = matrix->getElements();
207     const int * row = matrix->getIndices();
208     const CoinBigIndex * columnStart = matrix->getVectorStarts();
209     const int * columnLength = matrix->getVectorLengths();
210     double direction = model->optimizationDirection();
211     const double * dual = dualSolution_ + numberColumns;
212#if PSEUDO==2
213     double * activeWeight = new double [numberRows];
214     const double * rowLower = model->rowLower();
215     const double * rowUpper = model->rowUpper();
216     const double * rowActivity = model->primalRowSolution();
217     double tolerance = 1.0e-6;
218     for (int iRow = 0; iRow < numberRows; iRow++) {
219          // could use pi to see if active or activity
220          if (rowActivity[iRow] > rowUpper[iRow] - tolerance
221                    || rowActivity[iRow] < rowLower[iRow] + tolerance) {
222               activeWeight[iRow] = 0.0;
223          } else {
224               activeWeight[iRow] = -1.0;
225          }
226     }
227     for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
228          if (integerType[iColumn]) {
229               double value = solution[iColumn];
230               if (fabs(value - floor(value + 0.5)) > 1.0e-6) {
231                    CoinBigIndex start = columnStart[iColumn];
232                    CoinBigIndex end = start + columnLength[iColumn];
233                    for (CoinBigIndex j = start; j < end; j++) {
234                         int iRow = row[j];
235                         if (activeWeight[iRow] >= 0.0)
236                              activeWeight[iRow] += 1.0;
237                    }
238               }
239          }
240     }
241     for (int iRow = 0; iRow < numberRows; iRow++) {
242          if (activeWeight[iRow] > 0.0) {
243               // could use pi
244               activeWeight[iRow] = 1.0 / activeWeight[iRow];
245          } else {
246               activeWeight[iRow] = 0.0;
247          }
248     }
249#endif
250#endif
251     const double * downPseudo = stuff->downPseudo_;
252     const int * numberDown = stuff->numberDown_;
253     const int * numberDownInfeasible = stuff->numberDownInfeasible_;
254     const double * upPseudo = stuff->upPseudo_;
255     const int * priority = stuff->priority_;
256     const int * numberUp = stuff->numberUp_;
257     const int * numberUpInfeasible = stuff->numberUpInfeasible_;
258     int numberBeforeTrust = stuff->numberBeforeTrust_;
259     int stateOfSearch = stuff->stateOfSearch_;
260     int iInteger = 0;
261     // weight at 1.0 is max min (CbcBranch was 0.8,0.1) (ClpNode was 0.9,0.9)
262#define WEIGHT_AFTER 0.9
263#define WEIGHT_BEFORE 0.2
264     //Stolen from Constraint Integer Programming book (with epsilon change)
265#define WEIGHT_PRODUCT
266#ifdef WEIGHT_PRODUCT
267     double smallChange = stuff->smallChange_;
268#endif
269     for (iColumn = 0; iColumn < numberColumns; iColumn++) {
270          if (integerType[iColumn]) {
271               double value = solution[iColumn];
272               value = CoinMax(value, static_cast<double> (lower[iColumn]));
273               value = CoinMin(value, static_cast<double> (upper[iColumn]));
274               double nearest = floor(value + 0.5);
275               if (fabs(value - nearest) > integerTolerance) {
276                    numberInfeasibilities_++;
277                    sumInfeasibilities_ += fabs(value - nearest);
278#if PSEUDO==1 || PSEUDO ==2
279                    double upValue = 0.0;
280                    double downValue = 0.0;
281                    double value2 = direction * objective[iColumn];
282                    //double dj2=value2;
283                    if (value2) {
284                         if (value2 > 0.0)
285                              upValue += 1.5 * value2;
286                         else
287                              downValue -= 1.5 * value2;
288                    }
289                    CoinBigIndex start = columnStart[iColumn];
290                    CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
291                    for (CoinBigIndex j = start; j < end; j++) {
292                         int iRow = row[j];
293                         value2 = -dual[iRow];
294                         if (value2) {
295                              value2 *= element[j];
296                              //dj2 += value2;
297#if PSEUDO==2
298                              assert (activeWeight[iRow] > 0.0 || fabs(dual[iRow]) < 1.0e-6);
299                              value2 *= activeWeight[iRow];
300#endif
301                              if (value2 > 0.0)
302                                   upValue += value2;
303                              else
304                                   downValue -= value2;
305                         }
306                    }
307                    //assert (fabs(dj2)<1.0e-4);
308                    int nUp = numberUp[iInteger];
309                    double upValue2 = (upPseudo[iInteger] / (1.0 + nUp));
310                    // Extra for infeasible branches
311                    if (nUp) {
312                         double ratio = 1.0 + static_cast<double>(numberUpInfeasible[iInteger]) /
313                                        static_cast<double>(nUp);
314                         upValue2 *= ratio;
315                    }
316                    int nDown = numberDown[iInteger];
317                    double downValue2 = (downPseudo[iInteger] / (1.0 + nDown));
318                    if (nDown) {
319                         double ratio = 1.0 + static_cast<double>(numberDownInfeasible[iInteger]) /
320                                        static_cast<double>(nDown);
321                         downValue2 *= ratio;
322                    }
323                    //printf("col %d - downPi %g up %g, downPs %g up %g\n",
324                    //     iColumn,upValue,downValue,upValue2,downValue2);
325                    upValue = CoinMax(0.1 * upValue, upValue2);
326                    downValue = CoinMax(0.1 * downValue, downValue2);
327                    //upValue = CoinMax(upValue,1.0e-8);
328                    //downValue = CoinMax(downValue,1.0e-8);
329                    upValue *= ceil(value) - value;
330                    downValue *= value - floor(value);
331                    double infeasibility;
332                    //if (depth>1000)
333                    //infeasibility = CoinMax(upValue,downValue)+integerTolerance;
334                    //else
335                    if (stateOfSearch <= 2) {
336                         // no solution
337                         infeasibility = (1.0 - WEIGHT_BEFORE) * CoinMax(upValue, downValue) +
338                                         WEIGHT_BEFORE * CoinMin(upValue, downValue) + integerTolerance;
339                    } else {
340#ifndef WEIGHT_PRODUCT
341                         infeasibility = (1.0 - WEIGHT_AFTER) * CoinMax(upValue, downValue) +
342                                         WEIGHT_AFTER * CoinMin(upValue, downValue) + integerTolerance;
343#else
344                         infeasibility = CoinMax(CoinMax(upValue, downValue), smallChange) *
345                                         CoinMax(CoinMin(upValue, downValue), smallChange);
346#endif
347                    }
348                    estimatedSolution_ += CoinMin(upValue2, downValue2);
349#elif PSEUDO==3
350                    int nUp = numberUp[iInteger];
351                    int nDown = numberDown[iInteger];
352                    // Extra 100% for infeasible branches
353                    double upValue = (ceil(value) - value) * (upPseudo[iInteger] /
354                                     (1.0 + nUp));
355                    if (nUp) {
356                         double ratio = 1.0 + static_cast<double>(numberUpInfeasible[iInteger]) /
357                                        static_cast<double>(nUp);
358                         upValue *= ratio;
359                    }
360                    double downValue = (value - floor(value)) * (downPseudo[iInteger] /
361                                       (1.0 + nDown));
362                    if (nDown) {
363                         double ratio = 1.0 + static_cast<double>(numberDownInfeasible[iInteger]) /
364                                        static_cast<double>(nDown);
365                         downValue *= ratio;
366                    }
367                    if (nUp < numberBeforeTrust || nDown < numberBeforeTrust) {
368                         upValue *= 10.0;
369                         downValue *= 10.0;
370                    }
371
372                    double infeasibility;
373                    //if (depth>1000)
374                    //infeasibility = CoinMax(upValue,downValue)+integerTolerance;
375                    //else
376                    if (stateOfSearch <= 2) {
377                         // no solution
378                         infeasibility = (1.0 - WEIGHT_BEFORE) * CoinMax(upValue, downValue) +
379                                         WEIGHT_BEFORE * CoinMin(upValue, downValue) + integerTolerance;
380                    } else {
381#ifndef WEIGHT_PRODUCT
382                         infeasibility = (1.0 - WEIGHT_AFTER) * CoinMax(upValue, downValue) +
383                                         WEIGHT_AFTER * CoinMin(upValue, downValue) + integerTolerance;
384#else
385                    infeasibility = CoinMax(CoinMax(upValue, downValue), smallChange) *
386                                    CoinMax(CoinMin(upValue, downValue), smallChange);
387                    //infeasibility += CoinMin(upValue,downValue)*smallChange;
388#endif
389                    }
390                    //infeasibility = 0.1*CoinMax(upValue,downValue)+
391                    //0.9*CoinMin(upValue,downValue) + integerTolerance;
392                    estimatedSolution_ += CoinMin(upValue, downValue);
393#else
394                    double infeasibility = fabs(value - nearest);
395#endif
396                    assert (infeasibility > 0.0);
397                    if (priority[iInteger] < bestPriority) {
398                         mostAway = 0.0;
399                         bestPriority = priority[iInteger];
400                    } else if (priority[iInteger] > bestPriority) {
401                         infeasibility = 0.0;
402                    }
403                    if (infeasibility > mostAway) {
404                         mostAway = infeasibility;
405                         sequence_ = iColumn;
406                         branchingValue_ = value;
407                         branchState_.branch = 0;
408#if PSEUDO>0
409                         if (upValue <= downValue)
410                              branchState_.firstBranch = 1; // up
411                         else
412                              branchState_.firstBranch = 0; // down
413#else
414                         if (value <= nearest)
415                              branchState_.firstBranch = 1; // up
416                         else
417                              branchState_.firstBranch = 0; // down
418#endif
419                    }
420               } else if (model->getColumnStatus(iColumn) == ClpSimplex::atLowerBound) {
421                    bool fix = false;
422                    if (columnScale) {
423                         if (dualSolution_[iColumn] > gap * columnScale[iColumn])
424                              fix = true;
425                    } else {
426                         if (dualSolution_[iColumn] > gap)
427                              fix = true;
428                    }
429                    if (fix) {
430                         nFix++;
431                         //printf("fixed %d to zero gap %g dj %g %g\n",iColumn,
432                         // gap,dualSolution_[iColumn], columnScale ? columnScale[iColumn]:1.0);
433                         model->setColumnStatus(iColumn, ClpSimplex::isFixed);
434                    }
435               } else if (model->getColumnStatus(iColumn) == ClpSimplex::atUpperBound) {
436                    bool fix = false;
437                    if (columnScale) {
438                         if (-dualSolution_[iColumn] > gap * columnScale[iColumn])
439                              fix = true;
440                    } else {
441                         if (-dualSolution_[iColumn] > gap)
442                              fix = true;
443                    }
444                    if (fix) {
445                         nFix++;
446                         //printf("fixed %d to one gap %g dj %g %g\n",iColumn,
447                         // gap,dualSolution_[iColumn], columnScale ? columnScale[iColumn]:1.0);
448                         model->setColumnStatus(iColumn, ClpSimplex::isFixed);
449                    }
450               }
451               iInteger++;
452          }
453     }
454     //printf("Choosing %d inf %g pri %d\n",
455     // sequence_,mostAway,bestPriority);
456#if PSEUDO == 2
457     delete [] activeWeight;
458#endif
459     if (lower_) {
460          // save bounds
461          if (iInteger > maximumIntegers_) {
462               delete [] lower_;
463               delete [] upper_;
464               maximumIntegers_ = iInteger;
465               lower_ = new int [maximumIntegers_];
466               upper_ = new int [maximumIntegers_];
467          }
468          iInteger = 0;
469          for (iColumn = 0; iColumn < numberColumns; iColumn++) {
470               if (integerType[iColumn]) {
471                    lower_[iInteger] = static_cast<int> (lower[iColumn]);
472                    upper_[iInteger] = static_cast<int> (upper[iColumn]);
473                    iInteger++;
474               }
475          }
476     }
477     // Could omit save of fixed if doing full save of bounds
478     if (sequence_ >= 0 && nFix) {
479          if (nFix > maximumFixed_) {
480               delete [] fixed_;
481               fixed_ = new int [nFix];
482               maximumFixed_ = nFix;
483          }
484          numberFixed_ = 0;
485          unsigned char * status = model->statusArray();
486          for (iColumn = 0; iColumn < numberColumns; iColumn++) {
487               if (status[iColumn] != status_[iColumn]) {
488                    if (solution[iColumn] <= lower[iColumn] + 2.0 * integerTolerance) {
489                         model->setColumnUpper(iColumn, lower[iColumn]);
490                         fixed_[numberFixed_++] = iColumn;
491                    } else {
492                         assert (solution[iColumn] >= upper[iColumn] - 2.0 * integerTolerance);
493                         model->setColumnLower(iColumn, upper[iColumn]);
494                         fixed_[numberFixed_++] = iColumn | 0x10000000;
495                    }
496               }
497          }
498          //printf("%d fixed\n",numberFixed_);
499     }
500}
501
502//-------------------------------------------------------------------
503// Copy constructor
504//-------------------------------------------------------------------
505ClpNode::ClpNode (const ClpNode & )
506{
507     printf("ClpNode copy not implemented\n");
508     abort();
509}
510
511//-------------------------------------------------------------------
512// Destructor
513//-------------------------------------------------------------------
514ClpNode::~ClpNode ()
515{
516     delete factorization_;
517     delete weights_;
518     delete [] status_;
519     delete [] primalSolution_;
520     delete [] dualSolution_;
521     delete [] lower_;
522     delete [] upper_;
523     delete [] pivotVariables_;
524     delete [] fixed_;
525}
526
527//----------------------------------------------------------------
528// Assignment operator
529//-------------------------------------------------------------------
530ClpNode &
531ClpNode::operator=(const ClpNode& rhs)
532{
533     if (this != &rhs) {
534          printf("ClpNode = not implemented\n");
535          abort();
536     }
537     return *this;
538}
539// Create odd arrays
540void
541ClpNode::createArrays(ClpSimplex * model)
542{
543     int numberColumns = model->numberColumns();
544     const char * integerType = model->integerInformation();
545     int iColumn;
546     int numberIntegers = 0;
547     for (iColumn = 0; iColumn < numberColumns; iColumn++) {
548          if (integerType[iColumn])
549               numberIntegers++;
550     }
551     if (numberIntegers > maximumIntegers_ || !lower_) {
552          delete [] lower_;
553          delete [] upper_;
554          maximumIntegers_ = numberIntegers;
555          lower_ = new int [numberIntegers];
556          upper_ = new int [numberIntegers];
557     }
558}
559// Clean up as crunch is different model
560void
561ClpNode::cleanUpForCrunch()
562{
563     delete weights_;
564     weights_ = NULL;
565}
566/* Applies node to model
567   0 - just tree bounds
568   1 - tree bounds and basis etc
569   2 - saved bounds and basis etc
570*/
571void
572ClpNode::applyNode(ClpSimplex * model, int doBoundsEtc )
573{
574     int numberColumns = model->numberColumns();
575     const double * lower = model->columnLower();
576     const double * upper = model->columnUpper();
577     if (doBoundsEtc < 2) {
578          // current bound
579          int way = branchState_.firstBranch;
580          if (branchState_.branch > 0)
581               way = 1 - way;
582          if (!way) {
583               // This should also do underlying internal bound
584               model->setColumnUpper(sequence_, floor(branchingValue_));
585          } else {
586               // This should also do underlying internal bound
587               model->setColumnLower(sequence_, ceil(branchingValue_));
588          }
589          // apply dj fixings
590          for (int i = 0; i < numberFixed_; i++) {
591               int iColumn = fixed_[i];
592               if ((iColumn & 0x10000000) != 0) {
593                    iColumn &= 0xfffffff;
594                    model->setColumnLower(iColumn, upper[iColumn]);
595               } else {
596                    model->setColumnUpper(iColumn, lower[iColumn]);
597               }
598          }
599     } else {
600          // restore bounds
601          assert (lower_);
602          int iInteger = -1;
603          const char * integerType = model->integerInformation();
604          for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
605               if (integerType[iColumn]) {
606                    iInteger++;
607                    if (lower_[iInteger] != static_cast<int> (lower[iColumn]))
608                         model->setColumnLower(iColumn, lower_[iInteger]);
609                    if (upper_[iInteger] != static_cast<int> (upper[iColumn]))
610                         model->setColumnUpper(iColumn, upper_[iInteger]);
611               }
612          }
613     }
614     if (doBoundsEtc && doBoundsEtc < 3) {
615          //model->copyFactorization(*factorization_);
616          model->copyFactorization(*factorization_);
617          ClpDualRowSteepest* pivot =
618               dynamic_cast< ClpDualRowSteepest*>(model->dualRowPivot());
619          if (pivot && weights_) {
620               pivot->fill(*weights_);
621          }
622          int numberRows = model->numberRows();
623          int numberTotal = numberRows + numberColumns;
624          CoinMemcpyN(status_, numberTotal, model->statusArray());
625          if (doBoundsEtc < 2) {
626               CoinMemcpyN(primalSolution_, numberTotal, model->solutionRegion());
627               CoinMemcpyN(dualSolution_, numberTotal, model->djRegion());
628               CoinMemcpyN(pivotVariables_, numberRows, model->pivotVariable());
629               CoinMemcpyN(dualSolution_ + numberColumns, numberRows, model->dualRowSolution());
630          } else {
631               CoinMemcpyN(primalSolution_, numberColumns, model->primalColumnSolution());
632               CoinMemcpyN(dualSolution_, numberColumns, model->dualColumnSolution());
633               CoinMemcpyN(dualSolution_ + numberColumns, numberRows, model->dualRowSolution());
634               if (model->columnScale()) {
635                    // See if just primal will work
636                    double * solution = model->primalColumnSolution();
637                    const double * columnScale = model->columnScale();
638                    int i;
639                    for (i = 0; i < numberColumns; i++) {
640                         solution[i] *= columnScale[i];
641                    }
642               }
643          }
644          model->setObjectiveValue(objectiveValue_);
645     }
646}
647// Choose a new variable
648void
649ClpNode::chooseVariable(ClpSimplex * , ClpNodeStuff * /*info*/)
650{
651#if 0
652     int way = branchState_.firstBranch;
653     if (branchState_.branch > 0)
654          way = 1 - way;
655     assert (!branchState_.branch);
656     // We need to use pseudo costs to choose a variable
657     int numberColumns = model->numberColumns();
658#endif
659}
660// Fix on reduced costs
661int
662ClpNode::fixOnReducedCosts(ClpSimplex * )
663{
664
665     return 0;
666}
667/* Way for integer variable -1 down , +1 up */
668int
669ClpNode::way() const
670{
671     int way = branchState_.firstBranch;
672     if (branchState_.branch > 0)
673          way = 1 - way;
674     return way == 0 ? -1 : +1;
675}
676// Return true if branch exhausted
677bool
678ClpNode::fathomed() const
679{
680     return branchState_.branch >= 1
681            ;
682}
683// Change state of variable i.e. go other way
684void
685ClpNode::changeState()
686{
687     branchState_.branch++;
688     assert (branchState_.branch <= 2);
689}
690//#############################################################################
691// Constructors / Destructor / Assignment
692//#############################################################################
693
694//-------------------------------------------------------------------
695// Default Constructor
696//-------------------------------------------------------------------
697ClpNodeStuff::ClpNodeStuff () :
698     integerTolerance_(1.0e-7),
699     integerIncrement_(1.0e-8),
700     smallChange_(1.0e-8),
701     downPseudo_(NULL),
702     upPseudo_(NULL),
703     priority_(NULL),
704     numberDown_(NULL),
705     numberUp_(NULL),
706     numberDownInfeasible_(NULL),
707     numberUpInfeasible_(NULL),
708     saveCosts_(NULL),
709     nodeInfo_(NULL),
710     large_(NULL),
711     whichRow_(NULL),
712     whichColumn_(NULL),
713     nBound_(0),
714     saveOptions_(0),
715     solverOptions_(0),
716     maximumNodes_(0),
717     numberBeforeTrust_(0),
718     stateOfSearch_(0),
719     nDepth_(-1),
720     nNodes_(0),
721     numberNodesExplored_(0),
722     numberIterations_(0),
723     presolveType_(0)
724{
725
726}
727
728//-------------------------------------------------------------------
729// Copy constructor
730//-------------------------------------------------------------------
731ClpNodeStuff::ClpNodeStuff (const ClpNodeStuff & rhs)
732     : integerTolerance_(rhs.integerTolerance_),
733       integerIncrement_(rhs.integerIncrement_),
734       smallChange_(rhs.smallChange_),
735       downPseudo_(NULL),
736       upPseudo_(NULL),
737       priority_(NULL),
738       numberDown_(NULL),
739       numberUp_(NULL),
740       numberDownInfeasible_(NULL),
741       numberUpInfeasible_(NULL),
742       saveCosts_(NULL),
743       nodeInfo_(NULL),
744       large_(NULL),
745       whichRow_(NULL),
746       whichColumn_(NULL),
747       nBound_(0),
748       saveOptions_(rhs.saveOptions_),
749       solverOptions_(rhs.solverOptions_),
750       maximumNodes_(rhs.maximumNodes_),
751       numberBeforeTrust_(rhs.numberBeforeTrust_),
752       stateOfSearch_(rhs.stateOfSearch_),
753       nDepth_(rhs.nDepth_),
754       nNodes_(rhs.nNodes_),
755       numberNodesExplored_(rhs.numberNodesExplored_),
756       numberIterations_(rhs.numberIterations_),
757       presolveType_(rhs.presolveType_)
758{
759}
760//----------------------------------------------------------------
761// Assignment operator
762//-------------------------------------------------------------------
763ClpNodeStuff &
764ClpNodeStuff::operator=(const ClpNodeStuff& rhs)
765{
766     if (this != &rhs) {
767          integerTolerance_ = rhs.integerTolerance_;
768          integerIncrement_ = rhs.integerIncrement_;
769          smallChange_ = rhs.smallChange_;
770          downPseudo_ = NULL;
771          upPseudo_ = NULL;
772          priority_ = NULL;
773          numberDown_ = NULL;
774          numberUp_ = NULL;
775          numberDownInfeasible_ = NULL;
776          numberUpInfeasible_ = NULL;
777          saveCosts_ = NULL;
778          nodeInfo_ = NULL;
779          large_ = NULL;
780          whichRow_ = NULL;
781          whichColumn_ = NULL;
782          nBound_ = 0;
783          saveOptions_ = rhs.saveOptions_;
784          solverOptions_ = rhs.solverOptions_;
785          maximumNodes_ = rhs.maximumNodes_;
786          numberBeforeTrust_ = rhs.numberBeforeTrust_;
787          stateOfSearch_ = rhs.stateOfSearch_;
788          int n = maximumNodes();
789          if (n) {
790               for (int i = 0; i < n; i++)
791                    delete nodeInfo_[i];
792          }
793          delete [] nodeInfo_;
794          nodeInfo_ = NULL;
795          nDepth_ = rhs.nDepth_;
796          nNodes_ = rhs.nNodes_;
797          numberNodesExplored_ = rhs.numberNodesExplored_;
798          numberIterations_ = rhs.numberIterations_;
799          presolveType_ = rhs.presolveType_;
800     }
801     return *this;
802}
803// Zaps stuff 1 - arrays, 2 ints, 3 both
804void
805ClpNodeStuff::zap(int type)
806{
807     if ((type & 1) != 0) {
808          downPseudo_ = NULL;
809          upPseudo_ = NULL;
810          priority_ = NULL;
811          numberDown_ = NULL;
812          numberUp_ = NULL;
813          numberDownInfeasible_ = NULL;
814          numberUpInfeasible_ = NULL;
815          saveCosts_ = NULL;
816          nodeInfo_ = NULL;
817          large_ = NULL;
818          whichRow_ = NULL;
819          whichColumn_ = NULL;
820     }
821     if ((type & 2) != 0) {
822          nBound_ = 0;
823          saveOptions_ = 0;
824          solverOptions_ = 0;
825          maximumNodes_ = 0;
826          numberBeforeTrust_ = 0;
827          stateOfSearch_ = 0;
828          nDepth_ = -1;
829          nNodes_ = 0;
830          presolveType_ = 0;
831          numberNodesExplored_ = 0;
832          numberIterations_ = 0;
833     }
834}
835
836//-------------------------------------------------------------------
837// Destructor
838//-------------------------------------------------------------------
839ClpNodeStuff::~ClpNodeStuff ()
840{
841     delete [] downPseudo_;
842     delete [] upPseudo_;
843     delete [] priority_;
844     delete [] numberDown_;
845     delete [] numberUp_;
846     delete [] numberDownInfeasible_;
847     delete [] numberUpInfeasible_;
848     int n = maximumNodes();
849     if (n) {
850          for (int i = 0; i < n; i++)
851               delete nodeInfo_[i];
852     }
853     delete [] nodeInfo_;
854#ifdef CLP_INVESTIGATE
855     // Should be NULL - find out why not?
856     assert (!saveCosts_);
857#endif
858     delete [] saveCosts_;
859}
860// Return maximum number of nodes
861int
862ClpNodeStuff::maximumNodes() const
863{
864     int n = 0;
865#if 0
866     if (nDepth_ != -1) {
867          if ((solverOptions_ & 32) == 0)
868               n = (1 << nDepth_);
869          else if (nDepth_)
870               n = 1;
871     }
872     assert (n == maximumNodes_ - (1 + nDepth_) || n == 0);
873#else
874     if (nDepth_ != -1) {
875          n = maximumNodes_ - (1 + nDepth_);
876          assert (n > 0);
877     }
878#endif
879     return n;
880}
881// Return maximum space for nodes
882int
883ClpNodeStuff::maximumSpace() const
884{
885     return maximumNodes_;
886}
887/* Fill with pseudocosts */
888void
889ClpNodeStuff::fillPseudoCosts(const double * down, const double * up,
890                              const int * priority,
891                              const int * numberDown, const int * numberUp,
892                              const int * numberDownInfeasible,
893                              const int * numberUpInfeasible,
894                              int number)
895{
896     delete [] downPseudo_;
897     delete [] upPseudo_;
898     delete [] priority_;
899     delete [] numberDown_;
900     delete [] numberUp_;
901     delete [] numberDownInfeasible_;
902     delete [] numberUpInfeasible_;
903     downPseudo_ = CoinCopyOfArray(down, number);
904     upPseudo_ = CoinCopyOfArray(up, number);
905     priority_ = CoinCopyOfArray(priority, number);
906     numberDown_ = CoinCopyOfArray(numberDown, number);
907     numberUp_ = CoinCopyOfArray(numberUp, number);
908     numberDownInfeasible_ = CoinCopyOfArray(numberDownInfeasible, number);
909     numberUpInfeasible_ = CoinCopyOfArray(numberUpInfeasible, number);
910     // scale
911     for (int i = 0; i < number; i++) {
912          int n;
913          n = numberDown_[i];
914          if (n)
915               downPseudo_[i] *= n;
916          n = numberUp_[i];
917          if (n)
918               upPseudo_[i] *= n;
919     }
920}
921// Update pseudo costs
922void
923ClpNodeStuff::update(int way, int sequence, double change, bool feasible)
924{
925     assert (numberDown_[sequence] >= numberDownInfeasible_[sequence]);
926     assert (numberUp_[sequence] >= numberUpInfeasible_[sequence]);
927     if (way < 0) {
928          numberDown_[sequence]++;
929          if (!feasible)
930               numberDownInfeasible_[sequence]++;
931          downPseudo_[sequence] += CoinMax(change, 1.0e-12);
932     } else {
933          numberUp_[sequence]++;
934          if (!feasible)
935               numberUpInfeasible_[sequence]++;
936          upPseudo_[sequence] += CoinMax(change, 1.0e-12);
937     }
938}
939//#############################################################################
940// Constructors / Destructor / Assignment
941//#############################################################################
942
943//-------------------------------------------------------------------
944// Default Constructor
945//-------------------------------------------------------------------
946ClpHashValue::ClpHashValue () :
947     hash_(NULL),
948     numberHash_(0),
949     maxHash_(0),
950     lastUsed_(-1)
951{
952}
953//-------------------------------------------------------------------
954// Useful Constructor from model
955//-------------------------------------------------------------------
956ClpHashValue::ClpHashValue (ClpSimplex * model) :
957     hash_(NULL),
958     numberHash_(0),
959     maxHash_(0),
960     lastUsed_(-1)
961{
962     maxHash_ = 1000;
963     int numberColumns = model->numberColumns();
964     const double * columnLower = model->columnLower();
965     const double * columnUpper = model->columnUpper();
966     int numberRows = model->numberRows();
967     const double * rowLower = model->rowLower();
968     const double * rowUpper = model->rowUpper();
969     const double * objective = model->objective();
970     CoinPackedMatrix * matrix = model->matrix();
971     const int * columnLength = matrix->getVectorLengths();
972     const CoinBigIndex * columnStart = matrix->getVectorStarts();
973     const double * elementByColumn = matrix->getElements();
974     int i;
975     int ipos;
976
977     hash_ = new CoinHashLink[maxHash_];
978
979     for ( i = 0; i < maxHash_; i++ ) {
980          hash_[i].value = -1.0e-100;
981          hash_[i].index = -1;
982          hash_[i].next = -1;
983     }
984     // Put in +0
985     hash_[0].value = 0.0;
986     hash_[0].index = 0;
987     numberHash_ = 1;
988     /*
989      * Initialize the hash table.  Only the index of the first value that
990      * hashes to a value is entered in the table; subsequent values that
991      * collide with it are not entered.
992      */
993     for ( i = 0; i < numberColumns; i++ ) {
994          int length = columnLength[i];
995          CoinBigIndex start = columnStart[i];
996          for (CoinBigIndex i = start; i < start + length; i++) {
997               double value = elementByColumn[i];
998               ipos = hash ( value);
999               if ( hash_[ipos].index == -1 ) {
1000                    hash_[ipos].index = numberHash_;
1001                    numberHash_++;
1002                    hash_[ipos].value = elementByColumn[i];
1003               }
1004          }
1005     }
1006
1007     /*
1008      * Now take care of the values that collided in the preceding loop,
1009      * Also do other stuff
1010      */
1011     for ( i = 0; i < numberRows; i++ ) {
1012          if (numberHash_ * 2 > maxHash_)
1013               resize(true);
1014          double value;
1015          value = rowLower[i];
1016          ipos = index(value);
1017          if (ipos < 0)
1018               addValue(value);
1019          value = rowUpper[i];
1020          ipos = index(value);
1021          if (ipos < 0)
1022               addValue(value);
1023     }
1024     for ( i = 0; i < numberColumns; i++ ) {
1025          int length = columnLength[i];
1026          CoinBigIndex start = columnStart[i];
1027          if (numberHash_ * 2 > maxHash_)
1028               resize(true);
1029          double value;
1030          value = objective[i];
1031          ipos = index(value);
1032          if (ipos < 0)
1033               addValue(value);
1034          value = columnLower[i];
1035          ipos = index(value);
1036          if (ipos < 0)
1037               addValue(value);
1038          value = columnUpper[i];
1039          ipos = index(value);
1040          if (ipos < 0)
1041               addValue(value);
1042          for (CoinBigIndex j = start; j < start + length; j++) {
1043               if (numberHash_ * 2 > maxHash_)
1044                    resize(true);
1045               value = elementByColumn[j];
1046               ipos = index(value);
1047               if (ipos < 0)
1048                    addValue(value);
1049          }
1050     }
1051     resize(false);
1052}
1053
1054//-------------------------------------------------------------------
1055// Copy constructor
1056//-------------------------------------------------------------------
1057ClpHashValue::ClpHashValue (const ClpHashValue & rhs) :
1058     hash_(NULL),
1059     numberHash_(rhs.numberHash_),
1060     maxHash_(rhs.maxHash_),
1061     lastUsed_(rhs.lastUsed_)
1062{
1063     if (maxHash_) {
1064          CoinHashLink * newHash = new CoinHashLink[maxHash_];
1065          int i;
1066          for ( i = 0; i < maxHash_; i++ ) {
1067               newHash[i].value = rhs.hash_[i].value;
1068               newHash[i].index = rhs.hash_[i].index;
1069               newHash[i].next = rhs.hash_[i].next;
1070          }
1071     }
1072}
1073
1074//-------------------------------------------------------------------
1075// Destructor
1076//-------------------------------------------------------------------
1077ClpHashValue::~ClpHashValue ()
1078{
1079     delete [] hash_;
1080}
1081
1082//----------------------------------------------------------------
1083// Assignment operator
1084//-------------------------------------------------------------------
1085ClpHashValue &
1086ClpHashValue::operator=(const ClpHashValue& rhs)
1087{
1088     if (this != &rhs) {
1089          numberHash_ = rhs.numberHash_;
1090          maxHash_ = rhs.maxHash_;
1091          lastUsed_ = rhs.lastUsed_;
1092          delete [] hash_;
1093          if (maxHash_) {
1094               CoinHashLink * newHash = new CoinHashLink[maxHash_];
1095               int i;
1096               for ( i = 0; i < maxHash_; i++ ) {
1097                    newHash[i].value = rhs.hash_[i].value;
1098                    newHash[i].index = rhs.hash_[i].index;
1099                    newHash[i].next = rhs.hash_[i].next;
1100               }
1101          } else {
1102               hash_ = NULL;
1103          }
1104     }
1105     return *this;
1106}
1107// Return index or -1 if not found
1108int
1109ClpHashValue::index(double value) const
1110{
1111     if (!value)
1112          return 0;
1113     int ipos = hash ( value);
1114     int returnCode = -1;
1115     while ( hash_[ipos].index >= 0 ) {
1116          if (value == hash_[ipos].value) {
1117               returnCode = hash_[ipos].index;
1118               break;
1119          } else {
1120               int k = hash_[ipos].next;
1121               if ( k == -1 ) {
1122                    break;
1123               } else {
1124                    ipos = k;
1125               }
1126          }
1127     }
1128     return returnCode;
1129}
1130// Add value to list and return index
1131int
1132ClpHashValue::addValue(double value)
1133{
1134     int ipos = hash ( value);
1135
1136     assert (value != hash_[ipos].value);
1137     if (hash_[ipos].index == -1) {
1138          // can put in here
1139          hash_[ipos].index = numberHash_;
1140          numberHash_++;
1141          hash_[ipos].value = value;
1142          return numberHash_ - 1;
1143     }
1144     int k = hash_[ipos].next;
1145     while (k != -1) {
1146          ipos = k;
1147          k = hash_[ipos].next;
1148     }
1149     while ( true ) {
1150          ++lastUsed_;
1151          assert (lastUsed_ <= maxHash_);
1152          if ( hash_[lastUsed_].index == -1 ) {
1153               break;
1154          }
1155     }
1156     hash_[ipos].next = lastUsed_;
1157     hash_[lastUsed_].index = numberHash_;
1158     numberHash_++;
1159     hash_[lastUsed_].value = value;
1160     return numberHash_ - 1;
1161}
1162
1163namespace {
1164  /*
1165    Originally a local static variable in ClpHashValue::hash.
1166        Local static variables are a problem when building DLLs on Windows, but
1167        file-local constants seem to be ok.  -- lh, 101016 --
1168  */
1169  const int mmult_for_hash[] = {
1170          262139, 259459, 256889, 254291, 251701, 249133, 246709, 244247,
1171          241667, 239179, 236609, 233983, 231289, 228859, 226357, 223829,
1172          221281, 218849, 216319, 213721, 211093, 208673, 206263, 203773,
1173          201233, 198637, 196159, 193603, 191161, 188701, 186149, 183761,
1174          181303, 178873, 176389, 173897, 171469, 169049, 166471, 163871,
1175          161387, 158941, 156437, 153949, 151531, 149159, 146749, 144299,
1176          141709, 139369, 136889, 134591, 132169, 129641, 127343, 124853,
1177          122477, 120163, 117757, 115361, 112979, 110567, 108179, 105727,
1178          103387, 101021, 98639, 96179, 93911, 91583, 89317, 86939, 84521,
1179          82183, 79939, 77587, 75307, 72959, 70793, 68447, 66103
1180     };
1181}
1182int
1183ClpHashValue::hash ( double value) const
1184{
1185     
1186     union {
1187          double d;
1188          char c[8];
1189     } v1;
1190     assert (sizeof(double) == 8);
1191     v1.d = value;
1192     int n = 0;
1193     int j;
1194
1195     for ( j = 0; j < 8; ++j ) {
1196          int ichar = v1.c[j];
1197          n += mmult_for_hash[j] * ichar;
1198     }
1199     return ( abs ( n ) % maxHash_ );   /* integer abs */
1200}
1201void
1202ClpHashValue::resize(bool increaseMax)
1203{
1204     int newSize = increaseMax ? ((3 * maxHash_) >> 1) + 1000 : maxHash_;
1205     CoinHashLink * newHash = new CoinHashLink[newSize];
1206     int i;
1207     for ( i = 0; i < newSize; i++ ) {
1208          newHash[i].value = -1.0e-100;
1209          newHash[i].index = -1;
1210          newHash[i].next = -1;
1211     }
1212     // swap
1213     CoinHashLink * oldHash = hash_;
1214     hash_ = newHash;
1215     int oldSize = maxHash_;
1216     maxHash_ = newSize;
1217     /*
1218      * Initialize the hash table.  Only the index of the first value that
1219      * hashes to a value is entered in the table; subsequent values that
1220      * collide with it are not entered.
1221      */
1222     int ipos;
1223     int n = 0;
1224     for ( i = 0; i < oldSize; i++ ) {
1225          if (oldHash[i].index >= 0) {
1226               ipos = hash ( oldHash[i].value);
1227               if ( hash_[ipos].index == -1 ) {
1228                    hash_[ipos].index = n;
1229                    n++;
1230                    hash_[ipos].value = oldHash[i].value;
1231                    // unmark
1232                    oldHash[i].index = -1;
1233               }
1234          }
1235     }
1236     /*
1237      * Now take care of the values that collided in the preceding loop,
1238      * by finding some other entry in the table for them.
1239      * Since there are as many entries in the table as there are values,
1240      * there must be room for them.
1241      */
1242     lastUsed_ = -1;
1243     for ( i = 0; i < oldSize; ++i ) {
1244          if (oldHash[i].index >= 0) {
1245               double value = oldHash[i].value;
1246               ipos = hash ( value);
1247               int k;
1248               while ( true ) {
1249                    assert (value != hash_[ipos].value);
1250                    k = hash_[ipos].next;
1251                    if ( k == -1 ) {
1252                         while ( true ) {
1253                              ++lastUsed_;
1254                              assert (lastUsed_ <= maxHash_);
1255                              if ( hash_[lastUsed_].index == -1 ) {
1256                                   break;
1257                              }
1258                         }
1259                         hash_[ipos].next = lastUsed_;
1260                         hash_[lastUsed_].index = n;
1261                         n++;
1262                         hash_[lastUsed_].value = value;
1263                         break;
1264                    } else {
1265                         ipos = k;
1266                    }
1267               }
1268          }
1269     }
1270     assert (n == numberHash_);
1271     delete [] oldHash;
1272}
Note: See TracBrowser for help on using the repository browser.