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

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

changes for factorization and aux region

  • Property svn:executable set to *
File size: 32.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  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());
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=integerTolerance;
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,(double) lower[iColumn]);
254      value = CoinMin(value,(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        double upValue = (ceil(value)-value)*(upPseudo[iInteger]/(1.0+numberUp[iInteger]+numberUpInfeasible[iInteger]));
298        double downValue = (value-floor(value))*(downPseudo[iInteger]/(1.0+numberDown[iInteger]+numberDownInfeasible[iInteger]));
299        double infeasibility;
300        if (depth>1000)
301          infeasibility = CoinMax(upValue,downValue)+integerTolerance;
302        else
303          infeasibility = 0.1*CoinMax(upValue,downValue)+
304            0.9*CoinMin(upValue,downValue) + integerTolerance;
305        estimatedSolution_ += CoinMin(upValue,downValue);
306#else
307        double infeasibility = fabs(value-nearest);
308#endif
309        if (infeasibility>mostAway) {
310          mostAway=infeasibility;
311          sequence_=iColumn;
312          branchingValue_=value;
313          branchState_.branch=0;
314#if PSEUDO>0
315          if (upValue<=downValue)
316            branchState_.firstBranch=1; // up
317          else
318            branchState_.firstBranch=0; // down
319#else
320          if (value<=nearest)
321            branchState_.firstBranch=1; // up
322          else
323            branchState_.firstBranch=0; // down
324#endif
325        }
326      } else if (model->getColumnStatus(iColumn)==ClpSimplex::atLowerBound) {
327        bool fix=false;
328        if (columnScale) {
329          if (dualSolution_[iColumn]>gap*columnScale[iColumn]) 
330            fix=true;
331        } else {
332          if (dualSolution_[iColumn]>gap) 
333            fix=true;
334        }
335        if (fix) {
336          nFix++;
337          //printf("fixed %d to zero gap %g dj %g %g\n",iColumn,
338          // gap,dualSolution_[iColumn], columnScale ? columnScale[iColumn]:1.0);
339          model->setColumnStatus(iColumn,ClpSimplex::isFixed);
340        }
341      } else if (model->getColumnStatus(iColumn)==ClpSimplex::atUpperBound) {
342        bool fix=false;
343        if (columnScale) {
344          if (-dualSolution_[iColumn]>gap*columnScale[iColumn]) 
345            fix=true;
346        } else {
347          if (-dualSolution_[iColumn]>gap) 
348            fix=true;
349        }
350        if (fix) {
351          nFix++;
352          //printf("fixed %d to one gap %g dj %g %g\n",iColumn,
353          // gap,dualSolution_[iColumn], columnScale ? columnScale[iColumn]:1.0);
354          model->setColumnStatus(iColumn,ClpSimplex::isFixed);
355        }
356      }
357      iInteger++;
358    }
359  }
360#if PSEUDO == 2
361  delete [] activeWeight;
362#endif
363  if (lower_) {
364    // save bounds
365    if (iInteger>maximumIntegers_) {
366      delete [] lower_;
367      delete [] upper_;
368      maximumIntegers_ = iInteger;
369      lower_ = new int [maximumIntegers_];
370      upper_ = new int [maximumIntegers_];
371    }
372    iInteger=0;
373    for (iColumn=0;iColumn<numberColumns;iColumn++) {
374      if (integerType[iColumn]) {
375        lower_[iInteger]=(int) lower[iColumn];
376        upper_[iInteger]=(int) upper[iColumn];
377        iInteger++;
378      }
379    }
380  }
381  // Could omit save of fixed if doing full save of bounds
382  if (sequence_>=0&&nFix) {
383    if (nFix>maximumFixed_) {
384      delete [] fixed_;
385      fixed_ = new int [nFix];
386      maximumFixed_=nFix;
387    }
388    numberFixed_=0;
389    unsigned char * status = model->statusArray();
390    for (iColumn=0;iColumn<numberColumns;iColumn++) {
391      if (status[iColumn]!=status_[iColumn]) {
392        if (solution[iColumn]<=lower[iColumn]+2.0*integerTolerance) {
393          model->setColumnUpper(iColumn,lower[iColumn]);
394          fixed_[numberFixed_++]=iColumn;
395        } else {
396          assert (solution[iColumn]>=upper[iColumn]-2.0*integerTolerance);
397          model->setColumnLower(iColumn,upper[iColumn]);
398          fixed_[numberFixed_++]=iColumn|0x10000000;
399        }
400      }
401    }
402    //printf("%d fixed\n",numberFixed_);
403  }
404}
405
406//-------------------------------------------------------------------
407// Copy constructor
408//-------------------------------------------------------------------
409ClpNode::ClpNode (const ClpNode & source) 
410{ 
411  printf("ClpNode copy not implemented\n");
412  abort();
413}
414
415//-------------------------------------------------------------------
416// Destructor
417//-------------------------------------------------------------------
418ClpNode::~ClpNode ()
419{
420  delete factorization_;
421  delete weights_;
422  delete [] status_;
423  delete [] primalSolution_;
424  delete [] dualSolution_;
425  delete [] lower_;
426  delete [] upper_;
427  delete [] pivotVariables_;
428  delete [] fixed_;
429}
430
431//----------------------------------------------------------------
432// Assignment operator
433//-------------------------------------------------------------------
434ClpNode &
435ClpNode::operator=(const ClpNode& rhs)
436{
437  if (this != &rhs) {
438    printf("ClpNode = not implemented\n");
439    abort();
440  }
441  return *this;
442}
443// Create odd arrays
444void 
445ClpNode::createArrays(ClpSimplex * model)
446{
447  int numberColumns = model->numberColumns();
448  const char * integerType = model->integerInformation();
449  int iColumn;
450  int numberIntegers=0;
451  for (iColumn=0;iColumn<numberColumns;iColumn++) {
452    if (integerType[iColumn]) 
453      numberIntegers++;
454  }
455  if (numberIntegers>maximumIntegers_||!lower_) {
456    delete [] lower_;
457    delete [] upper_;
458    maximumIntegers_ = numberIntegers;
459    lower_ = new int [numberIntegers];
460    upper_ = new int [numberIntegers];
461  }
462}
463// Clean up as crunch is different model
464void 
465ClpNode::cleanUpForCrunch()
466{
467  delete weights_;
468  weights_ = NULL;
469}
470/* Applies node to model
471   0 - just tree bounds
472   1 - tree bounds and basis etc
473   2 - saved bounds and basis etc
474*/
475void 
476ClpNode::applyNode(ClpSimplex * model, int doBoundsEtc )
477{
478  int numberColumns = model->numberColumns();
479  const double * lower = model->columnLower();
480  const double * upper = model->columnUpper();
481  if (doBoundsEtc<2) {
482    // current bound
483    int way=branchState_.firstBranch;
484    if (branchState_.branch>0)
485      way=1-way;
486    if (!way) {
487      // This should also do underlying internal bound
488      model->setColumnUpper(sequence_,floor(branchingValue_));
489    } else {
490      // This should also do underlying internal bound
491      model->setColumnLower(sequence_,ceil(branchingValue_));
492    }
493    // apply dj fixings
494    for (int i=0;i<numberFixed_;i++) {
495      int iColumn = fixed_[i];
496      if ((iColumn&0x10000000)!=0) {
497        iColumn &= 0xfffffff;
498        model->setColumnLower(iColumn,upper[iColumn]);
499    } else {
500        model->setColumnUpper(iColumn,lower[iColumn]);
501      }
502    }
503  } else {
504    // restore bounds
505    assert (lower_);
506    int iInteger=-1;
507    const char * integerType = model->integerInformation();
508    for (int iColumn=0;iColumn<numberColumns;iColumn++) {
509      if (integerType[iColumn]) {
510        iInteger++;
511        if (lower_[iInteger]!=(int) lower[iColumn])
512          model->setColumnLower(iColumn,lower_[iInteger]);
513        if (upper_[iInteger]!=(int) upper[iColumn])
514          model->setColumnUpper(iColumn,upper_[iInteger]);
515      }
516    }
517  }
518  if (doBoundsEtc&&doBoundsEtc<3) {
519    //model->copyFactorization(*factorization_);
520    model->copyFactorization(*factorization_);
521    ClpDualRowSteepest* pivot =
522      dynamic_cast< ClpDualRowSteepest*>(model->dualRowPivot());
523    if (pivot&&weights_) {
524      pivot->fill(*weights_); 
525    }
526    int numberRows = model->numberRows();
527    int numberTotal = numberRows+numberColumns;
528    CoinMemcpyN(status_,numberTotal,model->statusArray());
529    if (doBoundsEtc<2) {
530      CoinMemcpyN(primalSolution_,numberTotal,model->solutionRegion());
531      CoinMemcpyN(dualSolution_,numberTotal,model->djRegion());
532      CoinMemcpyN(pivotVariables_,numberRows,model->pivotVariable());
533      CoinMemcpyN(dualSolution_+numberColumns,numberRows,model->dualRowSolution());
534    } else {
535      CoinMemcpyN(primalSolution_,numberColumns,model->primalColumnSolution());
536      CoinMemcpyN(dualSolution_,numberColumns,model->dualColumnSolution());
537      CoinMemcpyN(dualSolution_+numberColumns,numberRows,model->dualRowSolution());
538      if (model->columnScale()) {
539        // See if just primal will work
540        double * solution = model->primalColumnSolution();
541        const double * columnScale = model->columnScale();
542        int i;
543        for (i=0;i<numberColumns;i++) {
544          solution[i] *= columnScale[i];
545        }
546      }
547    }
548    model->setObjectiveValue(objectiveValue_);
549  }
550}
551// Choose a new variable
552void 
553ClpNode::chooseVariable(ClpSimplex * model, ClpNodeStuff * info)
554{
555#if 0
556  int way=branchState_.firstBranch;
557  if (branchState_.branch>0)
558    way=1-way;
559  assert (!branchState_.branch);
560  // We need to use pseudo costs to choose a variable
561  int numberColumns = model->numberColumns();
562#endif
563}
564// Fix on reduced costs
565int 
566ClpNode::fixOnReducedCosts(ClpSimplex * model)
567{
568 
569  return 0;
570}
571/* Way for integer variable -1 down , +1 up */
572int 
573ClpNode::way() const
574{
575  int way=branchState_.firstBranch;
576  if (branchState_.branch>0)
577    way=1-way;
578  return way==0 ? -1 : +1;
579}
580// Return true if branch exhausted
581bool 
582ClpNode::fathomed() const
583{ 
584  return branchState_.branch>=1
585;
586}
587// Change state of variable i.e. go other way
588void 
589ClpNode::changeState()
590{
591  branchState_.branch++;
592  assert (branchState_.branch<=2);
593}
594//#############################################################################
595// Constructors / Destructor / Assignment
596//#############################################################################
597
598//-------------------------------------------------------------------
599// Default Constructor
600//-------------------------------------------------------------------
601ClpNodeStuff::ClpNodeStuff () :
602  integerTolerance_(1.0e-7),
603  integerIncrement_(1.0e-8),
604  downPseudo_(NULL),
605  upPseudo_(NULL),
606  numberDown_(NULL),
607  numberUp_(NULL),
608  numberDownInfeasible_(NULL),
609  numberUpInfeasible_(NULL),
610  saveCosts_(NULL),
611  nodeInfo_(NULL),
612  large_(NULL),
613  whichRow_(NULL),
614  whichColumn_(NULL),
615  nBound_(0),
616  saveOptions_(0),
617  solverOptions_(0),
618  nDepth_(-1),
619  nNodes_(0),
620  numberNodesExplored_(0),
621  numberIterations_(0),
622  presolveType_(0)
623{
624
625}
626
627//-------------------------------------------------------------------
628// Copy constructor
629//-------------------------------------------------------------------
630ClpNodeStuff::ClpNodeStuff (const ClpNodeStuff & rhs) 
631  : integerTolerance_(rhs.integerTolerance_),
632    integerIncrement_(rhs.integerIncrement_),
633    downPseudo_(NULL),
634    upPseudo_(NULL),
635    numberDown_(NULL),
636    numberUp_(NULL),
637    numberDownInfeasible_(NULL),
638    numberUpInfeasible_(NULL),
639    saveCosts_(NULL),
640    nodeInfo_(NULL),
641    large_(NULL),
642    whichRow_(NULL),
643    whichColumn_(NULL),
644    nBound_(0),
645    saveOptions_(rhs.saveOptions_),
646    solverOptions_(rhs.solverOptions_),
647    nDepth_(rhs.nDepth_),
648    nNodes_(rhs.nNodes_),
649    numberNodesExplored_(rhs.numberNodesExplored_),
650    numberIterations_(rhs.numberIterations_),
651    presolveType_(rhs.presolveType_)
652{ 
653}
654//----------------------------------------------------------------
655// Assignment operator
656//-------------------------------------------------------------------
657ClpNodeStuff &
658ClpNodeStuff::operator=(const ClpNodeStuff& rhs)
659{
660  if (this != &rhs) {
661    integerTolerance_ = rhs.integerTolerance_;
662    integerIncrement_ = rhs.integerIncrement_;
663    downPseudo_ = NULL;
664    upPseudo_ = NULL;
665    numberDown_ = NULL;
666    numberUp_ = NULL;
667    numberDownInfeasible_ = NULL;
668    numberUpInfeasible_ = NULL;
669    saveCosts_ = NULL;
670    nodeInfo_ = NULL;
671    large_ = NULL;
672    whichRow_ = NULL;
673    whichColumn_ = NULL;
674    nBound_ = 0;
675    saveOptions_ = rhs.saveOptions_;
676    solverOptions_ = rhs.solverOptions_;
677    int n = maximumNodes();
678    if (n) {
679      for (int i=0;i<n;i++) 
680        delete nodeInfo_[i];
681      delete [] nodeInfo_;
682    }
683    nDepth_ = rhs.nDepth_;
684    nNodes_ = rhs.nNodes_;
685    numberNodesExplored_ = rhs.numberNodesExplored_;
686    numberIterations_ = rhs.numberIterations_;
687    presolveType_ = rhs.presolveType_;
688  }
689  return *this;
690}
691// Zaps stuff 1 - arrays, 2 ints, 3 both
692void 
693ClpNodeStuff::zap(int type)
694{
695  if ((type&1)!=0) {
696    downPseudo_ = NULL;
697    upPseudo_ = NULL;
698    numberDown_ = NULL;
699    numberUp_ = NULL;
700    numberDownInfeasible_ = NULL;
701    numberUpInfeasible_ = NULL;
702    saveCosts_ = NULL;
703    nodeInfo_ = NULL;
704    large_ = NULL;
705    whichRow_ = NULL;
706    whichColumn_ = NULL;
707  }
708  if ((type&2)!=0) {
709    nBound_ = 0;
710    saveOptions_ = 0;
711    solverOptions_ = 0;
712    nDepth_ = -1;
713    nNodes_ = 0;
714    presolveType_ = 0;
715    numberNodesExplored_ = 0;
716    numberIterations_ = 0;
717  }
718}
719
720//-------------------------------------------------------------------
721// Destructor
722//-------------------------------------------------------------------
723ClpNodeStuff::~ClpNodeStuff ()
724{
725  delete [] downPseudo_;
726  delete [] upPseudo_;
727  delete [] numberDown_;
728  delete [] numberUp_;
729  delete [] numberDownInfeasible_;
730  delete [] numberUpInfeasible_;
731  int n = maximumNodes();
732  if (n) {
733    for (int i=0;i<n;i++) 
734      delete nodeInfo_[i];
735    delete [] nodeInfo_;
736  }
737}
738// Return maximum number of nodes
739int 
740ClpNodeStuff::maximumNodes() const
741{
742  int n;
743  if ((solverOptions_&32)==0) 
744    n = (1<<nDepth_)+1+nDepth_;
745  else if (nDepth_) 
746    n = 1+1+nDepth_;
747  else
748    n = 0;
749  return n;
750}
751/* Fill with pseudocosts */
752void 
753ClpNodeStuff::fillPseudoCosts(const double * down, const double * up, 
754                              const int * numberDown, const int * numberUp,
755                              const int * numberDownInfeasible, 
756                              const int * numberUpInfeasible,
757                              int number)
758{
759  delete [] downPseudo_;
760  delete [] upPseudo_;
761  delete [] numberDown_;
762  delete [] numberUp_;
763  delete [] numberDownInfeasible_;
764  delete [] numberUpInfeasible_;
765  downPseudo_ = CoinCopyOfArray(down,number);
766  upPseudo_ = CoinCopyOfArray(up,number);
767  numberDown_ = CoinCopyOfArray(numberDown,number);
768  numberUp_ = CoinCopyOfArray(numberUp,number);
769  numberDownInfeasible_ = CoinCopyOfArray(numberDownInfeasible,number);
770  numberUpInfeasible_ = CoinCopyOfArray(numberUpInfeasible,number);
771  // scale
772  for (int i=0;i<number;i++) {
773    int n;
774    n = numberDown_[i]+numberDownInfeasible_[i];
775    if (n)
776      downPseudo_[i] *= n;
777    n = numberUp_[i]+numberUpInfeasible_[i];
778    if (n)
779      upPseudo_[i] *= n;
780  }
781}
782// Update pseudo costs
783void 
784ClpNodeStuff::update(int way,int sequence,double change,bool feasible)
785{
786  if (way<0) {
787    if (feasible)
788      numberDown_[sequence]++;
789    else
790      numberDownInfeasible_[sequence]++;
791    downPseudo_[sequence] += CoinMax(change,1.0e-12);
792  } else {
793    if (feasible)
794      numberUp_[sequence]++;
795    else
796      numberUpInfeasible_[sequence]++;
797    upPseudo_[sequence] += CoinMax(change,1.0e-12);
798  }
799}
800//#############################################################################
801// Constructors / Destructor / Assignment
802//#############################################################################
803
804//-------------------------------------------------------------------
805// Default Constructor
806//-------------------------------------------------------------------
807ClpHashValue::ClpHashValue () :
808  hash_(NULL),
809  numberHash_(0),
810  maxHash_(0),
811  lastUsed_(-1)
812{
813}
814//-------------------------------------------------------------------
815// Useful Constructor from model
816//-------------------------------------------------------------------
817ClpHashValue::ClpHashValue (ClpSimplex * model) :
818  hash_(NULL),
819  numberHash_(0),
820  maxHash_(0),
821  lastUsed_(-1)
822{
823  maxHash_ = 1000;
824  int numberColumns = model->numberColumns();
825  const double * columnLower = model->columnLower();
826  const double * columnUpper = model->columnUpper();
827  int numberRows = model->numberRows();
828  const double * rowLower = model->rowLower();
829  const double * rowUpper = model->rowUpper();
830  const double * objective = model->objective();
831  CoinPackedMatrix * matrix = model->matrix();
832  const int * columnLength = matrix->getVectorLengths();
833  const CoinBigIndex * columnStart = matrix->getVectorStarts();
834  const double * elementByColumn = matrix->getElements();
835  int i;
836  int ipos;
837
838  hash_ = new CoinHashLink[maxHash_];
839 
840  for ( i = 0; i < maxHash_; i++ ) {
841    hash_[i].value = -1.0e-100;
842    hash_[i].index = -1;
843    hash_[i].next = -1;
844  }
845  // Put in +0
846  hash_[0].value=0.0;
847  hash_[0].index=0;
848  numberHash_=1;
849  /*
850   * Initialize the hash table.  Only the index of the first value that
851   * hashes to a value is entered in the table; subsequent values that
852   * collide with it are not entered.
853   */
854  for ( i = 0; i < numberColumns; i++ ) {
855    int length=columnLength[i];
856    CoinBigIndex start = columnStart[i];
857    for (CoinBigIndex i=start;i<start+length;i++) {
858      double value = elementByColumn[i];
859      ipos = hash ( value);
860      if ( hash_[ipos].index == -1 ) {
861        hash_[ipos].index = numberHash_;
862        numberHash_++;
863        hash_[ipos].value = elementByColumn[i];
864      }
865    }
866  }
867 
868  /*
869   * Now take care of the values that collided in the preceding loop,
870   * Also do other stuff
871   */
872  for ( i = 0; i < numberRows; i++ ) {
873    if (numberHash_*2>maxHash_)
874      resize(true);
875    double value;
876    value = rowLower[i];
877    ipos = index(value);
878    if (ipos<0) 
879      addValue(value);
880    value = rowUpper[i];
881    ipos = index(value);
882    if (ipos<0) 
883      addValue(value);
884  }
885  for ( i = 0; i < numberColumns; i++ ) {
886    int length=columnLength[i];
887    CoinBigIndex start = columnStart[i];
888    if (numberHash_*2>maxHash_)
889      resize(true);
890    double value;
891    value = objective[i];
892    ipos = index(value);
893    if (ipos<0) 
894      addValue(value);
895    value = columnLower[i];
896    ipos = index(value);
897    if (ipos<0) 
898      addValue(value);
899    value = columnUpper[i];
900    ipos = index(value);
901    if (ipos<0) 
902      addValue(value);
903    for (CoinBigIndex j=start;j<start+length;j++) {
904      if (numberHash_*2>maxHash_)
905        resize(true);
906      value = elementByColumn[j];
907      ipos = index(value);
908      if (ipos<0) 
909        addValue(value);
910    }
911  }
912  resize(false);
913}
914
915//-------------------------------------------------------------------
916// Copy constructor
917//-------------------------------------------------------------------
918ClpHashValue::ClpHashValue (const ClpHashValue & rhs) :
919  hash_(NULL),
920  numberHash_(rhs.numberHash_),
921  maxHash_(rhs.maxHash_),
922  lastUsed_(rhs.lastUsed_)
923{ 
924  if (maxHash_) {
925    CoinHashLink * newHash = new CoinHashLink[maxHash_];
926    int i;
927    for ( i = 0; i < maxHash_; i++ ) {
928      newHash[i].value = rhs.hash_[i].value;
929      newHash[i].index = rhs.hash_[i].index;
930      newHash[i].next = rhs.hash_[i].next;
931    }
932  }
933}
934
935//-------------------------------------------------------------------
936// Destructor
937//-------------------------------------------------------------------
938ClpHashValue::~ClpHashValue ()
939{
940  delete [] hash_;
941}
942
943//----------------------------------------------------------------
944// Assignment operator
945//-------------------------------------------------------------------
946ClpHashValue &
947ClpHashValue::operator=(const ClpHashValue& rhs)
948{
949  if (this != &rhs) {
950    numberHash_ = rhs.numberHash_;
951    maxHash_ = rhs.maxHash_;
952    lastUsed_ = rhs.lastUsed_;
953    delete [] hash_;
954    if (maxHash_) {
955      CoinHashLink * newHash = new CoinHashLink[maxHash_];
956      int i;
957      for ( i = 0; i < maxHash_; i++ ) {
958        newHash[i].value = rhs.hash_[i].value;
959        newHash[i].index = rhs.hash_[i].index;
960        newHash[i].next = rhs.hash_[i].next;
961      }
962    } else {
963      hash_ = NULL;
964    }
965  }
966  return *this;
967}
968// Return index or -1 if not found
969int 
970ClpHashValue::index(double value) const
971{
972  if (!value)
973    return 0;
974  int ipos = hash ( value);
975  int returnCode=-1;
976  while ( hash_[ipos].index>=0 ) {
977    if (value==hash_[ipos].value) {
978      returnCode = hash_[ipos].index;
979      break;
980    } else { 
981      int k = hash_[ipos].next;
982      if ( k == -1 ) {
983        break;
984      } else {
985        ipos = k;
986      }
987    }
988  }
989  return returnCode;
990}
991// Add value to list and return index
992int 
993ClpHashValue::addValue(double value) 
994{
995  int ipos = hash ( value);
996 
997  assert (value!=hash_[ipos].value);
998  if (hash_[ipos].index==-1) {
999    // can put in here
1000    hash_[ipos].index = numberHash_;
1001    numberHash_++;
1002    hash_[ipos].value = value;
1003    return numberHash_-1;
1004  }
1005  int k = hash_[ipos].next;
1006  while (k!=-1) {
1007    ipos = k;
1008    k = hash_[ipos].next;
1009  }
1010  while ( true ) {
1011    ++lastUsed_;
1012    assert (lastUsed_<=maxHash_);
1013    if ( hash_[lastUsed_].index == -1 ) {
1014      break;
1015    }
1016  }
1017  hash_[ipos].next = lastUsed_;
1018  hash_[lastUsed_].index = numberHash_;
1019  numberHash_++;
1020  hash_[lastUsed_].value = value;
1021  return numberHash_-1;
1022}
1023
1024int
1025ClpHashValue::hash ( double value) const
1026{
1027  static int mmult[] = {
1028    262139, 259459, 256889, 254291, 251701, 249133, 246709, 244247,
1029    241667, 239179, 236609, 233983, 231289, 228859, 226357, 223829,
1030    221281, 218849, 216319, 213721, 211093, 208673, 206263, 203773,
1031    201233, 198637, 196159, 193603, 191161, 188701, 186149, 183761,
1032    181303, 178873, 176389, 173897, 171469, 169049, 166471, 163871,
1033    161387, 158941, 156437, 153949, 151531, 149159, 146749, 144299,
1034    141709, 139369, 136889, 134591, 132169, 129641, 127343, 124853,
1035    122477, 120163, 117757, 115361, 112979, 110567, 108179, 105727,
1036    103387, 101021, 98639, 96179, 93911, 91583, 89317, 86939, 84521,
1037    82183, 79939, 77587, 75307, 72959, 70793, 68447, 66103
1038  };
1039  union { double d; char c[8]; } v1;
1040  assert (sizeof(double)==8);
1041  v1.d = value;
1042  int n = 0;
1043  int j;
1044
1045  for ( j = 0; j < 8; ++j ) {
1046    int ichar = v1.c[j];
1047    n += mmult[j] * ichar;
1048  }
1049  return ( abs ( n ) % maxHash_ );      /* integer abs */
1050}
1051void
1052ClpHashValue::resize(bool increaseMax)
1053{
1054  int newSize = increaseMax ? ((3*maxHash_)>>1)+1000 : maxHash_;
1055  CoinHashLink * newHash = new CoinHashLink[newSize];
1056  int i;
1057  for ( i = 0; i < newSize; i++ ) {
1058    newHash[i].value = -1.0e-100;
1059    newHash[i].index = -1;
1060    newHash[i].next = -1;
1061  }
1062  // swap
1063  CoinHashLink * oldHash = hash_;
1064  hash_ = newHash;
1065  int oldSize = maxHash_;
1066  maxHash_ = newSize;
1067  /*
1068   * Initialize the hash table.  Only the index of the first value that
1069   * hashes to a value is entered in the table; subsequent values that
1070   * collide with it are not entered.
1071   */
1072  int ipos;
1073  int n=0;
1074  for ( i = 0; i < oldSize; i++ ) {
1075    if (oldHash[i].index>=0) {
1076      ipos = hash ( oldHash[i].value);
1077      if ( hash_[ipos].index == -1 ) {
1078        hash_[ipos].index = n;
1079        n++;
1080        hash_[ipos].value = oldHash[i].value;
1081        // unmark
1082        oldHash[i].index=-1;
1083      }
1084    }
1085  }
1086  /*
1087   * Now take care of the values that collided in the preceding loop,
1088   * by finding some other entry in the table for them.
1089   * Since there are as many entries in the table as there are values,
1090   * there must be room for them.
1091   */
1092  lastUsed_ = -1;
1093  for ( i = 0; i < oldSize; ++i ) {
1094    if (oldHash[i].index>=0) {
1095      double value = oldHash[i].value;
1096      ipos = hash ( value);
1097      int k;
1098      while ( true ) {
1099        assert (value!=hash_[ipos].value);
1100        k = hash_[ipos].next;
1101        if ( k == -1 ) {
1102          while ( true ) {
1103            ++lastUsed_;
1104            assert (lastUsed_<=maxHash_);
1105            if ( hash_[lastUsed_].index == -1 ) {
1106              break;
1107            }
1108          }
1109          hash_[ipos].next = lastUsed_;
1110          hash_[lastUsed_].index = n;
1111          n++;
1112          hash_[lastUsed_].value = value;
1113          break;
1114        } else {
1115          ipos = k;
1116        }
1117      }
1118    }
1119  }
1120  assert (n==numberHash_);
1121  delete [] oldHash;
1122}
Note: See TracBrowser for help on using the repository browser.