Ignore:
File:
1 edited

Legend:

Unmodified
Added
Removed
  • stable/2.8/Cbc/src/CbcHeuristicDive.cpp

    • Property svn:eol-style set to native
    • Property svn:keywords set to Author Date Id Revision
    r1675 r1902  
    1 /* $Id: CbcHeuristicDive.cpp 1240 2009-10-02 18:41:44Z forrest $ */
     1/* $Id$ */
    22// Copyright (C) 2008, International Business Machines
    33// Corporation and others.  All Rights Reserved.
     
    1111#include "CbcHeuristicDive.hpp"
    1212#include "CbcStrategy.hpp"
     13#include "CbcModel.hpp"
     14#include "CbcSubProblem.hpp"
    1315#include "OsiAuxInfo.hpp"
    1416#include  "CoinTime.hpp"
     
    190192}
    191193
    192 // See if diving will give better solution
    193 // Sets value of solution
    194 // Returns 1 if solution, 0 if not
    195 int
    196 CbcHeuristicDive::solution(double & solutionValue,
    197                            double * betterSolution)
    198 {
    199     int nodeCount = model_->getNodeCount();
    200     if (feasibilityPumpOptions_>0 && (nodeCount % feasibilityPumpOptions_) != 0)
    201         return 0;
    202 #ifdef DIVE_DEBUG
    203     std::cout << "solutionValue = " << solutionValue << std::endl;
    204 #endif
    205     ++numCouldRun_;
    206 
    207     // test if the heuristic can run
    208     if (!canHeuristicRun())
    209         return 0;
    210 
    211 #ifdef JJF_ZERO
    212     // See if to do
    213     if (!when() || (when() % 10 == 1 && model_->phase() != 1) ||
    214             (when() % 10 == 2 && (model_->phase() != 2 && model_->phase() != 3)))
    215         return 0; // switched off
    216 #endif
    217 
     194// inner part of dive
     195int
     196CbcHeuristicDive::solution(double & solutionValue, int & numberNodes,
     197                           int & numberCuts, OsiRowCut ** cuts,
     198                           CbcSubProblem ** & nodes,
     199                           double * newSolution)
     200{
    218201#ifdef DIVE_DEBUG
    219202    int nRoundInfeasible = 0;
    220203    int nRoundFeasible = 0;
     204#endif
    221205    int reasonToStop = 0;
    222 #endif
    223206    double time1 = CoinCpuTime();
    224207    int numberSimplexIterations = 0;
    225208    int maxSimplexIterations = (model_->getNodeCount()) ? maxSimplexIterations_
    226209                               : maxSimplexIterationsAtRoot_;
    227 
     210    // but can't be exactly coin_int_max
     211    maxSimplexIterations = CoinMin(maxSimplexIterations,COIN_INT_MAX>>3);
    228212    OsiSolverInterface * solver = cloneBut(6); // was model_->solver()->clone();
    229213# ifdef COIN_HAS_CLP
     
    231215    = dynamic_cast<OsiClpSolverInterface *> (solver);
    232216    if (clpSolver) {
     217      ClpSimplex * clpSimplex = clpSolver->getModelPtr();
     218      int oneSolveIts = clpSimplex->maximumIterations();
     219      oneSolveIts = CoinMin(1000+2*(clpSimplex->numberRows()+clpSimplex->numberColumns()),oneSolveIts);
     220      clpSimplex->setMaximumIterations(oneSolveIts);
     221      if (!nodes) {
    233222        // say give up easily
    234         ClpSimplex * clpSimplex = clpSolver->getModelPtr();
    235223        clpSimplex->setMoreSpecialOptions(clpSimplex->moreSpecialOptions() | 64);
     224      } else {
     225        // get ray
     226        int specialOptions = clpSimplex->specialOptions();
     227        specialOptions &= ~0x3100000;
     228        specialOptions |= 32;
     229        clpSimplex->setSpecialOptions(specialOptions);
     230        clpSolver->setSpecialOptions(clpSolver->specialOptions() | 1048576);
     231        if ((model_->moreSpecialOptions()&16777216)!=0) {
     232          // cutoff is constraint
     233          clpSolver->setDblParam(OsiDualObjectiveLimit, COIN_DBL_MAX);
     234        }
     235      }
    236236    }
    237237# endif
     
    268268    // Get solution array for heuristic solution
    269269    int numberColumns = solver->getNumCols();
    270     double * newSolution = new double [numberColumns];
    271270    memcpy(newSolution, solution, numberColumns*sizeof(double));
    272271
    273272    // vectors to store the latest variables fixed at their bounds
    274273    int* columnFixed = new int [numberIntegers];
    275     double* originalBound = new double [numberIntegers];
     274    double* originalBound = new double [numberIntegers+2*numberColumns];
     275    double * lowerBefore = originalBound+numberIntegers;
     276    double * upperBefore = lowerBefore+numberColumns;
     277    memcpy(lowerBefore,lower,numberColumns*sizeof(double));
     278    memcpy(upperBefore,upper,numberColumns*sizeof(double));
     279    double * lastDjs=newSolution+numberColumns;
    276280    bool * fixedAtLowerBound = new bool [numberIntegers];
    277281    PseudoReducedCost * candidate = new PseudoReducedCost [numberIntegers];
     
    279283
    280284    int maxNumberAtBoundToFix = static_cast<int> (floor(percentageToFix_ * numberIntegers));
     285    assert (!maxNumberAtBoundToFix||!nodes);
    281286
    282287    // count how many fractional variables
     
    308313        bool canRound = selectVariableToBranch(solver, newSolution,
    309314                                               bestColumn, bestRound);
    310 
    311315        // if the solution is not trivially roundable, we don't try to round;
    312316        // if the solution is trivially roundable, we try to round. However,
     
    340344                nRoundFeasible++;
    341345#endif
    342                 // Round all the fractional variables
    343                 for (int i = 0; i < numberIntegers; i++) {
     346                if (!nodes||bestColumn<0) {
     347                  // Round all the fractional variables
     348                  for (int i = 0; i < numberIntegers; i++) {
    344349                    int iColumn = integerVariable[i];
    345350                    double value = newSolution[iColumn];
    346351                    if (fabs(floor(value + 0.5) - value) > integerTolerance) {
    347                         assert(downLocks_[i] == 0 || upLocks_[i] == 0);
    348                         if (downLocks_[i] == 0 && upLocks_[i] == 0) {
    349                             if (direction * objective[iColumn] >= 0.0)
    350                                 newSolution[iColumn] = floor(value);
    351                             else
    352                                 newSolution[iColumn] = ceil(value);
    353                         } else if (downLocks_[i] == 0)
    354                             newSolution[iColumn] = floor(value);
    355                         else
    356                             newSolution[iColumn] = ceil(value);
     352                      assert(downLocks_[i] == 0 || upLocks_[i] == 0);
     353                      if (downLocks_[i] == 0 && upLocks_[i] == 0) {
     354                        if (direction * objective[iColumn] >= 0.0)
     355                          newSolution[iColumn] = floor(value);
     356                        else
     357                          newSolution[iColumn] = ceil(value);
     358                      } else if (downLocks_[i] == 0)
     359                        newSolution[iColumn] = floor(value);
     360                      else
     361                        newSolution[iColumn] = ceil(value);
    357362                    }
    358                 }
    359                 break;
    360             }
     363                  }
     364                  break;
     365                } else {
     366                  // can't round if going to use in branching
     367                  int i;
     368                  for (i = 0; i < numberIntegers; i++) {
     369                    int iColumn = integerVariable[i];
     370                    double value = newSolution[bestColumn];
     371                    if (fabs(floor(value + 0.5) - value) > integerTolerance) {
     372                      if (iColumn==bestColumn) {
     373                        assert(downLocks_[i] == 0 || upLocks_[i] == 0);
     374                        double obj = objective[bestColumn];
     375                        if (downLocks_[i] == 0 && upLocks_[i] == 0) {
     376                          if (direction * obj >= 0.0)
     377                            bestRound=-1;
     378                          else
     379                            bestRound=1;
     380                        } else if (downLocks_[i] == 0)
     381                          bestRound=-1;
     382                        else
     383                          bestRound=1;
     384                        break;
     385                      }
     386                    }
     387                  }
     388                }
     389            }
    361390#ifdef DIVE_DEBUG
    362391            else
     
    549578
    550579        double originalBoundBestColumn;
     580        double bestColumnValue;
     581        int whichWay;
    551582        if (bestColumn >= 0) {
     583            bestColumnValue = newSolution[bestColumn];
    552584            if (bestRound < 0) {
    553585                originalBoundBestColumn = upper[bestColumn];
    554                 solver->setColUpper(bestColumn, floor(newSolution[bestColumn]));
     586                solver->setColUpper(bestColumn, floor(bestColumnValue));
     587                whichWay=0;
    555588            } else {
    556589                originalBoundBestColumn = lower[bestColumn];
    557                 solver->setColLower(bestColumn, ceil(newSolution[bestColumn]));
     590                solver->setColLower(bestColumn, ceil(bestColumnValue));
     591                whichWay=1;
    558592            }
    559593        } else {
     
    562596        int originalBestRound = bestRound;
    563597        int saveModelOptions = model_->specialOptions();
     598       
    564599        while (1) {
    565600
     
    567602            solver->resolve();
    568603            model_->setSpecialOptions(saveModelOptions);
    569             if (!solver->isAbandoned()) {
     604            if (!solver->isAbandoned()&&!solver->isIterationLimitReached()) {
    570605                numberSimplexIterations += solver->getIterationCount();
    571606            } else {
    572607                numberSimplexIterations = maxSimplexIterations + 1;
     608                reasonToStop += 100;
    573609                break;
    574610            }
    575611
    576612            if (!solver->isProvenOptimal()) {
     613                if (nodes) {
     614                  if (solver->isProvenPrimalInfeasible()) {
     615                    if (maxSimplexIterationsAtRoot_!=COIN_INT_MAX) {
     616                      // stop now
     617                      printf("stopping on first infeasibility\n");
     618                      break;
     619                    } else if (cuts) {
     620                      // can do conflict cut
     621                      printf("could do intermediate conflict cut\n");
     622                      bool localCut;
     623                      OsiRowCut * cut = model_->conflictCut(solver,localCut);
     624                      if (cut) {
     625                        if (!localCut) {
     626                          model_->makePartialCut(cut,solver);
     627                          cuts[numberCuts++]=cut;
     628                        } else {
     629                          delete cut;
     630                        }
     631                      }
     632                    }
     633                  } else {
     634                    reasonToStop += 10;
     635                    break;
     636                  }
     637                }
    577638                if (numberAtBoundFixed > 0) {
    578639                    // Remove the bound fix for variables that were at bounds
     
    587648                } else if (bestRound == originalBestRound) {
    588649                    bestRound *= (-1);
     650                    whichWay |=2;
    589651                    if (bestRound < 0) {
    590652                        solver->setColLower(bestColumn, originalBoundBestColumn);
    591                         solver->setColUpper(bestColumn, floor(newSolution[bestColumn]));
     653                        solver->setColUpper(bestColumn, floor(bestColumnValue));
    592654                    } else {
    593                         solver->setColLower(bestColumn, ceil(newSolution[bestColumn]));
     655                        solver->setColLower(bestColumn, ceil(bestColumnValue));
    594656                        solver->setColUpper(bestColumn, originalBoundBestColumn);
    595657                    }
     
    602664        if (!solver->isProvenOptimal() ||
    603665                direction*solver->getObjValue() >= solutionValue) {
    604 #ifdef DIVE_DEBUG
    605             reasonToStop = 1;
    606 #endif
    607             break;
    608         }
    609 
    610         if (iteration > maxIterations_) {
    611 #ifdef DIVE_DEBUG
    612             reasonToStop = 2;
    613 #endif
    614             break;
    615         }
    616 
    617         if (CoinCpuTime() - time1 > maxTime_) {
    618 #ifdef DIVE_DEBUG
    619             reasonToStop = 3;
    620 #endif
    621             break;
    622         }
    623 
    624         if (numberSimplexIterations > maxSimplexIterations) {
    625 #ifdef DIVE_DEBUG
    626             reasonToStop = 4;
    627 #endif
     666            reasonToStop += 1;
     667        } else if (iteration > maxIterations_) {
     668            reasonToStop += 2;
     669        } else if (CoinCpuTime() - time1 > maxTime_) {
     670            reasonToStop += 3;
     671        } else if (numberSimplexIterations > maxSimplexIterations) {
     672            reasonToStop += 4;
    628673            // also switch off
    629674#ifdef CLP_INVESTIGATE
     
    632677#endif
    633678            when_ = 0;
    634             break;
    635         }
    636 
    637         if (solver->getIterationCount() > 1000 && iteration > 3) {
    638 #ifdef DIVE_DEBUG
    639             reasonToStop = 5;
    640 #endif
     679        } else if (solver->getIterationCount() > 1000 && iteration > 3 && !nodes) {
     680            reasonToStop += 5;
    641681            // also switch off
    642682#ifdef CLP_INVESTIGATE
     
    645685#endif
    646686            when_ = 0;
    647             break;
    648687        }
    649688
    650689        memcpy(newSolution, solution, numberColumns*sizeof(double));
    651690        numberFractionalVariables = 0;
     691        double sumFractionalVariables=0.0;
    652692        for (int i = 0; i < numberIntegers; i++) {
    653693            int iColumn = integerVariable[i];
    654694            double value = newSolution[iColumn];
    655             if (fabs(floor(value + 0.5) - value) > integerTolerance) {
     695            double away = fabs(floor(value + 0.5) - value);
     696            if (away > integerTolerance) {
    656697                numberFractionalVariables++;
    657             }
    658         }
    659 
    660     }
    661 
    662 
    663     double * rowActivity = new double[numberRows];
    664     memset(rowActivity, 0, numberRows*sizeof(double));
    665 
     698                sumFractionalVariables += away;
     699            }
     700        }
     701        if (nodes) {
     702          // save information
     703          //branchValues[numberNodes]=bestColumnValue;
     704          //statuses[numberNodes]=whichWay+(bestColumn<<2);
     705          //bases[numberNodes]=solver->getWarmStart();
     706          ClpSimplex * simplex = clpSolver->getModelPtr();
     707          CbcSubProblem * sub =
     708            new CbcSubProblem(clpSolver,lowerBefore,upperBefore,
     709                          simplex->statusArray(),numberNodes);
     710          nodes[numberNodes]=sub;
     711          // other stuff
     712          sub->branchValue_=bestColumnValue;
     713          sub->problemStatus_=whichWay;
     714          sub->branchVariable_=bestColumn;
     715          sub->objectiveValue_ = simplex->objectiveValue();
     716          sub->sumInfeasibilities_ = sumFractionalVariables;
     717          sub->numberInfeasibilities_ = numberFractionalVariables;
     718          printf("DiveNode %d column %d way %d bvalue %g obj %g\n",
     719                 numberNodes,sub->branchVariable_,sub->problemStatus_,
     720                 sub->branchValue_,sub->objectiveValue_);
     721          numberNodes++;
     722          if (solver->isProvenOptimal()) {
     723            memcpy(lastDjs,solver->getReducedCost(),numberColumns*sizeof(double));
     724            memcpy(lowerBefore,lower,numberColumns*sizeof(double));
     725            memcpy(upperBefore,upper,numberColumns*sizeof(double));
     726          }
     727        }
     728        if (!numberFractionalVariables||reasonToStop)
     729          break;
     730    }
     731    if (nodes) {
     732      printf("Exiting dive for reason %d\n",reasonToStop);
     733      if (reasonToStop>1) {
     734        printf("problems in diving\n");
     735        int whichWay=nodes[numberNodes-1]->problemStatus_;
     736        CbcSubProblem * sub;
     737        if ((whichWay&2)==0) {
     738          // leave both ways
     739          sub = new CbcSubProblem(*nodes[numberNodes-1]);
     740          nodes[numberNodes++]=sub;
     741        } else {
     742          sub = nodes[numberNodes-1];
     743        }
     744        if ((whichWay&1)==0)
     745          sub->problemStatus_=whichWay|1;
     746        else
     747          sub->problemStatus_=whichWay&~1;
     748      }
     749      if (!numberNodes) {
     750        // was good at start! - create fake
     751        clpSolver->resolve();
     752        ClpSimplex * simplex = clpSolver->getModelPtr();
     753        CbcSubProblem * sub =
     754          new CbcSubProblem(clpSolver,lowerBefore,upperBefore,
     755                            simplex->statusArray(),numberNodes);
     756        nodes[numberNodes]=sub;
     757        // other stuff
     758        sub->branchValue_=0.0;
     759        sub->problemStatus_=0;
     760        sub->branchVariable_=-1;
     761        sub->objectiveValue_ = simplex->objectiveValue();
     762        sub->sumInfeasibilities_ = 0.0;
     763        sub->numberInfeasibilities_ = 0;
     764        printf("DiveNode %d column %d way %d bvalue %g obj %g\n",
     765               numberNodes,sub->branchVariable_,sub->problemStatus_,
     766               sub->branchValue_,sub->objectiveValue_);
     767        numberNodes++;
     768        assert (solver->isProvenOptimal());
     769      }
     770      nodes[numberNodes-1]->problemStatus_ |= 256*reasonToStop;
     771      // use djs as well
     772      if (solver->isProvenPrimalInfeasible()&&cuts) {
     773        // can do conflict cut and re-order
     774        printf("could do final conflict cut\n");
     775        bool localCut;
     776        OsiRowCut * cut = model_->conflictCut(solver,localCut);
     777        if (cut) {
     778          printf("cut - need to use conflict and previous djs\n");
     779          if (!localCut) {
     780            model_->makePartialCut(cut,solver);
     781            cuts[numberCuts++]=cut;
     782          } else {
     783            delete cut;
     784          }
     785        } else {
     786          printf("bad conflict - just use previous djs\n");
     787        }
     788      }
     789    }
     790   
    666791    // re-compute new solution value
    667792    double objOffset = 0.0;
     
    669794    newSolutionValue = -objOffset;
    670795    for (int i = 0 ; i < numberColumns ; i++ )
    671         newSolutionValue += objective[i] * newSolution[i];
     796      newSolutionValue += objective[i] * newSolution[i];
    672797    newSolutionValue *= direction;
    673798    //printf("new solution value %g %g\n",newSolutionValue,solutionValue);
    674     if (newSolutionValue < solutionValue) {
    675         // paranoid check
    676         memset(rowActivity, 0, numberRows*sizeof(double));
    677         for (int i = 0; i < numberColumns; i++) {
    678             int j;
    679             double value = newSolution[i];
    680             if (value) {
    681                 for (j = columnStart[i];
    682                         j < columnStart[i] + columnLength[i]; j++) {
    683                     int iRow = row[j];
    684                     rowActivity[iRow] += value * element[j];
    685                 }
    686             }
    687         }
    688         // check was approximately feasible
    689         bool feasible = true;
    690         for (int i = 0; i < numberRows; i++) {
    691             if (rowActivity[i] < rowLower[i]) {
    692                 if (rowActivity[i] < rowLower[i] - 1000.0*primalTolerance)
    693                     feasible = false;
    694             } else if (rowActivity[i] > rowUpper[i]) {
    695                 if (rowActivity[i] > rowUpper[i] + 1000.0*primalTolerance)
    696                     feasible = false;
    697             }
    698         }
    699         for (int i = 0; i < numberIntegers; i++) {
    700             int iColumn = integerVariable[i];
    701             double value = newSolution[iColumn];
    702             if (fabs(floor(value + 0.5) - value) > integerTolerance) {
    703                 feasible = false;
    704                 break;
    705             }
    706         }
    707         if (feasible) {
    708             // new solution
    709             memcpy(betterSolution, newSolution, numberColumns*sizeof(double));
    710             solutionValue = newSolutionValue;
    711             //printf("** Solution of %g found by CbcHeuristicDive\n",newSolutionValue);
    712             returnCode = 1;
    713         } else {
    714             // Can easily happen
    715             //printf("Debug CbcHeuristicDive giving bad solution\n");
    716         }
     799    if (newSolutionValue < solutionValue && !reasonToStop) {
     800      double * rowActivity = new double[numberRows];
     801      memset(rowActivity, 0, numberRows*sizeof(double));
     802      // paranoid check
     803      memset(rowActivity, 0, numberRows*sizeof(double));
     804      for (int i = 0; i < numberColumns; i++) {
     805        int j;
     806        double value = newSolution[i];
     807        if (value) {
     808          for (j = columnStart[i];
     809               j < columnStart[i] + columnLength[i]; j++) {
     810            int iRow = row[j];
     811            rowActivity[iRow] += value * element[j];
     812          }
     813        }
     814      }
     815      // check was approximately feasible
     816      bool feasible = true;
     817      for (int i = 0; i < numberRows; i++) {
     818        if (rowActivity[i] < rowLower[i]) {
     819          if (rowActivity[i] < rowLower[i] - 1000.0*primalTolerance)
     820            feasible = false;
     821        } else if (rowActivity[i] > rowUpper[i]) {
     822          if (rowActivity[i] > rowUpper[i] + 1000.0*primalTolerance)
     823            feasible = false;
     824        }
     825      }
     826      for (int i = 0; i < numberIntegers; i++) {
     827        int iColumn = integerVariable[i];
     828        double value = newSolution[iColumn];
     829        if (fabs(floor(value + 0.5) - value) > integerTolerance) {
     830          feasible = false;
     831          break;
     832        }
     833      }
     834      if (feasible) {
     835        // new solution
     836        solutionValue = newSolutionValue;
     837        //printf("** Solution of %g found by CbcHeuristicDive\n",newSolutionValue);
     838        //if (cuts)
     839        //clpSolver->getModelPtr()->writeMps("good8.mps", 2);
     840        returnCode = 1;
     841      } else {
     842        // Can easily happen
     843        printf("Debug CbcHeuristicDive giving bad solution\n");
     844      }
     845      delete [] rowActivity;
    717846    }
    718847
     
    726855#endif
    727856
    728     delete [] newSolution;
    729857    delete [] columnFixed;
    730858    delete [] originalBound;
    731859    delete [] fixedAtLowerBound;
    732860    delete [] candidate;
    733     delete [] rowActivity;
    734861    delete [] random;
    735862    delete [] downArray_;
     
    739866    delete solver;
    740867    return returnCode;
     868}
     869// See if diving will give better solution
     870// Sets value of solution
     871// Returns 1 if solution, 0 if not
     872int
     873CbcHeuristicDive::solution(double & solutionValue,
     874                           double * betterSolution)
     875{
     876    int nodeCount = model_->getNodeCount();
     877    if (feasibilityPumpOptions_>0 && (nodeCount % feasibilityPumpOptions_) != 0)
     878        return 0;
     879#ifdef DIVE_DEBUG
     880    std::cout << "solutionValue = " << solutionValue << std::endl;
     881#endif
     882    ++numCouldRun_;
     883
     884    // test if the heuristic can run
     885    if (!canHeuristicRun())
     886        return 0;
     887
     888#ifdef JJF_ZERO
     889    // See if to do
     890    if (!when() || (when() % 10 == 1 && model_->phase() != 1) ||
     891            (when() % 10 == 2 && (model_->phase() != 2 && model_->phase() != 3)))
     892        return 0; // switched off
     893#endif
     894    // Get solution array for heuristic solution
     895    int numberColumns = model_->solver()->getNumCols();
     896    double * newSolution = new double [numberColumns];
     897    int numberCuts=0;
     898    int numberNodes=-1;
     899    CbcSubProblem ** nodes=NULL;
     900    int returnCode=solution(solutionValue,numberNodes,numberCuts,
     901                            NULL,nodes,
     902                            newSolution);
     903  if (returnCode==1)
     904    memcpy(betterSolution, newSolution, numberColumns*sizeof(double));
     905
     906    delete [] newSolution;
     907    return returnCode;
     908}
     909/* returns 0 if no solution, 1 if valid solution
     910   with better objective value than one passed in
     911   also returns list of nodes
     912   This does Fractional Diving
     913*/
     914int
     915CbcHeuristicDive::fathom(CbcModel * model, int & numberNodes,
     916                         CbcSubProblem ** & nodes)
     917{
     918  double solutionValue = model->getCutoff();
     919  numberNodes=0;
     920  // Get solution array for heuristic solution
     921  int numberColumns = model_->solver()->getNumCols();
     922  double * newSolution = new double [4*numberColumns];
     923  double * lastDjs = newSolution+numberColumns;
     924  double * originalLower = lastDjs+numberColumns;
     925  double * originalUpper = originalLower+numberColumns;
     926  memcpy(originalLower,model_->solver()->getColLower(),
     927         numberColumns*sizeof(double));
     928  memcpy(originalUpper,model_->solver()->getColUpper(),
     929         numberColumns*sizeof(double));
     930  int numberCuts=0;
     931  OsiRowCut ** cuts = NULL; //new OsiRowCut * [maxIterations_];
     932  nodes=new CbcSubProblem * [maxIterations_+2];
     933  int returnCode=solution(solutionValue,numberNodes,numberCuts,
     934                          cuts,nodes,
     935                          newSolution);
     936
     937  if (returnCode==1) {
     938    // copy to best solution ? or put in solver
     939    printf("Solution from heuristic fathom\n");
     940  }
     941  int numberFeasibleNodes=numberNodes;
     942  if (returnCode!=1)
     943    numberFeasibleNodes--;
     944  if (numberFeasibleNodes>0) {
     945    CoinWarmStartBasis * basis = nodes[numberFeasibleNodes-1]->status_;
     946    //double * sort = new double [numberFeasibleNodes];
     947    //int * whichNode = new int [numberFeasibleNodes];
     948    //int numberNodesNew=0;
     949    // use djs on previous unless feasible
     950    for (int iNode=0;iNode<numberFeasibleNodes;iNode++) {
     951      CbcSubProblem * sub = nodes[iNode];
     952      double branchValue = sub->branchValue_;
     953      int iStatus=sub->problemStatus_;
     954      int iColumn = sub->branchVariable_;
     955      bool secondBranch = (iStatus&2)!=0;
     956      bool branchUp;
     957      if (!secondBranch)
     958        branchUp = (iStatus&1)!=0;
     959      else
     960        branchUp = (iStatus&1)==0;
     961      double djValue=lastDjs[iColumn];
     962      sub->djValue_=fabs(djValue);
     963      if (!branchUp&&floor(branchValue)==originalLower[iColumn]
     964          &&basis->getStructStatus(iColumn) == CoinWarmStartBasis::atLowerBound) {
     965        if (djValue>0.0) {
     966          // naturally goes to LB
     967          printf("ignoring branch down on %d (node %d) from value of %g - branch was %s - dj %g\n",
     968                 iColumn,iNode,branchValue,secondBranch ? "second" : "first",
     969                 djValue);
     970          sub->problemStatus_ |= 4;
     971          //} else {
     972          // put on list
     973          //sort[numberNodesNew]=djValue;
     974          //whichNode[numberNodesNew++]=iNode;
     975        }
     976      } else if (branchUp&&ceil(branchValue)==originalUpper[iColumn]
     977          &&basis->getStructStatus(iColumn) == CoinWarmStartBasis::atUpperBound) {
     978        if (djValue<0.0) {
     979          // naturally goes to UB
     980          printf("ignoring branch up on %d (node %d) from value of %g - branch was %s - dj %g\n",
     981                 iColumn,iNode,branchValue,secondBranch ? "second" : "first",
     982                 djValue);
     983          sub->problemStatus_ |= 4;
     984          //} else {
     985          // put on list
     986          //sort[numberNodesNew]=-djValue;
     987          //whichNode[numberNodesNew++]=iNode;
     988        }
     989      }
     990    }
     991    // use conflict to order nodes
     992    for (int iCut=0;iCut<numberCuts;iCut++) {
     993    }
     994    //CoinSort_2(sort,sort+numberNodesNew,whichNode);
     995    // create nodes
     996    // last node will have one way already done
     997  }
     998  for (int iCut=0;iCut<numberCuts;iCut++) {
     999    delete cuts[iCut];
     1000  }
     1001  delete [] cuts;
     1002  delete [] newSolution;
     1003  return returnCode;
    7411004}
    7421005
Note: See TracChangeset for help on using the changeset viewer.