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

Last change on this file since 310 was 310, checked in by andreasw, 13 years ago

first commit for autotools conversion to be able to move more files

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