source: branches/devel/Cbc/src/CbcHeuristicLocal.cpp @ 491

Last change on this file since 491 was 439, checked in by forrest, 13 years ago

towards common use with other solvers

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 15.6 KB
Line 
1// Copyright (C) 2002, International Business Machines
2// Corporation and others.  All Rights Reserved.
3#if defined(_MSC_VER)
4// Turn off compiler warning about long names
5#  pragma warning(disable:4786)
6#endif
7#include <cassert>
8#include <cmath>
9#include <cfloat>
10
11#include "OsiSolverInterface.hpp"
12#include "CbcModel.hpp"
13#include "CbcMessage.hpp"
14#include "CbcHeuristicLocal.hpp"
15#include "CbcBranchActual.hpp"
16#include "CbcStrategy.hpp"
17#include "CglPreProcess.hpp"
18
19// Default Constructor
20CbcHeuristicLocal::CbcHeuristicLocal() 
21  :CbcHeuristic()
22{
23  numberSolutions_=0;
24  swap_=0;
25  used_=NULL;
26}
27
28// Constructor with model - assumed before cuts
29
30CbcHeuristicLocal::CbcHeuristicLocal(CbcModel & model)
31  :CbcHeuristic(model)
32{
33  numberSolutions_=0;
34  swap_=0;
35  // Get a copy of original matrix
36  assert(model.solver());
37  matrix_ = *model.solver()->getMatrixByCol();
38  int numberColumns = model.solver()->getNumCols();
39  used_ = new char[numberColumns];
40  memset(used_,0,numberColumns);
41}
42
43// Destructor
44CbcHeuristicLocal::~CbcHeuristicLocal ()
45{
46  delete [] used_;
47}
48
49// Clone
50CbcHeuristic *
51CbcHeuristicLocal::clone() const
52{
53  return new CbcHeuristicLocal(*this);
54}
55// Create C++ lines to get to current state
56void 
57CbcHeuristicLocal::generateCpp( FILE * fp) 
58{
59  CbcHeuristicLocal other;
60  fprintf(fp,"0#include \"CbcHeuristicLocal.hpp\"\n");
61  fprintf(fp,"3  CbcHeuristicLocal heuristicLocal(*cbcModel);\n");
62  if (swap_!=other.swap_)
63    fprintf(fp,"3  heuristicLocal.setSearchType(%d);\n",swap_);
64  else
65    fprintf(fp,"4  heuristicLocal.setSearchType(%d);\n",swap_);
66  fprintf(fp,"3  cbcModel->addHeuristic(&heuristicLocal);\n");
67}
68
69// Copy constructor
70CbcHeuristicLocal::CbcHeuristicLocal(const CbcHeuristicLocal & rhs)
71:
72  CbcHeuristic(rhs),
73  matrix_(rhs.matrix_),
74  numberSolutions_(rhs.numberSolutions_),
75  swap_(rhs.swap_)
76{
77  setWhen(rhs.when());
78  if (model_&&rhs.used_) {
79    int numberColumns = model_->solver()->getNumCols();
80    used_ = new char[numberColumns];
81    memcpy(used_,rhs.used_,numberColumns);
82  } else {
83    used_=NULL;
84  }
85}
86// Resets stuff if model changes
87void 
88CbcHeuristicLocal::resetModel(CbcModel * model)
89{
90  //CbcHeuristic::resetModel(model);
91  delete [] used_;
92  if (model_&&used_) {
93    int numberColumns = model_->solver()->getNumCols();
94    used_ = new char[numberColumns];
95    memset(used_,0,numberColumns);
96  } else {
97    used_=NULL;
98  }
99}
100// This version fixes stuff and does IP
101int 
102CbcHeuristicLocal::solutionFix(double & objectiveValue,
103                            double * newSolution,
104                            const int * keep)
105{
106  // See if to do
107  if (!when()||(when()==1&&model_->phase()!=1))
108    return 0; // switched off
109  // Don't do if it was this heuristic which found solution!
110  if (this==model_->lastHeuristic())
111    return 0;
112  OsiSolverInterface * newSolver = model_->continuousSolver()->clone();
113  const double * colLower = newSolver->getColLower();
114  //const double * colUpper = newSolver->getColUpper();
115
116  int numberIntegers = model_->numberIntegers();
117  const int * integerVariable = model_->integerVariable();
118 
119  int i;
120  int nFix=0;
121  for (i=0;i<numberIntegers;i++) {
122    int iColumn=integerVariable[i];
123    const OsiObject * object = model_->object(i);
124    // get original bounds
125    double originalLower;
126    double originalUpper;
127    getIntegerInformation( object,originalLower, originalUpper); 
128    newSolver->setColLower(iColumn,CoinMax(colLower[iColumn],originalLower));
129    if (!used_[iColumn]) {
130      newSolver->setColUpper(iColumn,colLower[iColumn]);
131      nFix++;
132    }
133  }
134  int returnCode = smallBranchAndBound(newSolver,200,newSolution,objectiveValue,
135                                         objectiveValue,"CbcHeuristicLocal");
136
137  delete newSolver;
138  return returnCode;
139}
140/*
141  First tries setting a variable to better value.  If feasible then
142  tries setting others.  If not feasible then tries swaps
143  Returns 1 if solution, 0 if not */
144int
145CbcHeuristicLocal::solution(double & solutionValue,
146                         double * betterSolution)
147{
148
149  if (numberSolutions_==model_->getSolutionCount())
150    return 0;
151  numberSolutions_=model_->getSolutionCount();
152  if (model_->getNumCols()>100000&&model_->getNumCols()>
153      10*model_->getNumRows()||numberSolutions_<=1)
154    return 0; // probably not worth it
155  // worth trying
156
157  OsiSolverInterface * solver = model_->solver();
158  const double * rowLower = solver->getRowLower();
159  const double * rowUpper = solver->getRowUpper();
160  const double * solution = model_->bestSolution();
161  if (!solution)
162    return 0; // No solution found yet
163  const double * objective = solver->getObjCoefficients();
164  double primalTolerance;
165  solver->getDblParam(OsiPrimalTolerance,primalTolerance);
166
167  int numberRows = matrix_.getNumRows();
168
169  int numberIntegers = model_->numberIntegers();
170  const int * integerVariable = model_->integerVariable();
171 
172  int i;
173  double direction = solver->getObjSense();
174  double newSolutionValue = model_->getObjValue()*direction;
175  int returnCode = 0;
176
177  // Column copy
178  const double * element = matrix_.getElements();
179  const int * row = matrix_.getIndices();
180  const CoinBigIndex * columnStart = matrix_.getVectorStarts();
181  const int * columnLength = matrix_.getVectorLengths();
182
183  // Get solution array for heuristic solution
184  int numberColumns = solver->getNumCols();
185  double * newSolution = new double [numberColumns];
186  memcpy(newSolution,solution,numberColumns*sizeof(double));
187
188  // way is 1 if down possible, 2 if up possible, 3 if both possible
189  char * way = new char[numberIntegers];
190  // corrected costs
191  double * cost = new double[numberIntegers];
192  // for array to mark infeasible rows after iColumn branch
193  char * mark = new char[numberRows];
194  memset(mark,0,numberRows);
195  // space to save values so we don't introduce rounding errors
196  double * save = new double[numberRows];
197
198  // clean solution
199  for (i=0;i<numberIntegers;i++) {
200    int iColumn = integerVariable[i];
201    const OsiObject * object = model_->object(i);
202    // get original bounds
203    double originalLower;
204    double originalUpper;
205    getIntegerInformation( object,originalLower, originalUpper); 
206    double value=newSolution[iColumn];
207    if (value<originalLower) {
208      value=originalLower;
209      newSolution[iColumn]=value;
210    } else if (value>originalUpper) {
211      value=originalUpper;
212      newSolution[iColumn]=value;
213    }
214    double nearest=floor(value+0.5);
215    //assert(fabs(value-nearest)<10.0*primalTolerance);
216    value=nearest;
217    newSolution[iColumn]=nearest;
218    // if away from lower bound mark that fact
219    if (nearest>originalLower) {
220      used_[iColumn]=1;
221    }
222    cost[i] = direction*objective[iColumn];
223    int iway=0;
224   
225    if (value>originalLower+0.5) 
226      iway = 1;
227    if (value<originalUpper-0.5) 
228      iway |= 2;
229    way[i]=iway;
230  }
231  // get row activities
232  double * rowActivity = new double[numberRows];
233  memset(rowActivity,0,numberRows*sizeof(double));
234
235  for (i=0;i<numberColumns;i++) {
236    int j;
237    double value = newSolution[i];
238    if (value) {
239      for (j=columnStart[i];
240           j<columnStart[i]+columnLength[i];j++) {
241        int iRow=row[j];
242        rowActivity[iRow] += value*element[j];
243      }
244    }
245  }
246  // check was feasible - if not adjust (cleaning may move)
247  // if very infeasible then give up
248  bool tryHeuristic=true;
249  for (i=0;i<numberRows;i++) {
250    if(rowActivity[i]<rowLower[i]) {
251      if (rowActivity[i]<rowLower[i]-10.0*primalTolerance)
252        tryHeuristic=false;
253      rowActivity[i]=rowLower[i];
254    } else if(rowActivity[i]>rowUpper[i]) {
255      if (rowActivity[i]<rowUpper[i]+10.0*primalTolerance)
256        tryHeuristic=false;
257      rowActivity[i]=rowUpper[i];
258    }
259  }
260  // Switch off if may take too long
261  if (model_->getNumCols()>10000&&model_->getNumCols()>
262      10*model_->getNumRows())
263    tryHeuristic=false;
264  if (tryHeuristic) {
265   
266    // best change in objective
267    double bestChange=0.0;
268   
269    for (i=0;i<numberIntegers;i++) {
270      int iColumn = integerVariable[i];
271     
272      double objectiveCoefficient = cost[i];
273      int k;
274      int j;
275      int goodK=-1;
276      int wayK=-1,wayI=-1;
277      if ((way[i]&1)!=0) {
278        int numberInfeasible=0;
279        // save row activities and adjust
280        for (j=columnStart[iColumn];
281             j<columnStart[iColumn]+columnLength[iColumn];j++) {
282          int iRow = row[j];
283          save[iRow]=rowActivity[iRow];
284          rowActivity[iRow] -= element[j];
285          if(rowActivity[iRow]<rowLower[iRow]-primalTolerance||
286             rowActivity[iRow]>rowUpper[iRow]+primalTolerance) {
287            // mark row
288            mark[iRow]=1;
289            numberInfeasible++;
290          }
291        }
292        // try down
293        for (k=i+1;k<numberIntegers;k++) {
294          if ((way[k]&1)!=0) {
295            // try down
296            if (-objectiveCoefficient-cost[k]<bestChange) {
297              // see if feasible down
298              bool good=true;
299              int numberMarked=0;
300              int kColumn = integerVariable[k];
301              for (j=columnStart[kColumn];
302                   j<columnStart[kColumn]+columnLength[kColumn];j++) {
303                int iRow = row[j];
304                double newValue = rowActivity[iRow] - element[j];
305                if(newValue<rowLower[iRow]-primalTolerance||
306                   newValue>rowUpper[iRow]+primalTolerance) {
307                  good=false;
308                  break;
309                } else if (mark[iRow]) {
310                  // made feasible
311                  numberMarked++;
312                }
313              }
314              if (good&&numberMarked==numberInfeasible) {
315                // better solution
316                goodK=k;
317                wayK=-1;
318                wayI=-1;
319                bestChange = -objectiveCoefficient-cost[k];
320              }
321            }
322          }
323          if ((way[k]&2)!=0) {
324            // try up
325            if (-objectiveCoefficient+cost[k]<bestChange) {
326              // see if feasible up
327              bool good=true;
328              int numberMarked=0;
329              int kColumn = integerVariable[k];
330              for (j=columnStart[kColumn];
331                   j<columnStart[kColumn]+columnLength[kColumn];j++) {
332                int iRow = row[j];
333                double newValue = rowActivity[iRow] + element[j];
334                if(newValue<rowLower[iRow]-primalTolerance||
335                   newValue>rowUpper[iRow]+primalTolerance) {
336                  good=false;
337                  break;
338                } else if (mark[iRow]) {
339                  // made feasible
340                  numberMarked++;
341                }
342              }
343              if (good&&numberMarked==numberInfeasible) {
344                // better solution
345                goodK=k;
346                wayK=1;
347                wayI=-1;
348                bestChange = -objectiveCoefficient+cost[k];
349              }
350            }
351          }
352        }
353        // restore row activities
354        for (j=columnStart[iColumn];
355             j<columnStart[iColumn]+columnLength[iColumn];j++) {
356          int iRow = row[j];
357          rowActivity[iRow] = save[iRow];
358          mark[iRow]=0;
359        }
360      }
361      if ((way[i]&2)!=0) {
362        int numberInfeasible=0;
363        // save row activities and adjust
364        for (j=columnStart[iColumn];
365             j<columnStart[iColumn]+columnLength[iColumn];j++) {
366          int iRow = row[j];
367          save[iRow]=rowActivity[iRow];
368          rowActivity[iRow] += element[j];
369          if(rowActivity[iRow]<rowLower[iRow]-primalTolerance||
370             rowActivity[iRow]>rowUpper[iRow]+primalTolerance) {
371            // mark row
372            mark[iRow]=1;
373            numberInfeasible++;
374          }
375        }
376        // try up
377        for (k=i+1;k<numberIntegers;k++) {
378          if ((way[k]&1)!=0) {
379            // try down
380            if (objectiveCoefficient-cost[k]<bestChange) {
381              // see if feasible down
382              bool good=true;
383              int numberMarked=0;
384              int kColumn = integerVariable[k];
385              for (j=columnStart[kColumn];
386                   j<columnStart[kColumn]+columnLength[kColumn];j++) {
387                int iRow = row[j];
388                double newValue = rowActivity[iRow] - element[j];
389                if(newValue<rowLower[iRow]-primalTolerance||
390                   newValue>rowUpper[iRow]+primalTolerance) {
391                  good=false;
392                  break;
393                } else if (mark[iRow]) {
394                  // made feasible
395                  numberMarked++;
396                }
397              }
398              if (good&&numberMarked==numberInfeasible) {
399                // better solution
400                goodK=k;
401                wayK=-1;
402                wayI=1;
403                bestChange = objectiveCoefficient-cost[k];
404              }
405            }
406          }
407          if ((way[k]&2)!=0) {
408            // try up
409            if (objectiveCoefficient+cost[k]<bestChange) {
410              // see if feasible up
411              bool good=true;
412              int numberMarked=0;
413              int kColumn = integerVariable[k];
414              for (j=columnStart[kColumn];
415                   j<columnStart[kColumn]+columnLength[kColumn];j++) {
416                int iRow = row[j];
417                double newValue = rowActivity[iRow] + element[j];
418                if(newValue<rowLower[iRow]-primalTolerance||
419                   newValue>rowUpper[iRow]+primalTolerance) {
420                  good=false;
421                  break;
422                } else if (mark[iRow]) {
423                  // made feasible
424                  numberMarked++;
425                }
426              }
427              if (good&&numberMarked==numberInfeasible) {
428                // better solution
429                goodK=k;
430                wayK=1;
431                wayI=1;
432                bestChange = objectiveCoefficient+cost[k];
433              }
434            }
435          }
436        }
437        // restore row activities
438        for (j=columnStart[iColumn];
439             j<columnStart[iColumn]+columnLength[iColumn];j++) {
440          int iRow = row[j];
441          rowActivity[iRow] = save[iRow];
442          mark[iRow]=0;
443        }
444      }
445      if (goodK>=0) {
446        // we found something - update solution
447        for (j=columnStart[iColumn];
448             j<columnStart[iColumn]+columnLength[iColumn];j++) {
449          int iRow = row[j];
450          rowActivity[iRow]  += wayI * element[j];
451        }
452        newSolution[iColumn] += wayI;
453        int kColumn = integerVariable[goodK];
454        for (j=columnStart[kColumn];
455             j<columnStart[kColumn]+columnLength[kColumn];j++) {
456          int iRow = row[j];
457          rowActivity[iRow]  += wayK * element[j];
458        }
459        newSolution[kColumn] += wayK;
460        // See if k can go further ?
461        const OsiObject * object = model_->object(goodK);
462        // get original bounds
463        double originalLower;
464        double originalUpper;
465        getIntegerInformation( object,originalLower, originalUpper); 
466       
467        double value=newSolution[kColumn];
468        int iway=0;
469       
470        if (value>originalLower+0.5) 
471          iway = 1;
472        if (value<originalUpper-0.5) 
473          iway |= 2;
474        way[goodK]=iway;
475      }
476    }
477    if (bestChange+newSolutionValue<solutionValue) {
478      // paranoid check
479      memset(rowActivity,0,numberRows*sizeof(double));
480     
481      for (i=0;i<numberColumns;i++) {
482        int j;
483        double value = newSolution[i];
484        if (value) {
485          for (j=columnStart[i];
486               j<columnStart[i]+columnLength[i];j++) {
487            int iRow=row[j];
488            rowActivity[iRow] += value*element[j];
489          }
490        }
491      }
492      int numberBad=0;
493      double sumBad=0.0;
494      // check was approximately feasible
495      for (i=0;i<numberRows;i++) {
496        if(rowActivity[i]<rowLower[i]) {
497          sumBad += rowLower[i]-rowActivity[i];
498          if (rowActivity[i]<rowLower[i]-10.0*primalTolerance)
499            numberBad++;
500        } else if(rowActivity[i]>rowUpper[i]) {
501          sumBad += rowUpper[i]-rowActivity[i];
502          if (rowActivity[i]>rowUpper[i]+10.0*primalTolerance)
503            numberBad++;
504        }
505      }
506      if (!numberBad) {
507        for (i=0;i<numberIntegers;i++) {
508          int iColumn = integerVariable[i];
509          const OsiObject * object = model_->object(i);
510          // get original bounds
511          double originalLower;
512          double originalUpper;
513          getIntegerInformation( object,originalLower, originalUpper); 
514         
515          double value=newSolution[iColumn];
516          // if away from lower bound mark that fact
517          if (value>originalLower) {
518            used_[iColumn]=1;
519          }
520        }
521        // new solution
522        memcpy(betterSolution,newSolution,numberColumns*sizeof(double));
523        returnCode=1;
524        solutionValue = newSolutionValue + bestChange;
525        if (bestChange<-1.0e-1)
526          model_->messageHandler()->message(CBC_HEURISTIC_SOLUTION,model_->messages())
527            << solutionValue
528            << "CbcHeuristicLocal"<<model_->getCurrentSeconds()
529            <<CoinMessageEol;
530      } else {
531        // bad solution - should not happen so debug if see message
532        printf("Local search got bad solution with %d infeasibilities summing to %g\n",
533               numberBad,sumBad);
534      }
535    }
536  }
537  delete [] newSolution;
538  delete [] rowActivity;
539  delete [] way;
540  delete [] cost;
541  delete [] save;
542  delete [] mark;
543  if (numberSolutions_>1&&swap_==1) {
544    // try merge
545    int returnCode2=solutionFix( solutionValue, betterSolution,NULL);
546    if (returnCode2)
547      returnCode=1;
548  }
549  return returnCode;
550}
551// update model
552void CbcHeuristicLocal::setModel(CbcModel * model)
553{
554  model_ = model;
555  // Get a copy of original matrix
556  assert(model_->solver());
557  matrix_ = *model_->solver()->getMatrixByCol();
558  delete [] used_;
559  int numberColumns = model->solver()->getNumCols();
560  used_ = new char[numberColumns];
561  memset(used_,0,numberColumns);
562}
563
564 
Note: See TracBrowser for help on using the repository browser.