Ignore:
Timestamp:
Jun 26, 2007 5:17:15 AM (13 years ago)
Author:
forrest
Message:

trunk from branches/devel

Location:
trunk
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • trunk

    • Property svn:externals
      •  

        old new  
        1 MSVisualStudio   https://projects.coin-or.org/svn/MSVisualStudio/trunk/ExternalsDirs/Cbc
        2 BuildTools    https://projects.coin-or.org/svn/BuildTools/stable/0.5
         1MSVisualStudio   https://projects.coin-or.org/svn/MSVisualStudio/branches/devel/ExternalsDirs/Cbc
         2BuildTools    https://projects.coin-or.org/svn/BuildTools/trunk
        33ThirdParty/ASL https://projects.coin-or.org/svn/BuildTools/ThirdParty/ASL/stable/1.0
         4ThirdParty/Blas https://projects.coin-or.org/svn/BuildTools/ThirdParty/Blas/stable/1.0
         5ThirdParty/Lapack https://projects.coin-or.org/svn/BuildTools/ThirdParty/Lapack/stable/1.0
        46Data/Netlib   https://projects.coin-or.org/svn/Data/stable/1.0/Netlib
        57Data/Sample   https://projects.coin-or.org/svn/Data/stable/1.0/Sample
  • trunk/Cbc/src/CbcHeuristicFPump.cpp

    r626 r640  
    2323   startTime_(0.0),
    2424   maximumTime_(0.0),
     25   fakeCutoff_(COIN_DBL_MAX),
     26   absoluteIncrement_(0.0),
     27   relativeIncrement_(0.0),
     28   defaultRounding_(0.49999),
     29   initialWeight_(0.0),
     30   weightFactor_(0.1),
    2531   maximumPasses_(100),
    26    downValue_(0.5),
     32   maximumRetries_(1),
     33   accumulate_(0),
    2734   roundExpensive_(false)
    2835{
     
    3643   startTime_(0.0),
    3744   maximumTime_(0.0),
     45   fakeCutoff_(COIN_DBL_MAX),
     46   absoluteIncrement_(0.0),
     47   relativeIncrement_(0.0),
     48   defaultRounding_(downValue),
     49   initialWeight_(0.0),
     50   weightFactor_(0.1),
    3851   maximumPasses_(100),
    39    downValue_(downValue),
     52   maximumRetries_(1),
     53   accumulate_(0),
    4054   roundExpensive_(roundExpensive)
    4155{
     
    6175  fprintf(fp,"0#include \"CbcHeuristicFPump.hpp\"\n");
    6276  fprintf(fp,"3  CbcHeuristicFPump heuristicFPump(*cbcModel);\n");
     77  CbcHeuristic::generateCpp(fp,"heuristicFPump");
    6378  if (maximumPasses_!=other.maximumPasses_)
    6479    fprintf(fp,"3  heuristicFPump.setMaximumPasses(%d);\n",maximumPasses_);
    6580  else
    6681    fprintf(fp,"4  heuristicFPump.setMaximumPasses(%d);\n",maximumPasses_);
     82  if (maximumRetries_!=other.maximumRetries_)
     83    fprintf(fp,"3  heuristicFPump.setMaximumRetries(%d);\n",maximumRetries_);
     84  else
     85    fprintf(fp,"4  heuristicFPump.setMaximumRetries(%d);\n",maximumRetries_);
     86  if (accumulate_!=other.accumulate_)
     87    fprintf(fp,"3  heuristicFPump.setAccumulate(%d);\n",accumulate_);
     88  else
     89    fprintf(fp,"4  heuristicFPump.setAccumulate(%d);\n",accumulate_);
    6790  if (maximumTime_!=other.maximumTime_)
    6891    fprintf(fp,"3  heuristicFPump.setMaximumTime(%g);\n",maximumTime_);
    6992  else
    7093    fprintf(fp,"4  heuristicFPump.setMaximumTime(%g);\n",maximumTime_);
     94  if (fakeCutoff_!=other.fakeCutoff_)
     95    fprintf(fp,"3  heuristicFPump.setFakeCutoff(%g);\n",fakeCutoff_);
     96  else
     97    fprintf(fp,"4  heuristicFPump.setFakeCutoff(%g);\n",fakeCutoff_);
     98  if (absoluteIncrement_!=other.absoluteIncrement_)
     99    fprintf(fp,"3  heuristicFPump.setAbsoluteIncrement(%g);\n",absoluteIncrement_);
     100  else
     101    fprintf(fp,"4  heuristicFPump.setAbsoluteIncrement(%g);\n",absoluteIncrement_);
     102  if (relativeIncrement_!=other.relativeIncrement_)
     103    fprintf(fp,"3  heuristicFPump.setRelativeIncrement(%g);\n",relativeIncrement_);
     104  else
     105    fprintf(fp,"4  heuristicFPump.setRelativeIncrement(%g);\n",relativeIncrement_);
     106  if (defaultRounding_!=other.defaultRounding_)
     107    fprintf(fp,"3  heuristicFPump.setDefaultRounding(%g);\n",defaultRounding_);
     108  else
     109    fprintf(fp,"4  heuristicFPump.setDefaultRounding(%g);\n",defaultRounding_);
    71110  fprintf(fp,"3  cbcModel->addHeuristic(&heuristicFPump);\n");
     111  if (initialWeight_!=other.initialWeight_)
     112    fprintf(fp,"3  heuristicFPump.setInitialWeight(%g);\n",initialWeight_);
     113  else
     114    fprintf(fp,"4  heuristicFPump.setInitialWeight(%g);\n",initialWeight_);
     115  if (weightFactor_!=other.weightFactor_)
     116    fprintf(fp,"3  heuristicFPump.setWeightFactor(%g);\n",weightFactor_);
     117  else
     118    fprintf(fp,"4  heuristicFPump.setWeightFactor(%g);\n",weightFactor_);
    72119}
    73120
     
    78125  startTime_(rhs.startTime_),
    79126  maximumTime_(rhs.maximumTime_),
     127  fakeCutoff_(rhs.fakeCutoff_),
     128  absoluteIncrement_(rhs.absoluteIncrement_),
     129  relativeIncrement_(rhs.relativeIncrement_),
     130  defaultRounding_(rhs.defaultRounding_),
     131  initialWeight_(rhs.initialWeight_),
     132  weightFactor_(rhs.weightFactor_),
    80133  maximumPasses_(rhs.maximumPasses_),
    81   downValue_(rhs.downValue_),
     134  maximumRetries_(rhs.maximumRetries_),
     135  accumulate_(rhs.accumulate_),
    82136  roundExpensive_(rhs.roundExpensive_)
    83137{
    84   setWhen(rhs.when());
    85 }
     138}
     139
     140// Assignment operator
     141CbcHeuristicFPump &
     142CbcHeuristicFPump::operator=( const CbcHeuristicFPump& rhs)
     143{
     144  if (this!=&rhs) {
     145    CbcHeuristic::operator=(rhs);
     146    startTime_ = rhs.startTime_;
     147    maximumTime_ = rhs.maximumTime_;
     148    fakeCutoff_ = rhs.fakeCutoff_;
     149    absoluteIncrement_ = rhs.absoluteIncrement_;
     150    relativeIncrement_ = rhs.relativeIncrement_;
     151    defaultRounding_ = rhs.defaultRounding_;
     152    initialWeight_ = rhs.initialWeight_;
     153    weightFactor_ = rhs.weightFactor_;
     154    maximumPasses_ = rhs.maximumPasses_;
     155    maximumRetries_ = rhs.maximumRetries_;
     156    accumulate_ = rhs.accumulate_;
     157    roundExpensive_ = rhs.roundExpensive_;
     158  }
     159  return *this;
     160}
     161
    86162// Resets stuff if model changes
    87163void
     
    109185  // probably a good idea
    110186  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  
     187  // loop round doing repeated pumps
     188  double cutoff;
     189  model_->solver()->getDblParam(OsiDualObjectiveLimit,cutoff);
     190  double direction = model_->solver()->getObjSense();
     191  cutoff *= direction;
     192  double firstCutoff = fabs(cutoff);
     193  cutoff = CoinMin(cutoff,solutionValue);
     194  // check plausible and space for rounded solution
    121195  int numberColumns = model_->getNumCols();
    122196  int numberIntegers = model_->numberIntegers();
    123   const int * integerVariable = model_->integerVariable();
    124 
    125 // 1. initially check 0-1
     197  const int * integerVariableOrig = model_->integerVariable();
     198
     199  // 1. initially check 0-1
    126200  int i,j;
    127201  int general=0;
     202  int * integerVariable = new int[numberIntegers];
     203  const double * lower = model_->solver()->getColLower();
     204  const double * upper = model_->solver()->getColUpper();
     205  bool doGeneral = (accumulate_&(4+8))!=0;
     206  j=0;
    128207  for (i=0;i<numberIntegers;i++) {
    129     int iColumn = integerVariable[i];
     208    int iColumn = integerVariableOrig[i];
    130209#ifndef NDEBUG
    131     const CbcObject * object = model_->object(i);
     210    const OsiObject * object = model_->object(i);
    132211    const CbcSimpleInteger * integerObject =
    133212      dynamic_cast<const  CbcSimpleInteger *> (object);
    134     assert(integerObject);
     213    const OsiSimpleInteger * integerObject2 =
     214      dynamic_cast<const  OsiSimpleInteger *> (object);
     215    assert(integerObject||integerObject2);
    135216#endif
    136217    if (upper[iColumn]-lower[iColumn]>1.000001) {
    137218      general++;
     219      if (doGeneral)
     220        integerVariable[j++]=iColumn;
     221    } else {
     222      integerVariable[j++]=iColumn;
     223    }
     224  }
     225  if (general*3>2*numberIntegers&&!doGeneral) {
     226    delete [] integerVariable;
     227    return 0;
     228  } else if ((accumulate_&8)==0) {
     229    doGeneral=false;
     230  }
     231  if (!general)
     232    doGeneral=false;
     233  // For solution closest to feasible if none found
     234  int * closestSolution = general ? NULL : new int[numberIntegers];
     235  double closestObjectiveValue = COIN_DBL_MAX;
     236 
     237  int numberIntegersOrig = numberIntegers;
     238  numberIntegers = j;
     239  double * newSolution = new double [numberColumns];
     240  double newSolutionValue=COIN_DBL_MAX;
     241  bool solutionFound=false;
     242  char * usedColumn = NULL;
     243  double * lastSolution=NULL;
     244  int fixContinuous=0;
     245  bool fixInternal=false;
     246  bool allSlack=false;
     247  if (when_>=21&&when_<=25) {
     248    when_ -= 10;
     249    allSlack=true;
     250  }
     251  if (when_>=11&&when_<=15) {
     252    fixInternal = when_ >11&&when_<15;
     253    if (when_<13)
     254      fixContinuous = 0;
     255    else if (when_!=14)
     256      fixContinuous=1;
     257    else
     258      fixContinuous=2;
     259    when_=1;
     260    usedColumn = new char [numberColumns];
     261    memset(usedColumn,0,numberColumns);
     262    model_->solver()->resolve();
     263    assert (model_->solver()->isProvenOptimal());
     264    lastSolution = CoinCopyOfArray(model_->solver()->getColSolution(),numberColumns);
     265  }
     266  int finalReturnCode=0;
     267  int totalNumberPasses=0;
     268  int numberTries=0;
     269  while (true) {
     270    int numberPasses=0;
     271    numberTries++;
     272    // Clone solver - otherwise annoys root node computations
     273    OsiSolverInterface * solver = model_->solver()->clone();
     274    // if cutoff exists then add constraint
     275    bool useCutoff = (fabs(cutoff)<1.0e20&&(fakeCutoff_!=COIN_DBL_MAX||numberTries>1));
     276    // but there may be a close one
     277    if (firstCutoff<2.0*solutionValue&&numberTries==1)
     278      useCutoff=true;
     279    if (useCutoff) {
     280      cutoff = CoinMin(cutoff,fakeCutoff_);
     281      const double * objective = solver->getObjCoefficients();
     282      int numberColumns = solver->getNumCols();
     283      int * which = new int[numberColumns];
     284      double * els = new double[numberColumns];
     285      int nel=0;
     286      for (int i=0;i<numberColumns;i++) {
     287        double value = objective[i];
     288        if (value) {
     289          which[nel]=i;
     290          els[nel++] = direction*value;
     291        }
     292      }
     293      solver->addRow(nel,which,els,-COIN_DBL_MAX,cutoff);
     294      delete [] which;
     295      delete [] els;
     296      bool takeHint;
     297      OsiHintStrength strength;
     298      solver->getHintParam(OsiDoDualInResolve,takeHint,strength);
     299      solver->setHintParam(OsiDoDualInResolve,true,OsiHintDo);
     300      solver->resolve();
     301      solver->setHintParam(OsiDoDualInResolve,takeHint,strength);
     302    }
     303    solver->setDblParam(OsiDualObjectiveLimit,1.0e50);
     304    solver->resolve();
     305    // Solver may not be feasible
     306    if (!solver->isProvenOptimal()) {
    138307      break;
    139308    }
    140   }
    141   if (general*3>numberIntegers) {
    142     delete solver;
    143     return 0;
    144   }
    145 
    146 // 2. space for rounded solution
    147   double * newSolution = new double [numberColumns];
    148   // space for last rounded solutions
     309    const double * lower = solver->getColLower();
     310    const double * upper = solver->getColUpper();
     311    const double * solution = solver->getColSolution();
     312    if (lastSolution)
     313      memcpy(lastSolution,solution,numberColumns*sizeof(double));
     314    double primalTolerance;
     315    solver->getDblParam(OsiPrimalTolerance,primalTolerance);
     316   
     317   
     318    // 2 space for last rounded solutions
    149319#define NUMBER_OLD 4
    150   double ** oldSolution = new double * [NUMBER_OLD];
    151   for (j=0;j<NUMBER_OLD;j++) {
    152     oldSolution[j]= new double[numberColumns];
    153     for (i=0;i<numberColumns;i++) oldSolution[j][i]=-COIN_DBL_MAX;
    154   }
    155 
    156 // 3. Replace objective with an initial 0-valued objective
    157   double * saveObjective = new double [numberColumns];
    158   memcpy(saveObjective,solver->getObjCoefficients(),numberColumns*sizeof(double));
    159   for (i=0;i<numberColumns;i++) {
    160     solver->setObjCoeff(i,0.0);
    161   }
    162   bool finished=false;
    163   double direction = solver->getObjSense();
    164   int returnCode=0;
    165   bool takeHint;
    166   OsiHintStrength strength;
    167   solver->getHintParam(OsiDoDualInResolve,takeHint,strength);
    168   solver->setHintParam(OsiDoDualInResolve,false);
    169   solver->messageHandler()->setLogLevel(0);
    170 
    171 // 4. Save objective offset so we can see progress
    172   double saveOffset;
    173   solver->getDblParam(OsiObjOffset,saveOffset);
    174 
    175 // 5. MAIN WHILE LOOP
    176   int numberPasses=0;
    177   bool newLineNeeded=false;
    178   while (!finished) {
    179     returnCode=0;
    180     if (numberPasses>=maximumPasses_) {
    181       break;
    182     }
    183     if (maximumTime_>0.0&&CoinCpuTime()>=startTime_+maximumTime_) break;
    184     numberPasses++;
    185     memcpy(newSolution,solution,numberColumns*sizeof(double));
    186     int flip;
    187     returnCode = rounds(solver,newSolution,saveObjective,roundExpensive_,downValue_,&flip);
    188     if (returnCode) {
    189       // SOLUTION IS INTEGER
    190       // Put back correct objective
    191       for (i=0;i<numberColumns;i++)
    192         solver->setObjCoeff(i,saveObjective[i]);
    193       // solution - but may not be better
    194       // Compute using dot product
    195       solver->setDblParam(OsiObjOffset,saveOffset);
    196       double newSolutionValue = direction*solver->OsiSolverInterface::getObjValue();
    197       if (model_->logLevel())
    198         printf(" - solution found\n");
    199       newLineNeeded=false;
    200       if (newSolutionValue<solutionValue) {
    201         if (general) {
    202           int numberLeft=0;
    203           for (i=0;i<numberIntegers;i++) {
    204             int iColumn = integerVariable[i];
    205             double value = floor(newSolution[iColumn]+0.5);
    206             if(solver->isBinary(iColumn)) {
    207               solver->setColLower(iColumn,value);
    208               solver->setColUpper(iColumn,value);
    209             } else {
    210               if (fabs(value-newSolution[iColumn])>1.0e-7)
    211                 numberLeft++;
    212             }
    213           }
    214           if (numberLeft) {
    215             returnCode = smallBranchAndBound(solver,200,newSolution,newSolutionValue,
    216                                          solutionValue,"CbcHeuristicFpump");
    217           }
    218         }
    219         if (returnCode) {
    220           memcpy(betterSolution,newSolution,numberColumns*sizeof(double));
    221           solutionValue=newSolutionValue;
    222         }
     320    double ** oldSolution = new double * [NUMBER_OLD];
     321    for (j=0;j<NUMBER_OLD;j++) {
     322      oldSolution[j]= new double[numberColumns];
     323      for (i=0;i<numberColumns;i++) oldSolution[j][i]=-COIN_DBL_MAX;
     324    }
     325   
     326    // 3. Replace objective with an initial 0-valued objective
     327    double * saveObjective = new double [numberColumns];
     328    memcpy(saveObjective,solver->getObjCoefficients(),numberColumns*sizeof(double));
     329    for (i=0;i<numberColumns;i++) {
     330      solver->setObjCoeff(i,0.0);
     331    }
     332    bool finished=false;
     333    double direction = solver->getObjSense();
     334    int returnCode=0;
     335    bool takeHint;
     336    OsiHintStrength strength;
     337    solver->getHintParam(OsiDoDualInResolve,takeHint,strength);
     338    solver->setHintParam(OsiDoDualInResolve,false);
     339    solver->messageHandler()->setLogLevel(0);
     340   
     341    // 4. Save objective offset so we can see progress
     342    double saveOffset;
     343    solver->getDblParam(OsiObjOffset,saveOffset);
     344    // Get amount for original objective
     345    double scaleFactor = 0.0;
     346    for (i=0;i<numberColumns;i++)
     347      scaleFactor += saveObjective[i]*saveObjective[i];
     348    if (scaleFactor)
     349      scaleFactor = (initialWeight_*sqrt((double) numberIntegers))/sqrt(scaleFactor);
     350    // 5. MAIN WHILE LOOP
     351    bool newLineNeeded=false;
     352    while (!finished) {
     353      returnCode=0;
     354      // see what changed
     355      if (usedColumn) {
     356        for (i=0;i<numberColumns;i++) {
     357          if (fabs(solution[i]-lastSolution[i])>1.0e-8)
     358            usedColumn[i]=1;
     359          lastSolution[i]=solution[i];
     360        }
     361      }
     362      if (numberPasses>=maximumPasses_) {
     363        break;
     364      }
     365      if (maximumTime_>0.0&&CoinCpuTime()>=startTime_+maximumTime_) break;
     366      numberPasses++;
     367      memcpy(newSolution,solution,numberColumns*sizeof(double));
     368      int flip;
     369      returnCode = rounds(solver,newSolution,saveObjective,numberIntegers,integerVariable,
     370                          roundExpensive_,defaultRounding_,&flip);
     371      if (returnCode) {
     372        // SOLUTION IS INTEGER
     373        // Put back correct objective
     374        for (i=0;i<numberColumns;i++)
     375          solver->setObjCoeff(i,saveObjective[i]);
     376
     377        // solution - but may not be better
     378        // Compute using dot product
     379        solver->setDblParam(OsiObjOffset,saveOffset);
     380        newSolutionValue = -saveOffset;
     381        for (  i=0 ; i<numberColumns ; i++ )
     382          newSolutionValue += saveObjective[i]*newSolution[i];
     383        newSolutionValue *= direction;
     384        if (model_->logLevel())
     385          printf(" - solution found of %g",newSolutionValue);
     386        newLineNeeded=false;
     387        if (newSolutionValue<solutionValue) {
     388          double saveValue = newSolutionValue;
     389          if (!doGeneral) {
     390            int numberLeft=0;
     391            for (i=0;i<numberIntegersOrig;i++) {
     392              int iColumn = integerVariableOrig[i];
     393              double value = floor(newSolution[iColumn]+0.5);
     394              if(solver->isBinary(iColumn)) {
     395                solver->setColLower(iColumn,value);
     396                solver->setColUpper(iColumn,value);
     397              } else {
     398                if (fabs(value-newSolution[iColumn])>1.0e-7)
     399                  numberLeft++;
     400              }
     401            }
     402            if (numberLeft) {
     403              returnCode = smallBranchAndBound(solver,numberNodes_,newSolution,newSolutionValue,
     404                                               solutionValue,"CbcHeuristicFpump");
     405              if ((returnCode&2)!=0) {
     406                // could add cut
     407                returnCode &= ~2;
     408              }
     409            }
     410          }
     411          if (returnCode) {
     412            memcpy(betterSolution,newSolution,numberColumns*sizeof(double));
     413            if ((accumulate_&1)!=0)
     414              model_->incrementUsed(betterSolution); // for local search
     415            solutionValue=newSolutionValue;
     416            solutionFound=true;
     417            if (general&&saveValue!=newSolutionValue)
     418              printf(" - cleaned solution of %g\n",solutionValue);
     419            else
     420              printf("\n");
     421          } else {
     422          if (model_->logLevel())
     423            printf(" - not good enough after small branch and bound\n");
     424          }
     425        } else {
     426          if (model_->logLevel())
     427            printf(" - not good enough\n");
     428          newLineNeeded=false;
     429          returnCode=0;
     430        }     
     431        break;
    223432      } else {
    224         returnCode=0;
    225       }     
    226       break;
    227     } else {
    228       // SOLUTION IS not INTEGER
    229       // 1. check for loop
    230       bool matched;
    231       for (int k = NUMBER_OLD-1; k > 0; k--) {
     433        // SOLUTION IS not INTEGER
     434        // 1. check for loop
     435        bool matched;
     436        for (int k = NUMBER_OLD-1; k > 0; k--) {
    232437          double * b = oldSolution[k];
    233438          matched = true;
    234439          for (i = 0; i <numberIntegers; i++) {
     440            int iColumn = integerVariable[i];
     441            if (newSolution[iColumn]!=b[iColumn]) {
     442              matched=false;
     443              break;
     444            }
     445          }
     446          if (matched) break;
     447        }
     448        if (matched || numberPasses%100 == 0) {
     449          // perturbation
     450          if (model_->logLevel())
     451            printf("Perturbation applied");
     452          newLineNeeded=true;
     453          for (i=0;i<numberIntegers;i++) {
     454            int iColumn = integerVariable[i];
     455            double value = max(0.0,CoinDrand48()-0.3);
     456            double difference = fabs(solution[iColumn]-newSolution[iColumn]);
     457            if (difference+value>0.5) {
     458              if (newSolution[iColumn]<lower[iColumn]+primalTolerance) {
     459                newSolution[iColumn] += 1.0;
     460              } else if (newSolution[iColumn]>upper[iColumn]-primalTolerance) {
     461                newSolution[iColumn] -= 1.0;
     462              } else {
     463                // general integer
     464                if (difference+value>0.75)
     465                  newSolution[iColumn] += 1.0;
     466                else
     467                  newSolution[iColumn] -= 1.0;
     468              }
     469            }
     470          }
     471        } else {
     472          for (j=NUMBER_OLD-1;j>0;j--) {
     473            for (i = 0; i < numberColumns; i++) oldSolution[j][i]=oldSolution[j-1][i];
     474          }
     475          for (j = 0; j < numberColumns; j++) oldSolution[0][j] = newSolution[j];
     476        }
     477       
     478        // 2. update the objective function based on the new rounded solution
     479        double offset=0.0;
     480        double costValue = (1.0-scaleFactor)*solver->getObjSense();
     481       
     482        for (i=0;i<numberColumns;i++) {
     483          // below so we can keep original code and allow for objective
     484          int iColumn = i;
     485          if(!solver->isBinary(iColumn)&&!doGeneral)
     486            continue;
     487          // deal with fixed variables (i.e., upper=lower)
     488          if (fabs(lower[iColumn]-upper[iColumn]) < primalTolerance) {
     489            //if (lower[iColumn] > 1. - primalTolerance) solver->setObjCoeff(iColumn,-costValue);
     490            //else                                       solver->setObjCoeff(iColumn,costValue);
     491            continue;
     492          }
     493          if (newSolution[iColumn]<lower[iColumn]+primalTolerance) {
     494            solver->setObjCoeff(iColumn,costValue+scaleFactor*saveObjective[iColumn]);
     495          } else {
     496            if (newSolution[iColumn]>upper[iColumn]-primalTolerance) {
     497              solver->setObjCoeff(iColumn,-costValue+scaleFactor*saveObjective[iColumn]);
     498            } else {
     499              solver->setObjCoeff(iColumn,0.0);
     500            }
     501          }
     502          offset += costValue*newSolution[iColumn];
     503        }
     504        solver->setDblParam(OsiObjOffset,-offset);
     505        if (!general&&false) {
     506          // Solve in two goes - first keep satisfied ones fixed
     507          double * saveLower = new double [numberIntegers];
     508          double * saveUpper = new double [numberIntegers];
     509          for (i=0;i<numberIntegers;i++) {
     510            int iColumn = integerVariable[i];
     511            saveLower[i]=COIN_DBL_MAX;
     512            saveUpper[i]=-COIN_DBL_MAX;
     513            if (solution[iColumn]<lower[iColumn]+primalTolerance) {
     514              saveUpper[i]=upper[iColumn];
     515              solver->setColUpper(iColumn,lower[iColumn]);
     516            } else if (solution[iColumn]>upper[iColumn]-primalTolerance) {
     517              saveLower[i]=lower[iColumn];
     518              solver->setColLower(iColumn,upper[iColumn]);
     519            }
     520          }
     521          solver->resolve();
     522          assert (solver->isProvenOptimal());
     523          for (i=0;i<numberIntegers;i++) {
     524            int iColumn = integerVariable[i];
     525            if (saveLower[i]!=COIN_DBL_MAX)
     526              solver->setColLower(iColumn,saveLower[i]);
     527            if (saveUpper[i]!=-COIN_DBL_MAX)
     528              solver->setColUpper(iColumn,saveUpper[i]);
     529            saveUpper[i]=-COIN_DBL_MAX;
     530          }
     531          memcpy(newSolution,solution,numberColumns*sizeof(double));
     532          int flip;
     533          returnCode = rounds(solver,newSolution,saveObjective,numberIntegers,integerVariable,
     534                              roundExpensive_,defaultRounding_,&flip);
     535          if (returnCode) {
     536            // solution - but may not be better
     537            // Compute using dot product
     538            double newSolutionValue = -saveOffset;
     539            for (  i=0 ; i<numberColumns ; i++ )
     540              newSolutionValue += saveObjective[i]*newSolution[i];
     541            newSolutionValue *= direction;
     542            if (model_->logLevel())
     543              printf(" - intermediate solution found of %g",newSolutionValue);
     544            if (newSolutionValue<solutionValue) {
     545              memcpy(betterSolution,newSolution,numberColumns*sizeof(double));
     546              if ((accumulate_&1)!=0)
     547                model_->incrementUsed(betterSolution); // for local search
     548              solutionValue=newSolutionValue;
     549              solutionFound=true;
     550            } else {
     551              returnCode=0;
     552            }
     553          }     
     554        }
     555        if (!doGeneral) {
     556          // faster to do from all slack!!!!
     557          if (allSlack) {
     558            CoinWarmStartBasis dummy;
     559            solver->setWarmStart(&dummy);
     560          }
     561          solver->resolve();
     562          assert (solver->isProvenOptimal());
     563          // in case very dubious solver
     564          lower = solver->getColLower();
     565          upper = solver->getColUpper();
     566          solution = solver->getColSolution();
     567        } else {
     568          int * addStart = new int[2*general+1];
     569          int * addIndex = new int[4*general];
     570          double * addElement = new double[4*general];
     571          double * addLower = new double[2*general];
     572          double * addUpper = new double[2*general];
     573          double * obj = new double[general];
     574          int nAdd=0;
     575          for (i=0;i<numberIntegers;i++) {
     576            int iColumn = integerVariable[i];
     577            if (newSolution[iColumn]>lower[iColumn]+primalTolerance&&
     578                newSolution[iColumn]<upper[iColumn]-primalTolerance) {
     579              obj[nAdd]=1.0;
     580              addLower[nAdd]=0.0;
     581              addUpper[nAdd]=COIN_DBL_MAX;
     582              nAdd++;
     583            }
     584          }
     585          OsiSolverInterface * solver2 = solver;
     586          if (nAdd) {
     587            CoinZeroN(addStart,nAdd+1);
     588            solver2 = solver->clone();
     589            solver2->addCols(nAdd,addStart,NULL,NULL,addLower,addUpper,obj);
     590            // feasible solution
     591            double * sol = new double[nAdd+numberColumns];
     592            memcpy(sol,solution,numberColumns*sizeof(double));
     593            // now rows
     594            int nAdd=0;
     595            int nEl=0;
     596            int nAddRow=0;
     597            for (i=0;i<numberIntegers;i++) {
    235598              int iColumn = integerVariable[i];
    236               if(!solver->isBinary(iColumn))
    237                 continue;
    238               if (newSolution[iColumn]!=b[iColumn]) {
    239                 matched=false;
    240                 break;
     599              if (newSolution[iColumn]>lower[iColumn]+primalTolerance&&
     600                  newSolution[iColumn]<upper[iColumn]-primalTolerance) {
     601                addLower[nAddRow]=-newSolution[iColumn];;
     602                addUpper[nAddRow]=COIN_DBL_MAX;
     603                addIndex[nEl] = iColumn;
     604                addElement[nEl++]=-1.0;
     605                addIndex[nEl] = numberColumns+nAdd;
     606                addElement[nEl++]=1.0;
     607                nAddRow++;
     608                addStart[nAddRow]=nEl;
     609                addLower[nAddRow]=newSolution[iColumn];;
     610                addUpper[nAddRow]=COIN_DBL_MAX;
     611                addIndex[nEl] = iColumn;
     612                addElement[nEl++]=1.0;
     613                addIndex[nEl] = numberColumns+nAdd;
     614                addElement[nEl++]=1.0;
     615                nAddRow++;
     616                addStart[nAddRow]=nEl;
     617                sol[nAdd+numberColumns] = fabs(sol[iColumn]-newSolution[iColumn]);
     618                nAdd++;
    241619              }
    242           }
    243           if (matched) break;
     620            }
     621            solver2->setColSolution(sol);
     622            delete [] sol;
     623            solver2->addRows(nAddRow,addStart,addIndex,addElement,addLower,addUpper);
     624          }
     625          delete [] addStart;
     626          delete [] addIndex;
     627          delete [] addElement;
     628          delete [] addLower;
     629          delete [] addUpper;
     630          delete [] obj;
     631          solver2->resolve();
     632          assert (solver2->isProvenOptimal());
     633          if (nAdd) {
     634            solver->setColSolution(solver2->getColSolution());
     635            delete solver2;
     636          }
     637        }
     638        if (model_->logLevel())
     639          printf("\npass %3d: obj. %10.5f --> ", numberPasses+totalNumberPasses,solver->getObjValue());
     640        if (closestSolution&&solver->getObjValue()<closestObjectiveValue) {
     641          int i;
     642          const double * objective = solver->getObjCoefficients();
     643          for (i=0;i<numberIntegersOrig;i++) {
     644            int iColumn=integerVariableOrig[i];
     645            if (objective[iColumn]>0.0)
     646              closestSolution[i]=0;
     647            else
     648              closestSolution[i]=1;
     649          }
     650          closestObjectiveValue = solver->getObjValue();
     651        }
     652        newLineNeeded=true;
     653       
    244654      }
    245       if (matched || numberPasses%100 == 0) {
    246          // perturbation
    247         if (model_->logLevel())
    248           printf("Perturbation applied");
    249          newLineNeeded=true;
    250          for (i=0;i<numberIntegers;i++) {
    251              int iColumn = integerVariable[i];
    252              if(!solver->isBinary(iColumn))
    253                continue;
    254              double value = max(0.0,CoinDrand48()-0.3);
    255              double difference = fabs(solution[iColumn]-newSolution[iColumn]);
    256              if (difference+value>0.5) {
    257                 if (newSolution[iColumn]<lower[iColumn]+primalTolerance) newSolution[iColumn] += 1.0;
    258              else if (newSolution[iColumn]>upper[iColumn]-primalTolerance) newSolution[iColumn] -= 1.0;
    259                   else abort();
    260              }
    261          }
     655      // reduce scale factor
     656      scaleFactor *= weightFactor_;
     657    } // END WHILE
     658    if (newLineNeeded&&model_->logLevel())
     659      printf(" - no solution found\n");
     660    delete solver;
     661    for ( j=0;j<NUMBER_OLD;j++)
     662      delete [] oldSolution[j];
     663    delete [] oldSolution;
     664    delete [] saveObjective;
     665    if (usedColumn) {
     666      OsiSolverInterface * newSolver = model_->continuousSolver()->clone();
     667      const double * colLower = newSolver->getColLower();
     668      const double * colUpper = newSolver->getColUpper();
     669      int i;
     670      int nFix=0;
     671      int nFixI=0;
     672      int nFixC=0;
     673      int nFixC2=0;
     674      for (i=0;i<numberIntegersOrig;i++) {
     675        int iColumn=integerVariableOrig[i];
     676        //const OsiObject * object = model_->object(i);
     677        //double originalLower;
     678        //double originalUpper;
     679        //getIntegerInformation( object,originalLower, originalUpper);
     680        //assert(colLower[iColumn]==originalLower);
     681        //newSolver->setColLower(iColumn,CoinMax(colLower[iColumn],originalLower));
     682        newSolver->setColLower(iColumn,colLower[iColumn]);
     683        //assert(colUpper[iColumn]==originalUpper);
     684        //newSolver->setColUpper(iColumn,CoinMin(colUpper[iColumn],originalUpper));
     685        newSolver->setColUpper(iColumn,colUpper[iColumn]);
     686        if (!usedColumn[iColumn]) {
     687          double value=lastSolution[iColumn];
     688          double nearest = floor(value+0.5);
     689          if (fabs(value-nearest)<1.0e-7) {
     690            if (nearest==colLower[iColumn]) {
     691              newSolver->setColUpper(iColumn,colLower[iColumn]);
     692              nFix++;
     693            } else if (nearest==colUpper[iColumn]) {
     694              newSolver->setColLower(iColumn,colUpper[iColumn]);
     695              nFix++;
     696            } else if (fixInternal) {
     697              newSolver->setColLower(iColumn,nearest);
     698              newSolver->setColUpper(iColumn,nearest);
     699              nFix++;
     700              nFixI++;
     701            }
     702          }
     703        }
     704      }
     705      if (fixContinuous) {
     706        for (int iColumn=0;iColumn<numberColumns;iColumn++) {
     707          if (!newSolver->isInteger(iColumn)&&!usedColumn[iColumn]) {
     708            double value=lastSolution[iColumn];
     709            if (value<colLower[iColumn]+1.0e-8) {
     710              newSolver->setColUpper(iColumn,colLower[iColumn]);
     711              nFixC++;
     712            } else if (value>colUpper[iColumn]-1.0e-8) {
     713              newSolver->setColLower(iColumn,colUpper[iColumn]);
     714              nFixC++;
     715            } else if (fixContinuous==2) {
     716              newSolver->setColLower(iColumn,value);
     717              newSolver->setColUpper(iColumn,value);
     718              nFixC++;
     719              nFixC2++;
     720            }
     721          }
     722        }
     723      }
     724      newSolver->initialSolve();
     725      if (!newSolver->isProvenOptimal()) {
     726        newSolver->writeMps("bad.mps");
     727        assert (newSolver->isProvenOptimal());
     728        break;
     729      }
     730      printf("%d integers at bound fixed and %d continuous",
     731             nFix,nFixC);
     732      if (nFixC2+nFixI==0)
     733        printf("\n");
     734      else
     735        printf(" of which %d were internal integer and %d internal continuous\n",
     736             nFixI,nFixC2);
     737      double saveValue = newSolutionValue;
     738      returnCode = smallBranchAndBound(newSolver,numberNodes_,newSolution,newSolutionValue,
     739                                       newSolutionValue,"CbcHeuristicLocalAfterFPump");
     740      if ((returnCode&2)!=0) {
     741        // could add cut
     742        returnCode &= ~2;
     743      }
     744      if (returnCode) {
     745        printf("old sol of %g new of %g\n",saveValue,newSolutionValue);
     746        memcpy(betterSolution,newSolution,numberColumns*sizeof(double));
     747        if (fixContinuous) {
     748          // may be able to do even better
     749          const double * lower = model_->solver()->getColLower();
     750          const double * upper = model_->solver()->getColUpper();
     751          for (int iColumn=0;iColumn<numberColumns;iColumn++) {
     752            if (newSolver->isInteger(iColumn)) {
     753              double value=floor(newSolution[iColumn]+0.5);
     754              newSolver->setColLower(iColumn,value);
     755              newSolver->setColUpper(iColumn,value);
     756            } else {
     757              newSolver->setColLower(iColumn,lower[iColumn]);
     758              newSolver->setColUpper(iColumn,upper[iColumn]);
     759            }
     760          }
     761          newSolver->initialSolve();
     762          if (newSolver->isProvenOptimal()) {
     763            double value = newSolver->getObjValue()*newSolver->getObjSense();
     764            if (value<newSolutionValue) {
     765              printf("freeing continuous gives a solution of %g\n", value);
     766              newSolutionValue=value;
     767              memcpy(betterSolution,newSolver->getColSolution(),numberColumns*sizeof(double));
     768            }
     769          } else {
     770            newSolver->writeMps("bad3.mps");
     771          }
     772        }
     773        if ((accumulate_&1)!=0)
     774          model_->incrementUsed(betterSolution); // for local search
     775        solutionValue=newSolutionValue;
     776        solutionFound=true;
     777      }
     778      delete newSolver;
     779    }
     780    if (solutionFound) finalReturnCode=1;
     781    cutoff = CoinMin(cutoff,solutionValue);
     782    if (numberTries>=maximumRetries_||!solutionFound) {
     783      break;
     784    } else if (absoluteIncrement_>0.0||relativeIncrement_>0.0) {
     785      solutionFound=false;
     786      double gap = relativeIncrement_*fabs(solutionValue);
     787      cutoff -= CoinMax(CoinMax(gap,absoluteIncrement_),model_->getCutoffIncrement());
     788      printf("round again with cutoff of %g\n",cutoff);
     789      if ((accumulate_&3)<2&&usedColumn)
     790        memset(usedColumn,0,numberColumns);
     791      totalNumberPasses += numberPasses;
     792    } else {
     793      break;
     794    }
     795  }
     796  if (!finalReturnCode&&closestSolution&&closestObjectiveValue <= 10.0&&usedColumn) {
     797    // try a bit of branch and bound
     798    OsiSolverInterface * newSolver = model_->continuousSolver()->clone();
     799    const double * colLower = newSolver->getColLower();
     800    const double * colUpper = newSolver->getColUpper();
     801    int i;
     802    double rhs = 0.0;
     803    for (i=0;i<numberIntegersOrig;i++) {
     804      int iColumn=integerVariableOrig[i];
     805      int direction = closestSolution[i];
     806      closestSolution[i]=iColumn;
     807      if (direction==0) {
     808        // keep close to LB
     809        rhs += colLower[iColumn];
     810        lastSolution[i]=1.0;
    262811      } else {
    263          for (j=NUMBER_OLD-1;j>0;j--) {
    264              for (i = 0; i < numberColumns; i++) oldSolution[j][i]=oldSolution[j-1][i];
    265          }
    266          for (j = 0; j < numberColumns; j++) oldSolution[0][j] = newSolution[j];
     812        // keep close to UB
     813        rhs -= colUpper[iColumn];
     814        lastSolution[i]=-1.0;
    267815      }
    268 
    269       // 2. update the objective function based on the new rounded solution
    270       double offset=0.0;
    271       for (i=0;i<numberIntegers;i++) {
    272         int iColumn = integerVariable[i];
    273         if(!solver->isBinary(iColumn))
    274           continue;
    275         double costValue = 1.0;
    276         // deal with fixed variables (i.e., upper=lower)
    277         if (fabs(lower[iColumn]-upper[iColumn]) < primalTolerance) {
    278            if (lower[iColumn] > 1. - primalTolerance) solver->setObjCoeff(iColumn,-costValue);
    279            else                                       solver->setObjCoeff(iColumn,costValue);
    280            continue;
    281         }
    282         if (newSolution[iColumn]<lower[iColumn]+primalTolerance) {
    283           solver->setObjCoeff(iColumn,costValue);
    284         } else {
    285           if (newSolution[iColumn]>upper[iColumn]-primalTolerance) {
    286             solver->setObjCoeff(iColumn,-costValue);
    287           } else {
    288             solver->setObjCoeff(iColumn,0.0);
    289           }
    290         }
    291         offset += costValue*newSolution[iColumn];
    292       }
    293       solver->setDblParam(OsiObjOffset,-offset);
    294       solver->resolve();
    295       if (model_->logLevel())
    296         printf("\npass %3d: obj. %10.5lf --> ", numberPasses,solver->getObjValue());
    297       newLineNeeded=true;
    298 
    299     }
    300   } // END WHILE
    301   if (newLineNeeded&&model_->logLevel())
    302     printf(" - no solution found\n");
    303   delete solver;
     816    }
     817    newSolver->addRow(numberIntegersOrig,closestSolution,
     818                      lastSolution,-COIN_DBL_MAX,rhs+10.0);
     819    double saveValue = newSolutionValue;
     820    //newSolver->writeMps("sub");
     821    int returnCode = smallBranchAndBound(newSolver,numberNodes_,newSolution,newSolutionValue,
     822                                     newSolutionValue,"CbcHeuristicLocalAfterFPump");
     823    if ((returnCode&2)!=0) {
     824      // could add cut
     825      returnCode &= ~2;
     826    }
     827    if (returnCode) {
     828      //printf("old sol of %g new of %g\n",saveValue,newSolutionValue);
     829      memcpy(betterSolution,newSolution,numberColumns*sizeof(double));
     830      //abort();
     831      solutionValue=newSolutionValue;
     832      solutionFound=true;
     833    }
     834    delete newSolver;
     835  }
     836  delete [] usedColumn;
     837  delete [] lastSolution;
    304838  delete [] newSolution;
    305   for ( j=0;j<NUMBER_OLD;j++)
    306     delete [] oldSolution[j];
    307   delete [] oldSolution;
    308   delete [] saveObjective;
    309   return returnCode;
     839  delete [] closestSolution;
     840  delete [] integerVariable;
     841  return finalReturnCode;
    310842}
    311843
     
    322854*/
    323855int
    324 CbcHeuristicFPump::rounds(OsiSolverInterface * solver, double * solution,
     856CbcHeuristicFPump::rounds(OsiSolverInterface * solver,double * solution,
    325857                          const double * objective,
     858                          int numberIntegers, const int * integerVariable,
    326859                          bool roundExpensive, double downValue, int *flip)
    327860{
     
    329862  double primalTolerance ;
    330863  solver->getDblParam(OsiPrimalTolerance,primalTolerance) ;
    331   int numberIntegers = model_->numberIntegers();
    332   const int * integerVariable = model_->integerVariable();
    333864
    334865  int i;
     
    351882  for (i=0;i<numberIntegers;i++) {
    352883    int iColumn = integerVariable[i];
    353     if(!solver->isBinary(iColumn))
    354       continue;
    355 #ifndef NDEBUG
    356     const CbcObject * object = model_->object(i);
    357     const CbcSimpleInteger * integerObject =
    358       dynamic_cast<const  CbcSimpleInteger *> (object);
    359     assert(integerObject);
    360 #endif
    361884    double value=solution[iColumn];
    362885    double round = floor(value+primalTolerance);
    363     if (value-round > .5) round += 1.;
     886    if (value-round > downValue) round += 1.;
    364887    if (round < integerTolerance && tmp[iColumn] < -1. + integerTolerance) flip_down++;
    365888    if (round > 1. - integerTolerance && tmp[iColumn] > 1. - integerTolerance) flip_up++;
     
    387910  delete [] tmp;
    388911
     912  const double * columnLower = solver->getColLower();
     913  const double * columnUpper = solver->getColUpper();
    389914  if (*flip == 0 && iter != 0) {
    390915    if(model_->logLevel())
    391916      printf(" -- rand = %4d (%4d) ", nnv, nn);
    392      for (i = 0; i < nnv; i++) solution[list[i]] = 1. - solution[list[i]];
     917     for (i = 0; i < nnv; i++) {
     918       // was solution[list[i]] = 1. - solution[list[i]]; but does that work for 7>=x>=6
     919       int index = list[i];
     920       double value = solution[index];
     921       if (value<=1.0) {
     922         solution[index] = 1.0-value;
     923       } else if (value<columnLower[index]+integerTolerance) {
     924         solution[index] = value+1.0;
     925       } else if (value>columnUpper[index]-integerTolerance) {
     926         solution[index] = value-1.0;
     927       } else {
     928         solution[index] = value-1.0;
     929       }
     930     }
    393931     *flip = nnv;
    394932  } else if (model_->logLevel()) {
Note: See TracChangeset for help on using the changeset viewer.