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

Last change on this file since 2470 was 2385, checked in by unxusr, 9 months ago

formatting

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