/* $Id$ */
// Copyright (C) 2008, International Business Machines
// Corporation and others.  All Rights Reserved.
// This code is licensed under the terms of the Eclipse Public License (EPL).

#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)
  , depth_(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)
  , depth_(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_;
  depth_ = depth;
  // 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_;
    assert(!factorization_);
    factorization_ = new ClpFactorization(*model->factorization(), numberRows);
    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);
    ClpDualRowSteepest *pivot = dynamic_cast< ClpDualRowSteepest * >(model->dualRowPivot());
    if (pivot) {
      assert(!weights_);
      weights_ = new ClpDualRowSteepest(*pivot);
    }
  } 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?
        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);
          }
        }
        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?
        ClpDualRowSteepest *pivot = dynamic_cast< ClpDualRowSteepest * >(model->dualRowPivot());
        if (pivot) {
          assert(!weights_);
          weights_ = new ClpDualRowSteepest(*pivot);
        }
      } 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() && model->numberRows() == numberRows)
        CoinMemcpyN(model->pivotVariable(), numberRows, pivotVariables_);
      else
        CoinFillN(pivotVariables_, numberRows, -1);
    }
  }
  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 = 0.0;
  int bestPriority = COIN_INT_MAX;
  sumInfeasibilities_ = 0.0;
  numberInfeasibilities_ = 0;
  int nFix = 0;
  double gap = CoinMax(model->dualObjectiveLimit() - objectiveValue_, 1.0e-4);
#define PSEUDO 3
#if PSEUDO == 1 || PSEUDO == 2
  // Column copy of matrix
  ClpPackedMatrix *matrix = model->clpScaledMatrix();
  const double *objective = model->costRegion();
  if (!objective) {
    objective = model->objective();
    //if (!matrix)
    matrix = dynamic_cast< ClpPackedMatrix * >(model->clpMatrix());
  } else 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();
  double direction = model->optimizationDirection();
  const double *dual = dualSolution_ + numberColumns;
#if PSEUDO == 2
  double *activeWeight = new double[numberRows];
  const double *rowLower = model->rowLower();
  const double *rowUpper = model->rowUpper();
  const double *rowActivity = model->primalRowSolution();
  double tolerance = 1.0e-6;
  for (int iRow = 0; iRow < numberRows; iRow++) {
    // could use pi to see if active or activity
    if (rowActivity[iRow] > rowUpper[iRow] - tolerance
      || rowActivity[iRow] < rowLower[iRow] + tolerance) {
      activeWeight[iRow] = 0.0;
    } else {
      activeWeight[iRow] = -1.0;
    }
  }
  for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
    if (integerType[iColumn]) {
      double value = solution[iColumn];
      if (fabs(value - floor(value + 0.5)) > 1.0e-6) {
        CoinBigIndex start = columnStart[iColumn];
        CoinBigIndex end = start + columnLength[iColumn];
        for (CoinBigIndex j = start; j < end; j++) {
          int iRow = row[j];
          if (activeWeight[iRow] >= 0.0)
            activeWeight[iRow] += 1.0;
        }
      }
    }
  }
  for (int iRow = 0; iRow < numberRows; iRow++) {
    if (activeWeight[iRow] > 0.0) {
      // could use pi
      activeWeight[iRow] = 1.0 / activeWeight[iRow];
    } else {
      activeWeight[iRow] = 0.0;
    }
  }
#endif
#endif
  const double *downPseudo = stuff->downPseudo_;
  const int *numberDown = stuff->numberDown_;
  const int *numberDownInfeasible = stuff->numberDownInfeasible_;
  const double *upPseudo = stuff->upPseudo_;
  const int *priority = stuff->priority_;
  const int *numberUp = stuff->numberUp_;
  const int *numberUpInfeasible = stuff->numberUpInfeasible_;
  int numberBeforeTrust = stuff->numberBeforeTrust_;
  int stateOfSearch = stuff->stateOfSearch_;
  int iInteger = 0;
  // weight at 1.0 is max min (CbcBranch was 0.8,0.1) (ClpNode was 0.9,0.9)
#define WEIGHT_AFTER 0.9
#define WEIGHT_BEFORE 0.2
  //Stolen from Constraint Integer Programming book (with epsilon change)
#define WEIGHT_PRODUCT
#ifdef WEIGHT_PRODUCT
  double smallChange = stuff->smallChange_;
#endif
#ifndef INFEAS_MULTIPLIER
#define INFEAS_MULTIPLIER 1.0
#endif
  for (iColumn = 0; iColumn < numberColumns; iColumn++) {
    if (integerType[iColumn]) {
      double value = solution[iColumn];
      value = CoinMax(value, static_cast< double >(lower[iColumn]));
      value = CoinMin(value, static_cast< double >(upper[iColumn]));
      double nearest = floor(value + 0.5);
      if (fabs(value - nearest) > integerTolerance) {
        numberInfeasibilities_++;
        sumInfeasibilities_ += fabs(value - nearest);
#if PSEUDO == 1 || PSEUDO == 2
        double upValue = 0.0;
        double downValue = 0.0;
        double value2 = direction * objective[iColumn];
        //double dj2=value2;
        if (value2) {
          if (value2 > 0.0)
            upValue += 1.5 * value2;
          else
            downValue -= 1.5 * 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];
            //dj2 += value2;
#if PSEUDO == 2
            assert(activeWeight[iRow] > 0.0 || fabs(dual[iRow]) < 1.0e-6);
            value2 *= activeWeight[iRow];
#endif
            if (value2 > 0.0)
              upValue += value2;
            else
              downValue -= value2;
          }
        }
        //assert (fabs(dj2)<1.0e-4);
        int nUp = numberUp[iInteger];
        double upValue2 = (upPseudo[iInteger] / (1.0 + nUp));
        // Extra for infeasible branches
        if (nUp) {
          double ratio = 1.0 + INFEAS_MULTIPLIER * static_cast< double >(numberUpInfeasible[iInteger]) / static_cast< double >(nUp);
          upValue2 *= ratio;
        }
        int nDown = numberDown[iInteger];
        double downValue2 = (downPseudo[iInteger] / (1.0 + nDown));
        if (nDown) {
          double ratio = 1.0 + INFEAS_MULTIPLIER * static_cast< double >(numberDownInfeasible[iInteger]) / static_cast< double >(nDown);
          downValue2 *= ratio;
        }
        //printf("col %d - downPi %g up %g, downPs %g up %g\n",
        //     iColumn,upValue,downValue,upValue2,downValue2);
        upValue = CoinMax(0.1 * upValue, upValue2);
        downValue = CoinMax(0.1 * downValue, downValue2);
        //upValue = CoinMax(upValue,1.0e-8);
        //downValue = CoinMax(downValue,1.0e-8);
        upValue *= ceil(value) - value;
        downValue *= value - floor(value);
        double infeasibility;
        //if (depth>1000)
        //infeasibility = CoinMax(upValue,downValue)+integerTolerance;
        //else
        if (stateOfSearch <= 2) {
          // no solution
          infeasibility = (1.0 - WEIGHT_BEFORE) * CoinMax(upValue, downValue) + WEIGHT_BEFORE * CoinMin(upValue, downValue) + integerTolerance;
        } else {
#ifndef WEIGHT_PRODUCT
          infeasibility = (1.0 - WEIGHT_AFTER) * CoinMax(upValue, downValue) + WEIGHT_AFTER * CoinMin(upValue, downValue) + integerTolerance;
#else
          infeasibility = CoinMax(CoinMax(upValue, downValue), smallChange) * CoinMax(CoinMin(upValue, downValue), smallChange);
#endif
        }
        estimatedSolution_ += CoinMin(upValue2, downValue2);
#elif PSEUDO == 3
        int nUp = numberUp[iInteger];
        int nDown = numberDown[iInteger];
        // Extra 100% for infeasible branches
        double upValue = (ceil(value) - value) * (upPseudo[iInteger] / (1.0 + nUp));
        if (nUp) {
          double ratio = 1.0 + INFEAS_MULTIPLIER * static_cast< double >(numberUpInfeasible[iInteger]) / static_cast< double >(nUp);
          upValue *= ratio;
        }
        double downValue = (value - floor(value)) * (downPseudo[iInteger] / (1.0 + nDown));
        if (nDown) {
          double ratio = 1.0 + INFEAS_MULTIPLIER * static_cast< double >(numberDownInfeasible[iInteger]) / static_cast< double >(nDown);
          downValue *= ratio;
        }
        if (nUp < numberBeforeTrust || nDown < numberBeforeTrust) {
          upValue *= 10.0;
          downValue *= 10.0;
        }

        double infeasibility;
        //if (depth>1000)
        //infeasibility = CoinMax(upValue,downValue)+integerTolerance;
        //else
        if (stateOfSearch <= 2) {
          // no solution
          infeasibility = (1.0 - WEIGHT_BEFORE) * CoinMax(upValue, downValue) + WEIGHT_BEFORE * CoinMin(upValue, downValue) + integerTolerance;
        } else {
#ifndef WEIGHT_PRODUCT
          infeasibility = (1.0 - WEIGHT_AFTER) * CoinMax(upValue, downValue) + WEIGHT_AFTER * CoinMin(upValue, downValue) + integerTolerance;
#else
          infeasibility = CoinMax(CoinMax(upValue, downValue), smallChange) * CoinMax(CoinMin(upValue, downValue), smallChange);
          //infeasibility += CoinMin(upValue,downValue)*smallChange;
#endif
        }
        //infeasibility = 0.1*CoinMax(upValue,downValue)+
        //0.9*CoinMin(upValue,downValue) + integerTolerance;
        estimatedSolution_ += CoinMin(upValue, downValue);
#else
        double infeasibility = fabs(value - nearest);
#endif
        assert(infeasibility > 0.0);
        if (priority[iInteger] < bestPriority) {
          mostAway = 0.0;
          bestPriority = priority[iInteger];
        } else if (priority[iInteger] > bestPriority) {
          infeasibility = 0.0;
        }
        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++;
    }
  }
  //printf("Choosing %d inf %g pri %d\n",
  // sequence_,mostAway,bestPriority);
#if PSEUDO == 2
  delete[] activeWeight;
#endif
  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] = static_cast< int >(lower[iColumn]);
        upper_[iInteger] = static_cast< 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 &)
{
  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] != static_cast< int >(lower[iColumn]))
          model->setColumnLower(iColumn, lower_[iInteger]);
        if (upper_[iInteger] != static_cast< int >(upper[iColumn]))
          model->setColumnUpper(iColumn, upper_[iInteger]);
      }
    }
  }
  if (doBoundsEtc && doBoundsEtc < 3) {
    //model->copyFactorization(*factorization_);
    model->copyFactorization(*factorization_);
    ClpDualRowSteepest *pivot = dynamic_cast< ClpDualRowSteepest * >(model->dualRowPivot());
    if (pivot && weights_) {
      pivot->fill(*weights_);
    }
    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_);
  }
}
// Choose a new variable
void ClpNode::chooseVariable(ClpSimplex *, ClpNodeStuff * /*info*/)
{
#if 0
     int way = branchState_.firstBranch;
     if (branchState_.branch > 0)
          way = 1 - way;
     assert (!branchState_.branch);
     // We need to use pseudo costs to choose a variable
     int numberColumns = model->numberColumns();
#endif
}
// Fix on reduced costs
int ClpNode::fixOnReducedCosts(ClpSimplex *)
{

  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)
  , smallChange_(1.0e-8)
  , downPseudo_(NULL)
  , upPseudo_(NULL)
  , priority_(NULL)
  , numberDown_(NULL)
  , numberUp_(NULL)
  , numberDownInfeasible_(NULL)
  , numberUpInfeasible_(NULL)
  , saveCosts_(NULL)
  , nodeInfo_(NULL)
  , large_(NULL)
  , whichRow_(NULL)
  , whichColumn_(NULL)
  ,
#ifndef NO_FATHOM_PRINT
  handler_(NULL)
  ,
#endif
  nBound_(0)
  , saveOptions_(0)
  , solverOptions_(0)
  , maximumNodes_(0)
  , numberBeforeTrust_(0)
  , stateOfSearch_(0)
  , nDepth_(-1)
  , nNodes_(0)
  , numberNodesExplored_(0)
  , numberIterations_(0)
  , presolveType_(0)
#ifndef NO_FATHOM_PRINT
  , startingDepth_(-1)
  , nodeCalled_(-1)
#endif
{
}

//-------------------------------------------------------------------
// Copy constructor
//-------------------------------------------------------------------
ClpNodeStuff::ClpNodeStuff(const ClpNodeStuff &rhs)
  : integerTolerance_(rhs.integerTolerance_)
  , integerIncrement_(rhs.integerIncrement_)
  , smallChange_(rhs.smallChange_)
  , downPseudo_(NULL)
  , upPseudo_(NULL)
  , priority_(NULL)
  , numberDown_(NULL)
  , numberUp_(NULL)
  , numberDownInfeasible_(NULL)
  , numberUpInfeasible_(NULL)
  , saveCosts_(NULL)
  , nodeInfo_(NULL)
  , large_(NULL)
  , whichRow_(NULL)
  , whichColumn_(NULL)
  ,
#ifndef NO_FATHOM_PRINT
  handler_(rhs.handler_)
  ,
#endif
  nBound_(0)
  , saveOptions_(rhs.saveOptions_)
  , solverOptions_(rhs.solverOptions_)
  , maximumNodes_(rhs.maximumNodes_)
  , numberBeforeTrust_(rhs.numberBeforeTrust_)
  , stateOfSearch_(rhs.stateOfSearch_)
  , nDepth_(rhs.nDepth_)
  , nNodes_(rhs.nNodes_)
  , numberNodesExplored_(rhs.numberNodesExplored_)
  , numberIterations_(rhs.numberIterations_)
  , presolveType_(rhs.presolveType_)
#ifndef NO_FATHOM_PRINT
  , startingDepth_(rhs.startingDepth_)
  , nodeCalled_(rhs.nodeCalled_)
#endif
{
}
//----------------------------------------------------------------
// Assignment operator
//-------------------------------------------------------------------
ClpNodeStuff &
ClpNodeStuff::operator=(const ClpNodeStuff &rhs)
{
  if (this != &rhs) {
    integerTolerance_ = rhs.integerTolerance_;
    integerIncrement_ = rhs.integerIncrement_;
    smallChange_ = rhs.smallChange_;
    downPseudo_ = NULL;
    upPseudo_ = NULL;
    priority_ = 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_;
    maximumNodes_ = rhs.maximumNodes_;
    numberBeforeTrust_ = rhs.numberBeforeTrust_;
    stateOfSearch_ = rhs.stateOfSearch_;
    int n = maximumNodes();
    if (n) {
      for (int i = 0; i < n; i++)
        delete nodeInfo_[i];
    }
    delete[] nodeInfo_;
    nodeInfo_ = NULL;
    nDepth_ = rhs.nDepth_;
    nNodes_ = rhs.nNodes_;
    numberNodesExplored_ = rhs.numberNodesExplored_;
    numberIterations_ = rhs.numberIterations_;
    presolveType_ = rhs.presolveType_;
#ifndef NO_FATHOM_PRINT
    handler_ = rhs.handler_;
    startingDepth_ = rhs.startingDepth_;
    nodeCalled_ = rhs.nodeCalled_;
#endif
  }
  return *this;
}
// Zaps stuff 1 - arrays, 2 ints, 3 both
void ClpNodeStuff::zap(int type)
{
  if ((type & 1) != 0) {
    downPseudo_ = NULL;
    upPseudo_ = NULL;
    priority_ = 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;
    maximumNodes_ = 0;
    numberBeforeTrust_ = 0;
    stateOfSearch_ = 0;
    nDepth_ = -1;
    nNodes_ = 0;
    presolveType_ = 0;
    numberNodesExplored_ = 0;
    numberIterations_ = 0;
  }
}

//-------------------------------------------------------------------
// Destructor
//-------------------------------------------------------------------
ClpNodeStuff::~ClpNodeStuff()
{
  delete[] downPseudo_;
  delete[] upPseudo_;
  delete[] priority_;
  delete[] numberDown_;
  delete[] numberUp_;
  delete[] numberDownInfeasible_;
  delete[] numberUpInfeasible_;
  int n = maximumNodes();
  if (n) {
    for (int i = 0; i < n; i++)
      delete nodeInfo_[i];
  }
  delete[] nodeInfo_;
#ifdef CLP_INVESTIGATE
  // Should be NULL - find out why not?
  assert(!saveCosts_);
#endif
  delete[] saveCosts_;
}
// Return maximum number of nodes
int ClpNodeStuff::maximumNodes() const
{
  int n = 0;
#if 0
     if (nDepth_ != -1) {
          if ((solverOptions_ & 32) == 0)
               n = (1 << nDepth_);
          else if (nDepth_)
               n = 1;
     }
     assert (n == maximumNodes_ - (1 + nDepth_) || n == 0);
#else
  if (nDepth_ != -1) {
    n = maximumNodes_ - (1 + nDepth_);
    assert(n > 0);
  }
#endif
  return n;
}
// Return maximum space for nodes
int ClpNodeStuff::maximumSpace() const
{
  return maximumNodes_;
}
/* Fill with pseudocosts */
void ClpNodeStuff::fillPseudoCosts(const double *down, const double *up,
  const int *priority,
  const int *numberDown, const int *numberUp,
  const int *numberDownInfeasible,
  const int *numberUpInfeasible,
  int number)
{
  delete[] downPseudo_;
  delete[] upPseudo_;
  delete[] priority_;
  delete[] numberDown_;
  delete[] numberUp_;
  delete[] numberDownInfeasible_;
  delete[] numberUpInfeasible_;
  downPseudo_ = CoinCopyOfArray(down, number);
  upPseudo_ = CoinCopyOfArray(up, number);
  priority_ = CoinCopyOfArray(priority, 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];
    if (n)
      downPseudo_[i] *= n;
    n = numberUp_[i];
    if (n)
      upPseudo_[i] *= n;
  }
}
// Update pseudo costs
void ClpNodeStuff::update(int way, int sequence, double change, bool feasible)
{
  assert(numberDown_[sequence] >= numberDownInfeasible_[sequence]);
  assert(numberUp_[sequence] >= numberUpInfeasible_[sequence]);
  if (way < 0) {
    numberDown_[sequence]++;
    if (!feasible)
      numberDownInfeasible_[sequence]++;
    downPseudo_[sequence] += CoinMax(change, 1.0e-12);
  } else {
    numberUp_[sequence]++;
    if (!feasible)
      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;
}

namespace {
/*
    Originally a local static variable in ClpHashValue::hash.
	Local static variables are a problem when building DLLs on Windows, but
	file-local constants seem to be ok.  -- lh, 101016 --
  */
const int mmult_for_hash[] = {
  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
};
}
int ClpHashValue::hash(double value) const
{

  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_for_hash[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;
}

/* vi: softtabstop=2 shiftwidth=2 expandtab tabstop=2
*/