source: trunk/Samples/CbcHeuristicUserB.cpp @ 124

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

more examples

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