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

Last change on this file since 417 was 356, checked in by ladanyi, 14 years ago

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