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

Last change on this file since 1202 was 1202, checked in by forrest, 12 years ago

for some of my experiments - not normally turned on

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