/* $Id: ClpNode.cpp 1910 2013-01-27 02:00:13Z stefan $ */ // 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 (lower[iColumn])); value = CoinMin(value, static_cast (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(numberUpInfeasible[iInteger]) / static_cast(nUp); upValue2 *= ratio; } int nDown = numberDown[iInteger]; double downValue2 = (downPseudo[iInteger] / (1.0 + nDown)); if (nDown) { double ratio = 1.0 + INFEAS_MULTIPLIER*static_cast(numberDownInfeasible[iInteger]) / static_cast(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(numberUpInfeasible[iInteger]) / static_cast(nUp); upValue *= ratio; } double downValue = (value - floor(value)) * (downPseudo[iInteger] / (1.0 + nDown)); if (nDown) { double ratio = 1.0 + INFEAS_MULTIPLIER*static_cast(numberDownInfeasible[iInteger]) / static_cast(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 (lower[iColumn]); upper_[iInteger] = static_cast (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 (lower[iColumn])) model->setColumnLower(iColumn, lower_[iInteger]); if (upper_[iInteger] != static_cast (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; }