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

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

wrong #define

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