source: stable/1.7/Clp/src/ClpNode.cpp @ 1213

Last change on this file since 1213 was 1213, checked in by forrest, 13 years ago

for Lou Hafer to update Osi stable

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