source: trunk/Samples/CbcHeuristicUser.cpp @ 71

Last change on this file since 71 was 71, checked in by forrest, 16 years ago

when_

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