source: stable/2.0/Cbc/src/CbcHeuristicLocal.cpp @ 905

Last change on this file since 905 was 905, checked in by ladanyi, 13 years ago

include cstdlib before cmath to get things to compile on AIX with xlC (same as changeset 904 in trunk)

  • 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.