source: trunk/Samples/CbcHeuristicUser.cpp @ 155

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

return if no solution

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 16.9 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  if (model_->getNumCols()>1000&&model_->getNumCols()>
193      10*model_->getNumRows())
194    return 0; // probably not worth it
195  // worth trying
196  numberSolutions_=model_->getSolutionCount();
197
198  OsiSolverInterface * solver = model_->solver();
199  const double * rowLower = solver->getRowLower();
200  const double * rowUpper = solver->getRowUpper();
201  const double * solution = model_->bestSolution();
202  if (!solution)
203    return 0; // No solution found yet
204  const double * objective = solver->getObjCoefficients();
205  double primalTolerance;
206  solver->getDblParam(OsiPrimalTolerance,primalTolerance);
207
208  int numberRows = matrix_.getNumRows();
209
210  int numberIntegers = model_->numberIntegers();
211  const int * integerVariable = model_->integerVariable();
212 
213  int i;
214  double direction = solver->getObjSense();
215  double newSolutionValue = model_->getObjValue()*direction;
216  int returnCode = 0;
217
218  // Column copy
219  const double * element = matrix_.getElements();
220  const int * row = matrix_.getIndices();
221  const CoinBigIndex * columnStart = matrix_.getVectorStarts();
222  const int * columnLength = matrix_.getVectorLengths();
223
224  // Get solution array for heuristic solution
225  int numberColumns = solver->getNumCols();
226  double * newSolution = new double [numberColumns];
227  memcpy(newSolution,solution,numberColumns*sizeof(double));
228
229  // way is 1 if down possible, 2 if up possible, 3 if both possible
230  char * way = new char[numberIntegers];
231  // corrected costs
232  double * cost = new double[numberIntegers];
233  // for array to mark infeasible rows after iColumn branch
234  char * mark = new char[numberRows];
235  memset(mark,0,numberRows);
236  // space to save values so we don't introduce rounding errors
237  double * save = new double[numberRows];
238
239  // clean solution
240  for (i=0;i<numberIntegers;i++) {
241    int iColumn = integerVariable[i];
242    const CbcObject * object = model_->object(i);
243    const CbcSimpleInteger * integerObject = 
244      dynamic_cast<const  CbcSimpleInteger *> (object);
245    assert(integerObject);
246    // get original bounds
247    double originalLower = integerObject->originalLowerBound();
248    double originalUpper = integerObject->originalUpperBound();
249
250    double value=newSolution[iColumn];
251    if (value<originalLower) {
252      value=originalLower;
253      newSolution[iColumn]=value;
254    } else if (value>originalUpper) {
255      value=originalUpper;
256      newSolution[iColumn]=value;
257    }
258    double nearest=floor(value+0.5);
259    //assert(fabs(value-nearest)<10.0*primalTolerance);
260    value=nearest;
261    newSolution[iColumn]=nearest;
262    // if away from lower bound mark that fact
263    if (nearest>originalLower) {
264      used_[iColumn]=1;
265    }
266    cost[i] = direction*objective[iColumn];
267    int iway=0;
268   
269    if (value>originalLower+0.5) 
270      iway = 1;
271    if (value<originalUpper-0.5) 
272      iway |= 2;
273    way[i]=iway;
274  }
275  // get row activities
276  double * rowActivity = new double[numberRows];
277  memset(rowActivity,0,numberRows*sizeof(double));
278
279  for (i=0;i<numberColumns;i++) {
280    int j;
281    double value = newSolution[i];
282    if (value) {
283      for (j=columnStart[i];
284           j<columnStart[i]+columnLength[i];j++) {
285        int iRow=row[j];
286        rowActivity[iRow] += value*element[j];
287      }
288    }
289  }
290  // check was feasible - if not adjust (cleaning may move)
291  // if very infeasible then give up
292  bool tryHeuristic=true;
293  for (i=0;i<numberRows;i++) {
294    if(rowActivity[i]<rowLower[i]) {
295      if (rowActivity[i]<rowLower[i]-10.0*primalTolerance)
296        tryHeuristic=false;
297      rowActivity[i]=rowLower[i];
298    } else if(rowActivity[i]>rowUpper[i]) {
299      if (rowActivity[i]<rowUpper[i]+10.0*primalTolerance)
300        tryHeuristic=false;
301      rowActivity[i]=rowUpper[i];
302    }
303  }
304  if (tryHeuristic) {
305   
306    // best change in objective
307    double bestChange=0.0;
308   
309    for (i=0;i<numberIntegers;i++) {
310      int iColumn = integerVariable[i];
311     
312      double objectiveCoefficient = cost[i];
313      int k;
314      int j;
315      int goodK=-1;
316      int wayK=-1,wayI=-1;
317      if ((way[i]&1)!=0) {
318        int numberInfeasible=0;
319        // save row activities and adjust
320        for (j=columnStart[iColumn];
321             j<columnStart[iColumn]+columnLength[iColumn];j++) {
322          int iRow = row[j];
323          save[iRow]=rowActivity[iRow];
324          rowActivity[iRow] -= element[j];
325          if(rowActivity[iRow]<rowLower[iRow]-primalTolerance||
326             rowActivity[iRow]>rowUpper[iRow]+primalTolerance) {
327            // mark row
328            mark[iRow]=1;
329            numberInfeasible++;
330          }
331        }
332        // try down
333        for (k=i+1;k<numberIntegers;k++) {
334          if ((way[k]&1)!=0) {
335            // try down
336            if (-objectiveCoefficient-cost[k]<bestChange) {
337              // see if feasible down
338              bool good=true;
339              int numberMarked=0;
340              int kColumn = integerVariable[k];
341              for (j=columnStart[kColumn];
342                   j<columnStart[kColumn]+columnLength[kColumn];j++) {
343                int iRow = row[j];
344                double newValue = rowActivity[iRow] - element[j];
345                if(newValue<rowLower[iRow]-primalTolerance||
346                   newValue>rowUpper[iRow]+primalTolerance) {
347                  good=false;
348                  break;
349                } else if (mark[iRow]) {
350                  // made feasible
351                  numberMarked++;
352                }
353              }
354              if (good&&numberMarked==numberInfeasible) {
355                // better solution
356                goodK=k;
357                wayK=-1;
358                wayI=-1;
359                bestChange = -objectiveCoefficient-cost[k];
360              }
361            }
362          }
363          if ((way[k]&2)!=0) {
364            // try up
365            if (-objectiveCoefficient+cost[k]<bestChange) {
366              // see if feasible up
367              bool good=true;
368              int numberMarked=0;
369              int kColumn = integerVariable[k];
370              for (j=columnStart[kColumn];
371                   j<columnStart[kColumn]+columnLength[kColumn];j++) {
372                int iRow = row[j];
373                double newValue = rowActivity[iRow] + element[j];
374                if(newValue<rowLower[iRow]-primalTolerance||
375                   newValue>rowUpper[iRow]+primalTolerance) {
376                  good=false;
377                  break;
378                } else if (mark[iRow]) {
379                  // made feasible
380                  numberMarked++;
381                }
382              }
383              if (good&&numberMarked==numberInfeasible) {
384                // better solution
385                goodK=k;
386                wayK=1;
387                wayI=-1;
388                bestChange = -objectiveCoefficient+cost[k];
389              }
390            }
391          }
392        }
393        // restore row activities
394        for (j=columnStart[iColumn];
395             j<columnStart[iColumn]+columnLength[iColumn];j++) {
396          int iRow = row[j];
397          rowActivity[iRow] = save[iRow];
398          mark[iRow]=0;
399        }
400      }
401      if ((way[i]&2)!=0) {
402        int numberInfeasible=0;
403        // save row activities and adjust
404        for (j=columnStart[iColumn];
405             j<columnStart[iColumn]+columnLength[iColumn];j++) {
406          int iRow = row[j];
407          save[iRow]=rowActivity[iRow];
408          rowActivity[iRow] += element[j];
409          if(rowActivity[iRow]<rowLower[iRow]-primalTolerance||
410             rowActivity[iRow]>rowUpper[iRow]+primalTolerance) {
411            // mark row
412            mark[iRow]=1;
413            numberInfeasible++;
414          }
415        }
416        // try up
417        for (k=i+1;k<numberIntegers;k++) {
418          if ((way[k]&1)!=0) {
419            // try down
420            if (objectiveCoefficient-cost[k]<bestChange) {
421              // see if feasible down
422              bool good=true;
423              int numberMarked=0;
424              int kColumn = integerVariable[k];
425              for (j=columnStart[kColumn];
426                   j<columnStart[kColumn]+columnLength[kColumn];j++) {
427                int iRow = row[j];
428                double newValue = rowActivity[iRow] - element[j];
429                if(newValue<rowLower[iRow]-primalTolerance||
430                   newValue>rowUpper[iRow]+primalTolerance) {
431                  good=false;
432                  break;
433                } else if (mark[iRow]) {
434                  // made feasible
435                  numberMarked++;
436                }
437              }
438              if (good&&numberMarked==numberInfeasible) {
439                // better solution
440                goodK=k;
441                wayK=-1;
442                wayI=1;
443                bestChange = objectiveCoefficient-cost[k];
444              }
445            }
446          }
447          if ((way[k]&2)!=0) {
448            // try up
449            if (objectiveCoefficient+cost[k]<bestChange) {
450              // see if feasible up
451              bool good=true;
452              int numberMarked=0;
453              int kColumn = integerVariable[k];
454              for (j=columnStart[kColumn];
455                   j<columnStart[kColumn]+columnLength[kColumn];j++) {
456                int iRow = row[j];
457                double newValue = rowActivity[iRow] + element[j];
458                if(newValue<rowLower[iRow]-primalTolerance||
459                   newValue>rowUpper[iRow]+primalTolerance) {
460                  good=false;
461                  break;
462                } else if (mark[iRow]) {
463                  // made feasible
464                  numberMarked++;
465                }
466              }
467              if (good&&numberMarked==numberInfeasible) {
468                // better solution
469                goodK=k;
470                wayK=1;
471                wayI=1;
472                bestChange = objectiveCoefficient+cost[k];
473              }
474            }
475          }
476        }
477        // restore row activities
478        for (j=columnStart[iColumn];
479             j<columnStart[iColumn]+columnLength[iColumn];j++) {
480          int iRow = row[j];
481          rowActivity[iRow] = save[iRow];
482          mark[iRow]=0;
483        }
484      }
485      if (goodK>=0) {
486        // we found something - update solution
487        for (j=columnStart[iColumn];
488             j<columnStart[iColumn]+columnLength[iColumn];j++) {
489          int iRow = row[j];
490          rowActivity[iRow]  += wayI * element[j];
491        }
492        newSolution[iColumn] += wayI;
493        int kColumn = integerVariable[goodK];
494        for (j=columnStart[kColumn];
495             j<columnStart[kColumn]+columnLength[kColumn];j++) {
496          int iRow = row[j];
497          rowActivity[iRow]  += wayK * element[j];
498        }
499        newSolution[kColumn] += wayK;
500        // See if k can go further ?
501        const CbcObject * object = model_->object(goodK);
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[kColumn];
509        int iway=0;
510       
511        if (value>originalLower+0.5) 
512          iway = 1;
513        if (value<originalUpper-0.5) 
514          iway |= 2;
515        way[goodK]=iway;
516      }
517    }
518    if (bestChange+newSolutionValue<solutionValue) {
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      int numberBad=0;
534      double sumBad=0.0;
535      // check was approximately feasible
536      for (i=0;i<numberRows;i++) {
537        if(rowActivity[i]<rowLower[i]) {
538          sumBad += rowLower[i]-rowActivity[i];
539          if (rowActivity[i]<rowLower[i]-10.0*primalTolerance)
540            numberBad++;
541        } else if(rowActivity[i]>rowUpper[i]) {
542          sumBad += rowUpper[i]-rowActivity[i];
543          if (rowActivity[i]>rowUpper[i]+10.0*primalTolerance)
544            numberBad++;
545        }
546      }
547      if (!numberBad) {
548        for (i=0;i<numberIntegers;i++) {
549          int iColumn = integerVariable[i];
550          const CbcObject * object = model_->object(i);
551          const CbcSimpleInteger * integerObject = 
552            dynamic_cast<const  CbcSimpleInteger *> (object);
553          // get original bounds
554          double originalLower = integerObject->originalLowerBound();
555          //double originalUpper = integerObject->originalUpperBound();
556         
557          double value=newSolution[iColumn];
558          // if away from lower bound mark that fact
559          if (value>originalLower) {
560            used_[iColumn]=1;
561          }
562        }
563        // new solution
564        memcpy(betterSolution,newSolution,numberColumns*sizeof(double));
565        returnCode=1;
566        solutionValue = newSolutionValue + bestChange;
567        if (bestChange>1.0e-12)
568          printf("Local search heuristic improved solution by %g\n",
569                 -bestChange);
570      } else {
571        // bad solution - should not happen so debug if see message
572        printf("Local search got bad solution with %d infeasibilities summing to %g\n",
573               numberBad,sumBad);
574      }
575    }
576  }
577  delete [] newSolution;
578  delete [] rowActivity;
579  delete [] way;
580  delete [] cost;
581  delete [] save;
582  delete [] mark;
583  if (numberSolutions_>1&&swap_==1) {
584    // try merge
585    int returnCode2=solutionFix( solutionValue, betterSolution,NULL);
586    if (returnCode2)
587      returnCode=1;
588  }
589  return returnCode;
590}
591// update model
592void CbcLocalSearch::setModel(CbcModel * model)
593{
594  model_ = model;
595  // Get a copy of original matrix
596  assert(model_->solver());
597  matrix_ = *model_->solver()->getMatrixByCol();
598  delete [] used_;
599  int numberColumns = model->solver()->getNumCols();
600  used_ = new char[numberColumns];
601  memset(used_,0,numberColumns);
602}
603
604 
Note: See TracBrowser for help on using the repository browser.