Changeset 150


Ignore:
Timestamp:
Jun 26, 2005 5:56:08 AM (15 years ago)
Author:
forrest
Message:

add statistics

Location:
trunk
Files:
2 added
3 edited

Legend:

Unmodified
Added
Removed
  • trunk/CbcModel.cpp

    r148 r150  
    2121#include "CbcHeuristic.hpp"
    2222#include "CbcModel.hpp"
     23#include "CbcStatistics.hpp"
    2324#include "CbcStrategy.hpp"
    2425#include "CbcMessage.hpp"
     
    376377*/
    377378
    378 void CbcModel::branchAndBound()
     379void CbcModel::branchAndBound(int doStatistics)
    379380
    380381{
     
    565566  usedInSolution_ = new int[numberColumns];
    566567  CoinZeroN(usedInSolution_,numberColumns);
     568/*
     569  For printing totals and for CbcNode (numberNodes_)
     570*/
     571  numberIterations_ = 0 ;
     572  numberNodes_ = 0 ;
     573  // For printing in case node at n000 cutoff
     574  int lastPrinted=0;
     575  int lastRedoTree=0;
     576  int maximumStatistics=0;
     577  CbcStatistics ** statistics = NULL;
     578  // Do on switch
     579  if (doStatistics) {
     580    maximumStatistics=10000;
     581    statistics = new CbcStatistics * [maximumStatistics];
     582    memset(statistics,0,maximumStatistics*sizeof(CbcStatistics *));
     583  }
     584
    567585  { int iObject ;
    568586    int preferredWay ;
     
    575593          object_[iObject]->infeasibility(preferredWay) ;
    576594      if (infeasibility ) numberUnsatisfied++ ; }
    577     if (numberUnsatisfied)
    578     { feasible = solveWithCuts(cuts,maximumCutPassesAtRoot_,
     595    if (numberUnsatisfied) {
     596      feasible = solveWithCuts(cuts,maximumCutPassesAtRoot_,
    579597                               NULL,numberOldActiveCuts,numberNewCuts,
    580                                maximumWhich, whichGenerator) ; } }
     598                               maximumWhich, whichGenerator) ;
     599    }
     600  }
    581601  // make cut generators less aggressive
    582602  for (iCutGenerator = 0;iCutGenerator<numberCutGenerators_;iCutGenerator++) {
     
    602622  bool resolved = false ;
    603623  CbcNode *newNode = NULL ;
    604 /*
    605   For printing totals and for CbcNode (numberNodes_)
    606 */
    607   numberIterations_ = 0 ;
    608   numberNodes_ = 0 ;
    609 
    610624  if (feasible)
    611625  { newNode = new CbcNode ;
     
    726740      newNode->initializeInfo() ;
    727741      tree_->push(newNode) ;
     742      if (statistics) {
     743        if (numberNodes_==maximumStatistics) {
     744          maximumStatistics = 2*maximumStatistics;
     745          CbcStatistics ** temp = new CbcStatistics * [maximumStatistics];
     746          memset(temp,0,maximumStatistics*sizeof(CbcStatistics *));
     747          memcpy(temp,statistics,numberNodes_*sizeof(CbcStatistics *));
     748          delete [] statistics;
     749          statistics=temp;
     750        }
     751        assert (!statistics[numberNodes_]);
     752        statistics[numberNodes_]=new CbcStatistics(newNode);
     753      }
     754      numberNodes_++;
    728755#     ifdef CHECK_NODE
    729756      printf("Node %x on tree\n",newNode) ;
     
    785812    + print a summary line to let the user know we're working
    786813*/
    787     if ((numberNodes_%1000) == 0) {
    788       bool redoTree=nodeCompare_->every1000Nodes(this, numberNodes_) ;
     814    if (numberNodes_==lastRedoTree+1000||numberNodes_==lastRedoTree+1001) {
     815      lastRedoTree=numberNodes_&(~1);
     816      bool redoTree=nodeCompare_->every1000Nodes(this, lastRedoTree) ;
    789817      // redo tree if wanted
    790818      if (redoTree)
    791819        tree_->setComparison(*nodeCompare_) ;
    792820    }
    793     if ((numberNodes_%printFrequency_) == 0) {
     821    if (numberNodes_==lastPrinted+printFrequency_||
     822         numberNodes_==lastPrinted+printFrequency_+1) {
     823      lastPrinted=numberNodes_&(~1);
    794824      int j ;
    795825      int nNodes = tree_->size() ;
     
    801831      }
    802832      messageHandler()->message(CBC_STATUS,messages())
    803         << numberNodes_<< nNodes<< bestObjective_<< bestPossibleObjective_
     833        << lastPrinted<< nNodes<< bestObjective_<< bestPossibleObjective_
    804834        << CoinMessageEol ;
    805835      if (eventHandler) {
     
    886916      bool deleteNode ;
    887917      if (node->branch())
    888       { tree_->push(node) ;
     918      {
     919        // set nodenumber correctly
     920        node->nodeInfo()->setNodeNumber(numberNodes_);
     921        tree_->push(node) ;
     922        if (statistics) {
     923          if (numberNodes_==maximumStatistics) {
     924            maximumStatistics = 2*maximumStatistics;
     925            CbcStatistics ** temp = new CbcStatistics * [maximumStatistics];
     926            memset(temp,0,maximumStatistics*sizeof(CbcStatistics *));
     927            memcpy(temp,statistics,numberNodes_*sizeof(CbcStatistics *));
     928            delete [] statistics;
     929            statistics=temp;
     930          }
     931          assert (!statistics[numberNodes_]);
     932          statistics[numberNodes_]=new CbcStatistics(node);
     933        }
     934        numberNodes_++;
    889935        nodeOnTree=true; // back on tree
    890936        deleteNode = false ;
     
    922968      cuts = OsiCuts() ;
    923969      currentNumberCuts = solver_->getNumRows()-numberRowsAtContinuous_ ;
     970      int saveNumber = numberIterations_;
    924971      feasible = solveWithCuts(cuts,maximumCutPasses_,node,
    925972                               numberOldActiveCuts,numberNewCuts,
    926973                               maximumWhich,whichGenerator) ;
     974      if (statistics) {
     975        assert (numberNodes_);
     976        assert (statistics[numberNodes_-1]);
     977        assert (statistics[numberNodes_-1]->node()==numberNodes_-1);
     978        statistics[numberNodes_-1]->endOfBranch(numberIterations_-saveNumber,
     979                                               feasible ? solver_->getObjValue()
     980                                               : COIN_DBL_MAX);
     981    }
    927982/*
    928983  Check for abort on limits: node count, solution count, time, integrality gap.
     
    10621117        }
    10631118        assert (!newNode || newNode->objectiveValue() <= getCutoff()) ;
     1119        if (statistics) {
     1120          assert (numberNodes_);
     1121          assert (statistics[numberNodes_-1]);
     1122          assert (statistics[numberNodes_-1]->node()==numberNodes_-1);
     1123          if (newNode)
     1124            statistics[numberNodes_-1]->updateInfeasibility(newNode->numberUnsatisfied());
     1125          else
     1126            statistics[numberNodes_-1]->sayInfeasible();
     1127        }
    10641128        if (newNode)
    10651129        { if (newNode->variable() >= 0)
     
    11031167            newNode->setGuessedObjectiveValue(estValue) ;
    11041168            tree_->push(newNode) ;
     1169            if (statistics) {
     1170              if (numberNodes_==maximumStatistics) {
     1171                maximumStatistics = 2*maximumStatistics;
     1172                CbcStatistics ** temp = new CbcStatistics * [maximumStatistics];
     1173                memset(temp,0,maximumStatistics*sizeof(CbcStatistics *));
     1174                memcpy(temp,statistics,numberNodes_*sizeof(CbcStatistics *));
     1175                delete [] statistics;
     1176                statistics=temp;
     1177              }
     1178              assert (!statistics[numberNodes_]);
     1179              statistics[numberNodes_]=new CbcStatistics(newNode);
     1180            }
     1181            numberNodes_++;
    11051182#           ifdef CHECK_NODE
    11061183            printf("Node %x pushed on tree c\n",newNode) ;
     
    12041281      << CoinMessageEol ;
    12051282  }
     1283  if (statistics) {
     1284    // report in some way
     1285    int * lookup = new int[numberObjects_];
     1286    int i;
     1287    for (i=0;i<numberObjects_;i++)
     1288      lookup[i]=-1;
     1289    bool goodIds=true;
     1290    for (i=0;i<numberObjects_;i++) {
     1291      int id = object_[i]->id();
     1292      int iColumn = object_[i]->columnNumber();
     1293      if (iColumn<0)
     1294        iColumn = id+numberColumns;
     1295      if(id>=0&&id<numberObjects_) {
     1296        if (lookup[id]==-1) {
     1297          lookup[id]=iColumn;
     1298        } else {
     1299          goodIds=false;
     1300          break;
     1301        }
     1302      } else {
     1303        goodIds=false;
     1304        break;
     1305      }
     1306    }
     1307    if (!goodIds) {
     1308      delete [] lookup;
     1309      lookup=NULL;
     1310    }
     1311    if (doStatistics==3) {
     1312      printf("  node parent depth column   value                    obj      inf\n");
     1313      for ( i=0;i<numberNodes_;i++) {
     1314        statistics[i]->print(lookup);
     1315      }
     1316    }
     1317    if (doStatistics>1) {
     1318      // Find last solution
     1319      int k;
     1320      for (k=numberNodes_-1;k>=0;k--) {
     1321        if (statistics[k]->endingObjective()!=COIN_DBL_MAX&&
     1322            !statistics[k]->endingInfeasibility())
     1323          break;
     1324      }
     1325      if (k>=0) {
     1326        int depth=statistics[k]->depth();
     1327        int * which = new int[depth+1];
     1328        for (i=depth;i>=0;i--) {
     1329          which[i]=k;
     1330          k=statistics[k]->parentNode();
     1331        }
     1332        printf("  node parent depth column   value                    obj      inf\n");
     1333        for (i=0;i<=depth;i++) {
     1334          statistics[which[i]]->print(lookup);
     1335        }
     1336        delete [] which;
     1337      }
     1338    }
     1339    // now summary
     1340    int maxDepth=0;
     1341    double averageSolutionDepth=0.0;
     1342    int numberSolutions=0;
     1343    double averageCutoffDepth=0.0;
     1344    double averageSolvedDepth=0.0;
     1345    int numberCutoff=0;
     1346    int numberDown=0;
     1347    int numberFirstDown=0;
     1348    double averageInfDown=0.0;
     1349    double averageObjDown=0.0;
     1350    int numberCutoffDown=0;
     1351    int numberUp=0;
     1352    int numberFirstUp=0;
     1353    double averageInfUp=0.0;
     1354    double averageObjUp=0.0;
     1355    int numberCutoffUp=0;
     1356    double averageNumberIterations1=0.0;
     1357    double averageValue=0.0;
     1358    for ( i=0;i<numberNodes_;i++) {
     1359      int depth =  statistics[i]->depth();
     1360      int way =  statistics[i]->way();
     1361      double value = statistics[i]->value();
     1362      double startingObjective =  statistics[i]->startingObjective();
     1363      int startingInfeasibility = statistics[i]->startingInfeasibility();
     1364      double endingObjective = statistics[i]->endingObjective();
     1365      int endingInfeasibility = statistics[i]->endingInfeasibility();
     1366      maxDepth = CoinMax(depth,maxDepth);
     1367      // Only for completed
     1368      averageNumberIterations1 += statistics[i]->numberIterations();
     1369      averageValue += value;
     1370      if (endingObjective!=COIN_DBL_MAX&&!endingInfeasibility) {
     1371        numberSolutions++;
     1372        averageSolutionDepth += depth;
     1373      }
     1374      if (endingObjective==COIN_DBL_MAX) {
     1375        numberCutoff++;
     1376        averageCutoffDepth += depth;
     1377        if (way<0) {
     1378          numberDown++;
     1379          numberCutoffDown++;
     1380          if (way==-1)
     1381            numberFirstDown++;
     1382        } else {
     1383          numberUp++;
     1384          numberCutoffUp++;
     1385          if (way==1)
     1386            numberFirstUp++;
     1387        }
     1388      } else {
     1389        averageSolvedDepth += depth;
     1390        if (way<0) {
     1391          numberDown++;
     1392          averageInfDown += startingInfeasibility-endingInfeasibility;
     1393          averageObjDown += endingObjective-startingObjective;
     1394          if (way==-1)
     1395            numberFirstDown++;
     1396        } else {
     1397          numberUp++;
     1398          averageInfUp += startingInfeasibility-endingInfeasibility;
     1399          averageObjUp += endingObjective-startingObjective;
     1400          if (way==1)
     1401            numberFirstUp++;
     1402        }
     1403      }
     1404    }
     1405    // Now print
     1406    if (numberSolutions)
     1407      averageSolutionDepth /= (double) numberSolutions;
     1408    int numberSolved = numberNodes_-numberCutoff;
     1409    double averageNumberIterations2=numberIterations_-averageNumberIterations1;
     1410    if(numberCutoff) {
     1411      averageCutoffDepth /= (double) numberCutoff;
     1412      averageNumberIterations2 /= (double) numberCutoff;
     1413    }
     1414    if (numberNodes_)
     1415      averageValue /= (double) numberNodes_;
     1416    if (numberSolved) {
     1417      averageNumberIterations1 /= (double) numberSolved;
     1418      averageSolvedDepth /= (double) numberSolved;
     1419    }
     1420    printf("%d solution(s) were found (by branching) at an average depth of %g\n",
     1421           numberSolutions,averageSolutionDepth);
     1422    printf("average value of variable being branched on was %g\n",
     1423           averageValue);
     1424    printf("%d nodes were cutoff at an average depth of %g with iteration count of %g\n",
     1425           numberCutoff,averageCutoffDepth,averageNumberIterations2);
     1426    printf("%d nodes were solved at an average depth of %g with iteration count of %g\n",
     1427           numberSolved,averageSolvedDepth,averageNumberIterations1);
     1428    if (numberDown) {
     1429      averageInfDown /= (double) numberDown;
     1430      averageObjDown /= (double) numberDown;
     1431    }
     1432    printf("Down %d nodes (%d first, %d second) - %d cutoff, rest decrease numinf %g increase obj %g\n",
     1433           numberDown,numberFirstDown,numberDown-numberFirstDown,numberCutoffDown,
     1434           averageInfDown,averageObjDown);
     1435    if (numberUp) {
     1436      averageInfUp /= (double) numberUp;
     1437      averageObjUp /= (double) numberUp;
     1438    }
     1439    printf("Up %d nodes (%d first, %d second) - %d cutoff, rest decrease numinf %g increase obj %g\n",
     1440           numberUp,numberFirstUp,numberUp-numberFirstUp,numberCutoffUp,
     1441           averageInfUp,averageObjUp);
     1442    for ( i=0;i<numberNodes_;i++)
     1443      delete statistics[i];
     1444    delete [] statistics;
     1445    delete [] lookup;
     1446  }
    12061447/*
    12071448  If we think we have a solution, restore and confirm it with a call to
     
    12481489*/
    12491490  globalCuts_= OsiCuts() ;
    1250  
    12511491  return ; }
    12521492
     
    23282568      delete addCuts[i];
    23292569    delete [] addCuts;
    2330     numberNodes_++;
     2570    //numberNodes_++;
    23312571    return 0;
    23322572  }
     
    31863426        OsiClpSolverInterface * clpSolver
    31873427          = dynamic_cast<OsiClpSolverInterface *> (solver_);
    3188         if (clpSolver&&0) {
     3428        if (clpSolver) {
    31893429          // Maybe solver might like to know only column bounds will change
    31903430          int options = clpSolver->specialOptions();
  • trunk/CbcNode.cpp

    r146 r150  
    779779  we should see at least one bound change as a consequence of processing this
    780780  subproblem. Different types of branching objects could break this assertion.
    781   Yes - cuts will break
     781  Not true at all - we have not applied current branch - JJF.
    782782*/
    783783    double *boundChanges = new double [2*numberColumns] ;
     
    805805    printf("%d changed bounds\n",numberChangedBounds) ;
    806806#endif
    807     if (lastNode->branchingObject()->boundBranch())
    808       assert (numberChangedBounds);
     807    //if (lastNode->branchingObject()->boundBranch())
     808    //assert (numberChangedBounds);
    809809/*
    810810  Hand the lot over to the CbcPartialNodeInfo constructor, then clean up and
     
    19281928          } else {
    19291929            int iColumn = dynamicObject->columnNumber();
    1930             sort[numberToDo]=saveSolution[iColumn]-floor(saveSolution[iColumn]);
    1931             if (fabs(sort[numberToDo])>bestNot) {
     1930            double part =saveSolution[iColumn]-floor(saveSolution[iColumn]);
     1931            sort[numberToDo]=part;
     1932            if (1.0-fabs(part-0.5)>bestNot) {
    19321933              iBestNot=numberToDo;
    1933               bestNot = fabs(sort[numberToDo]);
     1934              bestNot = 1.0-fabs(part-0.5);
    19341935            }
    19351936            objectMark[i]=numberThisDown+numberThisUp;
     
    20152016      branch_=object->createBranch(preferredWay);
    20162017      branch_->way(preferredWay);
     2018      delete ws;
     2019      ws=NULL;
    20172020      break;
    20182021    } else {
     
    20382041          sort[iObj] += add*average;
    20392042        else
    2040           sort[iObj] = -CoinMax(sort[iObj]*averageDown,(1.0-sort[iObj])*averageUp);
     2043          sort[iObj] = -(sort[iObj]*averageDown+(1.0-sort[iObj])*averageUp);
    20412044      }
    20422045      // Sort
     
    20782081        int put=0;
    20792082        const double * downCost=osiclp->upRange();
    2080         const double * upCost=osiclp->downRange();;
     2083        const double * upCost=osiclp->downRange();
     2084        // Bug - so switch off for now
     2085        double distanceToCutoff=COIN_DBL_MAX;
    20812086        for ( i=0;i<numberToDo;i++) {
    20822087          int iObject = whichObject[i];
     
    21092114                //printf("%d infeas both penalty %g %g\n",iObject,upPenalty,downPenalty);
    21102115                anyAction=-2;
     2116                delete ws;
     2117                ws=NULL;
    21112118                break;
    21122119              } else if (iAction==2) {
     
    21412148            // just estimate
    21422149            double add = objectMark[iObject];
    2143             sort[iObj] += add*average;
    21442150            double trueEstimate = sort[i] - add*average;
    21452151            sort[put]=trueEstimate;
     
    23772383              bestChoice = choice.objectNumber;
    23782384              whichChoice = iDo;
    2379               if (numberStrong<=1)
     2385              if (numberStrong<=1) {
     2386                delete ws;
     2387                ws=NULL;
    23802388                break;
     2389              }
    23812390            } else {
    23822391              delete choice.possibleBranch;
    2383               if (iDo>=2*numberStrong)
     2392              if (iDo>=2*numberStrong) {
     2393                delete ws;
     2394                ws=NULL;
    23842395                break;
     2396              }
    23852397              if (!dynamicObject||dynamicObject->numberTimesUp()>1) {
    23862398                if (iDo-whichChoice>=numberStrong)
    23872399                  break; // give up
    23882400              } else {
    2389                 if (iDo-whichChoice>=2*numberStrong)
     2401                if (iDo-whichChoice>=2*numberStrong) {
     2402                  delete ws;
     2403                  ws=NULL;
    23902404                  break; // give up
     2405                }
    23912406              }
    23922407            }
     
    24002415              choice.possibleBranch->branch();
    24012416              delete choice.possibleBranch;
     2417              delete ws;
     2418              ws=NULL;
    24022419              break;
    24032420            } else {
     
    24162433              choice.possibleBranch->branch();
    24172434              delete choice.possibleBranch;
     2435              delete ws;
     2436              ws=NULL;
    24182437              break;
    24192438            } else {
     
    24272446            //printf("Both infeasible for choice %d sequence %d\n",i,
    24282447            // model->object(choice.objectNumber)->columnNumber());
     2448            delete ws;
     2449            ws=NULL;
    24292450            break;
    24302451          }
     
    24462467          }
    24472468          anyAction=0;
     2469          delete ws;
     2470          ws=NULL;
    24482471          break;
    24492472        }
     
    26602683  CbcNodeInfo * parent = nodeInfo_->parent();
    26612684  int parentNodeNumber = -1;
     2685  //CbcBranchingObject * object1 = branch_->object_;
     2686  //CbcObject * object = object1->
     2687  //int sequence = object->columnNumber);
     2688  int id=-1;
     2689  double value=0.0;
     2690  if (branch_) {
     2691    id = branch_->variable();
     2692    value = branch_->value();
     2693  }
     2694  printf("id %d value %g objvalue %g\n",id,value,objectiveValue_);
    26622695  if (parent)
    26632696    parentNodeNumber = parent->nodeNumber();
  • trunk/include/CbcModel.hpp

    r143 r150  
    154154      branch & cut search tree. The search ends when the tree is exhausted or
    155155      one of several execution limits is reached.
     156      If doStatistics is 1 summary statistics are printed
     157      if 2 then also the path to best solution (if found by branching)
     158      if 3 then also one line per node
    156159    */
    157      void branchAndBound();
     160     void branchAndBound(int doStatistics=0);
    158161
    159162    /** \brief create a clean model from partially fixed problem
Note: See TracChangeset for help on using the changeset viewer.