Ignore:
Timestamp:
Jan 16, 2013 1:41:25 PM (7 years ago)
Author:
forrest
Message:

multiple root solvers, stronger strong branching and cutoff as constraint

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/Cbc/src/CbcModel.cpp

    r1832 r1839  
    6262#include "CbcFeasibilityBase.hpp"
    6363#include "CbcFathom.hpp"
     64#include "CbcFullNodeInfo.hpp"
    6465// include Probing
    6566#include "CglProbing.hpp"
     
    16931694                // try and redo debugger
    16941695                OsiRowCutDebugger * debugger = const_cast<OsiRowCutDebugger *> (solver_->getRowCutDebuggerAlways());
    1695                 if (debugger)
     1696                if (debugger) {
     1697                  if (numberColumns<=debugger->numberColumns())
    16961698                    debugger->redoSolution(numberColumns, originalColumns);
     1699                  else
     1700                    debugger=NULL; // no idea how to handle (SOS?)
     1701                }
    16971702                // User-provided solution might have been best. Synchronise.
    16981703                if (bestSolution_) {
     
    22182223    continuousSolver_ = solver_->clone() ;
    22192224
     2225    // add cutoff as constraint if wanted
     2226    if ((moreSpecialOptions_&16777216)!=0) {
     2227      if (!parentModel_) {
     2228        int numberColumns=solver_->getNumCols();
     2229        double * obj = CoinCopyOfArray(solver_->getObjCoefficients(),numberColumns);
     2230        int * indices = new int [numberColumns];
     2231        int n=0;
     2232        for (int i=0;i<numberColumns;i++) {
     2233          if (obj[i]) {
     2234            indices[n]=i;
     2235            obj[n++]=obj[i];
     2236          }
     2237        }
     2238        if (n) {
     2239          double cutoff=getCutoff();
     2240          double offset;
     2241          solver_->getDblParam(OsiObjOffset, offset);
     2242          solver_->addRow(n,indices,obj,-COIN_DBL_MAX,CoinMin(cutoff,1.0e25)+offset);
     2243        } else {
     2244          // no objective!
     2245          moreSpecialOptions_ &= ~16777216;
     2246        }
     2247        delete [] indices;
     2248        delete [] obj;
     2249      } else {
     2250        // switch off
     2251        moreSpecialOptions_ &= ~16777216;
     2252      }
     2253    }
    22202254    numberRowsAtContinuous_ = getNumRows() ;
    22212255    solver_->saveBaseModel();
     
    23462380        CglCutGenerator * generator = generator_[iCutGenerator]->generator();
    23472381        generator->setAggressiveness(generator->getAggressiveness() + 100);
     2382        if (!generator->canDoGlobalCuts())
     2383          generator->setGlobalCuts(false);
    23482384    }
    23492385    OsiCuts cuts ;
     
    24142450      will be removed from the heuristics list by doHeuristicsAtRoot.
    24152451    */
     2452    // See if multiple runs wanted
     2453    CbcModel ** rootModels=NULL;
     2454    if (!parentModel_&&multipleRootTries_%100) {
     2455      double rootTimeCpu=CoinCpuTime();
     2456      double startTimeRoot=CoinGetTimeOfDay();
     2457      int numberRootThreads=1;
     2458      /* undocumented fine tuning
     2459         aabbcc where cc is number of tries
     2460         bb if nonzero is number of threads
     2461         aa if nonzero just do heuristics
     2462      */
     2463      int numberModels = multipleRootTries_%100;
     2464#ifdef CBC_THREAD
     2465      numberRootThreads = (multipleRootTries_/100)%100;
     2466      if (!numberRootThreads)
     2467        numberRootThreads=numberModels;
     2468#endif
     2469      int otherOptions = (multipleRootTries_/10000)%100;
     2470      rootModels = new CbcModel * [numberModels];
     2471      unsigned int newSeed = randomSeed_;
     2472      if (newSeed==0) {
     2473        double time = fabs(CoinGetTimeOfDay());
     2474        while (time>=COIN_INT_MAX)
     2475          time *= 0.5;
     2476        newSeed = static_cast<unsigned int>(time);
     2477      } else if (newSeed<0) {
     2478        newSeed = 123456789;
     2479#ifdef COIN_HAS_CLP
     2480        OsiClpSolverInterface * clpSolver
     2481          = dynamic_cast<OsiClpSolverInterface *> (solver_);
     2482        if (clpSolver) {
     2483          newSeed += clpSolver->getModelPtr()->randomNumberGenerator()->getSeed();
     2484        }
     2485#endif
     2486      }
     2487      CoinWarmStartBasis * basis = dynamic_cast<CoinWarmStartBasis *> (solver_->getEmptyWarmStart());
     2488      for (int i=0;i<numberModels;i++) {
     2489        rootModels[i]=new CbcModel(*this);
     2490        rootModels[i]->setNumberThreads(0);
     2491        rootModels[i]->setMaximumNodes(otherOptions ? -1 : 0);
     2492        rootModels[i]->setRandomSeed(newSeed+10000000*i);
     2493        rootModels[i]->randomNumberGenerator()->setSeed(newSeed+50000000*i);
     2494        rootModels[i]->setMultipleRootTries(0);
     2495        // use seed
     2496        rootModels[i]->setSpecialOptions(specialOptions_ |(4194304|8388608));
     2497        rootModels[i]->solver_->setWarmStart(basis);
     2498#ifdef COIN_HAS_CLP
     2499        OsiClpSolverInterface * clpSolver
     2500          = dynamic_cast<OsiClpSolverInterface *> (rootModels[i]->solver_);
     2501        if (clpSolver) {
     2502          clpSolver->getModelPtr()->setRandomSeed(newSeed+20000000*i);
     2503          clpSolver->getModelPtr()->allSlackBasis();
     2504        }
     2505#endif
     2506      }
     2507      delete basis;
     2508      void * doRootCbcThread(void * voidInfo);
     2509#ifdef CBC_THREAD
     2510      if (numberRootThreads==1) {
     2511#endif
     2512        for (int iModel=0;iModel<numberModels;iModel++) {
     2513          doRootCbcThread(rootModels[iModel]);
     2514          // see if solved at root node
     2515          if (rootModels[iModel]->getMaximumNodes()) {
     2516            feasible=false;
     2517            break;
     2518          }
     2519        }
     2520#ifdef CBC_THREAD
     2521      } else {
     2522        Coin_pthread_t * threadId = new Coin_pthread_t [numberRootThreads];
     2523        for (int kModel=0;kModel<numberModels;kModel+=numberRootThreads) {
     2524          bool finished=false;
     2525          for (int iModel=kModel;iModel<CoinMin(numberModels,kModel+numberRootThreads);iModel++) {
     2526            pthread_create(&(threadId[iModel-kModel].thr), NULL,
     2527                           doRootCbcThread,
     2528                           rootModels[iModel]);
     2529          }
     2530          // wait
     2531          for (int iModel=kModel;iModel<CoinMin(numberModels,kModel+numberRootThreads);iModel++) {
     2532            pthread_join(threadId[iModel-kModel].thr, NULL);
     2533          }
     2534          // see if solved at root node
     2535          for (int iModel=kModel;iModel<CoinMin(numberModels,kModel+numberRootThreads);iModel++) {
     2536            if (rootModels[iModel]->getMaximumNodes())
     2537              finished=true;
     2538          }
     2539          if (finished) {
     2540            feasible=false;
     2541            break;
     2542          }
     2543        }
     2544        delete [] threadId;
     2545      }
     2546#endif
     2547      // sort solutions
     2548      int * which = new int [numberModels];
     2549      double * value = new double [numberModels];
     2550      int numberSolutions=0;
     2551      for (int iModel=0;iModel<numberModels;iModel++) {
     2552        if (rootModels[iModel]->bestSolution()) {
     2553          which[numberSolutions]=iModel;
     2554          value[numberSolutions++]=
     2555            -rootModels[iModel]->getMinimizationObjValue();
     2556        }
     2557      }
     2558      char general[100];
     2559      rootTimeCpu=CoinCpuTime()-rootTimeCpu;
     2560      if (numberRootThreads==1)
     2561        sprintf(general,"Multiple root solvers took a total of %.2f seconds\n",
     2562                rootTimeCpu);
     2563      else
     2564        sprintf(general,"Multiple root solvers took a total of %.2f seconds (%.2f elapsed)\n",
     2565                rootTimeCpu,CoinGetTimeOfDay()-startTimeRoot);
     2566      messageHandler()->message(CBC_GENERAL,
     2567                                messages())
     2568        << general << CoinMessageEol ;
     2569      CoinSort_2(value,value+numberSolutions,which);
     2570      // to get name
     2571      CbcHeuristicRINS dummyHeuristic;
     2572      dummyHeuristic.setHeuristicName("Multiple root solvers");
     2573      lastHeuristic_=&dummyHeuristic;
     2574      for (int i=0;i<numberSolutions;i++) {
     2575        double objValue = - value[i];
     2576        if (objValue<getCutoff()) {
     2577          int iModel=which[i];
     2578          setBestSolution(CBC_ROUNDING,objValue,
     2579                          rootModels[iModel]->bestSolution());
     2580        }
     2581      }
     2582      lastHeuristic_=NULL;
     2583      delete [] which;
     2584      delete [] value;
     2585    }
    24162586    // Do heuristics
    2417     if (numberObjects_)
     2587    if (numberObjects_&&!rootModels)
    24182588        doHeuristicsAtRoot();
    24192589    if (solverCharacteristics_->solutionAddsCuts()) {
     
    24962666            // User event
    24972667            if (!eventHappened_ && feasible) {
     2668                if (rootModels) {
     2669                  int numberModels = multipleRootTries_%100;
     2670                  const OsiSolverInterface ** solvers = new
     2671                    const OsiSolverInterface * [numberModels];
     2672                  int numberRows=continuousSolver_->getNumRows();
     2673                  int maxCuts=0;
     2674                  for (int i=0;i<numberModels;i++) {
     2675                    solvers[i]=rootModels[i]->solver();
     2676                    int numberRows2=solvers[i]->getNumRows();
     2677                    assert (numberRows2>=numberRows);
     2678                    maxCuts += numberRows2-numberRows;
     2679                    // accumulate statistics
     2680                    for (int j=0;j<numberCutGenerators_;j++) {
     2681                      generator_[j]->addStatistics(rootModels[i]->cutGenerator(j));
     2682                    }
     2683                  }
     2684                  for (int j=0;j<numberCutGenerators_;j++) {
     2685                    generator_[j]->scaleBackStatistics(numberModels);
     2686                  }
     2687                  CbcRowCuts rowCut(maxCuts);
     2688                  for (int iModel=0;iModel<numberModels;iModel++) {
     2689                    int numberRows2=solvers[iModel]->getNumRows();
     2690                    const CoinPackedMatrix * rowCopy = solvers[iModel]->getMatrixByRow();
     2691                    const int * rowLength = rowCopy->getVectorLengths();
     2692                    const double * elements = rowCopy->getElements();
     2693                    const int * column = rowCopy->getIndices();
     2694                    const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
     2695                    const double * rowLower = solvers[iModel]->getRowLower();
     2696                    const double * rowUpper = solvers[iModel]->getRowUpper();
     2697                    for (int iRow=numberRows;iRow<numberRows2;iRow++) {
     2698                      OsiRowCut rc;
     2699                      rc.setLb(rowLower[iRow]);
     2700                      rc.setUb(rowUpper[iRow]);
     2701                      CoinBigIndex start = rowStart[iRow];
     2702                      rc.setRow(rowLength[iRow],column+start,elements+start,false);
     2703                      rowCut.addCutIfNotDuplicate(rc);
     2704                    }
     2705                    //int cutsAdded=rowCut.numberCuts()-numberCuts;
     2706                    //numberCuts += cutsAdded;
     2707                    //printf("Model %d gave %d cuts (out of %d possible)\n",
     2708                    //     iModel,cutsAdded,numberRows2-numberRows);
     2709                  }
     2710                  rowCut.addCuts(cuts);
     2711                  char printBuffer[200];
     2712                  sprintf(printBuffer,"%d solvers added %d different cuts out of pool of %d",
     2713                          numberModels,cuts.sizeRowCuts(),maxCuts);
     2714                  messageHandler()->message(CBC_GENERAL, messages())
     2715                    << printBuffer << CoinMessageEol ;
     2716                  delete [] solvers;
     2717                }
    24982718                feasible = solveWithCuts(cuts, maximumCutPassesAtRoot_,
    24992719                                         NULL);
     
    25882808        }
    25892809    }
     2810    if (rootModels) {
     2811      int numberModels = multipleRootTries_%100;
     2812      for (int i=0;i<numberModels;i++)
     2813        delete rootModels[i];
     2814      delete [] rootModels;
     2815    }
    25902816    // check extra info on feasibility
    25912817    if (!solverCharacteristics_->mipFeasible())
    25922818        feasible = false;
    2593 
     2819    // If max nodes==0 - don't do strong branching
     2820    if (!getMaximumNodes()) {
     2821      if (feasible)
     2822        feasible=false;
     2823      else
     2824        setMaximumNodes(1); //allow to stop on success
     2825    }
     2826    topOfTree_=NULL;
    25942827#ifdef CLP_RESOLVE
    25952828    if ((moreSpecialOptions_&2097152)!=0&&!parentModel_&&feasible) {
     
    28263059    //solverCharacteristics_->setSolver(solver_);
    28273060    if (feasible) {
    2828         //#define HOTSTART -1
     3061#define HOTSTART -1
    28293062#if HOTSTART<0
    28303063        if (bestSolution_ && !parentModel_ && !hotstartSolution_ &&
     
    29313164                        changed = true;
    29323165                    } else {
    2933                         if (nForced + nZeroDj > 50 ||
    2934                                 (nForced + nZeroDj)*4 > numberIntegers_)
     3166                        if (nForced + nZeroDj > 5000 ||
     3167                                (nForced + nZeroDj)*2 > numberIntegers_)
    29353168                            possible = false;
    29363169                    }
     
    32103443                CglImplication implication(probingInfo_);
    32113444                addCutGenerator(&implication, 1, "ImplicationCuts", true, false, false, -200);
     3445                generator_[numberCutGenerators_-1]->setGlobalCuts(true);
    32123446            } else {
    32133447                delete probingInfo_;
     
    32973531            newNode->initializeInfo() ;
    32983532            tree_->push(newNode) ;
     3533            // save pointer to root node - so can pick up bounds
     3534            if (!topOfTree_)
     3535              topOfTree_ = dynamic_cast<CbcFullNodeInfo *>(newNode->nodeInfo()) ;
    32993536            if (statistics_) {
    33003537                if (numberNodes2_ == maximumStatistics_) {
     
    42074444      outside world will be able to obtain information about the solution using
    42084445      public methods.
     4446
     4447      Don't replace if we are trying to save cuts
    42094448    */
    4210     if (bestSolution_ && (solverCharacteristics_->solverType() < 2 || solverCharacteristics_->solverType() == 4)) {
     4449    if (bestSolution_ && (solverCharacteristics_->solverType() < 2 || solverCharacteristics_->solverType() == 4) &&
     4450        ((specialOptions_&8388608)==0||(specialOptions_&2048)!=0)) {
    42114451        setCutoff(1.0e50) ; // As best solution should be worse than cutoff
     4452        // change cutoff as constraint if wanted
     4453        if ((moreSpecialOptions_&16777216)!=0) {
     4454          if (solver_->getNumRows()>=numberRowsAtContinuous_)
     4455            solver_->setRowUpper(numberRowsAtContinuous_-1,1.0e50);
     4456        }
    42124457        // also in continuousSolver_
    42134458        if (continuousSolver_) {
     
    42684513      Destroy global cuts by replacing with an empty OsiCuts object.
    42694514    */
    4270     globalCuts_ = OsiCuts() ;
     4515    globalCuts_ = CbcRowCuts() ;
     4516    delete globalConflictCuts_;
     4517    globalConflictCuts_=NULL;
    42714518    if (!bestSolution_) {
    42724519        // make sure lp solver is infeasible
     
    44614708        currentSolution_(NULL),
    44624709        testSolution_(NULL),
     4710        globalConflictCuts_(NULL),
    44634711        minimumDrop_(1.0e-4),
    44644712        numberSolutions_(0),
     
    44984746        specialOptions_(0),
    44994747        moreSpecialOptions_(0),
     4748        topOfTree_(NULL),
    45004749        subTreeModel_(NULL),
     4750        heuristicModel_(NULL),
    45014751        numberStoppedSubTrees_(0),
    45024752        presolve_(0),
     
    45374787        maximumRows_(0),
    45384788        randomSeed_(-1),
     4789        multipleRootTries_(0),
    45394790        currentDepth_(0),
    45404791        whichGenerator_(NULL),
     
    45524803        numberNewCuts_(0),
    45534804        searchStrategy_(-1),
     4805        strongStrategy_(0),
    45544806        numberStrongIterations_(0),
    45554807        resolveAfterTakeOffCuts_(true),
     
    46234875        sumChangeObjective1_(0.0),
    46244876        sumChangeObjective2_(0.0),
     4877        globalConflictCuts_(NULL),
    46254878        minimumDrop_(1.0e-4),
    46264879        numberSolutions_(0),
     
    46564909        specialOptions_(0),
    46574910        moreSpecialOptions_(0),
     4911        topOfTree_(NULL),
    46584912        subTreeModel_(NULL),
     4913        heuristicModel_(NULL),
    46594914        numberStoppedSubTrees_(0),
    46604915        presolve_(0),
     
    46954950        maximumRows_(0),
    46964951        randomSeed_(-1),
     4952        multipleRootTries_(0),
    46974953        currentDepth_(0),
    46984954        whichGenerator_(NULL),
     
    47104966        numberNewCuts_(0),
    47114967        searchStrategy_(-1),
     4968        strongStrategy_(0),
    47124969        numberStrongIterations_(0),
    47134970        resolveAfterTakeOffCuts_(true),
     
    48975154        sumChangeObjective1_(rhs.sumChangeObjective1_),
    48985155        sumChangeObjective2_(rhs.sumChangeObjective2_),
     5156        globalConflictCuts_(NULL),
    48995157        minimumDrop_(rhs.minimumDrop_),
    49005158        numberSolutions_(rhs.numberSolutions_),
     
    49125170        specialOptions_(rhs.specialOptions_),
    49135171        moreSpecialOptions_(rhs.moreSpecialOptions_),
     5172        topOfTree_(NULL),
    49145173        subTreeModel_(rhs.subTreeModel_),
     5174        heuristicModel_(NULL),
    49155175        numberStoppedSubTrees_(rhs.numberStoppedSubTrees_),
    49165176        presolve_(rhs.presolve_),
     
    49405200        maximumRows_(0),
    49415201        randomSeed_(rhs.randomSeed_),
     5202        multipleRootTries_(rhs.multipleRootTries_),
    49425203        currentDepth_(0),
    49435204        whichGenerator_(NULL),
     
    49555216        numberNewCuts_(rhs.numberNewCuts_),
    49565217        searchStrategy_(rhs.searchStrategy_),
     5218        strongStrategy_(rhs.strongStrategy_),
    49575219        numberStrongIterations_(rhs.numberStrongIterations_),
    49585220        resolveAfterTakeOffCuts_(rhs.resolveAfterTakeOffCuts_),
     
    52625524        moreSpecialOptions_ = rhs.moreSpecialOptions_;
    52635525        subTreeModel_ = rhs.subTreeModel_;
     5526        heuristicModel_ = NULL;
    52645527        numberStoppedSubTrees_ = rhs.numberStoppedSubTrees_;
    52655528        presolve_ = rhs.presolve_;
     
    52855548        maximumCutPasses_ = rhs.maximumCutPasses_;
    52865549        randomSeed_ = rhs.randomSeed_;
     5550        multipleRootTries_ = rhs.multipleRootTries_;
    52875551        preferredWay_ = rhs.preferredWay_;
    52885552        currentPassNumber_ = rhs.currentPassNumber_;
     
    52905554        memcpy(dblParam_, rhs.dblParam_, sizeof(dblParam_));
    52915555        globalCuts_ = rhs.globalCuts_;
     5556        delete globalConflictCuts_;
     5557        globalConflictCuts_=NULL;
    52925558        int i;
    52935559        for (i = 0; i < numberCutGenerators_; i++) {
     
    53405606        masterThread_ = NULL;
    53415607        searchStrategy_ = rhs.searchStrategy_;
     5608        strongStrategy_ = rhs.strongStrategy_;
    53425609        numberStrongIterations_ = rhs.numberStrongIterations_;
    53435610        strongInfo_[0] = rhs.strongInfo_[0];
     
    55685835    delete cutModifier_;
    55695836    cutModifier_ = NULL;
     5837    topOfTree_ = NULL;
    55705838    resetModel();
    55715839}
     
    56485916        tree_->cleanTree(this, -1.0e100, bestPossibleObjective_) ;
    56495917    subTreeModel_ = NULL;
     5918    heuristicModel_ = NULL;
    56505919    numberStoppedSubTrees_ = 0;
    56515920    numberInfeasibleNodes_ = 0;
     
    56645933    numberNewCuts_ = 0;
    56655934    searchStrategy_ = -1;
     5935    strongStrategy_ = 0;
    56665936    numberStrongIterations_ = 0;
    56675937    // Parameters which need to be reset
     
    56715941    dblParam_[CbcCurrentObjectiveValue] = 1.0e100;
    56725942    dblParam_[CbcCurrentMinimizationObjectiveValue] = 1.0e100;
     5943    delete globalConflictCuts_;
     5944    globalConflictCuts_=NULL;
    56735945}
    56745946/* Most of copy constructor
     
    56925964    maximumCutPasses_ =  rhs.maximumCutPasses_;
    56935965    randomSeed_ = rhs.randomSeed_;
     5966    multipleRootTries_ = rhs.multipleRootTries_;
    56945967    preferredWay_ = rhs.preferredWay_;
    56955968    resolveAfterTakeOffCuts_ = rhs.resolveAfterTakeOffCuts_;
     
    66326905#endif
    66336906    bool feasible = true ;
    6634     int lastNumberCuts = 0 ;
    66356907    int violated = 0 ;
    66366908    int numberRowsAtStart = solver_->getNumRows() ;
     
    66406912
    66416913    numberOldActiveCuts_ = numberRowsAtStart - numberRowsAtContinuous_ ;
    6642     numberNewCuts_ = 0 ;
     6914    numberNewCuts_ = cuts.sizeRowCuts();
     6915    int lastNumberCuts = numberNewCuts_ ;
     6916    if (numberNewCuts_) {
     6917      // from multiple root passes
     6918      const OsiRowCut **addCuts = new const OsiRowCut* [numberNewCuts_];
     6919      for (int i = 0; i < numberNewCuts_; i++) {
     6920        addCuts[i] = cuts.rowCutPtr(i);
     6921      }
     6922      solver_->applyRowCuts(numberNewCuts_, addCuts);
     6923      delete [] addCuts;
     6924    }
    66436925
    66446926    bool onOptimalPath = false ;
     
    66806962        objectiveValue = node->objectiveValue();
    66816963    }
     6964    int save = moreSpecialOptions_;
     6965    if ((moreSpecialOptions_&4194304)!=0)
     6966      moreSpecialOptions_ |= 8388608;
    66826967    int returnCode = resolve(node ? node->nodeInfo() : NULL, 1);
     6968    moreSpecialOptions_=save;
     6969#ifdef CONFLICT_CUTS
     6970#ifdef COIN_HAS_CLP
     6971    // if infeasible conflict analysis
     6972    if (solver_->isProvenPrimalInfeasible()&&!parentModel_&&
     6973        (moreSpecialOptions_&4194304)!=0&&clpSolver) {
     6974      //printf("infeasible - do conflict analysis\n");
     6975      assert (topOfTree_);
     6976      int iType=1;
     6977      OsiRowCut * cut = clpSolver->modelCut(topOfTree_->lower(),
     6978                                            topOfTree_->upper(),
     6979                                            numberRowsAtContinuous_,whichGenerator_,iType);
     6980      if (cut) {
     6981        printf("XXXXXX cut\n");
     6982        //cut->print();
     6983        if (!iType) {
     6984          makeGlobalCut(cut) ;
     6985          if ((specialOptions_&1) != 0) {
     6986            debugger = continuousSolver_->getRowCutDebugger() ;
     6987            if (debugger)
     6988              if (debugger->invalidCut(*cut)) {
     6989                continuousSolver_->applyRowCuts(1,cut);
     6990                continuousSolver_->writeMps("bad");
     6991              }
     6992              CoinAssert (!debugger->invalidCut(*cut));
     6993          }
     6994        } else {
     6995          makePartialCut(cut);
     6996        }
     6997        delete cut;
     6998      }
     6999    }
     7000    if ((moreSpecialOptions_&4194304)!=0&&solver_->isProvenPrimalInfeasible()
     7001        &&clpSolver&&clpSolver->lastAlgorithm()==2&&
     7002        clpSolver->getModelPtr()->infeasibilityRay()&&
     7003        !parentModel_) {
     7004      printf("ray exists\n");
     7005    }
     7006#endif
     7007#endif   
    66837008#if COIN_DEVELOP>1
    66847009    //if (!solver_->getIterationCount()&&solver_->isProvenOptimal())
     
    70187343                (numberNodes_ % howOftenGlobalScan_) == 0 &&
    70197344                (doCutsNow(1) || true)) {
    7020             int numberCuts = globalCuts_.sizeColCuts() ;
    7021             int i;
     7345            // global column cuts now done in node at top of tree
     7346            int numberCuts = globalCuts_.sizeRowCuts() ;
    70227347            // possibly extend whichGenerator
    70237348            resizeWhichGenerator(numberViolated, numberViolated + numberCuts);
    7024             for ( i = 0 ; i < numberCuts ; i++) {
    7025                 OsiColCut *thisCut = globalCuts_.colCutPtr(i) ;
    7026                 if (thisCut->violated(cbcColSolution_) > primalTolerance ||
    7027                         thisCut->effectiveness() == COIN_DBL_MAX) {
    7028 #ifdef CLP_INVESTIGATE
    7029                     if (thisCut->violated(cbcColSolution_) > primalTolerance)
    7030                         printf("Global cut added - violation %g\n",
    7031                                thisCut->violated(cbcColSolution_)) ;
    7032 #endif
    7033                     whichGenerator_[numberViolated++] = -1;
    7034 #ifndef GLOBAL_CUTS_JUST_POINTERS
    7035                     theseCuts.insert(*thisCut) ;
    7036 #else
    7037                     theseCuts.insert(thisCut) ;
    7038 #endif
    7039                 }
    7040             }
    7041             numberCuts = globalCuts_.sizeRowCuts() ;
    7042             // possibly extend whichGenerator
    7043             resizeWhichGenerator(numberViolated, numberViolated + numberCuts);
    7044             for ( i = 0; i < numberCuts; i++) {
     7349            for (int i = 0; i < numberCuts; i++) {
    70457350                OsiRowCut * thisCut = globalCuts_.rowCutPtr(i) ;
    70467351                if (thisCut->violated(cbcColSolution_) > primalTolerance ||
    70477352                        thisCut->effectiveness() == COIN_DBL_MAX) {
    7048                     //printf("Global cut added - violation %g\n",
    7049                     // thisCut->violated(cbcColSolution_)) ;
     7353#ifndef NDEBUG
     7354                  printf("Global cut added - violation %g\n",
     7355                           thisCut->violated(cbcColSolution_)) ;
     7356#endif
    70507357                    whichGenerator_[numberViolated++] = -1;
    70517358#ifndef GLOBAL_CUTS_JUST_POINTERS
     
    76657972                        whichGenerator_[numberNewCuts_++] = -1;
    76667973#ifndef GLOBAL_CUTS_JUST_POINTERS
    7667                         cuts.insert(globalCuts_.rowCut(i)) ;
     7974                        cuts.insert(*globalCuts_.rowCutPtr(i)) ;
    76687975#else
    76697976                        OsiRowCut * rowCutPointer = globalCuts_.rowCutPtr(i);
     
    78308137        */
    78318138        if (!numberNodes_) {
     8139          //solver_->writeMps("second");
    78328140            if (numberRowsAdded)
    78338141                handler_->message(CBC_CUTS_STATS, messages_)
     
    80058313        for (i = 0; i < numberNewCuts_; i++) {
    80068314            int iGenerator = whichGenerator_[i];
     8315            if (iGenerator>=0)
     8316              iGenerator=iGenerator%10000;
    80078317            if (iGenerator >= 0 && iGenerator < numberCutGenerators_)
    80088318                count[iGenerator]++ ;
    80098319        }
     8320        // add in any active cuts if at root node (for multiple solvers)
     8321        if (!numberNodes_) {
     8322          for (i = 0; i < numberCutGenerators_; i++)
     8323            count[i] += generator_[i]->numberCutsActive();
     8324        }
    80108325        double totalCuts = 0.0 ;
    80118326        //#define JUST_ACTIVE
     
    83588673                for (i = 0; i < numberNewCuts_; i++) {
    83598674                    int iGenerator = whichGenerator_[i];
     8675                    if (iGenerator>=0)
     8676                      iGenerator=iGenerator%10000;
    83608677                    if (iGenerator >= 0)
    83618678                        generator_[iGenerator]->incrementNumberCutsActive();
     
    85758892            }
    85768893            whichGenerator_[numberBefore++] = i ;
     8894            if (!numberNodes_||generator_[i]->globalCuts())
     8895              whichGenerator_[numberBefore-1]=i+10000;
    85778896            if (thisCut->lb() > thisCut->ub())
    85788897                status = -1; // sub-problem is infeasible
     
    85828901                newCut.setGloballyValid(true);
    85838902                newCut.mutableRow().setTestForDuplicateIndex(false);
    8584                 globalCuts_.insert(newCut) ;
     8903                globalCuts_.addCutIfNotDuplicate(newCut) ;
     8904                whichGenerator_[numberBefore-1] = i+10000 ;
    85858905            }
    85868906        }
     
    85978917                const OsiRowCut * thisCut = theseCuts.rowCutPtr(j) ;
    85988918                whichGenerator_[numberBefore++] = i ;
     8919                if (!numberNodes_||generator_[i]->globalCuts())
     8920                  whichGenerator_[numberBefore-1]=i+10000;
    85998921                if (thisCut->globallyValid()) {
    86008922                    // add to global list
     
    86028924                    newCut.setGloballyValid(true);
    86038925                    newCut.mutableRow().setTestForDuplicateIndex(false);
    8604                     globalCuts_.insert(newCut) ;
     8926                    globalCuts_.addCutIfNotDuplicate(newCut) ;
     8927                    whichGenerator_[numberBefore-1]=i+10000;
    86058928                }
    86068929            }
     
    86108933            const OsiColCut * thisCut = theseCuts.colCutPtr(j) ;
    86118934            if (thisCut->globallyValid()) {
    8612                 // add to global list
    8613                 OsiColCut newCut(*thisCut);
    8614                 newCut.setGloballyValid(true);
    8615                 globalCuts_.insert(newCut) ;
     8935                // fix
     8936                makeGlobalCut(thisCut);
    86168937            }
    86178938        }
     
    86468967                    printf("Old cut added - violation %g\n",
    86478968                           thisCut->violated(cbcColSolution_)) ;
    8648                 whichGenerator_[numberOld++] = -1;
     8969                whichGenerator_[numberOld++] = -3;
    86498970                theseCuts.insert(*thisCut) ;
    86508971            }
     
    1089511216                            newCut.setGloballyValid(true);
    1089611217                            newCut.mutableRow().setTestForDuplicateIndex(false);
    10897                             globalCuts_.insert(newCut) ;
     11218                            globalCuts_.addCutIfNotDuplicate(newCut) ;
    1089811219                        } else {
    1089911220                            // obviously wrong
     
    1112811449            //setCutoff(cutoff*direction);
    1112911450            setCutoff(cutoff);
     11451            // change cutoff as constraint if wanted
     11452            if ((moreSpecialOptions_&16777216)!=0) {
     11453              if (solver_->getNumRows()>=numberRowsAtContinuous_) {
     11454                double offset;
     11455                solver_->getDblParam(OsiObjOffset, offset);
     11456                solver_->setRowUpper(numberRowsAtContinuous_-1,cutoff+offset);
     11457              }
     11458            }
    1113011459
    1113111460            if (how == CBC_ROUNDING)
     
    1119111520                            newCut.setGloballyValid(true);
    1119211521                            newCut.mutableRow().setTestForDuplicateIndex(false);
    11193                             globalCuts_.insert(newCut) ;
     11522                            globalCuts_.addCutIfNotDuplicate(newCut) ;
    1119411523                            generator_[i]->incrementNumberCutsInTotal();
    1119511524                        }
     
    1120111530                const OsiColCut * thisCut = theseCuts.colCutPtr(i);
    1120211531                if (thisCut->globallyValid()) {
    11203                     // add to global list
    11204                     OsiColCut newCut(*thisCut);
    11205                     newCut.setGloballyValid(true);
    11206                     globalCuts_.insert(newCut) ;
     11532                    // fix
     11533                    makeGlobalCut(thisCut);
    1120711534                }
    1120811535            }
     
    1130211629                //setCutoff(cutoff*direction);
    1130311630                setCutoff(cutoff);
     11631                // change cutoff as constraint if wanted
     11632                if ((moreSpecialOptions_&16777216)!=0) {
     11633                  if (solver_->getNumRows()>=numberRowsAtContinuous_) {
     11634                    double offset;
     11635                    solver_->getDblParam(OsiObjOffset, offset);
     11636                    solver_->setRowUpper(numberRowsAtContinuous_-1,cutoff+offset);
     11637                  }
     11638                }
    1130411639
    1130511640                numberSolutions_++;
     
    1192012255                thisCut.setRow(rowLength[iRow], column + start, elementByRow + start, false);
    1192112256                thisCut.setGloballyValid(true);
    11922                 globalCuts_.insert(thisCut) ;
     12257                globalCuts_.addCutIfNotDuplicate(thisCut) ;
    1192312258            }
    1192412259        }
     
    1193512270    newCut.setGloballyValidAsInteger(2);
    1193612271    newCut.mutableRow().setTestForDuplicateIndex(false);
    11937     globalCuts_.insert(newCut) ;
     12272    globalCuts_.addCutIfNotDuplicate(newCut) ;
    1193812273}
    1193912274// Make given cut into a global cut
     
    1194412279    newCut.setGloballyValid(true);
    1194512280    newCut.mutableRow().setTestForDuplicateIndex(false);
    11946     globalCuts_.insert(newCut) ;
     12281    globalCuts_.addCutIfNotDuplicate(newCut) ;
    1194712282}
    1194812283// Make given column cut into a global cut
     
    1195012285CbcModel::makeGlobalCut(const OsiColCut * cut)
    1195112286{
    11952     OsiColCut newCut(*cut);
    11953     newCut.setGloballyValidAsInteger(2);
    11954     globalCuts_.insert(newCut) ;
     12287  const double * lower;
     12288  const double * upper;
     12289  if (topOfTree_) {
     12290    lower = topOfTree_->lower();
     12291    upper = topOfTree_->upper();
     12292  } else {
     12293    lower = solver_->getColLower();
     12294    upper = solver_->getColUpper();
     12295  }
     12296  int nLower=cut->lbs().getNumElements();
     12297  const int * indexLower=cut->lbs().getIndices();
     12298  const double * boundLower=cut->lbs().getElements();
     12299  for (int i=0;i<nLower;i++) {
     12300    int iColumn=indexLower[i];
     12301    double newValue=CoinMax(lower[iColumn],boundLower[iColumn]);
     12302    if (topOfTree_)
     12303      topOfTree_->setColLower(iColumn,newValue);
     12304    else
     12305      solver_->setColLower(iColumn,newValue);
     12306  }
     12307  int nUpper=cut->ubs().getNumElements();
     12308  const int * indexUpper=cut->ubs().getIndices();
     12309  const double * boundUpper=cut->ubs().getElements();
     12310  for (int i=0;i<nUpper;i++) {
     12311    int iColumn=indexUpper[i];
     12312    double newValue=CoinMin(upper[iColumn],boundUpper[iColumn]);
     12313    if (topOfTree_)
     12314      topOfTree_->setColUpper(iColumn,newValue);
     12315    else
     12316      solver_->setColUpper(iColumn,newValue);
     12317  }
    1195512318}
    1195612319// Make given column cut into a global cut
     
    1195812321CbcModel::makeGlobalCut(const OsiColCut & cut)
    1195912322{
    11960     OsiColCut newCut(cut);
     12323  const double * lower;
     12324  const double * upper;
     12325  if (topOfTree_) {
     12326    lower = topOfTree_->lower();
     12327    upper = topOfTree_->upper();
     12328  } else {
     12329    lower = solver_->getColLower();
     12330    upper = solver_->getColUpper();
     12331  }
     12332  int nLower=cut.lbs().getNumElements();
     12333  const int * indexLower=cut.lbs().getIndices();
     12334  const double * boundLower=cut.lbs().getElements();
     12335  for (int i=0;i<nLower;i++) {
     12336    int iColumn=indexLower[i];
     12337    double newValue=CoinMax(lower[iColumn],boundLower[iColumn]);
     12338    if (topOfTree_)
     12339      topOfTree_->setColLower(iColumn,newValue);
     12340    else
     12341      solver_->setColLower(iColumn,newValue);
     12342  }
     12343  int nUpper=cut.ubs().getNumElements();
     12344  const int * indexUpper=cut.ubs().getIndices();
     12345  const double * boundUpper=cut.ubs().getElements();
     12346  for (int i=0;i<nUpper;i++) {
     12347    int iColumn=indexUpper[i];
     12348    double newValue=CoinMin(upper[iColumn],boundUpper[iColumn]);
     12349    if (topOfTree_)
     12350      topOfTree_->setColUpper(iColumn,newValue);
     12351    else
     12352      solver_->setColUpper(iColumn,newValue);
     12353  }
     12354}
     12355// Make partial cut into a global cut and save
     12356void
     12357CbcModel::makePartialCut(const OsiRowCut * partialCut)
     12358{
     12359  // get greedy cut
     12360  double bSum = partialCut->lb();
     12361  assert (bSum<0.0);
     12362  int nConflict = partialCut->row().getNumElements();
     12363  const int * column = partialCut->row().getIndices();
     12364  const double * element = partialCut->row().getElements();
     12365  double * originalLower = topOfTree_->mutableLower();
     12366  const double * columnLower = solver_->getColLower();
     12367  double * originalUpper = topOfTree_->mutableUpper();
     12368  const double * columnUpper = solver_->getColUpper();
     12369  int nC=nConflict;
     12370  while (nConflict) {
     12371    int iColumn = column[nConflict-1];
     12372    double farkasValue = element[nConflict-1];
     12373    double change;
     12374    if (farkasValue>0.0) {
     12375      change=farkasValue*(originalUpper[iColumn]-columnUpper[iColumn]);
     12376    } else {
     12377      change=farkasValue*(originalLower[iColumn]-columnLower[iColumn]);
     12378    }
     12379    if (bSum+change>-1.0e-4)
     12380      break;
     12381    nConflict--;
     12382    bSum += change;
     12383  }
     12384  OsiRowCut newCut;
     12385  newCut.setUb(COIN_DBL_MAX);
     12386  double lo=1.0;
     12387  double * values = new double[nConflict];
     12388  for (int i=0;i<nConflict;i++) {
     12389    int iColumn = column[i];
     12390    if (originalLower[iColumn]==columnLower[iColumn]) {
     12391      // must be at least one higher
     12392      values[i]=1.0;
     12393      lo += originalLower[iColumn];
     12394    } else {
     12395      // must be at least one lower
     12396      values[i]=-1.0;
     12397      lo -= originalUpper[iColumn];
     12398    }
     12399  }
     12400  newCut.setLb(lo);
     12401  newCut.setRow(nConflict,column,values);
     12402  printf("CUTa has %d (started at %d) - final bSum %g\n",nConflict,nC,bSum);
     12403  if (nConflict>1) {
    1196112404    newCut.setGloballyValidAsInteger(2);
    11962     globalCuts_.insert(newCut) ;
     12405    newCut.mutableRow().setTestForDuplicateIndex(false);
     12406    globalCuts_.addCutIfNotDuplicate(newCut) ;
     12407  } else {
     12408    // change bounds
     12409    int iColumn=column[0];
     12410    if (values[0]<0.0) {
     12411      // change upper bound
     12412      double newUpper = -lo;
     12413      assert (newUpper<originalUpper[iColumn]);
     12414      printf("Changing upper bound on %d from %g to %g\n",
     12415             iColumn,originalUpper[iColumn],newUpper);
     12416      originalUpper[iColumn]=newUpper;
     12417    } else {
     12418      // change lower bound
     12419      double newLower = lo;
     12420      assert (newLower>originalLower[iColumn]);
     12421      printf("Changing lower bound on %d from %g to %g\n",
     12422             iColumn,originalLower[iColumn],newLower);
     12423      originalLower[iColumn]=newLower;
     12424    }
     12425  }
     12426  // add to partial cuts
     12427  if (globalConflictCuts_) {
     12428    globalConflictCuts_->addCutIfNotDuplicateWhenGreedy(*partialCut,2);
     12429  }
     12430  delete [] values;
     12431}
     12432// Make partial cuts into global cuts
     12433void
     12434CbcModel::makeGlobalCuts()
     12435{
    1196312436}
    1196412437void
     
    1213012603        ClpSimplex * clpSimplex = clpSolver->getModelPtr();
    1213112604        int save = clpSimplex->specialOptions();
    12132         clpSimplex->setSpecialOptions(save | 0x11000000); // say is Cbc (and in branch and bound)
     12605        if ((moreSpecialOptions_&8388608)==0)
     12606          clpSimplex->setSpecialOptions(save | 0x11000000); // say is Cbc (and in branch and bound)
     12607        else
     12608          clpSimplex->setSpecialOptions(save | 0x11200000); // say is Cbc (and in branch and bound - but save ray)
    1213312609        int save2 = clpSolver->specialOptions();
    1213412610        if (false && (save2&2048) == 0) {
     
    1268613162            }
    1268713163#endif
     13164#ifdef COIN_HAS_CLP
     13165            int save;
     13166            OsiClpSolverInterface * clpSolver
     13167              = dynamic_cast<OsiClpSolverInterface *> (solver_);
     13168            if (clpSolver&&(moreSpecialOptions_&4194304)!=0) {
     13169              ClpSimplex * clpSimplex = clpSolver->getModelPtr();
     13170              save=clpSimplex->specialOptions();
     13171              clpSimplex->setSpecialOptions(save | 0x11200000); // say is Cbc (and in branch and bound - but save ray)
     13172            }
     13173#endif
    1268813174            if (numberBeforeTrust_ == 0 ) {
    1268913175                anyAction = newNode->chooseBranch(this, oldNode, numberPassesLeft) ;
     
    1269313179                    anyAction = newNode->chooseBranch(this, oldNode, numberPassesLeft) ; // dynamic did nothing
    1269413180            }
     13181#ifdef COIN_HAS_CLP
     13182            if (clpSolver&&(moreSpecialOptions_&4194304)!=0) {
     13183              ClpSimplex * clpSimplex = clpSolver->getModelPtr();
     13184              clpSimplex->setSpecialOptions(save);
     13185            }
     13186#endif
    1269513187            /*
    1269613188              We're on the new (Osi) side of the branching hierarchy.
     
    1307213564        cutoff = objectiveValue - increment ;
    1307313565        setCutoff(cutoff) ;
     13566        // change cutoff as constraint if wanted
     13567        if ((moreSpecialOptions_&16777216)!=0) {
     13568          if (solver_->getNumRows()>=numberRowsAtContinuous_) {
     13569            double offset;
     13570            solver_->getDblParam(OsiObjOffset, offset);
     13571            solver_->setRowUpper(numberRowsAtContinuous_-1,cutoff+offset);
     13572          }
     13573        }
    1307413574    }
    1307513575    int n = CoinMax(numberColumns, solver_->getNumCols());
     
    1377714277            = dynamic_cast<OsiClpSolverInterface *> (solver_);
    1377814278            if ((clpSolver || (specialOptions_&16384) != 0) && fastNodeDepth_ < -1
    13779                     && (specialOptions_&2048) == 0) {
     14279                && (specialOptions_&2048) == 0) {
    1378014280#define FATHOM_BIAS -2
    1378114281                if (numberNodes_ == 1) {
     
    1390414404                        numberExtraNodes_ += info->numberNodesExplored_;
    1390514405                        numberExtraIterations_ += info->numberIterations_;
     14406                        char general[200];
     14407                        int fathomStatus=info->nNodes_;
     14408                        if (feasible)
     14409                          fathomStatus=1;
     14410                        sprintf(general, "fathom took %d nodes, %d iterations - status %d",
     14411                                info->numberNodesExplored_,
     14412                                info->numberIterations_,fathomStatus);
     14413                        messageHandler()->message(CBC_FPUMP2,
     14414                                          messages())
     14415                          << general << CoinMessageEol ;
    1390614416                        if (info->numberNodesExplored_ > 10000 /* && !feasible */
    1390714417                            && (moreSpecialOptions_&524288) == 0  && info->nNodes_>=0) {
     
    1408514595        }
    1408614596        if ((specialOptions_&1) != 0 && onOptimalPath) {
    14087             if (!solver_->getRowCutDebugger() || !feasible) {
    14088                 // dj fix did something???
    14089                 solver_->writeMpsNative("infeas2.mps", NULL, NULL, 2);
    14090                 solver_->getRowCutDebuggerAlways()->printOptimalSolution(*solver_);
    14091                 assert (solver_->getRowCutDebugger()) ;
    14092                 assert (feasible);
     14597            if(solver_->getRowCutDebuggerAlways()->optimalValue()<getCutoff()) {
     14598                if (!solver_->getRowCutDebugger() || !feasible) {
     14599                  // dj fix did something???
     14600                  solver_->writeMpsNative("infeas2.mps", NULL, NULL, 2);
     14601                  solver_->getRowCutDebuggerAlways()->printOptimalSolution(*solver_);
     14602                  assert (solver_->getRowCutDebugger()) ;
     14603                  assert (feasible);
     14604                }
    1409314605            }
    1409414606        }
     
    1411714629            bool objLim = solver_->isDualObjectiveLimitReached() ;
    1411814630            if (!feasible && !objLim) {
     14631              if(solver_->getRowCutDebuggerAlways()->optimalValue()<getCutoff()) {
    1411914632                printf("infeas2\n");
    1412014633                solver_->getRowCutDebuggerAlways()->printOptimalSolution(*solver_);
     
    1412714640                solver_->initialSolve();
    1412814641                assert (!solver_->isProvenOptimal());
    14129             }
    14130             assert (feasible || objLim);
     14642                assert (feasible || objLim);
     14643              }
     14644            }
    1413114645        }
    1413214646        bool checkingNode = false;
    1413314647        if (feasible) {
    14134 #ifdef FUNNY_BRANCHING
     14648#ifdef FUNNY_BRANCHING2
    1413514649            // Far too clever
    1413614650            if ((numberThreads_ == -10 || true) && node->numberBranches() == 2) {
     
    1472915243    if (bestObjective > bestObjective_)
    1473015244        foundSolution = 2;
    14731     if (parallelMode() >= 0 && foundSolution) {
     15245    if (parallelMode() > 0 && foundSolution) {
    1473215246        lockThread();
    1473315247        // might as well mark all including continuous
     
    1473715251            usedInSolution_[i] = 0;
    1473815252        }
    14739         baseModel->numberSolutions_ = numberSolutions_;
     15253        baseModel->numberSolutions_++;
    1474015254        if (bestObjective_ < baseModel->bestObjective_ && bestObjective_ < baseModel->getCutoff()) {
    1474115255            baseModel->bestObjective_ = bestObjective_ ;
     
    1611816632    */
    1611916633    globalCuts_ = OsiCuts() ;
     16634    delete globalConflictCuts_;
     16635    globalConflictCuts_=NULL;
    1612016636    numberHeuristics_ = saveNumberHeuristics;
    1612116637
     
    1636316879}
    1636416880#endif
    16365 
     16881#ifdef CBC_THREAD
     16882void * doRootCbcThread(void * voidInfo)
     16883{
     16884    CbcModel * model = reinterpret_cast<CbcModel *> (voidInfo);
     16885#ifdef COIN_HAS_CLP
     16886    OsiClpSolverInterface * clpSolver
     16887      = dynamic_cast<OsiClpSolverInterface *> (model->solver());
     16888    char general[200];
     16889    if (clpSolver) {
     16890      clpSolver->getModelPtr()->dual(); // to get more randomness
     16891      clpSolver->setWarmStart(NULL);
     16892      sprintf(general,"Starting multiple root solver - %d iterations in initialSolve",
     16893              clpSolver->getIterationCount());
     16894    } else {
     16895#endif
     16896#ifdef COIN_HAS_CLP
     16897      model->initialSolve();
     16898      sprintf(general,"Solver did %d iterations in initialSolve\n",
     16899             model->solver()->getIterationCount());
     16900    }
     16901#endif
     16902    model->messageHandler()->message(CBC_GENERAL,
     16903                              model->messages())
     16904      << general << CoinMessageEol ;
     16905    model->branchAndBound();
     16906    sprintf(general,"Ending multiple root solver");
     16907    model->messageHandler()->message(CBC_GENERAL,
     16908                              model->messages())
     16909      << general << CoinMessageEol ;
     16910    return NULL;
     16911}
     16912#endif
     16913
     16914
Note: See TracChangeset for help on using the changeset viewer.