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

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

changes for fathoming

  • Property svn:executable set to *
  • Property svn:keywords set to Id
File size: 36.1 KB
Line 
1/* $Id: ClpNode.cpp 1430 2009-09-15 15:51:17Z 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}
854// Return maximum number of nodes
855int 
856ClpNodeStuff::maximumNodes() const
857{
858  int n=0;
859#if 0
860  if (nDepth_!=-1) {
861    if ((solverOptions_&32)==0)
862      n = (1<<nDepth_);
863    else if (nDepth_)
864      n = 1;
865  }
866  assert (n==maximumNodes_-(1+nDepth_)||n==0);
867#else
868  if (nDepth_!=-1) {
869    n = maximumNodes_-(1+nDepth_);
870    assert (n>0);
871  }
872#endif
873  return n;
874}
875// Return maximum space for nodes
876int 
877ClpNodeStuff::maximumSpace() const
878{
879  return maximumNodes_;
880}
881/* Fill with pseudocosts */
882void 
883ClpNodeStuff::fillPseudoCosts(const double * down, const double * up, 
884                              const int * priority,
885                              const int * numberDown, const int * numberUp,
886                              const int * numberDownInfeasible, 
887                              const int * numberUpInfeasible,
888                              int number)
889{
890  delete [] downPseudo_;
891  delete [] upPseudo_;
892  delete [] priority_;
893  delete [] numberDown_;
894  delete [] numberUp_;
895  delete [] numberDownInfeasible_;
896  delete [] numberUpInfeasible_;
897  downPseudo_ = CoinCopyOfArray(down,number);
898  upPseudo_ = CoinCopyOfArray(up,number);
899  priority_ = CoinCopyOfArray(priority,number);
900  numberDown_ = CoinCopyOfArray(numberDown,number);
901  numberUp_ = CoinCopyOfArray(numberUp,number);
902  numberDownInfeasible_ = CoinCopyOfArray(numberDownInfeasible,number);
903  numberUpInfeasible_ = CoinCopyOfArray(numberUpInfeasible,number);
904  // scale
905  for (int i=0;i<number;i++) {
906    int n;
907    n = numberDown_[i];
908    if (n)
909      downPseudo_[i] *= n;
910    n = numberUp_[i];
911    if (n)
912      upPseudo_[i] *= n;
913  }
914}
915// Update pseudo costs
916void 
917ClpNodeStuff::update(int way,int sequence,double change,bool feasible)
918{
919  assert (numberDown_[sequence]>=numberDownInfeasible_[sequence]);
920  assert (numberUp_[sequence]>=numberUpInfeasible_[sequence]);
921  if (way<0) {
922    numberDown_[sequence]++;
923    if (!feasible)
924      numberDownInfeasible_[sequence]++;
925    downPseudo_[sequence] += CoinMax(change,1.0e-12);
926  } else {
927    numberUp_[sequence]++;
928    if (!feasible)
929      numberUpInfeasible_[sequence]++;
930    upPseudo_[sequence] += CoinMax(change,1.0e-12);
931  }
932}
933//#############################################################################
934// Constructors / Destructor / Assignment
935//#############################################################################
936
937//-------------------------------------------------------------------
938// Default Constructor
939//-------------------------------------------------------------------
940ClpHashValue::ClpHashValue () :
941  hash_(NULL),
942  numberHash_(0),
943  maxHash_(0),
944  lastUsed_(-1)
945{
946}
947//-------------------------------------------------------------------
948// Useful Constructor from model
949//-------------------------------------------------------------------
950ClpHashValue::ClpHashValue (ClpSimplex * model) :
951  hash_(NULL),
952  numberHash_(0),
953  maxHash_(0),
954  lastUsed_(-1)
955{
956  maxHash_ = 1000;
957  int numberColumns = model->numberColumns();
958  const double * columnLower = model->columnLower();
959  const double * columnUpper = model->columnUpper();
960  int numberRows = model->numberRows();
961  const double * rowLower = model->rowLower();
962  const double * rowUpper = model->rowUpper();
963  const double * objective = model->objective();
964  CoinPackedMatrix * matrix = model->matrix();
965  const int * columnLength = matrix->getVectorLengths();
966  const CoinBigIndex * columnStart = matrix->getVectorStarts();
967  const double * elementByColumn = matrix->getElements();
968  int i;
969  int ipos;
970
971  hash_ = new CoinHashLink[maxHash_];
972 
973  for ( i = 0; i < maxHash_; i++ ) {
974    hash_[i].value = -1.0e-100;
975    hash_[i].index = -1;
976    hash_[i].next = -1;
977  }
978  // Put in +0
979  hash_[0].value=0.0;
980  hash_[0].index=0;
981  numberHash_=1;
982  /*
983   * Initialize the hash table.  Only the index of the first value that
984   * hashes to a value is entered in the table; subsequent values that
985   * collide with it are not entered.
986   */
987  for ( i = 0; i < numberColumns; i++ ) {
988    int length=columnLength[i];
989    CoinBigIndex start = columnStart[i];
990    for (CoinBigIndex i=start;i<start+length;i++) {
991      double value = elementByColumn[i];
992      ipos = hash ( value);
993      if ( hash_[ipos].index == -1 ) {
994        hash_[ipos].index = numberHash_;
995        numberHash_++;
996        hash_[ipos].value = elementByColumn[i];
997      }
998    }
999  }
1000 
1001  /*
1002   * Now take care of the values that collided in the preceding loop,
1003   * Also do other stuff
1004   */
1005  for ( i = 0; i < numberRows; i++ ) {
1006    if (numberHash_*2>maxHash_)
1007      resize(true);
1008    double value;
1009    value = rowLower[i];
1010    ipos = index(value);
1011    if (ipos<0) 
1012      addValue(value);
1013    value = rowUpper[i];
1014    ipos = index(value);
1015    if (ipos<0) 
1016      addValue(value);
1017  }
1018  for ( i = 0; i < numberColumns; i++ ) {
1019    int length=columnLength[i];
1020    CoinBigIndex start = columnStart[i];
1021    if (numberHash_*2>maxHash_)
1022      resize(true);
1023    double value;
1024    value = objective[i];
1025    ipos = index(value);
1026    if (ipos<0) 
1027      addValue(value);
1028    value = columnLower[i];
1029    ipos = index(value);
1030    if (ipos<0) 
1031      addValue(value);
1032    value = columnUpper[i];
1033    ipos = index(value);
1034    if (ipos<0) 
1035      addValue(value);
1036    for (CoinBigIndex j=start;j<start+length;j++) {
1037      if (numberHash_*2>maxHash_)
1038        resize(true);
1039      value = elementByColumn[j];
1040      ipos = index(value);
1041      if (ipos<0) 
1042        addValue(value);
1043    }
1044  }
1045  resize(false);
1046}
1047
1048//-------------------------------------------------------------------
1049// Copy constructor
1050//-------------------------------------------------------------------
1051ClpHashValue::ClpHashValue (const ClpHashValue & rhs) :
1052  hash_(NULL),
1053  numberHash_(rhs.numberHash_),
1054  maxHash_(rhs.maxHash_),
1055  lastUsed_(rhs.lastUsed_)
1056{ 
1057  if (maxHash_) {
1058    CoinHashLink * newHash = new CoinHashLink[maxHash_];
1059    int i;
1060    for ( i = 0; i < maxHash_; i++ ) {
1061      newHash[i].value = rhs.hash_[i].value;
1062      newHash[i].index = rhs.hash_[i].index;
1063      newHash[i].next = rhs.hash_[i].next;
1064    }
1065  }
1066}
1067
1068//-------------------------------------------------------------------
1069// Destructor
1070//-------------------------------------------------------------------
1071ClpHashValue::~ClpHashValue ()
1072{
1073  delete [] hash_;
1074}
1075
1076//----------------------------------------------------------------
1077// Assignment operator
1078//-------------------------------------------------------------------
1079ClpHashValue &
1080ClpHashValue::operator=(const ClpHashValue& rhs)
1081{
1082  if (this != &rhs) {
1083    numberHash_ = rhs.numberHash_;
1084    maxHash_ = rhs.maxHash_;
1085    lastUsed_ = rhs.lastUsed_;
1086    delete [] hash_;
1087    if (maxHash_) {
1088      CoinHashLink * newHash = new CoinHashLink[maxHash_];
1089      int i;
1090      for ( i = 0; i < maxHash_; i++ ) {
1091        newHash[i].value = rhs.hash_[i].value;
1092        newHash[i].index = rhs.hash_[i].index;
1093        newHash[i].next = rhs.hash_[i].next;
1094      }
1095    } else {
1096      hash_ = NULL;
1097    }
1098  }
1099  return *this;
1100}
1101// Return index or -1 if not found
1102int 
1103ClpHashValue::index(double value) const
1104{
1105  if (!value)
1106    return 0;
1107  int ipos = hash ( value);
1108  int returnCode=-1;
1109  while ( hash_[ipos].index>=0 ) {
1110    if (value==hash_[ipos].value) {
1111      returnCode = hash_[ipos].index;
1112      break;
1113    } else { 
1114      int k = hash_[ipos].next;
1115      if ( k == -1 ) {
1116        break;
1117      } else {
1118        ipos = k;
1119      }
1120    }
1121  }
1122  return returnCode;
1123}
1124// Add value to list and return index
1125int 
1126ClpHashValue::addValue(double value) 
1127{
1128  int ipos = hash ( value);
1129 
1130  assert (value!=hash_[ipos].value);
1131  if (hash_[ipos].index==-1) {
1132    // can put in here
1133    hash_[ipos].index = numberHash_;
1134    numberHash_++;
1135    hash_[ipos].value = value;
1136    return numberHash_-1;
1137  }
1138  int k = hash_[ipos].next;
1139  while (k!=-1) {
1140    ipos = k;
1141    k = hash_[ipos].next;
1142  }
1143  while ( true ) {
1144    ++lastUsed_;
1145    assert (lastUsed_<=maxHash_);
1146    if ( hash_[lastUsed_].index == -1 ) {
1147      break;
1148    }
1149  }
1150  hash_[ipos].next = lastUsed_;
1151  hash_[lastUsed_].index = numberHash_;
1152  numberHash_++;
1153  hash_[lastUsed_].value = value;
1154  return numberHash_-1;
1155}
1156
1157int
1158ClpHashValue::hash ( double value) const
1159{
1160  static int mmult[] = {
1161    262139, 259459, 256889, 254291, 251701, 249133, 246709, 244247,
1162    241667, 239179, 236609, 233983, 231289, 228859, 226357, 223829,
1163    221281, 218849, 216319, 213721, 211093, 208673, 206263, 203773,
1164    201233, 198637, 196159, 193603, 191161, 188701, 186149, 183761,
1165    181303, 178873, 176389, 173897, 171469, 169049, 166471, 163871,
1166    161387, 158941, 156437, 153949, 151531, 149159, 146749, 144299,
1167    141709, 139369, 136889, 134591, 132169, 129641, 127343, 124853,
1168    122477, 120163, 117757, 115361, 112979, 110567, 108179, 105727,
1169    103387, 101021, 98639, 96179, 93911, 91583, 89317, 86939, 84521,
1170    82183, 79939, 77587, 75307, 72959, 70793, 68447, 66103
1171  };
1172  union { double d; char c[8]; } v1;
1173  assert (sizeof(double)==8);
1174  v1.d = value;
1175  int n = 0;
1176  int j;
1177
1178  for ( j = 0; j < 8; ++j ) {
1179    int ichar = v1.c[j];
1180    n += mmult[j] * ichar;
1181  }
1182  return ( abs ( n ) % maxHash_ );      /* integer abs */
1183}
1184void
1185ClpHashValue::resize(bool increaseMax)
1186{
1187  int newSize = increaseMax ? ((3*maxHash_)>>1)+1000 : maxHash_;
1188  CoinHashLink * newHash = new CoinHashLink[newSize];
1189  int i;
1190  for ( i = 0; i < newSize; i++ ) {
1191    newHash[i].value = -1.0e-100;
1192    newHash[i].index = -1;
1193    newHash[i].next = -1;
1194  }
1195  // swap
1196  CoinHashLink * oldHash = hash_;
1197  hash_ = newHash;
1198  int oldSize = maxHash_;
1199  maxHash_ = newSize;
1200  /*
1201   * Initialize the hash table.  Only the index of the first value that
1202   * hashes to a value is entered in the table; subsequent values that
1203   * collide with it are not entered.
1204   */
1205  int ipos;
1206  int n=0;
1207  for ( i = 0; i < oldSize; i++ ) {
1208    if (oldHash[i].index>=0) {
1209      ipos = hash ( oldHash[i].value);
1210      if ( hash_[ipos].index == -1 ) {
1211        hash_[ipos].index = n;
1212        n++;
1213        hash_[ipos].value = oldHash[i].value;
1214        // unmark
1215        oldHash[i].index=-1;
1216      }
1217    }
1218  }
1219  /*
1220   * Now take care of the values that collided in the preceding loop,
1221   * by finding some other entry in the table for them.
1222   * Since there are as many entries in the table as there are values,
1223   * there must be room for them.
1224   */
1225  lastUsed_ = -1;
1226  for ( i = 0; i < oldSize; ++i ) {
1227    if (oldHash[i].index>=0) {
1228      double value = oldHash[i].value;
1229      ipos = hash ( value);
1230      int k;
1231      while ( true ) {
1232        assert (value!=hash_[ipos].value);
1233        k = hash_[ipos].next;
1234        if ( k == -1 ) {
1235          while ( true ) {
1236            ++lastUsed_;
1237            assert (lastUsed_<=maxHash_);
1238            if ( hash_[lastUsed_].index == -1 ) {
1239              break;
1240            }
1241          }
1242          hash_[ipos].next = lastUsed_;
1243          hash_[lastUsed_].index = n;
1244          n++;
1245          hash_[lastUsed_].value = value;
1246          break;
1247        } else {
1248          ipos = k;
1249        }
1250      }
1251    }
1252  }
1253  assert (n==numberHash_);
1254  delete [] oldHash;
1255}
Note: See TracBrowser for help on using the repository browser.