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

Last change on this file since 1484 was 1484, checked in by forrest, 10 years ago

memory leak

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