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

Last change on this file since 1370 was 1370, checked in by forrest, 11 years ago

add ids

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