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

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

changes to try and make faster

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