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

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

many changes to try and improve performance

  • Property svn:executable set to *
File size: 29.6 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 = 0.0;
763    hash_[i].index = -1;
764    hash_[i].next = -1;
765  }
766
767  /*
768   * Initialize the hash table.  Only the index of the first value that
769   * hashes to a value is entered in the table; subsequent values that
770   * collide with it are not entered.
771   */
772  for ( i = 0; i < numberColumns; i++ ) {
773    int length=columnLength[i];
774    CoinBigIndex start = columnStart[i];
775    for (CoinBigIndex i=start;i<start+length;i++) {
776      ipos = hash ( elementByColumn[i]);
777      if ( hash_[ipos].index == -1 ) {
778        hash_[ipos].index = numberHash_;
779        numberHash_++;
780        hash_[ipos].value = elementByColumn[i];
781      }
782    }
783  }
784 
785  /*
786   * Now take care of the values that collided in the preceding loop,
787   * Also do other stuff
788   */
789  for ( i = 0; i < numberRows; i++ ) {
790    if (numberHash_*2>maxHash_)
791      resize(true);
792    double value;
793    value = rowLower[i];
794    ipos = index(value);
795    if (ipos<0) 
796      addValue(value);
797    value = rowUpper[i];
798    ipos = index(value);
799    if (ipos<0) 
800      addValue(value);
801  }
802  for ( i = 0; i < numberColumns; i++ ) {
803    int length=columnLength[i];
804    CoinBigIndex start = columnStart[i];
805    if (numberHash_*2>maxHash_)
806      resize(true);
807    double value;
808    value = objective[i];
809    ipos = index(value);
810    if (ipos<0) 
811      addValue(value);
812    value = columnLower[i];
813    ipos = index(value);
814    if (ipos<0) 
815      addValue(value);
816    value = columnUpper[i];
817    ipos = index(value);
818    if (ipos<0) 
819      addValue(value);
820    for (CoinBigIndex j=start;j<start+length;j++) {
821      if (numberHash_*2>maxHash_)
822        resize(true);
823      value = elementByColumn[j];
824      ipos = index(value);
825      if (ipos<0) 
826        addValue(value);
827    }
828  }
829  resize(false);
830}
831
832//-------------------------------------------------------------------
833// Copy constructor
834//-------------------------------------------------------------------
835ClpHashValue::ClpHashValue (const ClpHashValue & rhs) :
836  hash_(NULL),
837  numberHash_(rhs.numberHash_),
838  maxHash_(rhs.maxHash_),
839  lastUsed_(rhs.lastUsed_)
840{ 
841  if (maxHash_) {
842    CoinHashLink * newHash = new CoinHashLink[maxHash_];
843    int i;
844    for ( i = 0; i < maxHash_; i++ ) {
845      newHash[i].value = rhs.hash_[i].value;
846      newHash[i].index = rhs.hash_[i].index;
847      newHash[i].next = rhs.hash_[i].next;
848    }
849  }
850}
851
852//-------------------------------------------------------------------
853// Destructor
854//-------------------------------------------------------------------
855ClpHashValue::~ClpHashValue ()
856{
857  delete [] hash_;
858}
859
860//----------------------------------------------------------------
861// Assignment operator
862//-------------------------------------------------------------------
863ClpHashValue &
864ClpHashValue::operator=(const ClpHashValue& rhs)
865{
866  if (this != &rhs) {
867    numberHash_ = rhs.numberHash_;
868    maxHash_ = rhs.maxHash_;
869    lastUsed_ = rhs.lastUsed_;
870    delete [] hash_;
871    if (maxHash_) {
872      CoinHashLink * newHash = new CoinHashLink[maxHash_];
873      int i;
874      for ( i = 0; i < maxHash_; i++ ) {
875        newHash[i].value = rhs.hash_[i].value;
876        newHash[i].index = rhs.hash_[i].index;
877        newHash[i].next = rhs.hash_[i].next;
878      }
879    } else {
880      hash_ = NULL;
881    }
882  }
883  return *this;
884}
885// Return index or -1 if not found
886int 
887ClpHashValue::index(double value) const
888{
889  int ipos = hash ( value);
890  int returnCode=-1;
891  while ( hash_[ipos].index>=0 ) {
892    if (value==hash_[ipos].value) {
893      returnCode = hash_[ipos].index;
894      break;
895    } else { 
896      int k = hash_[ipos].next;
897      if ( k == -1 ) {
898        break;
899      } else {
900        ipos = k;
901      }
902    }
903  }
904  return returnCode;
905}
906// Add value to list and return index
907int 
908ClpHashValue::addValue(double value) 
909{
910  int ipos = hash ( value);
911 
912  assert (value!=hash_[ipos].value);
913  int k = hash_[ipos].next;
914  while (k!=-1) {
915    ipos = k;
916    k = hash_[ipos].next;
917  }
918  while ( true ) {
919    ++lastUsed_;
920    assert (lastUsed_<=maxHash_);
921    if ( hash_[lastUsed_].index == -1 ) {
922      break;
923    }
924  }
925  hash_[ipos].next = lastUsed_;
926  hash_[lastUsed_].index = numberHash_;
927  numberHash_++;
928  hash_[lastUsed_].value = value;
929  return numberHash_-1;
930}
931
932int
933ClpHashValue::hash ( double value) const
934{
935  static int mmult[] = {
936    262139, 259459, 256889, 254291, 251701, 249133, 246709, 244247,
937    241667, 239179, 236609, 233983, 231289, 228859, 226357, 223829,
938    221281, 218849, 216319, 213721, 211093, 208673, 206263, 203773,
939    201233, 198637, 196159, 193603, 191161, 188701, 186149, 183761,
940    181303, 178873, 176389, 173897, 171469, 169049, 166471, 163871,
941    161387, 158941, 156437, 153949, 151531, 149159, 146749, 144299,
942    141709, 139369, 136889, 134591, 132169, 129641, 127343, 124853,
943    122477, 120163, 117757, 115361, 112979, 110567, 108179, 105727,
944    103387, 101021, 98639, 96179, 93911, 91583, 89317, 86939, 84521,
945    82183, 79939, 77587, 75307, 72959, 70793, 68447, 66103
946  };
947  union { double d; char c[8]; } v1;
948  assert (sizeof(double)==8);
949  v1.d = value;
950  int n = 0;
951  int j;
952
953  for ( j = 0; j < 8; ++j ) {
954    int ichar = v1.c[j];
955    n += mmult[j] * ichar;
956  }
957  return ( abs ( n ) % maxHash_ );      /* integer abs */
958}
959void
960ClpHashValue::resize(bool increaseMax)
961{
962  int newSize = increaseMax ? ((3*maxHash_)>>1)+1000 : maxHash_;
963  CoinHashLink * newHash = new CoinHashLink[newSize];
964  int i;
965  for ( i = 0; i < newSize; i++ ) {
966    newHash[i].value = 0.0;
967    newHash[i].index = -1;
968    newHash[i].next = -1;
969  }
970  // swap
971  CoinHashLink * oldHash = hash_;
972  hash_ = newHash;
973  int oldSize = maxHash_;
974  maxHash_ = newSize;
975  /*
976   * Initialize the hash table.  Only the index of the first value that
977   * hashes to a value is entered in the table; subsequent values that
978   * collide with it are not entered.
979   */
980  int ipos;
981  int n=0;
982  for ( i = 0; i < oldSize; i++ ) {
983    if (oldHash[i].index>=0) {
984      ipos = hash ( oldHash[i].value);
985      if ( hash_[ipos].index == -1 ) {
986        hash_[ipos].index = n;
987        n++;
988        hash_[ipos].value = hash_[i].value;
989      }
990    }
991  }
992  /*
993   * Now take care of the values that collided in the preceding loop,
994   * by finding some other entry in the table for them.
995   * Since there are as many entries in the table as there are values,
996   * there must be room for them.
997   */
998  lastUsed_ = -1;
999  for ( i = 0; i < oldSize; ++i ) {
1000    if (oldHash[i].index>=0) {
1001      double value = oldHash[i].value;
1002      ipos = hash ( value);
1003      int k;
1004      while ( true ) {
1005        if (value==hash_[ipos].value)
1006          break;
1007        else 
1008          k = hash_[ipos].next;
1009       
1010        if ( k == -1 ) {
1011          while ( true ) {
1012            ++lastUsed_;
1013            assert (lastUsed_<=maxHash_);
1014            if ( hash_[lastUsed_].index == -1 ) {
1015              break;
1016            }
1017          }
1018          hash_[ipos].next = lastUsed_;
1019          hash_[lastUsed_].index = n;
1020          n++;
1021          hash_[lastUsed_].value = value;
1022          break;
1023        } else {
1024          ipos = k;
1025        }
1026      }
1027    }
1028  }
1029  assert (n==numberHash_);
1030}
Note: See TracBrowser for help on using the repository browser.