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

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

include cstdlib before cmath to get things to compile on AIX with xlC

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