source: trunk/Samples/CbcHeuristicUser.cpp @ 63

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

mistake

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