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

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

update branches/devel for threads

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