Ignore:
File:
1 edited

Legend:

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

    r1912 r2101  
    1313#include "CbcModel.hpp"
    1414#include "CbcSubProblem.hpp"
     15#include "CbcSimpleInteger.hpp"
    1516#include "OsiAuxInfo.hpp"
    1617#include  "CoinTime.hpp"
     
    2223//#define DIVE_FIX_BINARY_VARIABLES
    2324//#define DIVE_DEBUG
     25#ifdef DIVE_DEBUG
     26#define DIVE_PRINT=2
     27#endif
    2428
    2529// Default Constructor
     
    3236    downArray_ = NULL;
    3337    upArray_ = NULL;
     38    priority_ = NULL;
    3439    percentageToFix_ = 0.2;
    3540    maxIterations_ = 100;
     
    3944    whereFrom_ = 255 - 2 - 16 + 256;
    4045    decayFactor_ = 1.0;
     46    smallObjective_ = 1.0e-10;
    4147}
    4248
     
    4955    downArray_ = NULL;
    5056    upArray_ = NULL;
     57    priority_ = NULL;
    5158    // Get a copy of original matrix
    5259    assert(model.solver());
     
    5966    }
    6067    percentageToFix_ = 0.2;
     68    maxTime_ = 600;
     69    smallObjective_ = 1.0e-10;
    6170    maxIterations_ = 100;
    6271    maxSimplexIterations_ = 10000;
    6372    maxSimplexIterationsAtRoot_ = 1000000;
    64     maxTime_ = 600;
    6573    whereFrom_ = 255 - 2 - 16 + 256;
    6674    decayFactor_ = 1.0;
     75    smallObjective_ = 1.0e-10;
     76    setPriorities();
    6777}
    6878
     
    7282    delete [] downLocks_;
    7383    delete [] upLocks_;
     84    delete [] priority_;
    7485    assert (!downArray_);
    7586}
     
    106117        matrixByRow_(rhs.matrixByRow_),
    107118        percentageToFix_(rhs.percentageToFix_),
     119        maxTime_(rhs.maxTime_),
     120        smallObjective_(rhs.smallObjective_),
    108121        maxIterations_(rhs.maxIterations_),
    109122        maxSimplexIterations_(rhs.maxSimplexIterations_),
    110         maxSimplexIterationsAtRoot_(rhs.maxSimplexIterationsAtRoot_),
    111         maxTime_(rhs.maxTime_)
     123        maxSimplexIterationsAtRoot_(rhs.maxSimplexIterationsAtRoot_)
    112124{
    113125    downArray_ = NULL;
     
    117129        downLocks_ = CoinCopyOfArray(rhs.downLocks_, numberIntegers);
    118130        upLocks_ = CoinCopyOfArray(rhs.upLocks_, numberIntegers);
     131        priority_ = CoinCopyOfArray(rhs.priority_, numberIntegers);
    119132    } else {
    120133        downLocks_ = NULL;
    121134        upLocks_ = NULL;
     135        priority_ = NULL;
    122136    }
    123137}
     
    136150        maxSimplexIterationsAtRoot_ = rhs.maxSimplexIterationsAtRoot_;
    137151        maxTime_ = rhs.maxTime_;
     152        smallObjective_ = rhs.smallObjective_;
    138153        delete [] downLocks_;
    139154        delete [] upLocks_;
     155        delete [] priority_;
    140156        if (rhs.downLocks_) {
    141157            int numberIntegers = model_->numberIntegers();
    142158            downLocks_ = CoinCopyOfArray(rhs.downLocks_, numberIntegers);
    143159            upLocks_ = CoinCopyOfArray(rhs.upLocks_, numberIntegers);
     160            priority_ = CoinCopyOfArray(rhs.priority_, numberIntegers);
    144161        } else {
    145162            downLocks_ = NULL;
    146163            upLocks_ = NULL;
     164            priority_ = NULL;
    147165        }
    148166    }
     
    164182        validate();
    165183    }
     184    setPriorities();
    166185}
    167186
    168187// update model
    169 void CbcHeuristicDive::setModel(CbcModel * model)
     188void
     189CbcHeuristicDive::setModel(CbcModel * model)
    170190{
    171191    model_ = model;
     
    179199        validate();
    180200    }
    181 }
     201    setPriorities();
     202}
     203
     204// Sets priorities if any
     205void
     206CbcHeuristicDive::setPriorities()
     207{
     208  delete [] priority_;
     209  assert (model_);
     210  priority_=NULL;
     211  if (!model_->objects())
     212    return;
     213  bool gotPriorities=false;
     214  int numberIntegers = model_->numberIntegers();
     215  int priority1=-COIN_INT_MAX;
     216  int priority2=COIN_INT_MAX;
     217  smallObjective_=0.0;
     218  const double * objective = model_->solver()->getObjCoefficients();
     219  for (int i = 0; i < numberIntegers; i++) {
     220    OsiObject * object = model_->modifiableObject(i);
     221    const CbcSimpleInteger * thisOne = dynamic_cast <const CbcSimpleInteger *> (object);
     222    assert (thisOne);
     223    int iColumn = thisOne->columnNumber();
     224    smallObjective_ += objective[iColumn];
     225    int level=thisOne->priority();
     226    priority1=CoinMax(priority1,level);
     227    priority2=CoinMin(priority2,level);
     228    if (thisOne->preferredWay()!=0)
     229      gotPriorities=true;
     230  }
     231  smallObjective_ = CoinMax(1.0e-10,1.0e-5*(smallObjective_/numberIntegers));
     232  if (gotPriorities || priority1>priority2) {
     233    priority_ = new PriorityType [numberIntegers];
     234    for (int i = 0; i < numberIntegers; i++) {
     235      OsiObject * object = model_->modifiableObject(i);
     236      const CbcSimpleInteger * thisOne = dynamic_cast <const CbcSimpleInteger *> (object);
     237      int level=thisOne->priority()-priority2;
     238      assert (level<(1<<29));
     239      priority_[i].priority=static_cast<unsigned int>(level);
     240      int direction=0;
     241      if (thisOne->preferredWay()<0)
     242        direction=1;
     243      else if (thisOne->preferredWay()>0)
     244        direction=1|1;
     245        // at present don't try other way is not used
     246      priority_[i].direction=static_cast<unsigned char>(direction);
     247    }
     248  }
     249}
     250
    182251
    183252bool CbcHeuristicDive::canHeuristicRun()
    184253{
     254    if (model_->bestSolution()||model_->getNodeCount()) {
     255      if (when_==3 || (when_==4 && numberSolutionsFound_) )
     256        return false;
     257    }
    185258    return shouldHeurRun_randomChoice();
    186259}
     
    199272                           double * newSolution)
    200273{
    201 #ifdef DIVE_DEBUG
     274#if DIVE_PRINT
    202275    int nRoundInfeasible = 0;
    203276    int nRoundFeasible = 0;
     277    printf("Entering %s - fix %.1f%% maxTime %.2f maxPasses %d - max iterations %d (at root %d) - when to do %d\n",
     278    heuristicName_.c_str(),percentageToFix_*100.0,maxTime_,maxIterations_,
     279    maxSimplexIterations_,maxSimplexIterationsAtRoot_,when());
    204280#endif
    205281    int reasonToStop = 0;
     
    208284    int maxSimplexIterations = (model_->getNodeCount()) ? maxSimplexIterations_
    209285                               : maxSimplexIterationsAtRoot_;
     286    int maxIterationsInOneSolve = (maxSimplexIterations<1000000) ? 1000 : 10000;
    210287    // but can't be exactly coin_int_max
    211288    maxSimplexIterations = CoinMin(maxSimplexIterations,COIN_INT_MAX>>3);
     289    bool fixGeneralIntegers=false;
     290    int maxIterations = maxIterations_;
     291    int saveSwitches = switches_;
     292    if ((maxIterations_%10)!=0) {
     293      int digit = maxIterations_%10;
     294      maxIterations -= digit;
     295      switches_ |= 65536;
     296      if ((digit&3)!=0)
     297        fixGeneralIntegers=true;
     298    }
     299
    212300    OsiSolverInterface * solver = cloneBut(6); // was model_->solver()->clone();
    213301# ifdef COIN_HAS_CLP
     
    218306      int oneSolveIts = clpSimplex->maximumIterations();
    219307      oneSolveIts = CoinMin(1000+2*(clpSimplex->numberRows()+clpSimplex->numberColumns()),oneSolveIts);
     308      if (maxSimplexIterations>1000000)
     309        maxIterationsInOneSolve=oneSolveIts;
    220310      clpSimplex->setMaximumIterations(oneSolveIts);
    221311      if (!nodes) {
     
    271361
    272362    // vectors to store the latest variables fixed at their bounds
    273     int* columnFixed = new int [numberIntegers];
     363    int* columnFixed = new int [numberIntegers+numberColumns];
     364    int * back = columnFixed + numberIntegers;
    274365    double* originalBound = new double [numberIntegers+2*numberColumns];
    275366    double * lowerBefore = originalBound+numberIntegers;
     
    287378    // count how many fractional variables
    288379    int numberFractionalVariables = 0;
     380    for (int i=0;i<numberColumns;i++)
     381      back[i]=-1;
    289382    for (int i = 0; i < numberIntegers; i++) {
    290383        random[i] = randomNumberGenerator_.randomDouble() + 0.3;
    291384        int iColumn = integerVariable[i];
     385        back[iColumn]=i;
    292386        double value = newSolution[iColumn];
     387        // clean
     388        value = CoinMin(value,upperBefore[iColumn]);
     389        value = CoinMax(value,lowerBefore[iColumn]);
     390        newSolution[iColumn] = value;
    293391        if (fabs(floor(value + 0.5) - value) > integerTolerance) {
    294392            numberFractionalVariables++;
     
    298396    const double* reducedCost = NULL;
    299397    // See if not NLP
    300     if (model_->solverCharacteristics()->reducedCostsAccurate())
     398    if (!model_->solverCharacteristics() ||
     399        model_->solverCharacteristics()->reducedCostsAccurate())
    301400        reducedCost = solver->getReducedCost();
    302401
    303402    int iteration = 0;
     403    int numberAtBoundFixed = 0;
     404    int numberGeneralFixed = 0; // fixed as satisfied but not at bound
     405    int numberReducedCostFixed = 0;
    304406    while (numberFractionalVariables) {
    305407        iteration++;
     
    341443            }
    342444            if (direction*(solver->getObjValue() + delta) < solutionValue) {
    343 #ifdef DIVE_DEBUG
     445#if DIVE_PRINT
    344446                nRoundFeasible++;
    345447#endif
     
    388490                }
    389491            }
    390 #ifdef DIVE_DEBUG
     492#if DIVE_PRINT
    391493            else
    392494                nRoundInfeasible++;
     
    395497
    396498        // do reduced cost fixing
    397 #ifdef DIVE_DEBUG
    398         int numberFixed = reducedCostFix(solver);
    399         std::cout << "numberReducedCostFixed = " << numberFixed << std::endl;
     499#if DIVE_PRINT>1
     500        numberReducedCostFixed = reducedCostFix(solver);
    400501#else
    401502        reducedCostFix(solver);
    402503#endif
    403504
    404         int numberAtBoundFixed = 0;
     505        numberAtBoundFixed = 0;
     506        numberGeneralFixed = 0; // fixed as satisfied but not at bound
    405507#ifdef DIVE_FIX_BINARY_VARIABLES
    406508        // fix binary variables based on pseudo reduced cost
     
    480582        double gap = 1.0e30;
    481583#endif
     584        int fixPriority=COIN_INT_MAX;
    482585        if (reducedCost && true) {
    483586#ifndef JJF_ONE
    484587            cnt = fixOtherVariables(solver, solution, candidate, random);
     588            if (priority_) {
     589              for (int i = 0; i < cnt; i++) {
     590                int iColumn = candidate[i].var;
     591                if (upper[iColumn] > lower[iColumn]) {
     592                  int j=back[iColumn];
     593                  fixPriority = CoinMin(fixPriority,static_cast<int>(priority_[j].priority));
     594                }
     595              }
     596            }
    485597#else
    486598#ifdef GAP
     
    504616                if (upper[iColumn] > lower[iColumn]) {
    505617                    numberFree++;
     618                    if (priority_) {
     619                      fixPriority = CoinMin(fixPriority,static_cast<int>(priority_[i].priority));
     620                    }
    506621                    double value = newSolution[iColumn];
    507622                    if (fabs(floor(value + 0.5) - value) <= integerTolerance) {
     
    539654                int iColumn = integerVariable[i];
    540655                if (upper[iColumn] > lower[iColumn]) {
     656                    if (priority_) {
     657                      fixPriority = CoinMin(fixPriority,static_cast<int>(priority_[i].priority));
     658                    }
    541659                    double value = newSolution[iColumn];
    542660                    if (fabs(floor(value + 0.5) - value) <= integerTolerance) {
     
    548666        }
    549667        std::sort(candidate, candidate + cnt, compareBinaryVars);
     668        // If getting on fix all
     669        if (iteration*3>maxIterations_*2)
     670          fixPriority=COIN_INT_MAX;
    550671        for (int i = 0; i < cnt; i++) {
    551672            int iColumn = candidate[i].var;
     
    555676                        numberAtBoundFixed < maxNumberAtBoundToFix) {
    556677                    // fix the variable at one of its bounds
    557                     if (fabs(lower[iColumn] - value) <= integerTolerance) {
    558                         columnFixed[numberAtBoundFixed] = iColumn;
    559                         originalBound[numberAtBoundFixed] = upper[iColumn];
    560                         fixedAtLowerBound[numberAtBoundFixed] = true;
    561                         solver->setColUpper(iColumn, lower[iColumn]);
     678                    if (fabs(lower[iColumn] - value) <= integerTolerance ||
     679                        fixGeneralIntegers) {
     680                        if (priority_) {
     681                          int j=back[iColumn];
     682                          if (priority_[j].priority>fixPriority)
     683                            continue; // skip - only fix ones at high priority
     684                          int thisRound=static_cast<int>(priority_[j].direction);
     685                          if ((thisRound&1)!=0) {
     686                            // for now force way
     687                            if((thisRound&2)!=0)
     688                              continue;
     689                          }
     690                        }
     691                        if (fabs(lower[iColumn] - value) <= integerTolerance) {
     692                          columnFixed[numberAtBoundFixed] = iColumn;
     693                          originalBound[numberAtBoundFixed] = upper[iColumn];
     694                          fixedAtLowerBound[numberAtBoundFixed] = true;
     695                          solver->setColUpper(iColumn, lower[iColumn]);
     696                        } else {
     697                          // fix to interior value
     698                          numberGeneralFixed++;
     699                          double fixValue = floor(value + 0.5);
     700                          columnFixed[numberAtBoundFixed] = iColumn;
     701                          originalBound[numberAtBoundFixed] = upper[iColumn];
     702                          fixedAtLowerBound[numberAtBoundFixed] = true;
     703                          solver->setColUpper(iColumn, fixValue);
     704                          numberAtBoundFixed++;
     705                          columnFixed[numberAtBoundFixed] = iColumn;
     706                          originalBound[numberAtBoundFixed] = lower[iColumn];
     707                          fixedAtLowerBound[numberAtBoundFixed] = false;
     708                          solver->setColLower(iColumn, fixValue);
     709                        }
     710                        //if (priority_)
     711                        //printf("fixing %d (priority %d) to lower bound of %g\n",
     712                        //       iColumn,priority_[back[iColumn]].priority,lower[iColumn]);
    562713                        numberAtBoundFixed++;
    563714                    } else if (fabs(upper[iColumn] - value) <= integerTolerance) {
     715                        if (priority_) {
     716                          int j=back[iColumn];
     717                          if (priority_[j].priority>fixPriority)
     718                            continue; // skip - only fix ones at high priority
     719                          int thisRound=static_cast<int>(priority_[j].direction);
     720                          if ((thisRound&1)!=0) {
     721                            // for now force way
     722                            if((thisRound&2)==0)
     723                              continue;
     724                          }
     725                        }
    564726                        columnFixed[numberAtBoundFixed] = iColumn;
    565727                        originalBound[numberAtBoundFixed] = lower[iColumn];
    566728                        fixedAtLowerBound[numberAtBoundFixed] = false;
    567729                        solver->setColLower(iColumn, upper[iColumn]);
     730                        //if (priority_)
     731                        //printf("fixing %d (priority %d) to upper bound of %g\n",
     732                        //       iColumn,priority_[back[iColumn]].priority,upper[iColumn]);
    568733                        numberAtBoundFixed++;
    569734                    }
     
    573738            }
    574739        }
    575 #ifdef DIVE_DEBUG
    576         std::cout << "numberAtBoundFixed = " << numberAtBoundFixed << std::endl;
    577 #endif
    578740
    579741        double originalBoundBestColumn;
     
    585747                originalBoundBestColumn = upper[bestColumn];
    586748                solver->setColUpper(bestColumn, floor(bestColumnValue));
     749#ifdef DIVE_DEBUG
     750                if (priority_) {
     751                  printf("setting %d (priority %d) upper bound to %g (%g)\n",
     752                         bestColumn,priority_[back[bestColumn]].priority,floor(bestColumnValue),bestColumnValue);
     753                }
     754#endif
    587755                whichWay=0;
    588756            } else {
    589757                originalBoundBestColumn = lower[bestColumn];
    590758                solver->setColLower(bestColumn, ceil(bestColumnValue));
     759#ifdef DIVE_DEBUG
     760                if (priority_) {
     761                  printf("setting %d (priority %d) lower bound to %g (%g)\n",
     762                         bestColumn,priority_[back[bestColumn]].priority,ceil(bestColumnValue),bestColumnValue);
     763                }
     764#endif
    591765                whichWay=1;
    592766            }
     
    596770        int originalBestRound = bestRound;
    597771        int saveModelOptions = model_->specialOptions();
    598        
    599772        while (1) {
    600773
    601774            model_->setSpecialOptions(saveModelOptions | 2048);
    602775            solver->resolve();
     776            numberSimplexIterations += solver->getIterationCount();
     777#if DIVE_PRINT>1
     778            int numberFractionalVariables = 0;
     779            double sumFractionalVariables=0.0;
     780            int numberFixed=0;
     781            for (int i = 0; i < numberIntegers; i++) {
     782              int iColumn = integerVariable[i];
     783              double value = newSolution[iColumn];
     784              double away = fabs(floor(value + 0.5) - value);
     785              if (away > integerTolerance) {
     786                numberFractionalVariables++;
     787                sumFractionalVariables += away;
     788              }
     789              if (upper[iColumn]==lower[iColumn])
     790                  numberFixed++;
     791            }
     792            printf("pass %d obj %g %s its %d total %d fixed %d +(%d,%d)",iteration,
     793                   solver->getObjValue(),solver->isProvenOptimal() ? "opt" : "infeasible",
     794                   solver->getIterationCount(),
     795                   solver->getIterationCount()+numberSimplexIterations,
     796                   numberFixed,
     797                   numberReducedCostFixed,
     798                   numberAtBoundFixed-numberGeneralFixed);
     799            if (solver->isProvenOptimal()) {
     800              printf(" - %d at bound, %d away (sum %g)\n",
     801                     numberIntegers-numberFixed-numberFractionalVariables,
     802                     numberFractionalVariables,sumFractionalVariables);
     803            } else {
     804              printf("\n");
     805              if (fixGeneralIntegers) {
     806                int digit = maxIterations_%10;
     807                if (digit==1) {
     808                  // switch off for now
     809                  switches_=saveSwitches;
     810                  fixGeneralIntegers=false;
     811                } else if (digit==2) {
     812                  // switch off always
     813                  switches_=saveSwitches;
     814                  fixGeneralIntegers=false;
     815                  maxIterations_ -= digit;
     816                }
     817              }
     818            }
     819#else
     820            if (!solver->isProvenOptimal()) {
     821              if (fixGeneralIntegers) {
     822                int digit = maxIterations_%10;
     823                if (digit==1) {
     824                  // switch off for now
     825                  switches_=saveSwitches;
     826                  fixGeneralIntegers=false;
     827                } else if (digit==2) {
     828                  // switch off always
     829                  switches_=saveSwitches;
     830                  fixGeneralIntegers=false;
     831                  maxIterations_ -= digit;
     832                }
     833              }
     834            }
     835#endif
    603836            model_->setSpecialOptions(saveModelOptions);
    604837            if (!solver->isAbandoned()&&!solver->isIterationLimitReached()) {
    605                 numberSimplexIterations += solver->getIterationCount();
     838              //numberSimplexIterations += solver->getIterationCount();
    606839            } else {
    607840                numberSimplexIterations = maxSimplexIterations + 1;
     
    672905            reasonToStop += 4;
    673906            // also switch off
    674 #ifdef CLP_INVESTIGATE
     907#if DIVE_PRINT
    675908            printf("switching off diving as too many iterations %d, %d allowed\n",
    676909                   numberSimplexIterations, maxSimplexIterations);
    677910#endif
    678911            when_ = 0;
    679         } else if (solver->getIterationCount() > 1000 && iteration > 3 && !nodes) {
     912        } else if (solver->getIterationCount() > maxIterationsInOneSolve && iteration > 3 && !nodes) {
    680913            reasonToStop += 5;
    681914            // also switch off
    682 #ifdef CLP_INVESTIGATE
     915#if DIVE_PRINT
    683916            printf("switching off diving one iteration took %d iterations (total %d)\n",
    684917                   solver->getIterationCount(), numberSimplexIterations);
     
    750983        // was good at start! - create fake
    751984        clpSolver->resolve();
     985        numberSimplexIterations += clpSolver->getIterationCount();
    752986        ClpSimplex * simplex = clpSolver->getModelPtr();
    753987        CbcSubProblem * sub =
     
    8461080    }
    8471081
    848 #ifdef DIVE_DEBUG
    849     std::cout << "nRoundInfeasible = " << nRoundInfeasible
     1082#if DIVE_PRINT
     1083    std::cout << heuristicName_
     1084              << " nRoundInfeasible = " << nRoundInfeasible
    8501085              << ", nRoundFeasible = " << nRoundFeasible
    8511086              << ", returnCode = " << returnCode
     
    8651100    upArray_ = NULL;
    8661101    delete solver;
     1102    switches_ = saveSwitches;
    8671103    return returnCode;
    8681104}
     
    8771113    if (feasibilityPumpOptions_>0 && (nodeCount % feasibilityPumpOptions_) != 0)
    8781114        return 0;
    879 #ifdef DIVE_DEBUG
    880     std::cout << "solutionValue = " << solutionValue << std::endl;
    881 #endif
    8821115    ++numCouldRun_;
    8831116
     
    8921125        return 0; // switched off
    8931126#endif
     1127#ifdef HEURISTIC_INFORM
     1128    printf("Entering heuristic %s - nRuns %d numCould %d when %d\n",
     1129           heuristicName(),numRuns_,numCouldRun_,when_);
     1130#endif
     1131#ifdef DIVE_DEBUG
     1132    std::cout << "solutionValue = " << solutionValue << std::endl;
     1133#endif
    8941134    // Get solution array for heuristic solution
    8951135    int numberColumns = model_->solver()->getNumCols();
    896     double * newSolution = new double [numberColumns];
     1136    double * newSolution = CoinCopyOfArray(model_->solver()->getColSolution(),
     1137                                           numberColumns);
    8971138    int numberCuts=0;
    8981139    int numberNodes=-1;
Note: See TracChangeset for help on using the changeset viewer.