source: trunk/Cbc/src/CbcHeuristicLocal.cpp @ 862

Last change on this file since 862 was 759, checked in by forrest, 12 years ago

start of work on new vub heuristic

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