source: stable/1.7/Clp/src/ClpNode.cpp @ 1213

// Copyright (C) 2008, International Business Machines
// Corporation and others.  All Rights Reserved.

#include "CoinPragma.hpp"
#include "ClpSimplex.hpp"
#include "ClpNode.hpp"
#include "ClpFactorization.hpp"
#include "ClpDualRowSteepest.hpp"

//#############################################################################
// Constructors / Destructor / Assignment
//#############################################################################

//-------------------------------------------------------------------
// Default Constructor 
//-------------------------------------------------------------------
ClpNode::ClpNode () :
  branchingValue_(0.5),
  objectiveValue_(0.0),
  sumInfeasibilities_(0.0),
  estimatedSolution_(0.0),
  factorization_(NULL),
  weights_(NULL),
  status_(NULL),
  primalSolution_(NULL),
  dualSolution_(NULL),
  lower_(NULL),
  upper_(NULL),
  pivotVariables_(NULL),
  fixed_(NULL),
  sequence_(1),
  numberInfeasibilities_(0),
  numberFixed_(0),
  flags_(0),
  maximumFixed_(0),
  maximumRows_(0),
  maximumColumns_(0),
  maximumIntegers_(0)
{
 branchState_.firstBranch=0;
 branchState_.branch=0;
}
//-------------------------------------------------------------------
// Useful Constructor from model
//-------------------------------------------------------------------
ClpNode::ClpNode (ClpSimplex * model, const ClpNodeStuff * stuff, int depth) :
  branchingValue_(0.5),
  objectiveValue_(0.0),
  sumInfeasibilities_(0.0),
  estimatedSolution_(0.0),
  factorization_(NULL),
  weights_(NULL),
  status_(NULL),
  primalSolution_(NULL),
  dualSolution_(NULL),
  lower_(NULL),
  upper_(NULL),
  pivotVariables_(NULL),
  fixed_(NULL),
  sequence_(1),
  numberInfeasibilities_(0),
  numberFixed_(0),
  flags_(0),
  maximumFixed_(0),
  maximumRows_(0),
  maximumColumns_(0),
  maximumIntegers_(0)
{
  branchState_.firstBranch=0;
  branchState_.branch=0;
  gutsOfConstructor(model,stuff,0,depth);
}

//-------------------------------------------------------------------
// Most of work of constructor from model
//-------------------------------------------------------------------
void
ClpNode::gutsOfConstructor (ClpSimplex * model, const ClpNodeStuff * stuff,
                            int arraysExist,int depth) 
{
  int numberRows = model->numberRows();
  int numberColumns = model->numberColumns();
  int numberTotal = numberRows+numberColumns;
  int maximumTotal = maximumRows_+maximumColumns_;
  // save stuff
  objectiveValue_ = model->objectiveValue()*model->optimizationDirection();
  estimatedSolution_ = objectiveValue_;
  flags_ = 1; //say scaled
  if (!arraysExist) {
    maximumRows_=CoinMax(maximumRows_,numberRows);
    maximumColumns_=CoinMax(maximumColumns_,numberColumns);
    maximumTotal = maximumRows_+maximumColumns_;
    factorization_ = new ClpFactorization(*model->factorization());
    status_ = CoinCopyOfArrayPartial(model->statusArray(),maximumTotal,numberTotal);
    primalSolution_ = CoinCopyOfArrayPartial(model->solutionRegion(),maximumTotal,numberTotal);
    dualSolution_ = CoinCopyOfArrayPartial(model->djRegion(),maximumTotal,numberTotal); //? has duals as well?
    pivotVariables_ = CoinCopyOfArrayPartial(model->pivotVariable(),maximumRows_,numberRows); 
    //#define DUAL_WEIGHTS
#ifdef DUAL_WEIGHTS
    ClpDualRowSteepest* pivot =
      dynamic_cast< ClpDualRowSteepest*>(model->dualRowPivot());
    if (pivot) {
      assert (!weights_);
      weights_ = new ClpDualRowSteepest(*pivot);
    }
#endif
  } else {
    if (arraysExist==2)
      assert(lower_);
    if (numberRows<=maximumRows_&&numberColumns<=maximumColumns_) {
      CoinMemcpyN(model->statusArray(),numberTotal,status_);
      if (arraysExist==1) {
        *factorization_ = *model->factorization();
        CoinMemcpyN(model->solutionRegion(),numberTotal,primalSolution_);
        CoinMemcpyN(model->djRegion(),numberTotal,dualSolution_); //? has duals as well?
#ifdef DUAL_WEIGHTS
        ClpDualRowSteepest* pivot =
          dynamic_cast< ClpDualRowSteepest*>(model->dualRowPivot());
        if (pivot) {
          if (weights_) {
            //if (weights_->numberRows()==pivot->numberRows()) {
              weights_->fill(*pivot);
              //} else {
              //delete weights_;
              //weights_ = new ClpDualRowSteepest(*pivot);
              //}
          } else {
            weights_ = new ClpDualRowSteepest(*pivot);
          }
        }
#endif
        CoinMemcpyN(model->pivotVariable(),numberRows,pivotVariables_); 
      } else {
        CoinMemcpyN(model->primalColumnSolution(),numberColumns,primalSolution_);
        CoinMemcpyN(model->dualColumnSolution(),numberColumns,dualSolution_);
        flags_ = 0;
        CoinMemcpyN(model->dualRowSolution(),numberRows,dualSolution_+numberColumns); 
      }
    } else {
      // size has changed
      maximumRows_=CoinMax(maximumRows_,numberRows);
      maximumColumns_=CoinMax(maximumColumns_,numberColumns);
      maximumTotal = maximumRows_+maximumColumns_;
      delete weights_;
      weights_ = NULL;
      delete [] status_;
      delete [] primalSolution_;
      delete [] dualSolution_;
      delete [] pivotVariables_;
      status_ = CoinCopyOfArrayPartial(model->statusArray(),maximumTotal,numberTotal);
      primalSolution_ = new double [maximumTotal*sizeof(double)];
      dualSolution_ = new double [maximumTotal*sizeof(double)];
      if (arraysExist==1) {
        *factorization_ = *model->factorization(); // I think this is OK
        CoinMemcpyN(model->solutionRegion(),numberTotal,primalSolution_);
        CoinMemcpyN(model->djRegion(),numberTotal,dualSolution_); //? has duals as well?
#ifdef DUAL_WEIGHTS
        ClpDualRowSteepest* pivot =
          dynamic_cast< ClpDualRowSteepest*>(model->dualRowPivot());
        if (pivot) {
          assert (!weights_);
          weights_ = new ClpDualRowSteepest(*pivot);
        }
#endif
      } else {
        CoinMemcpyN(model->primalColumnSolution(),numberColumns,primalSolution_);
        CoinMemcpyN(model->dualColumnSolution(),numberColumns,dualSolution_); 
        flags_ = 0;
        CoinMemcpyN(model->dualRowSolution(),numberRows,dualSolution_+numberColumns); 
      }
      pivotVariables_ = new int [maximumRows_];
      if (model->pivotVariable())
        CoinMemcpyN(model->pivotVariable(),numberRows,pivotVariables_); 
    }
  }
  numberFixed_=0;
  const double * lower = model->columnLower();
  const double * upper = model->columnUpper();
  const double * solution = model->primalColumnSolution();
  const char * integerType = model->integerInformation();
  const double * columnScale = model->columnScale();
  if (!flags_)
    columnScale=NULL; // as duals correct
  int iColumn;
  sequence_=-1;
  double integerTolerance = stuff->integerTolerance_;
  double mostAway=integerTolerance;
  sumInfeasibilities_ = 0.0;
  numberInfeasibilities_ = 0;
  int nFix=0;
  double gap = CoinMax(model->dualObjectiveLimit()-objectiveValue_,1.0e-4);
#define PSEUDO 2
#if PSEUDO==1
  // Column copy of matrix
  ClpPackedMatrix * matrix = model->clpScaledMatrix();
  if (!matrix)
    matrix = dynamic_cast< ClpPackedMatrix*>(model->clpMatrix());
  const double * element = matrix->getElements();
  const int * row = matrix->getIndices();
  const CoinBigIndex * columnStart = matrix->getVectorStarts();
  const int * columnLength = matrix->getVectorLengths();
  const double *objective = model->costRegion() ;
  double direction = model->optimizationDirection();
  const double * dual = dualSolution_+numberColumns;
#elif PSEUDO==2
  const double * downPseudo = stuff->downPseudo_;
  const int * numberDown = stuff->numberDown_;
  const int * numberDownInfeasible = stuff->numberDownInfeasible_;
  const double * upPseudo = stuff->upPseudo_;
  const int * numberUp = stuff->numberUp_;
  const int * numberUpInfeasible = stuff->numberUpInfeasible_;
#endif
  int iInteger=0;
  for (iColumn=0;iColumn<numberColumns;iColumn++) {
    if (integerType[iColumn]) {
      double value = solution[iColumn];
      value = max(value,(double) lower[iColumn]);
      value = min(value,(double) upper[iColumn]);
      double nearest = floor(value+0.5);
      if (fabs(value-nearest)>integerTolerance) {
        numberInfeasibilities_++;
        sumInfeasibilities_ += fabs(value-nearest);
#if PSEUDO==1
        double upValue = 0.0;
        double downValue = 0.0;
        double value2 = direction*objective[iColumn];
        if (value2) {
          if (value2>0.0)
            upValue += value2;
          else
            downValue -= value2;
        }
        CoinBigIndex start = columnStart[iColumn];
        CoinBigIndex end = columnStart[iColumn]+columnLength[iColumn];
        for (CoinBigIndex j=start;j<end;j++) {
          int iRow = row[j];
          value2 = -dual[iRow];
          if (value2) {
            value2 *= element[j];
            if (value2>0.0)
              upValue += value2;
            else
              downValue -= value2;
          }
        }
        upValue *= ceil(value)-value;
        downValue *= value-floor(value);
        double infeasibility;
        if (depth>1000)
          infeasibility = CoinMax(upValue,downValue)+integerTolerance;
        else
          infeasibility = 0.1*CoinMax(upValue,downValue)+
            0.9*CoinMin(upValue,downValue) + integerTolerance;
#elif PSEUDO==2
        double upValue = (ceil(value)-value)*(upPseudo[iInteger]/(1.0+numberUp[iInteger]+numberUpInfeasible[iInteger]));
        double downValue = (value-floor(value))*(downPseudo[iInteger]/(1.0+numberDown[iInteger]+numberDownInfeasible[iInteger]));
        double infeasibility;
        if (depth>1000)
          infeasibility = CoinMax(upValue,downValue)+integerTolerance;
        else
          infeasibility = 0.1*CoinMax(upValue,downValue)+
            0.9*CoinMin(upValue,downValue) + integerTolerance;
        estimatedSolution_ += CoinMin(upValue,downValue);
#else
        double infeasibility = fabs(value-nearest);
#endif
        if (infeasibility>mostAway) {
          mostAway=infeasibility;
          sequence_=iColumn;
          branchingValue_=value;
          branchState_.branch=0;
#if PSEUDO>0
          if (upValue<=downValue)
            branchState_.firstBranch=1; // up
          else
            branchState_.firstBranch=0; // down
#else
          if (value<=nearest)
            branchState_.firstBranch=1; // up
          else
            branchState_.firstBranch=0; // down
#endif
        }
      } else if (model->getColumnStatus(iColumn)==ClpSimplex::atLowerBound) {
        bool fix=false;
        if (columnScale) {
          if (dualSolution_[iColumn]>gap*columnScale[iColumn]) 
            fix=true;
        } else {
          if (dualSolution_[iColumn]>gap) 
            fix=true;
        }
        if (fix) {
          nFix++;
          //printf("fixed %d to zero gap %g dj %g %g\n",iColumn,
          // gap,dualSolution_[iColumn], columnScale ? columnScale[iColumn]:1.0);
          model->setColumnStatus(iColumn,ClpSimplex::isFixed);
        }
      } else if (model->getColumnStatus(iColumn)==ClpSimplex::atUpperBound) {
        bool fix=false;
        if (columnScale) {
          if (-dualSolution_[iColumn]>gap*columnScale[iColumn]) 
            fix=true;
        } else {
          if (-dualSolution_[iColumn]>gap) 
            fix=true;
        }
        if (fix) {
          nFix++;
          //printf("fixed %d to one gap %g dj %g %g\n",iColumn,
          // gap,dualSolution_[iColumn], columnScale ? columnScale[iColumn]:1.0);
          model->setColumnStatus(iColumn,ClpSimplex::isFixed);
        }
      }
      iInteger++;
    }
  }
  if (lower_) {
    // save bounds
    if (iInteger>maximumIntegers_) {
      delete [] lower_;
      delete [] upper_;
      maximumIntegers_ = iInteger;
      lower_ = new int [maximumIntegers_];
      upper_ = new int [maximumIntegers_];
    }
    iInteger=0;
    for (iColumn=0;iColumn<numberColumns;iColumn++) {
      if (integerType[iColumn]) {
        lower_[iInteger]=(int) lower[iColumn];
        upper_[iInteger]=(int) upper[iColumn];
        iInteger++;
      }
    }
  }
  // Could omit save of fixed if doing full save of bounds
  if (sequence_>=0&&nFix) {
    if (nFix>maximumFixed_) {
      delete [] fixed_;
      fixed_ = new int [nFix];
      maximumFixed_=nFix;
    }
    numberFixed_=0;
    unsigned char * status = model->statusArray();
    for (iColumn=0;iColumn<numberColumns;iColumn++) {
      if (status[iColumn]!=status_[iColumn]) {
        if (solution[iColumn]<=lower[iColumn]+2.0*integerTolerance) {
          model->setColumnUpper(iColumn,lower[iColumn]);
          fixed_[numberFixed_++]=iColumn;
        } else {
          assert (solution[iColumn]>=upper[iColumn]-2.0*integerTolerance);
          model->setColumnLower(iColumn,upper[iColumn]);
          fixed_[numberFixed_++]=iColumn|0x10000000;
        }
      }
    }
    //printf("%d fixed\n",numberFixed_);
  }
}

//-------------------------------------------------------------------
// Copy constructor 
//-------------------------------------------------------------------
ClpNode::ClpNode (const ClpNode & source) 
{  
  printf("ClpNode copy not implemented\n");
  abort();
}

//-------------------------------------------------------------------
// Destructor 
//-------------------------------------------------------------------
ClpNode::~ClpNode ()
{
  delete factorization_;
  delete weights_;
  delete [] status_;
  delete [] primalSolution_;
  delete [] dualSolution_;
  delete [] lower_;
  delete [] upper_;
  delete [] pivotVariables_;
  delete [] fixed_;
}

//----------------------------------------------------------------
// Assignment operator 
//-------------------------------------------------------------------
ClpNode &
ClpNode::operator=(const ClpNode& rhs)
{
  if (this != &rhs) {
    printf("ClpNode = not implemented\n");
    abort();
  }
  return *this;
}
// Create odd arrays
void 
ClpNode::createArrays(ClpSimplex * model)
{
  int numberColumns = model->numberColumns();
  const char * integerType = model->integerInformation();
  int iColumn;
  int numberIntegers=0;
  for (iColumn=0;iColumn<numberColumns;iColumn++) {
    if (integerType[iColumn]) 
      numberIntegers++;
  }
  if (numberIntegers>maximumIntegers_||!lower_) {
    delete [] lower_;
    delete [] upper_;
    maximumIntegers_ = numberIntegers;
    lower_ = new int [numberIntegers];
    upper_ = new int [numberIntegers];
  }
}
// Clean up as crunch is different model
void 
ClpNode::cleanUpForCrunch()
{
  delete weights_;
  weights_ = NULL;
}
/* Applies node to model
   0 - just tree bounds
   1 - tree bounds and basis etc
   2 - saved bounds and basis etc
*/
void 
ClpNode::applyNode(ClpSimplex * model, int doBoundsEtc )
{
  int numberColumns = model->numberColumns();
  const double * lower = model->columnLower();
  const double * upper = model->columnUpper();
  if (doBoundsEtc<2) {
    // current bound
    int way=branchState_.firstBranch;
    if (branchState_.branch>0)
      way=1-way;
    if (!way) {
      // This should also do underlying internal bound
      model->setColumnUpper(sequence_,floor(branchingValue_));
    } else {
      // This should also do underlying internal bound
      model->setColumnLower(sequence_,ceil(branchingValue_));
    }
    // apply dj fixings
    for (int i=0;i<numberFixed_;i++) {
      int iColumn = fixed_[i];
      if ((iColumn&0x10000000)!=0) {
        iColumn &= 0xfffffff;
        model->setColumnLower(iColumn,upper[iColumn]);
    } else {
        model->setColumnUpper(iColumn,lower[iColumn]);
      }
    }
  } else {
    // restore bounds
    assert (lower_);
    int iInteger=-1;
    const char * integerType = model->integerInformation();
    for (int iColumn=0;iColumn<numberColumns;iColumn++) {
      if (integerType[iColumn]) {
        iInteger++;
        if (lower_[iInteger]!=(int) lower[iColumn])
          model->setColumnLower(iColumn,lower_[iInteger]);
        if (upper_[iInteger]!=(int) upper[iColumn])
          model->setColumnUpper(iColumn,upper_[iInteger]);
      }
    }
  }
  if (doBoundsEtc&&doBoundsEtc<3) {
#ifdef DUAL_WEIGHTS
    model->copyFactorization(*factorization_);
    ClpDualRowSteepest* pivot =
      dynamic_cast< ClpDualRowSteepest*>(model->dualRowPivot());
    if (pivot&&weights_) {
      pivot->fill(*weights_); 
    }
#endif
    int numberRows = model->numberRows();
    int numberTotal = numberRows+numberColumns;
    CoinMemcpyN(status_,numberTotal,model->statusArray());
    if (doBoundsEtc<2) {
      CoinMemcpyN(primalSolution_,numberTotal,model->solutionRegion());
      CoinMemcpyN(dualSolution_,numberTotal,model->djRegion());
      CoinMemcpyN(pivotVariables_,numberRows,model->pivotVariable());
      CoinMemcpyN(dualSolution_+numberColumns,numberRows,model->dualRowSolution());
    } else {
      CoinMemcpyN(primalSolution_,numberColumns,model->primalColumnSolution());
      CoinMemcpyN(dualSolution_,numberColumns,model->dualColumnSolution());
      CoinMemcpyN(dualSolution_+numberColumns,numberRows,model->dualRowSolution());
      if (model->columnScale()) {
        // See if just primal will work
        double * solution = model->primalColumnSolution();
        const double * columnScale = model->columnScale();
        int i;
        for (i=0;i<numberColumns;i++) {
          solution[i] *= columnScale[i];
        }
      }
    }
    model->setObjectiveValue(objectiveValue_);
  }
}
// Fix on reduced costs
int 
ClpNode::fixOnReducedCosts(ClpSimplex * model)
{
  return 0;
}
/* Way for integer variable -1 down , +1 up */
int 
ClpNode::way() const
{
  int way=branchState_.firstBranch;
  if (branchState_.branch>0)
    way=1-way;
  return way==0 ? -1 : +1;
}
// Return true if branch exhausted
bool 
ClpNode::fathomed() const
{ 
  return branchState_.branch>=1
;
}
// Change state of variable i.e. go other way
void 
ClpNode::changeState()
{
  branchState_.branch++;
  assert (branchState_.branch<=2);
}
//#############################################################################
// Constructors / Destructor / Assignment
//#############################################################################

//-------------------------------------------------------------------
// Default Constructor 
//-------------------------------------------------------------------
ClpNodeStuff::ClpNodeStuff () :
  integerTolerance_(1.0e-7),
  integerIncrement_(1.0e-8),
  downPseudo_(NULL),
  upPseudo_(NULL),
  numberDown_(NULL),
  numberUp_(NULL),
  numberDownInfeasible_(NULL),
  numberUpInfeasible_(NULL),
  saveCosts_(NULL),
  nodeInfo_(NULL),
  large_(NULL),
  whichRow_(NULL),
  whichColumn_(NULL),
  nBound_(0),
  saveOptions_(0),
  solverOptions_(0),
  nDepth_(-1),
  nNodes_(0),
  numberNodesExplored_(0),
  numberIterations_(0),
  presolveType_(0)
{

}

//-------------------------------------------------------------------
// Copy constructor 
//-------------------------------------------------------------------
ClpNodeStuff::ClpNodeStuff (const ClpNodeStuff & rhs) 
  : integerTolerance_(rhs.integerTolerance_),
    integerIncrement_(rhs.integerIncrement_),
    downPseudo_(NULL),
    upPseudo_(NULL),
    numberDown_(NULL),
    numberUp_(NULL),
    numberDownInfeasible_(NULL),
    numberUpInfeasible_(NULL),
    saveCosts_(NULL),
    nodeInfo_(NULL),
    large_(NULL),
    whichRow_(NULL),
    whichColumn_(NULL),
    nBound_(0),
    saveOptions_(rhs.saveOptions_),
    solverOptions_(rhs.solverOptions_),
    nDepth_(rhs.nDepth_),
    nNodes_(rhs.nNodes_),
    numberNodesExplored_(rhs.numberNodesExplored_),
    numberIterations_(rhs.numberIterations_),
    presolveType_(rhs.presolveType_)
{  
}
//----------------------------------------------------------------
// Assignment operator 
//-------------------------------------------------------------------
ClpNodeStuff &
ClpNodeStuff::operator=(const ClpNodeStuff& rhs)
{
  if (this != &rhs) {
    integerTolerance_ = rhs.integerTolerance_;
    integerIncrement_ = rhs.integerIncrement_;
    downPseudo_ = NULL;
    upPseudo_ = NULL;
    numberDown_ = NULL;
    numberUp_ = NULL;
    numberDownInfeasible_ = NULL;
    numberUpInfeasible_ = NULL;
    saveCosts_ = NULL;
    nodeInfo_ = NULL;
    large_ = NULL;
    whichRow_ = NULL;
    whichColumn_ = NULL;
    nBound_ = 0;
    saveOptions_ = rhs.saveOptions_;
    solverOptions_ = rhs.solverOptions_;
    if (nDepth_>=0) {
      int n = (1<<nDepth_)+1+nDepth_;
      for (int i=0;i<n;i++) 
        delete nodeInfo_[i];
      delete [] nodeInfo_;
    }
    nDepth_ = rhs.nDepth_;
    nNodes_ = rhs.nNodes_;
    numberNodesExplored_ = rhs.numberNodesExplored_;
    numberIterations_ = rhs.numberIterations_;
    presolveType_ = rhs.presolveType_;
  }
  return *this;
}
// Zaps stuff 1 - arrays, 2 ints, 3 both
void 
ClpNodeStuff::zap(int type)
{
  if ((type&1)!=0) {
    downPseudo_ = NULL;
    upPseudo_ = NULL;
    numberDown_ = NULL;
    numberUp_ = NULL;
    numberDownInfeasible_ = NULL;
    numberUpInfeasible_ = NULL;
    saveCosts_ = NULL;
    nodeInfo_ = NULL;
    large_ = NULL;
    whichRow_ = NULL;
    whichColumn_ = NULL;
  }
  if ((type&2)!=0) {
    nBound_ = 0;
    saveOptions_ = 0;
    solverOptions_ = 0;
    nDepth_ = -1;
    nNodes_ = 0;
    presolveType_ = 0;
    numberNodesExplored_ = 0;
    numberIterations_ = 0;
  }
}

//-------------------------------------------------------------------
// Destructor 
//-------------------------------------------------------------------
ClpNodeStuff::~ClpNodeStuff ()
{
  delete [] downPseudo_;
  delete [] upPseudo_;
  delete [] numberDown_;
  delete [] numberUp_;
  delete [] numberDownInfeasible_;
  delete [] numberUpInfeasible_;
  if (nDepth_>=0) {
    int n = (1<<nDepth_)+1+nDepth_;
    for (int i=0;i<n;i++) 
      delete nodeInfo_[i];
    delete [] nodeInfo_;
  }
}
/* Fill with pseudocosts */
void 
ClpNodeStuff::fillPseudoCosts(const double * down, const double * up, 
                              const int * numberDown, const int * numberUp,
                              const int * numberDownInfeasible, 
                              const int * numberUpInfeasible,
                              int number)
{
  delete [] downPseudo_;
  delete [] upPseudo_;
  delete [] numberDown_;
  delete [] numberUp_;
  delete [] numberDownInfeasible_;
  delete [] numberUpInfeasible_;
  downPseudo_ = CoinCopyOfArray(down,number);
  upPseudo_ = CoinCopyOfArray(up,number);
  numberDown_ = CoinCopyOfArray(numberDown,number);
  numberUp_ = CoinCopyOfArray(numberUp,number);
  numberDownInfeasible_ = CoinCopyOfArray(numberDownInfeasible,number);
  numberUpInfeasible_ = CoinCopyOfArray(numberUpInfeasible,number);
  // scale
  for (int i=0;i<number;i++) {
    int n;
    n = numberDown_[i]+numberDownInfeasible_[i];
    if (n)
      downPseudo_[i] *= n;
    n = numberUp_[i]+numberUpInfeasible_[i];
    if (n)
      upPseudo_[i] *= n;
  }
}
// Update pseudo costs
void 
ClpNodeStuff::update(int way,int sequence,double change,bool feasible)
{
  if (way<0) {
    if (feasible)
      numberDown_[sequence]++;
    else
      numberDownInfeasible_[sequence]++;
    downPseudo_[sequence] += CoinMax(change,1.0e-12);
  } else {
    if (feasible)
      numberUp_[sequence]++;
    else
      numberUpInfeasible_[sequence]++;
    upPseudo_[sequence] += CoinMax(change,1.0e-12);
  }
}
//#############################################################################
// Constructors / Destructor / Assignment
//#############################################################################

//-------------------------------------------------------------------
// Default Constructor 
//-------------------------------------------------------------------
ClpHashValue::ClpHashValue () :
  hash_(NULL),
  numberHash_(0),
  maxHash_(0),
  lastUsed_(-1)
{
}
//-------------------------------------------------------------------
// Useful Constructor from model
//-------------------------------------------------------------------
ClpHashValue::ClpHashValue (ClpSimplex * model) :
  hash_(NULL),
  numberHash_(0),
  maxHash_(0),
  lastUsed_(-1)
{
  maxHash_ = 1000;
  int numberColumns = model->numberColumns();
  const double * columnLower = model->columnLower();
  const double * columnUpper = model->columnUpper();
  int numberRows = model->numberRows();
  const double * rowLower = model->rowLower();
  const double * rowUpper = model->rowUpper();
  const double * objective = model->objective();
  CoinPackedMatrix * matrix = model->matrix();
  const int * columnLength = matrix->getVectorLengths();
  const CoinBigIndex * columnStart = matrix->getVectorStarts();
  const double * elementByColumn = matrix->getElements();
  int i;
  int ipos;

  hash_ = new CoinHashLink[maxHash_];
  
  for ( i = 0; i < maxHash_; i++ ) {
    hash_[i].value = -1.0e-100;
    hash_[i].index = -1;
    hash_[i].next = -1;
  }
  // Put in +0
  hash_[0].value=0.0;
  hash_[0].index=0;
  numberHash_=1;
  /*
   * Initialize the hash table.  Only the index of the first value that
   * hashes to a value is entered in the table; subsequent values that
   * collide with it are not entered.
   */
  for ( i = 0; i < numberColumns; i++ ) {
    int length=columnLength[i];
    CoinBigIndex start = columnStart[i];
    for (CoinBigIndex i=start;i<start+length;i++) {
      double value = elementByColumn[i];
      ipos = hash ( value);
      if ( hash_[ipos].index == -1 ) {
        hash_[ipos].index = numberHash_;
        numberHash_++;
        hash_[ipos].value = elementByColumn[i];
      }
    }
  }
  
  /*
   * Now take care of the values that collided in the preceding loop,
   * Also do other stuff
   */
  for ( i = 0; i < numberRows; i++ ) {
    if (numberHash_*2>maxHash_)
      resize(true);
    double value;
    value = rowLower[i];
    ipos = index(value);
    if (ipos<0) 
      addValue(value);
    value = rowUpper[i];
    ipos = index(value);
    if (ipos<0) 
      addValue(value);
  }
  for ( i = 0; i < numberColumns; i++ ) {
    int length=columnLength[i];
    CoinBigIndex start = columnStart[i];
    if (numberHash_*2>maxHash_)
      resize(true);
    double value;
    value = objective[i];
    ipos = index(value);
    if (ipos<0) 
      addValue(value);
    value = columnLower[i];
    ipos = index(value);
    if (ipos<0) 
      addValue(value);
    value = columnUpper[i];
    ipos = index(value);
    if (ipos<0) 
      addValue(value);
    for (CoinBigIndex j=start;j<start+length;j++) {
      if (numberHash_*2>maxHash_)
        resize(true);
      value = elementByColumn[j];
      ipos = index(value);
      if (ipos<0) 
        addValue(value);
    }
  }
  resize(false);
}

//-------------------------------------------------------------------
// Copy constructor 
//-------------------------------------------------------------------
ClpHashValue::ClpHashValue (const ClpHashValue & rhs) :
  hash_(NULL),
  numberHash_(rhs.numberHash_),
  maxHash_(rhs.maxHash_),
  lastUsed_(rhs.lastUsed_)
{  
  if (maxHash_) {
    CoinHashLink * newHash = new CoinHashLink[maxHash_];
    int i;
    for ( i = 0; i < maxHash_; i++ ) {
      newHash[i].value = rhs.hash_[i].value;
      newHash[i].index = rhs.hash_[i].index;
      newHash[i].next = rhs.hash_[i].next;
    }
  }
}

//-------------------------------------------------------------------
// Destructor 
//-------------------------------------------------------------------
ClpHashValue::~ClpHashValue ()
{
  delete [] hash_;
}

//----------------------------------------------------------------
// Assignment operator 
//-------------------------------------------------------------------
ClpHashValue &
ClpHashValue::operator=(const ClpHashValue& rhs)
{
  if (this != &rhs) {
    numberHash_ = rhs.numberHash_;
    maxHash_ = rhs.maxHash_;
    lastUsed_ = rhs.lastUsed_;
    delete [] hash_;
    if (maxHash_) {
      CoinHashLink * newHash = new CoinHashLink[maxHash_];
      int i;
      for ( i = 0; i < maxHash_; i++ ) {
        newHash[i].value = rhs.hash_[i].value;
        newHash[i].index = rhs.hash_[i].index;
        newHash[i].next = rhs.hash_[i].next;
      }
    } else {
      hash_ = NULL;
    }
  }
  return *this;
}
// Return index or -1 if not found
int 
ClpHashValue::index(double value) const
{
  if (!value)
    return 0;
  int ipos = hash ( value);
  int returnCode=-1;
  while ( hash_[ipos].index>=0 ) {
    if (value==hash_[ipos].value) {
      returnCode = hash_[ipos].index;
      break;
    } else { 
      int k = hash_[ipos].next;
      if ( k == -1 ) {
        break;
      } else {
        ipos = k;
      }
    }
  }
  return returnCode;
}
// Add value to list and return index
int 
ClpHashValue::addValue(double value) 
{
  int ipos = hash ( value);
  
  assert (value!=hash_[ipos].value);
  if (hash_[ipos].index==-1) {
    // can put in here
    hash_[ipos].index = numberHash_;
    numberHash_++;
    hash_[ipos].value = value;
    return numberHash_-1;
  }
  int k = hash_[ipos].next;
  while (k!=-1) {
    ipos = k;
    k = hash_[ipos].next;
  }
  while ( true ) {
    ++lastUsed_;
    assert (lastUsed_<=maxHash_);
    if ( hash_[lastUsed_].index == -1 ) {
      break;
    }
  }
  hash_[ipos].next = lastUsed_;
  hash_[lastUsed_].index = numberHash_;
  numberHash_++;
  hash_[lastUsed_].value = value;
  return numberHash_-1;
}

int
ClpHashValue::hash ( double value) const
{
  static int mmult[] = {
    262139, 259459, 256889, 254291, 251701, 249133, 246709, 244247,
    241667, 239179, 236609, 233983, 231289, 228859, 226357, 223829,
    221281, 218849, 216319, 213721, 211093, 208673, 206263, 203773,
    201233, 198637, 196159, 193603, 191161, 188701, 186149, 183761,
    181303, 178873, 176389, 173897, 171469, 169049, 166471, 163871,
    161387, 158941, 156437, 153949, 151531, 149159, 146749, 144299,
    141709, 139369, 136889, 134591, 132169, 129641, 127343, 124853,
    122477, 120163, 117757, 115361, 112979, 110567, 108179, 105727,
    103387, 101021, 98639, 96179, 93911, 91583, 89317, 86939, 84521,
    82183, 79939, 77587, 75307, 72959, 70793, 68447, 66103
  };
  union { double d; char c[8]; } v1;
  assert (sizeof(double)==8);
  v1.d = value;
  int n = 0;
  int j;

  for ( j = 0; j < 8; ++j ) {
    int ichar = v1.c[j];
    n += mmult[j] * ichar;
  }
  return ( abs ( n ) % maxHash_ );      /* integer abs */
}
void
ClpHashValue::resize(bool increaseMax)
{
  int newSize = increaseMax ? ((3*maxHash_)>>1)+1000 : maxHash_;
  CoinHashLink * newHash = new CoinHashLink[newSize];
  int i;
  for ( i = 0; i < newSize; i++ ) {
    newHash[i].value = -1.0e-100;
    newHash[i].index = -1;
    newHash[i].next = -1;
  }
  // swap
  CoinHashLink * oldHash = hash_;
  hash_ = newHash;
  int oldSize = maxHash_;
  maxHash_ = newSize;
  /*
   * Initialize the hash table.  Only the index of the first value that
   * hashes to a value is entered in the table; subsequent values that
   * collide with it are not entered.
   */
  int ipos;
  int n=0;
  for ( i = 0; i < oldSize; i++ ) {
    if (oldHash[i].index>=0) {
      ipos = hash ( oldHash[i].value);
      if ( hash_[ipos].index == -1 ) {
        hash_[ipos].index = n;
        n++;
        hash_[ipos].value = oldHash[i].value;
        // unmark
        oldHash[i].index=-1;
      }
    }
  }
  /*
   * Now take care of the values that collided in the preceding loop,
   * by finding some other entry in the table for them.
   * Since there are as many entries in the table as there are values,
   * there must be room for them.
   */
  lastUsed_ = -1;
  for ( i = 0; i < oldSize; ++i ) {
    if (oldHash[i].index>=0) {
      double value = oldHash[i].value;
      ipos = hash ( value);
      int k;
      while ( true ) {
        assert (value!=hash_[ipos].value);
        k = hash_[ipos].next;
        if ( k == -1 ) {
          while ( true ) {
            ++lastUsed_;
            assert (lastUsed_<=maxHash_);
            if ( hash_[lastUsed_].index == -1 ) {
              break;
            }
          }
          hash_[ipos].next = lastUsed_;
          hash_[lastUsed_].index = n;
          n++;
          hash_[lastUsed_].value = value;
          break;
        } else {
          ipos = k;
        }
      }
    }
  }
  assert (n==numberHash_);
  delete [] oldHash;
}