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

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

out compiler warnings and stability improvements

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