Ignore:
Timestamp:
Oct 8, 2006 7:33:47 PM (13 years ago)
Author:
forrest
Message:

towards common use with other solvers

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/devel/Cbc/src/CbcHeuristicFPump.cpp

    r435 r439  
    2323   startTime_(0.0),
    2424   maximumTime_(0.0),
     25   downValue_(0.5),
     26   fakeCutoff_(COIN_DBL_MAX),
     27   absoluteIncrement_(0.0),
     28   relativeIncrement_(0.0),
    2529   maximumPasses_(100),
    26    downValue_(0.5),
     30   maximumRetries_(1),
    2731   roundExpensive_(false)
    2832{
     
    3640   startTime_(0.0),
    3741   maximumTime_(0.0),
     42   downValue_(downValue),
     43   fakeCutoff_(COIN_DBL_MAX),
     44   absoluteIncrement_(0.0),
     45   relativeIncrement_(0.0),
    3846   maximumPasses_(100),
    39    downValue_(downValue),
     47   maximumRetries_(1),
    4048   roundExpensive_(roundExpensive)
    4149{
     
    6573  else
    6674    fprintf(fp,"4  heuristicFPump.setMaximumPasses(%d);\n",maximumPasses_);
     75  if (maximumRetries_!=other.maximumRetries_)
     76    fprintf(fp,"3  heuristicFPump.setMaximumRetries(%d);\n",maximumRetries_);
     77  else
     78    fprintf(fp,"4  heuristicFPump.setMaximumRetries(%d);\n",maximumRetries_);
    6779  if (maximumTime_!=other.maximumTime_)
    6880    fprintf(fp,"3  heuristicFPump.setMaximumTime(%g);\n",maximumTime_);
    6981  else
    7082    fprintf(fp,"4  heuristicFPump.setMaximumTime(%g);\n",maximumTime_);
     83  if (fakeCutoff_!=other.fakeCutoff_)
     84    fprintf(fp,"3  heuristicFPump.setFakeCutoff(%g);\n",fakeCutoff_);
     85  else
     86    fprintf(fp,"4  heuristicFPump.setFakeCutoff(%g);\n",fakeCutoff_);
     87  if (absoluteIncrement_!=other.absoluteIncrement_)
     88    fprintf(fp,"3  heuristicFPump.setAbsoluteIncrement(%g);\n",absoluteIncrement_);
     89  else
     90    fprintf(fp,"4  heuristicFPump.setAbsoluteIncrement(%g);\n",absoluteIncrement_);
     91  if (relativeIncrement_!=other.relativeIncrement_)
     92    fprintf(fp,"3  heuristicFPump.setRelativeIncrement(%g);\n",relativeIncrement_);
     93  else
     94    fprintf(fp,"4  heuristicFPump.setRelativeIncrement(%g);\n",relativeIncrement_);
    7195  fprintf(fp,"3  cbcModel->addHeuristic(&heuristicFPump);\n");
    7296}
     
    78102  startTime_(rhs.startTime_),
    79103  maximumTime_(rhs.maximumTime_),
     104  downValue_(rhs.downValue_),
     105  fakeCutoff_(rhs.fakeCutoff_),
     106  absoluteIncrement_(rhs.absoluteIncrement_),
     107  relativeIncrement_(rhs.relativeIncrement_),
    80108  maximumPasses_(rhs.maximumPasses_),
    81   downValue_(rhs.downValue_),
     109  maximumRetries_(rhs.maximumRetries_),
    82110  roundExpensive_(rhs.roundExpensive_)
    83111{
     
    109137  // probably a good idea
    110138  if (model_->getSolutionCount()) return 0;
    111   // Clone solver - otherwise annoys root node computations
    112   OsiSolverInterface * solver = model_->solver()->clone();
    113   solver->setDblParam(OsiDualObjectiveLimit,1.0e50);
    114   solver->resolve();
    115   const double * lower = solver->getColLower();
    116   const double * upper = solver->getColUpper();
    117   const double * solution = solver->getColSolution();
    118   double primalTolerance;
    119   solver->getDblParam(OsiPrimalTolerance,primalTolerance);
    120  
     139  // loop round doing repeated pumps
     140  double cutoff;
     141  model_->solver()->getDblParam(OsiDualObjectiveLimit,cutoff);
     142  double direction = model_->solver()->getObjSense();
     143  cutoff *= direction;
     144  cutoff = CoinMin(cutoff,solutionValue);
     145  // space for rounded solution
    121146  int numberColumns = model_->getNumCols();
     147  double * newSolution = new double [numberColumns];
     148  double newSolutionValue=COIN_DBL_MAX;
     149  bool solutionFound=false;
    122150  char * usedColumn = NULL;
    123151  double * lastSolution=NULL;
     
    130158    usedColumn = new char [numberColumns];
    131159    memset(usedColumn,0,numberColumns);
    132     lastSolution = CoinCopyOfArray(solver->getColSolution(),numberColumns);
     160    model_->solver()->resolve();
     161    assert (model_->solver()->isProvenOptimal());
     162    lastSolution = CoinCopyOfArray(model_->solver()->getColSolution(),numberColumns);
    133163  }
    134   int numberIntegers = model_->numberIntegers();
    135   const int * integerVariable = model_->integerVariable();
    136 
    137 // 1. initially check 0-1
    138   int i,j;
    139   int general=0;
    140   for (i=0;i<numberIntegers;i++) {
    141     int iColumn = integerVariable[i];
     164  int finalReturnCode=0;
     165  int totalNumberPasses=0;
     166  int numberTries=0;
     167  while (true) {
     168    int numberPasses=0;
     169    numberTries++;
     170    // Clone solver - otherwise annoys root node computations
     171    OsiSolverInterface * solver = model_->solver()->clone();
     172    // if cutoff exists then add constraint
     173    if (fabs(cutoff)<1.0e20&&(fakeCutoff_!=COIN_DBL_MAX||numberTries>1)) {
     174      cutoff = CoinMin(cutoff,fakeCutoff_);
     175      const double * objective = solver->getObjCoefficients();
     176      int numberColumns = solver->getNumCols();
     177      int * which = new int[numberColumns];
     178      double * els = new double[numberColumns];
     179      int nel=0;
     180      for (int i=0;i<numberColumns;i++) {
     181        double value = objective[i];
     182        if (value) {
     183          which[nel]=i;
     184          els[nel++] = direction*value;
     185        }
     186      }
     187      solver->addRow(nel,which,els,-COIN_DBL_MAX,cutoff);
     188      delete [] which;
     189      delete [] els;
     190    }
     191    solver->setDblParam(OsiDualObjectiveLimit,1.0e50);
     192    solver->resolve();
     193    // Solver may not be feasible
     194    if (!solver->isProvenOptimal()) {
     195      break;
     196    }
     197    const double * lower = solver->getColLower();
     198    const double * upper = solver->getColUpper();
     199    const double * solution = solver->getColSolution();
     200    if (lastSolution)
     201      memcpy(lastSolution,solution,numberColumns*sizeof(double));
     202    double primalTolerance;
     203    solver->getDblParam(OsiPrimalTolerance,primalTolerance);
     204   
     205    int numberIntegers = model_->numberIntegers();
     206    const int * integerVariable = model_->integerVariable();
     207
     208    // 1. initially check 0-1
     209    int i,j;
     210    int general=0;
     211    for (i=0;i<numberIntegers;i++) {
     212      int iColumn = integerVariable[i];
    142213#ifndef NDEBUG
    143     const CbcObject * object = model_->object(i);
    144     const CbcSimpleInteger * integerObject =
    145       dynamic_cast<const  CbcSimpleInteger *> (object);
    146     assert(integerObject);
     214      const OsiObject * object = model_->object(i);
     215      const CbcSimpleInteger * integerObject =
     216        dynamic_cast<const  CbcSimpleInteger *> (object);
     217      const OsiSimpleInteger * integerObject2 =
     218        dynamic_cast<const  OsiSimpleInteger *> (object);
     219      assert(integerObject||integerObject2);
    147220#endif
    148     if (upper[iColumn]-lower[iColumn]>1.000001) {
    149       general++;
    150       break;
    151     }
    152   }
    153   if (general*3>numberIntegers) {
    154     delete solver;
    155     return 0;
    156   }
    157 
    158 // 2. space for rounded solution
    159   double * newSolution = new double [numberColumns];
    160   // space for last rounded solutions
     221      if (upper[iColumn]-lower[iColumn]>1.000001) {
     222        general++;
     223        break;
     224      }
     225    }
     226    if (general*3>numberIntegers) {
     227      delete solver;
     228      return 0;
     229    }
     230   
     231    // 2 space for last rounded solutions
    161232#define NUMBER_OLD 4
    162   double ** oldSolution = new double * [NUMBER_OLD];
    163   for (j=0;j<NUMBER_OLD;j++) {
    164     oldSolution[j]= new double[numberColumns];
    165     for (i=0;i<numberColumns;i++) oldSolution[j][i]=-COIN_DBL_MAX;
    166   }
    167 
    168 // 3. Replace objective with an initial 0-valued objective
    169   double * saveObjective = new double [numberColumns];
    170   memcpy(saveObjective,solver->getObjCoefficients(),numberColumns*sizeof(double));
    171   for (i=0;i<numberColumns;i++) {
    172     solver->setObjCoeff(i,0.0);
    173   }
    174   bool finished=false;
    175   double direction = solver->getObjSense();
    176   int returnCode=0;
    177   bool takeHint;
    178   OsiHintStrength strength;
    179   solver->getHintParam(OsiDoDualInResolve,takeHint,strength);
    180   solver->setHintParam(OsiDoDualInResolve,false);
    181   solver->messageHandler()->setLogLevel(0);
    182 
    183 // 4. Save objective offset so we can see progress
    184   double saveOffset;
    185   solver->getDblParam(OsiObjOffset,saveOffset);
    186 
    187 // 5. MAIN WHILE LOOP
    188   int numberPasses=0;
    189   double newSolutionValue=COIN_DBL_MAX;
    190   bool newLineNeeded=false;
    191   while (!finished) {
    192     returnCode=0;
    193     // see what changed
    194     if (usedColumn) {
    195       const double * solution = solver->getColSolution();
    196       for (i=0;i<numberColumns;i++) {
    197         if (fabs(solution[i]-lastSolution[i])>1.0e-8)
    198           usedColumn[i]=1;
    199         lastSolution[i]=solution[i];
    200       }
    201     }
    202     if (numberPasses>=maximumPasses_) {
    203       break;
    204     }
    205     if (maximumTime_>0.0&&CoinCpuTime()>=startTime_+maximumTime_) break;
    206     numberPasses++;
    207     memcpy(newSolution,solution,numberColumns*sizeof(double));
    208     int flip;
    209     returnCode = rounds(newSolution,saveObjective,roundExpensive_,downValue_,&flip);
    210     if (returnCode) {
    211       // SOLUTION IS INTEGER
    212       // Put back correct objective
    213       for (i=0;i<numberColumns;i++)
    214         solver->setObjCoeff(i,saveObjective[i]);
    215       // solution - but may not be better
    216       // Compute using dot product
    217       solver->setDblParam(OsiObjOffset,saveOffset);
    218       newSolutionValue = direction*solver->OsiSolverInterface::getObjValue();
    219       if (model_->logLevel())
    220         printf(" - solution found\n");
    221       newLineNeeded=false;
    222       if (newSolutionValue<solutionValue) {
    223         if (general) {
    224           int numberLeft=0;
    225           for (i=0;i<numberIntegers;i++) {
    226             int iColumn = integerVariable[i];
    227             double value = floor(newSolution[iColumn]+0.5);
    228             if(solver->isBinary(iColumn)) {
    229               solver->setColLower(iColumn,value);
    230               solver->setColUpper(iColumn,value);
    231             } else {
    232               if (fabs(value-newSolution[iColumn])>1.0e-7)
    233                 numberLeft++;
    234             }
    235           }
    236           if (numberLeft) {
    237             returnCode = smallBranchAndBound(solver,200,newSolution,newSolutionValue,
    238                                          solutionValue,"CbcHeuristicFpump");
    239           }
    240         }
    241         if (returnCode) {
    242           memcpy(betterSolution,newSolution,numberColumns*sizeof(double));
    243           solutionValue=newSolutionValue;
    244         }
     233    double ** oldSolution = new double * [NUMBER_OLD];
     234    for (j=0;j<NUMBER_OLD;j++) {
     235      oldSolution[j]= new double[numberColumns];
     236      for (i=0;i<numberColumns;i++) oldSolution[j][i]=-COIN_DBL_MAX;
     237    }
     238   
     239    // 3. Replace objective with an initial 0-valued objective
     240    double * saveObjective = new double [numberColumns];
     241    memcpy(saveObjective,solver->getObjCoefficients(),numberColumns*sizeof(double));
     242    for (i=0;i<numberColumns;i++) {
     243      solver->setObjCoeff(i,0.0);
     244    }
     245    bool finished=false;
     246    double direction = solver->getObjSense();
     247    int returnCode=0;
     248    bool takeHint;
     249    OsiHintStrength strength;
     250    solver->getHintParam(OsiDoDualInResolve,takeHint,strength);
     251    solver->setHintParam(OsiDoDualInResolve,false);
     252    solver->messageHandler()->setLogLevel(0);
     253   
     254    // 4. Save objective offset so we can see progress
     255    double saveOffset;
     256    solver->getDblParam(OsiObjOffset,saveOffset);
     257   
     258    // 5. MAIN WHILE LOOP
     259    bool newLineNeeded=false;
     260    while (!finished) {
     261      returnCode=0;
     262      // see what changed
     263      if (usedColumn) {
     264        for (i=0;i<numberColumns;i++) {
     265          if (fabs(solution[i]-lastSolution[i])>1.0e-8)
     266            usedColumn[i]=1;
     267          lastSolution[i]=solution[i];
     268        }
     269      }
     270      if (numberPasses>=maximumPasses_) {
     271        break;
     272      }
     273      if (maximumTime_>0.0&&CoinCpuTime()>=startTime_+maximumTime_) break;
     274      numberPasses++;
     275      memcpy(newSolution,solution,numberColumns*sizeof(double));
     276      int flip;
     277      returnCode = rounds(newSolution,saveObjective,roundExpensive_,downValue_,&flip);
     278      if (returnCode) {
     279        // SOLUTION IS INTEGER
     280        // Put back correct objective
     281        for (i=0;i<numberColumns;i++)
     282          solver->setObjCoeff(i,saveObjective[i]);
     283
     284        // solution - but may not be better
     285        // Compute using dot product
     286        solver->setDblParam(OsiObjOffset,saveOffset);
     287        newSolutionValue = -saveOffset;
     288        for (  i=0 ; i<numberColumns ; i++ )
     289          newSolutionValue += saveObjective[i]*newSolution[i];
     290        newSolutionValue *= direction;
     291        if (model_->logLevel())
     292          printf(" - solution found of %g",newSolutionValue);
     293        newLineNeeded=false;
     294        if (newSolutionValue<solutionValue) {
     295          double saveValue = newSolutionValue;
     296          if (general) {
     297            int numberLeft=0;
     298            for (i=0;i<numberIntegers;i++) {
     299              int iColumn = integerVariable[i];
     300              double value = floor(newSolution[iColumn]+0.5);
     301              if(solver->isBinary(iColumn)) {
     302                solver->setColLower(iColumn,value);
     303                solver->setColUpper(iColumn,value);
     304              } else {
     305                if (fabs(value-newSolution[iColumn])>1.0e-7)
     306                  numberLeft++;
     307              }
     308            }
     309            if (numberLeft) {
     310              returnCode = smallBranchAndBound(solver,200,newSolution,newSolutionValue,
     311                                               solutionValue,"CbcHeuristicFpump");
     312            }
     313          }
     314          if (returnCode) {
     315            memcpy(betterSolution,newSolution,numberColumns*sizeof(double));
     316            solutionValue=newSolutionValue;
     317            solutionFound=true;
     318            if (general&&saveValue!=newSolutionValue)
     319              printf(" - cleaned solution of %g\n",solutionValue);
     320            else
     321              printf("\n");
     322          } else {
     323          if (model_->logLevel())
     324            printf(" - not good enough after small branch and bound\n");
     325          }
     326        } else {
     327          if (model_->logLevel())
     328            printf(" - not good enough\n");
     329          newLineNeeded=false;
     330          returnCode=0;
     331        }     
     332        break;
    245333      } else {
    246         returnCode=0;
    247       }     
    248       break;
    249     } else {
    250       // SOLUTION IS not INTEGER
    251       // 1. check for loop
    252       bool matched;
    253       for (int k = NUMBER_OLD-1; k > 0; k--) {
     334        // SOLUTION IS not INTEGER
     335        // 1. check for loop
     336        bool matched;
     337        for (int k = NUMBER_OLD-1; k > 0; k--) {
    254338          double * b = oldSolution[k];
    255339          matched = true;
    256340          for (i = 0; i <numberIntegers; i++) {
    257               int iColumn = integerVariable[i];
    258               if(!solver->isBinary(iColumn))
    259                 continue;
    260               if (newSolution[iColumn]!=b[iColumn]) {
    261                 matched=false;
    262                 break;
    263               }
     341            int iColumn = integerVariable[i];
     342            if(!solver->isBinary(iColumn))
     343              continue;
     344            if (newSolution[iColumn]!=b[iColumn]) {
     345              matched=false;
     346              break;
     347            }
    264348          }
    265349          if (matched) break;
    266       }
    267       if (matched || numberPasses%100 == 0) {
    268          // perturbation
    269         if (model_->logLevel())
    270           printf("Perturbation applied");
    271          newLineNeeded=true;
    272          for (i=0;i<numberIntegers;i++) {
    273              int iColumn = integerVariable[i];
    274              if(!solver->isBinary(iColumn))
    275                continue;
    276              double value = max(0.0,CoinDrand48()-0.3);
    277              double difference = fabs(solution[iColumn]-newSolution[iColumn]);
    278              if (difference+value>0.5) {
    279                 if (newSolution[iColumn]<lower[iColumn]+primalTolerance) newSolution[iColumn] += 1.0;
    280              else if (newSolution[iColumn]>upper[iColumn]-primalTolerance) newSolution[iColumn] -= 1.0;
    281                   else abort();
    282              }
    283          }
    284       } else {
    285          for (j=NUMBER_OLD-1;j>0;j--) {
    286              for (i = 0; i < numberColumns; i++) oldSolution[j][i]=oldSolution[j-1][i];
    287          }
    288          for (j = 0; j < numberColumns; j++) oldSolution[0][j] = newSolution[j];
    289       }
    290 
    291       // 2. update the objective function based on the new rounded solution
    292       double offset=0.0;
     350        }
     351        if (matched || numberPasses%100 == 0) {
     352          // perturbation
     353          if (model_->logLevel())
     354            printf("Perturbation applied");
     355          newLineNeeded=true;
     356          for (i=0;i<numberIntegers;i++) {
     357            int iColumn = integerVariable[i];
     358            if(!solver->isBinary(iColumn))
     359              continue;
     360            double value = max(0.0,CoinDrand48()-0.3);
     361            double difference = fabs(solution[iColumn]-newSolution[iColumn]);
     362            if (difference+value>0.5) {
     363              if (newSolution[iColumn]<lower[iColumn]+primalTolerance) newSolution[iColumn] += 1.0;
     364              else if (newSolution[iColumn]>upper[iColumn]-primalTolerance) newSolution[iColumn] -= 1.0;
     365              else abort();
     366            }
     367          }
     368        } else {
     369          for (j=NUMBER_OLD-1;j>0;j--) {
     370            for (i = 0; i < numberColumns; i++) oldSolution[j][i]=oldSolution[j-1][i];
     371          }
     372          for (j = 0; j < numberColumns; j++) oldSolution[0][j] = newSolution[j];
     373        }
     374       
     375        // 2. update the objective function based on the new rounded solution
     376        double offset=0.0;
     377        for (i=0;i<numberIntegers;i++) {
     378          int iColumn = integerVariable[i];
     379          if(!solver->isBinary(iColumn))
     380            continue;
     381          double costValue = 1.0;
     382          // deal with fixed variables (i.e., upper=lower)
     383          if (fabs(lower[iColumn]-upper[iColumn]) < primalTolerance) {
     384            if (lower[iColumn] > 1. - primalTolerance) solver->setObjCoeff(iColumn,-costValue);
     385            else                                       solver->setObjCoeff(iColumn,costValue);
     386            continue;
     387          }
     388          if (newSolution[iColumn]<lower[iColumn]+primalTolerance) {
     389            solver->setObjCoeff(iColumn,costValue);
     390          } else {
     391            if (newSolution[iColumn]>upper[iColumn]-primalTolerance) {
     392              solver->setObjCoeff(iColumn,-costValue);
     393            } else {
     394              abort();
     395            }
     396          }
     397          offset += costValue*newSolution[iColumn];
     398        }
     399        solver->setDblParam(OsiObjOffset,-offset);
     400        if (!general&false) {
     401          // Solve in two goes - first keep satisfied ones fixed
     402          double * saveLower = new double [numberIntegers];
     403          double * saveUpper = new double [numberIntegers];
     404          for (i=0;i<numberIntegers;i++) {
     405            int iColumn = integerVariable[i];
     406            saveLower[i]=COIN_DBL_MAX;
     407            saveUpper[i]=-COIN_DBL_MAX;
     408            if (solution[iColumn]<lower[iColumn]+primalTolerance) {
     409              saveUpper[i]=upper[iColumn];
     410              solver->setColUpper(iColumn,lower[iColumn]);
     411            } else if (solution[iColumn]>upper[iColumn]-primalTolerance) {
     412              saveLower[i]=lower[iColumn];
     413              solver->setColLower(iColumn,upper[iColumn]);
     414            }
     415          }
     416          solver->resolve();
     417          assert (solver->isProvenOptimal());
     418          for (i=0;i<numberIntegers;i++) {
     419            int iColumn = integerVariable[i];
     420            if (saveLower[i]!=COIN_DBL_MAX)
     421              solver->setColLower(iColumn,saveLower[i]);
     422            if (saveUpper[i]!=-COIN_DBL_MAX)
     423              solver->setColUpper(iColumn,saveUpper[i]);
     424            saveUpper[i]=-COIN_DBL_MAX;
     425          }
     426          memcpy(newSolution,solution,numberColumns*sizeof(double));
     427          int flip;
     428          returnCode = rounds(newSolution,saveObjective,roundExpensive_,downValue_,&flip);
     429          if (returnCode) {
     430            // solution - but may not be better
     431            // Compute using dot product
     432            double newSolutionValue = -saveOffset;
     433            for (  i=0 ; i<numberColumns ; i++ )
     434              newSolutionValue += saveObjective[i]*newSolution[i];
     435            newSolutionValue *= direction;
     436            if (model_->logLevel())
     437              printf(" - intermediate solution found of %g",newSolutionValue);
     438            if (newSolutionValue<solutionValue) {
     439              memcpy(betterSolution,newSolution,numberColumns*sizeof(double));
     440              solutionValue=newSolutionValue;
     441              solutionFound=true;
     442            } else {
     443              returnCode=0;
     444            }
     445          }     
     446        }
     447        solver->resolve();
     448        assert (solver->isProvenOptimal());
     449        if (model_->logLevel())
     450          printf("\npass %3d: obj. %10.5f --> ", numberPasses+totalNumberPasses,solver->getObjValue());
     451        newLineNeeded=true;
     452       
     453      }
     454    } // END WHILE
     455    if (newLineNeeded&&model_->logLevel())
     456      printf(" - no solution found\n");
     457    delete solver;
     458    for ( j=0;j<NUMBER_OLD;j++)
     459      delete [] oldSolution[j];
     460    delete [] oldSolution;
     461    delete [] saveObjective;
     462    if (usedColumn) {
     463      OsiSolverInterface * newSolver = model_->continuousSolver()->clone();
     464      const double * colLower = newSolver->getColLower();
     465      const double * colUpper = newSolver->getColUpper();
     466      int i;
     467      int nFix=0;
     468      int nFixI=0;
     469      int nFixC=0;
    293470      for (i=0;i<numberIntegers;i++) {
    294         int iColumn = integerVariable[i];
    295         if(!solver->isBinary(iColumn))
    296           continue;
    297         double costValue = 1.0;
    298         // deal with fixed variables (i.e., upper=lower)
    299         if (fabs(lower[iColumn]-upper[iColumn]) < primalTolerance) {
    300            if (lower[iColumn] > 1. - primalTolerance) solver->setObjCoeff(iColumn,-costValue);
    301            else                                       solver->setObjCoeff(iColumn,costValue);
    302            continue;
    303         }
    304         if (newSolution[iColumn]<lower[iColumn]+primalTolerance) {
    305           solver->setObjCoeff(iColumn,costValue);
    306         } else {
    307           if (newSolution[iColumn]>upper[iColumn]-primalTolerance) {
    308             solver->setObjCoeff(iColumn,-costValue);
    309           } else {
    310             abort();
    311           }
    312         }
    313         offset += costValue*newSolution[iColumn];
    314       }
    315       solver->setDblParam(OsiObjOffset,-offset);
    316       solver->resolve();
    317       if (model_->logLevel())
    318         printf("\npass %3d: obj. %10.5f --> ", numberPasses,solver->getObjValue());
    319       newLineNeeded=true;
    320 
    321     }
    322   } // END WHILE
    323   if (newLineNeeded&&model_->logLevel())
    324     printf(" - no solution found\n");
    325   delete solver;
    326   for ( j=0;j<NUMBER_OLD;j++)
    327     delete [] oldSolution[j];
    328   delete [] oldSolution;
    329   delete [] saveObjective;
    330   if (usedColumn) {
    331     OsiSolverInterface * newSolver = model_->continuousSolver()->clone();
    332     const double * colLower = newSolver->getColLower();
    333     const double * colUpper = newSolver->getColUpper();
    334     int i;
    335     int nFix=0;
    336     int nFixI=0;
    337     int nFixC=0;
    338     for (i=0;i<numberIntegers;i++) {
    339       int iColumn=integerVariable[i];
    340       const CbcObject * object = model_->object(i);
    341       const CbcSimpleInteger * integerObject =
    342         dynamic_cast<const  CbcSimpleInteger *> (object);
    343       assert(integerObject);
    344       // get original bounds
    345       double originalLower = integerObject->originalLowerBound();
    346       assert(colLower[iColumn]==originalLower);
    347       newSolver->setColLower(iColumn,CoinMax(colLower[iColumn],originalLower));
    348       double originalUpper = integerObject->originalUpperBound();
    349       assert(colUpper[iColumn]==originalUpper);
    350       newSolver->setColUpper(iColumn,CoinMin(colUpper[iColumn],originalUpper));
    351       if (!usedColumn[iColumn]) {
    352         double value=lastSolution[iColumn];
    353         double nearest = floor(value+0.5);
    354         if (fabs(value-nearest)<1.0e-7) {
    355           if (nearest==colLower[iColumn]) {
    356             newSolver->setColUpper(iColumn,colLower[iColumn]);
    357             nFix++;
    358           } else if (nearest==colUpper[iColumn]) {
    359             newSolver->setColLower(iColumn,colUpper[iColumn]);
    360             nFix++;
    361           } else if (fixInternal) {
    362             newSolver->setColLower(iColumn,nearest);
    363             newSolver->setColUpper(iColumn,nearest);
    364             nFix++;
    365             nFixI++;
    366           }
    367         }
    368       }
    369     }
    370     if (fixContinuous) {
    371       for (int iColumn=0;iColumn<numberColumns;iColumn++) {
    372         if (!newSolver->isInteger(iColumn)&&!usedColumn[iColumn]) {
     471        int iColumn=integerVariable[i];
     472        const OsiObject * object = model_->object(i);
     473        double originalLower;
     474        double originalUpper;
     475        getIntegerInformation( object,originalLower, originalUpper);
     476        assert(colLower[iColumn]==originalLower);
     477        newSolver->setColLower(iColumn,CoinMax(colLower[iColumn],originalLower));
     478        assert(colUpper[iColumn]==originalUpper);
     479        newSolver->setColUpper(iColumn,CoinMin(colUpper[iColumn],originalUpper));
     480        if (!usedColumn[iColumn]) {
    373481          double value=lastSolution[iColumn];
    374           if (value<colLower[iColumn]+1.0e-8) {
    375             newSolver->setColUpper(iColumn,colLower[iColumn]);
    376             nFix++;
    377             nFixC++;
    378           } else if (value>colUpper[iColumn]-1.0e-8) {
    379             newSolver->setColLower(iColumn,colUpper[iColumn]);
    380             nFix++;
    381             nFixC++;
    382           }
    383         }
    384       }
    385     }
    386     newSolver->initialSolve();
    387     assert (newSolver->isProvenOptimal());
    388     printf("%d integers at bound fixed, %d internal and %d continuous\n",
    389            nFix,nFixI,nFixC);
    390     double saveValue = newSolutionValue;
    391     returnCode = smallBranchAndBound(newSolver,200,newSolution,newSolutionValue,
    392                                      newSolutionValue,"CbcHeuristicLocalAfterFPump");
    393     if (returnCode) {
    394       printf("old sol of %g new of %g\n",saveValue,newSolutionValue);
    395       memcpy(betterSolution,newSolution,numberColumns*sizeof(double));
    396       solutionValue=newSolutionValue;
    397     }
    398     delete newSolver;
    399     delete [] usedColumn;
    400     delete [] lastSolution;
     482          double nearest = floor(value+0.5);
     483          if (fabs(value-nearest)<1.0e-7) {
     484            if (nearest==colLower[iColumn]) {
     485              newSolver->setColUpper(iColumn,colLower[iColumn]);
     486              nFix++;
     487            } else if (nearest==colUpper[iColumn]) {
     488              newSolver->setColLower(iColumn,colUpper[iColumn]);
     489              nFix++;
     490            } else if (fixInternal) {
     491              newSolver->setColLower(iColumn,nearest);
     492              newSolver->setColUpper(iColumn,nearest);
     493              nFix++;
     494              nFixI++;
     495            }
     496          }
     497        }
     498      }
     499      if (fixContinuous) {
     500        for (int iColumn=0;iColumn<numberColumns;iColumn++) {
     501          if (!newSolver->isInteger(iColumn)&&!usedColumn[iColumn]) {
     502            double value=lastSolution[iColumn];
     503            if (value<colLower[iColumn]+1.0e-8) {
     504              newSolver->setColUpper(iColumn,colLower[iColumn]);
     505              nFixC++;
     506            } else if (value>colUpper[iColumn]-1.0e-8) {
     507              newSolver->setColLower(iColumn,colUpper[iColumn]);
     508              nFixC++;
     509            }
     510          }
     511        }
     512      }
     513      newSolver->initialSolve();
     514      if (!newSolver->isProvenOptimal()) {
     515        newSolver->writeMps("bad.mps");
     516        assert (newSolver->isProvenOptimal());
     517        break;
     518      }
     519      printf("%d integers at bound fixed (of which %d were internal) and %d continuous\n",
     520             nFix,nFixI,nFixC);
     521      double saveValue = newSolutionValue;
     522      returnCode = smallBranchAndBound(newSolver,200,newSolution,newSolutionValue,
     523                                       newSolutionValue,"CbcHeuristicLocalAfterFPump");
     524      if (returnCode) {
     525        printf("old sol of %g new of %g\n",saveValue,newSolutionValue);
     526        memcpy(betterSolution,newSolution,numberColumns*sizeof(double));
     527        solutionValue=newSolutionValue;
     528        solutionFound=true;
     529      }
     530      delete newSolver;
     531    }
     532    if (solutionFound) finalReturnCode=1;
     533    cutoff = CoinMin(cutoff,solutionValue);
     534    if (numberTries>=CoinAbs(maximumRetries_)||!solutionFound) {
     535      break;
     536    } else if (absoluteIncrement_>0.0||relativeIncrement_>0.0) {
     537      solutionFound=false;
     538      double gap = relativeIncrement_*fabs(solutionValue);
     539      cutoff -= CoinMax(CoinMax(gap,absoluteIncrement_),model_->getCutoffIncrement());
     540      printf("round again with cutoff of %g\n",cutoff);
     541      if (maximumRetries_<0)
     542        memset(usedColumn,0,numberColumns);
     543      totalNumberPasses += numberPasses;
     544    } else {
     545      break;
     546    }
    401547  }
     548  delete [] usedColumn;
     549  delete [] lastSolution;
    402550  delete [] newSolution;
    403   return returnCode;
     551  return finalReturnCode;
    404552}
    405553
     
    449597      continue;
    450598#ifndef NDEBUG
    451     const CbcObject * object = model_->object(i);
    452     const CbcSimpleInteger * integerObject =
    453       dynamic_cast<const  CbcSimpleInteger *> (object);
    454     assert(integerObject);
     599      const OsiObject * object = model_->object(i);
     600      const CbcSimpleInteger * integerObject =
     601        dynamic_cast<const  CbcSimpleInteger *> (object);
     602      const OsiSimpleInteger * integerObject2 =
     603        dynamic_cast<const  OsiSimpleInteger *> (object);
     604      assert(integerObject||integerObject2);
    455605#endif
    456606    double value=solution[iColumn];
Note: See TracChangeset for help on using the changeset viewer.