source: trunk/Clp/src/ClpNode.cpp @ 2030

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