source: stable/2.4/Cbc/src/CbcHeuristicLocal.cpp @ 1271

Last change on this file since 1271 was 1271, checked in by forrest, 10 years ago

Creating new stable branch 2.4 from trunk (rev 1270)

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 30.3 KB
Line 
1/* $Id: CbcHeuristicLocal.cpp 1271 2009-11-05 15:57:25Z forrest $ */
2// Copyright (C) 2002, International Business Machines
3// Corporation and others.  All Rights Reserved.
4#if defined(_MSC_VER)
5// Turn off compiler warning about long names
6#  pragma warning(disable:4786)
7#endif
8#include <cassert>
9#include <cstdlib>
10#include <cmath>
11#include <cfloat>
12
13#include "OsiSolverInterface.hpp"
14#include "CbcModel.hpp"
15#include "CbcMessage.hpp"
16#include "CbcHeuristicLocal.hpp"
17#include "CbcBranchActual.hpp"
18#include "CbcStrategy.hpp"
19#include "CglPreProcess.hpp"
20
21// Default Constructor
22CbcHeuristicLocal::CbcHeuristicLocal() 
23  :CbcHeuristic()
24{
25  numberSolutions_=0;
26  swap_=0;
27  used_=NULL;
28}
29
30// Constructor with model - assumed before cuts
31
32CbcHeuristicLocal::CbcHeuristicLocal(CbcModel & model)
33  :CbcHeuristic(model)
34{
35  numberSolutions_=0;
36  swap_=0;
37  // Get a copy of original matrix
38  assert(model.solver());
39  if (model.solver()->getNumRows()) {
40    matrix_ = *model.solver()->getMatrixByCol();
41  }
42  int numberColumns = model.solver()->getNumCols();
43  used_ = new int[numberColumns];
44  memset(used_,0,numberColumns*sizeof(int));
45}
46
47// Destructor
48CbcHeuristicLocal::~CbcHeuristicLocal ()
49{
50  delete [] used_;
51}
52
53// Clone
54CbcHeuristic *
55CbcHeuristicLocal::clone() const
56{
57  return new CbcHeuristicLocal(*this);
58}
59// Create C++ lines to get to current state
60void 
61CbcHeuristicLocal::generateCpp( FILE * fp) 
62{
63  CbcHeuristicLocal other;
64  fprintf(fp,"0#include \"CbcHeuristicLocal.hpp\"\n");
65  fprintf(fp,"3  CbcHeuristicLocal heuristicLocal(*cbcModel);\n");
66  CbcHeuristic::generateCpp(fp,"heuristicLocal");
67  if (swap_!=other.swap_)
68    fprintf(fp,"3  heuristicLocal.setSearchType(%d);\n",swap_);
69  else
70    fprintf(fp,"4  heuristicLocal.setSearchType(%d);\n",swap_);
71  fprintf(fp,"3  cbcModel->addHeuristic(&heuristicLocal);\n");
72}
73
74// Copy constructor
75CbcHeuristicLocal::CbcHeuristicLocal(const CbcHeuristicLocal & rhs)
76:
77  CbcHeuristic(rhs),
78  matrix_(rhs.matrix_),
79  numberSolutions_(rhs.numberSolutions_),
80  swap_(rhs.swap_)
81{
82  if (model_&&rhs.used_) {
83    int numberColumns = model_->solver()->getNumCols();
84    used_ = CoinCopyOfArray(rhs.used_,numberColumns);
85  } else {
86    used_=NULL;
87  }
88}
89
90// Assignment operator
91CbcHeuristicLocal & 
92CbcHeuristicLocal::operator=( const CbcHeuristicLocal& rhs)
93{
94  if (this!=&rhs) {
95    CbcHeuristic::operator=(rhs);
96    matrix_ = rhs.matrix_;
97    numberSolutions_ = rhs.numberSolutions_;
98    swap_ = rhs.swap_;
99    delete [] used_;
100    if (model_&&rhs.used_) {
101      int numberColumns = model_->solver()->getNumCols();
102      used_ = CoinCopyOfArray(rhs.used_,numberColumns);
103    } else {
104      used_=NULL;
105    }
106  }
107  return *this;
108}
109
110// Resets stuff if model changes
111void 
112CbcHeuristicLocal::resetModel(CbcModel * /*model*/)
113{
114  //CbcHeuristic::resetModel(model);
115  delete [] used_;
116  if (model_&&used_) {
117    int numberColumns = model_->solver()->getNumCols();
118    used_ = new int[numberColumns];
119    memset(used_,0,numberColumns*sizeof(int));
120  } else {
121    used_=NULL;
122  }
123}
124// This version fixes stuff and does IP
125int 
126CbcHeuristicLocal::solutionFix(double & objectiveValue,
127                            double * newSolution,
128                               const int * /*keep*/)
129{
130  numCouldRun_++;
131  // See if to do
132  if (!when()||(when()==1&&model_->phase()!=1))
133    return 0; // switched off
134  // Don't do if it was this heuristic which found solution!
135  if (this==model_->lastHeuristic())
136    return 0;
137  OsiSolverInterface * newSolver = model_->continuousSolver()->clone();
138  const double * colLower = newSolver->getColLower();
139  //const double * colUpper = newSolver->getColUpper();
140
141  int numberIntegers = model_->numberIntegers();
142  const int * integerVariable = model_->integerVariable();
143 
144  int i;
145  int nFix=0;
146  for (i=0;i<numberIntegers;i++) {
147    int iColumn=integerVariable[i];
148    const OsiObject * object = model_->object(i);
149    // get original bounds
150    double originalLower;
151    double originalUpper;
152    getIntegerInformation( object,originalLower, originalUpper); 
153    newSolver->setColLower(iColumn,CoinMax(colLower[iColumn],originalLower));
154    if (!used_[iColumn]) {
155      newSolver->setColUpper(iColumn,colLower[iColumn]);
156      nFix++;
157    }
158  }
159  int returnCode = 0;
160  if (nFix*10>numberIntegers) {
161    returnCode = smallBranchAndBound(newSolver,numberNodes_,newSolution,objectiveValue,
162                                     objectiveValue,"CbcHeuristicLocal");
163    if (returnCode<0) {
164      returnCode=0; // returned on size
165      int numberColumns=newSolver->getNumCols();
166      int numberContinuous = numberColumns-numberIntegers;
167      if (numberContinuous>2*numberIntegers&&
168          nFix*10<numberColumns) {
169#define LOCAL_FIX_CONTINUOUS
170#ifdef LOCAL_FIX_CONTINUOUS
171        //const double * colUpper = newSolver->getColUpper();
172        const double * colLower = newSolver->getColLower();
173        int nAtLb=0;
174        //double sumDj=0.0;
175        const double * dj = newSolver->getReducedCost();
176        double direction=newSolver->getObjSense();
177        for (int iColumn=0;iColumn<numberColumns;iColumn++) {
178          if (!newSolver->isInteger(iColumn)) {
179            if (!used_[iColumn]) {
180              //double djValue = dj[iColumn]*direction;
181              nAtLb++;
182              //sumDj += djValue;
183            }
184          }
185        }
186        if (nAtLb) {
187          // fix some continuous
188          double * sort = new double[nAtLb];
189          int * which = new int [nAtLb];
190          //double threshold = CoinMax((0.01*sumDj)/static_cast<double>(nAtLb),1.0e-6);
191          int nFix2=0;
192          for (int iColumn=0;iColumn<numberColumns;iColumn++) {
193            if (!newSolver->isInteger(iColumn)) {
194              if (!used_[iColumn]) {
195                double djValue = dj[iColumn]*direction;
196                if (djValue>1.0e-6) {
197                  sort[nFix2]=-djValue;
198                  which[nFix2++]=iColumn;
199                }
200              }
201            }
202          }
203          CoinSort_2(sort,sort+nFix2,which);
204          int divisor=2;
205          nFix2=CoinMin(nFix2,(numberColumns-nFix)/divisor);
206          for (int i=0;i<nFix2;i++) {
207            int iColumn = which[i];
208            newSolver->setColUpper(iColumn,colLower[iColumn]);
209          }
210          delete [] sort;
211          delete [] which;
212#ifdef CLP_INVESTIGATE2
213          printf("%d integers have zero value, and %d continuous fixed at lb\n",
214                 nFix,nFix2);
215#endif
216          returnCode = smallBranchAndBound(newSolver,
217                                           numberNodes_,newSolution,
218                                           objectiveValue,
219                                           objectiveValue,"CbcHeuristicLocal");
220          if (returnCode<0) 
221            returnCode=0; // returned on size
222        }
223#endif
224      }
225    }
226  }
227  if ((returnCode&2)!=0) {
228    // could add cut
229    returnCode &= ~2;
230  }
231
232  delete newSolver;
233  return returnCode;
234}
235/*
236  First tries setting a variable to better value.  If feasible then
237  tries setting others.  If not feasible then tries swaps
238  Returns 1 if solution, 0 if not */
239int
240CbcHeuristicLocal::solution(double & solutionValue,
241                         double * betterSolution)
242{
243
244  numCouldRun_++;
245  if (numberSolutions_==model_->getSolutionCount())
246    return 0;
247  numberSolutions_=model_->getSolutionCount();
248  if ((model_->getNumCols()>100000&&model_->getNumCols()>
249       10*model_->getNumRows())||numberSolutions_<=1)
250    return 0; // probably not worth it
251  // worth trying
252
253  OsiSolverInterface * solver = model_->solver();
254  const double * rowLower = solver->getRowLower();
255  const double * rowUpper = solver->getRowUpper();
256  const double * solution = model_->bestSolution();
257  if (!solution)
258    return 0; // No solution found yet
259  const double * objective = solver->getObjCoefficients();
260  double primalTolerance;
261  solver->getDblParam(OsiPrimalTolerance,primalTolerance);
262
263  int numberRows = matrix_.getNumRows();
264
265  int numberIntegers = model_->numberIntegers();
266  const int * integerVariable = model_->integerVariable();
267 
268  int i;
269  double direction = solver->getObjSense();
270  double newSolutionValue = model_->getObjValue()*direction;
271  int returnCode = 0;
272  numRuns_++;
273  // Column copy
274  const double * element = matrix_.getElements();
275  const int * row = matrix_.getIndices();
276  const CoinBigIndex * columnStart = matrix_.getVectorStarts();
277  const int * columnLength = matrix_.getVectorLengths();
278
279  // Get solution array for heuristic solution
280  int numberColumns = solver->getNumCols();
281  double * newSolution = new double [numberColumns];
282  memcpy(newSolution,solution,numberColumns*sizeof(double));
283#ifdef LOCAL_FIX_CONTINUOUS
284  // mark continuous used
285  const double * columnLower = solver->getColLower();
286  for (int iColumn=0;iColumn<numberColumns;iColumn++) {
287    if (!solver->isInteger(iColumn)) {
288      if (solution[iColumn]>columnLower[iColumn]+1.0e-8)
289        used_[iColumn]=numberSolutions_;
290    }
291  }
292#endif
293
294  // way is 1 if down possible, 2 if up possible, 3 if both possible
295  char * way = new char[numberIntegers];
296  // corrected costs
297  double * cost = new double[numberIntegers];
298  // for array to mark infeasible rows after iColumn branch
299  char * mark = new char[numberRows];
300  memset(mark,0,numberRows);
301  // space to save values so we don't introduce rounding errors
302  double * save = new double[numberRows];
303
304  // clean solution
305  for (i=0;i<numberIntegers;i++) {
306    int iColumn = integerVariable[i];
307    const OsiObject * object = model_->object(i);
308    // get original bounds
309    double originalLower;
310    double originalUpper;
311    getIntegerInformation( object,originalLower, originalUpper); 
312    double value=newSolution[iColumn];
313    if (value<originalLower) {
314      value=originalLower;
315      newSolution[iColumn]=value;
316    } else if (value>originalUpper) {
317      value=originalUpper;
318      newSolution[iColumn]=value;
319    }
320    double nearest=floor(value+0.5);
321    //assert(fabs(value-nearest)<10.0*primalTolerance);
322    value=nearest;
323    newSolution[iColumn]=nearest;
324    // if away from lower bound mark that fact
325    if (nearest>originalLower) {
326      used_[iColumn]=numberSolutions_;
327    }
328    cost[i] = direction*objective[iColumn];
329    int iway=0;
330   
331    if (value>originalLower+0.5) 
332      iway = 1;
333    if (value<originalUpper-0.5) 
334      iway |= 2;
335    way[i]=static_cast<char>(iway);
336  }
337  // get row activities
338  double * rowActivity = new double[numberRows];
339  memset(rowActivity,0,numberRows*sizeof(double));
340
341  for (i=0;i<numberColumns;i++) {
342    int j;
343    double value = newSolution[i];
344    if (value) {
345      for (j=columnStart[i];
346           j<columnStart[i]+columnLength[i];j++) {
347        int iRow=row[j];
348        rowActivity[iRow] += value*element[j];
349      }
350    }
351  }
352  // check was feasible - if not adjust (cleaning may move)
353  // if very infeasible then give up
354  bool tryHeuristic=true;
355  for (i=0;i<numberRows;i++) {
356    if(rowActivity[i]<rowLower[i]) {
357      if (rowActivity[i]<rowLower[i]-10.0*primalTolerance)
358        tryHeuristic=false;
359      rowActivity[i]=rowLower[i];
360    } else if(rowActivity[i]>rowUpper[i]) {
361      if (rowActivity[i]<rowUpper[i]+10.0*primalTolerance)
362        tryHeuristic=false;
363      rowActivity[i]=rowUpper[i];
364    }
365  }
366  // Switch off if may take too long
367  if (model_->getNumCols()>10000&&model_->getNumCols()>
368      10*model_->getNumRows())
369    tryHeuristic=false;
370  if (tryHeuristic) {
371   
372    // best change in objective
373    double bestChange=0.0;
374   
375    for (i=0;i<numberIntegers;i++) {
376      int iColumn = integerVariable[i];
377     
378      double objectiveCoefficient = cost[i];
379      int k;
380      int j;
381      int goodK=-1;
382      int wayK=-1,wayI=-1;
383      if ((way[i]&1)!=0) {
384        int numberInfeasible=0;
385        // save row activities and adjust
386        for (j=columnStart[iColumn];
387             j<columnStart[iColumn]+columnLength[iColumn];j++) {
388          int iRow = row[j];
389          save[iRow]=rowActivity[iRow];
390          rowActivity[iRow] -= element[j];
391          if(rowActivity[iRow]<rowLower[iRow]-primalTolerance||
392             rowActivity[iRow]>rowUpper[iRow]+primalTolerance) {
393            // mark row
394            mark[iRow]=1;
395            numberInfeasible++;
396          }
397        }
398        // try down
399        for (k=i+1;k<numberIntegers;k++) {
400          if ((way[k]&1)!=0) {
401            // try down
402            if (-objectiveCoefficient-cost[k]<bestChange) {
403              // see if feasible down
404              bool good=true;
405              int numberMarked=0;
406              int kColumn = integerVariable[k];
407              for (j=columnStart[kColumn];
408                   j<columnStart[kColumn]+columnLength[kColumn];j++) {
409                int iRow = row[j];
410                double newValue = rowActivity[iRow] - element[j];
411                if(newValue<rowLower[iRow]-primalTolerance||
412                   newValue>rowUpper[iRow]+primalTolerance) {
413                  good=false;
414                  break;
415                } else if (mark[iRow]) {
416                  // made feasible
417                  numberMarked++;
418                }
419              }
420              if (good&&numberMarked==numberInfeasible) {
421                // better solution
422                goodK=k;
423                wayK=-1;
424                wayI=-1;
425                bestChange = -objectiveCoefficient-cost[k];
426              }
427            }
428          }
429          if ((way[k]&2)!=0) {
430            // try up
431            if (-objectiveCoefficient+cost[k]<bestChange) {
432              // see if feasible up
433              bool good=true;
434              int numberMarked=0;
435              int kColumn = integerVariable[k];
436              for (j=columnStart[kColumn];
437                   j<columnStart[kColumn]+columnLength[kColumn];j++) {
438                int iRow = row[j];
439                double newValue = rowActivity[iRow] + element[j];
440                if(newValue<rowLower[iRow]-primalTolerance||
441                   newValue>rowUpper[iRow]+primalTolerance) {
442                  good=false;
443                  break;
444                } else if (mark[iRow]) {
445                  // made feasible
446                  numberMarked++;
447                }
448              }
449              if (good&&numberMarked==numberInfeasible) {
450                // better solution
451                goodK=k;
452                wayK=1;
453                wayI=-1;
454                bestChange = -objectiveCoefficient+cost[k];
455              }
456            }
457          }
458        }
459        // restore row activities
460        for (j=columnStart[iColumn];
461             j<columnStart[iColumn]+columnLength[iColumn];j++) {
462          int iRow = row[j];
463          rowActivity[iRow] = save[iRow];
464          mark[iRow]=0;
465        }
466      }
467      if ((way[i]&2)!=0) {
468        int numberInfeasible=0;
469        // save row activities and adjust
470        for (j=columnStart[iColumn];
471             j<columnStart[iColumn]+columnLength[iColumn];j++) {
472          int iRow = row[j];
473          save[iRow]=rowActivity[iRow];
474          rowActivity[iRow] += element[j];
475          if(rowActivity[iRow]<rowLower[iRow]-primalTolerance||
476             rowActivity[iRow]>rowUpper[iRow]+primalTolerance) {
477            // mark row
478            mark[iRow]=1;
479            numberInfeasible++;
480          }
481        }
482        // try up
483        for (k=i+1;k<numberIntegers;k++) {
484          if ((way[k]&1)!=0) {
485            // try down
486            if (objectiveCoefficient-cost[k]<bestChange) {
487              // see if feasible down
488              bool good=true;
489              int numberMarked=0;
490              int kColumn = integerVariable[k];
491              for (j=columnStart[kColumn];
492                   j<columnStart[kColumn]+columnLength[kColumn];j++) {
493                int iRow = row[j];
494                double newValue = rowActivity[iRow] - element[j];
495                if(newValue<rowLower[iRow]-primalTolerance||
496                   newValue>rowUpper[iRow]+primalTolerance) {
497                  good=false;
498                  break;
499                } else if (mark[iRow]) {
500                  // made feasible
501                  numberMarked++;
502                }
503              }
504              if (good&&numberMarked==numberInfeasible) {
505                // better solution
506                goodK=k;
507                wayK=-1;
508                wayI=1;
509                bestChange = objectiveCoefficient-cost[k];
510              }
511            }
512          }
513          if ((way[k]&2)!=0) {
514            // try up
515            if (objectiveCoefficient+cost[k]<bestChange) {
516              // see if feasible up
517              bool good=true;
518              int numberMarked=0;
519              int kColumn = integerVariable[k];
520              for (j=columnStart[kColumn];
521                   j<columnStart[kColumn]+columnLength[kColumn];j++) {
522                int iRow = row[j];
523                double newValue = rowActivity[iRow] + element[j];
524                if(newValue<rowLower[iRow]-primalTolerance||
525                   newValue>rowUpper[iRow]+primalTolerance) {
526                  good=false;
527                  break;
528                } else if (mark[iRow]) {
529                  // made feasible
530                  numberMarked++;
531                }
532              }
533              if (good&&numberMarked==numberInfeasible) {
534                // better solution
535                goodK=k;
536                wayK=1;
537                wayI=1;
538                bestChange = objectiveCoefficient+cost[k];
539              }
540            }
541          }
542        }
543        // restore row activities
544        for (j=columnStart[iColumn];
545             j<columnStart[iColumn]+columnLength[iColumn];j++) {
546          int iRow = row[j];
547          rowActivity[iRow] = save[iRow];
548          mark[iRow]=0;
549        }
550      }
551      if (goodK>=0) {
552        // we found something - update solution
553        for (j=columnStart[iColumn];
554             j<columnStart[iColumn]+columnLength[iColumn];j++) {
555          int iRow = row[j];
556          rowActivity[iRow]  += wayI * element[j];
557        }
558        newSolution[iColumn] += wayI;
559        int kColumn = integerVariable[goodK];
560        for (j=columnStart[kColumn];
561             j<columnStart[kColumn]+columnLength[kColumn];j++) {
562          int iRow = row[j];
563          rowActivity[iRow]  += wayK * element[j];
564        }
565        newSolution[kColumn] += wayK;
566        // See if k can go further ?
567        const OsiObject * object = model_->object(goodK);
568        // get original bounds
569        double originalLower;
570        double originalUpper;
571        getIntegerInformation( object,originalLower, originalUpper); 
572       
573        double value=newSolution[kColumn];
574        int iway=0;
575       
576        if (value>originalLower+0.5) 
577          iway = 1;
578        if (value<originalUpper-0.5) 
579          iway |= 2;
580        way[goodK]=static_cast<char>(iway);
581      }
582    }
583    if (bestChange+newSolutionValue<solutionValue) {
584      // paranoid check
585      memset(rowActivity,0,numberRows*sizeof(double));
586     
587      for (i=0;i<numberColumns;i++) {
588        int j;
589        double value = newSolution[i];
590        if (value) {
591          for (j=columnStart[i];
592               j<columnStart[i]+columnLength[i];j++) {
593            int iRow=row[j];
594            rowActivity[iRow] += value*element[j];
595          }
596        }
597      }
598      int numberBad=0;
599      double sumBad=0.0;
600      // check was approximately feasible
601      for (i=0;i<numberRows;i++) {
602        if(rowActivity[i]<rowLower[i]) {
603          sumBad += rowLower[i]-rowActivity[i];
604          if (rowActivity[i]<rowLower[i]-10.0*primalTolerance)
605            numberBad++;
606        } else if(rowActivity[i]>rowUpper[i]) {
607          sumBad += rowUpper[i]-rowActivity[i];
608          if (rowActivity[i]>rowUpper[i]+10.0*primalTolerance)
609            numberBad++;
610        }
611      }
612      if (!numberBad) {
613        for (i=0;i<numberIntegers;i++) {
614          int iColumn = integerVariable[i];
615          const OsiObject * object = model_->object(i);
616          // get original bounds
617          double originalLower;
618          double originalUpper;
619          getIntegerInformation( object,originalLower, originalUpper); 
620         
621          double value=newSolution[iColumn];
622          // if away from lower bound mark that fact
623          if (value>originalLower) {
624            used_[iColumn]=numberSolutions_;
625          }
626        }
627        // new solution
628        memcpy(betterSolution,newSolution,numberColumns*sizeof(double));
629        CoinWarmStartBasis * basis =
630          dynamic_cast<CoinWarmStartBasis *>(solver->getWarmStart()) ;
631        if (basis) {
632          model_->setBestSolutionBasis(* basis);
633          delete basis;
634        }
635        returnCode=1;
636        solutionValue = newSolutionValue + bestChange;
637      } else {
638        // bad solution - should not happen so debug if see message
639        printf("Local search got bad solution with %d infeasibilities summing to %g\n",
640               numberBad,sumBad);
641      }
642    }
643  }
644  delete [] newSolution;
645  delete [] rowActivity;
646  delete [] way;
647  delete [] cost;
648  delete [] save;
649  delete [] mark;
650  if (numberSolutions_>1&&swap_==1) {
651    // try merge
652    int returnCode2=solutionFix( solutionValue, betterSolution,NULL);
653    if (returnCode2)
654      returnCode=1;
655  }
656  return returnCode;
657}
658// update model
659void CbcHeuristicLocal::setModel(CbcModel * model)
660{
661  model_ = model;
662  // Get a copy of original matrix
663  assert(model_->solver());
664  if (model_->solver()->getNumRows()) {
665    matrix_ = *model_->solver()->getMatrixByCol();
666  }
667  delete [] used_;
668  int numberColumns = model->solver()->getNumCols();
669  used_ = new int[numberColumns];
670  memset(used_,0,numberColumns*sizeof(int));
671}
672// Default Constructor
673CbcHeuristicNaive::CbcHeuristicNaive() 
674  :CbcHeuristic()
675{
676  large_=1.0e6;
677}
678
679// Constructor with model - assumed before cuts
680
681CbcHeuristicNaive::CbcHeuristicNaive(CbcModel & model)
682  :CbcHeuristic(model)
683{
684  large_=1.0e6;
685}
686
687// Destructor
688CbcHeuristicNaive::~CbcHeuristicNaive ()
689{
690}
691
692// Clone
693CbcHeuristic *
694CbcHeuristicNaive::clone() const
695{
696  return new CbcHeuristicNaive(*this);
697}
698// Create C++ lines to get to current state
699void 
700CbcHeuristicNaive::generateCpp( FILE * fp) 
701{
702  CbcHeuristicNaive other;
703  fprintf(fp,"0#include \"CbcHeuristicLocal.hpp\"\n");
704  fprintf(fp,"3  CbcHeuristicNaive naive(*cbcModel);\n");
705  CbcHeuristic::generateCpp(fp,"naive");
706  if (large_!=other.large_)
707    fprintf(fp,"3  naive.setLarge(%g);\n",large_);
708  else
709    fprintf(fp,"4  naive.setLarge(%g);\n",large_);
710  fprintf(fp,"3  cbcModel->addHeuristic(&naive);\n");
711}
712
713// Copy constructor
714CbcHeuristicNaive::CbcHeuristicNaive(const CbcHeuristicNaive & rhs)
715:
716  CbcHeuristic(rhs),
717  large_(rhs.large_)
718{
719}
720
721// Assignment operator
722CbcHeuristicNaive & 
723CbcHeuristicNaive::operator=( const CbcHeuristicNaive& rhs)
724{
725  if (this!=&rhs) {
726    CbcHeuristic::operator=(rhs);
727    large_ = rhs.large_;
728  }
729  return *this;
730}
731
732// Resets stuff if model changes
733void 
734CbcHeuristicNaive::resetModel(CbcModel * model)
735{
736  CbcHeuristic::resetModel(model);
737}
738int
739CbcHeuristicNaive::solution(double & solutionValue,
740                         double * betterSolution)
741{
742  numCouldRun_++;
743  // See if to do
744  bool atRoot = model_->getNodeCount()==0;
745  int passNumber = model_->getCurrentPassNumber();
746  if (!when()||(when()==1&&model_->phase()!=1)||!atRoot||passNumber!=1)
747    return 0; // switched off
748  // Don't do if it was this heuristic which found solution!
749  if (this==model_->lastHeuristic())
750    return 0;
751  numRuns_++;
752  double cutoff;
753  model_->solver()->getDblParam(OsiDualObjectiveLimit,cutoff);
754  double direction = model_->solver()->getObjSense();
755  cutoff *= direction;
756  cutoff = CoinMin(cutoff,solutionValue);
757  OsiSolverInterface * solver = model_->continuousSolver();
758  if (!solver)
759    solver = model_->solver();
760  const double * colLower = solver->getColLower();
761  const double * colUpper = solver->getColUpper();
762  const double * objective = solver->getObjCoefficients();
763
764  int numberColumns = model_->getNumCols();
765  int numberIntegers = model_->numberIntegers();
766  const int * integerVariable = model_->integerVariable();
767 
768  int i;
769  bool solutionFound=false;
770  CoinWarmStartBasis saveBasis;
771  CoinWarmStartBasis * basis =
772    dynamic_cast<CoinWarmStartBasis *>(solver->getWarmStart()) ;
773  if (basis) {
774    saveBasis = * basis;
775    delete basis;
776  }
777  // First just fix all integers as close to zero as possible
778  OsiSolverInterface * newSolver = cloneBut(7); // wassolver->clone();
779  for (i=0;i<numberIntegers;i++) {
780    int iColumn=integerVariable[i];
781    double lower = colLower[iColumn];
782    double upper = colUpper[iColumn];
783    double value;
784    if (lower>0.0)
785      value=lower;
786    else if (upper<0.0)
787      value=upper;
788    else
789      value=0.0;
790    newSolver->setColLower(iColumn,value);
791    newSolver->setColUpper(iColumn,value);
792  }
793  newSolver->initialSolve();
794  if (newSolver->isProvenOptimal()) {
795    double solValue = newSolver->getObjValue()*direction ;
796    if (solValue<cutoff) {
797      // we have a solution
798      solutionFound=true;
799      solutionValue=solValue;
800      memcpy(betterSolution,newSolver->getColSolution(),
801             numberColumns*sizeof(double));
802      printf("Naive fixing close to zero gave solution of %g\n",solutionValue);
803      cutoff = solValue - model_->getCutoffIncrement();
804    }
805  }
806  // Now fix all integers as close to zero if zero or large cost
807  int nFix=0;
808  for (i=0;i<numberIntegers;i++) {
809    int iColumn=integerVariable[i];
810    double lower = colLower[iColumn];
811    double upper = colUpper[iColumn];
812    double value;
813    if (fabs(objective[i])>0.0&&fabs(objective[i])<large_) {
814      nFix++;
815      if (lower>0.0)
816        value=lower;
817      else if (upper<0.0)
818        value=upper;
819      else
820        value=0.0;
821      newSolver->setColLower(iColumn,value);
822      newSolver->setColUpper(iColumn,value);
823    } else {
824      // set back to original
825      newSolver->setColLower(iColumn,lower);
826      newSolver->setColUpper(iColumn,upper);
827    }
828  }
829  const double * solution = solver->getColSolution();
830  if (nFix) {
831    newSolver->setWarmStart(&saveBasis);
832    newSolver->setColSolution(solution);
833    newSolver->initialSolve();
834    if (newSolver->isProvenOptimal()) {
835      double solValue = newSolver->getObjValue()*direction ;
836      if (solValue<cutoff) {
837        // try branch and bound
838        double * newSolution = new double [numberColumns];
839        printf("%d fixed after fixing costs\n",nFix);
840        int returnCode = smallBranchAndBound(newSolver,
841                                             numberNodes_,newSolution,
842                                             solutionValue,
843                                             solutionValue,"CbcHeuristicNaive1");
844        if (returnCode<0)
845          returnCode=0; // returned on size
846        if ((returnCode&2)!=0) {
847          // could add cut
848          returnCode &= ~2;
849        }
850        if (returnCode==1) {
851          // solution
852          solutionFound=true;
853          memcpy(betterSolution,newSolution,
854                 numberColumns*sizeof(double));
855          printf("Naive fixing zeros gave solution of %g\n",solutionValue);
856          cutoff = solutionValue - model_->getCutoffIncrement();
857        }
858        delete [] newSolution;
859      }
860    }
861  }
862#if 1
863  newSolver->setObjSense(-direction); // maximize
864  newSolver->setWarmStart(&saveBasis);
865  newSolver->setColSolution(solution);
866  for (int iColumn=0;iColumn<numberColumns;iColumn++) {
867    double value = solution[iColumn];
868    double lower = colLower[iColumn];
869    double upper = colUpper[iColumn];
870    double newLower;
871    double newUpper;
872    if (newSolver->isInteger(iColumn)) {
873      newLower = CoinMax(lower,floor(value)-2.0);
874      newUpper = CoinMin(upper,ceil(value)+2.0);
875    } else {
876      newLower = CoinMax(lower,value-1.0e5);
877      newUpper = CoinMin(upper,value+1.0e-5);
878    }
879    newSolver->setColLower(iColumn,newLower);
880    newSolver->setColUpper(iColumn,newUpper);
881  }
882  newSolver->initialSolve();
883  if (newSolver->isProvenOptimal()) {
884    double solValue = newSolver->getObjValue()*direction ;
885    if (solValue<cutoff) {
886      nFix=0;
887      newSolver->setObjSense(direction); // correct direction
888      //const double * thisSolution = newSolver->getColSolution();
889      for (int iColumn=0;iColumn<numberColumns;iColumn++) {
890        double value = solution[iColumn];
891        double lower = colLower[iColumn];
892        double upper = colUpper[iColumn];
893        double newLower=lower;
894        double newUpper=upper;
895        if (newSolver->isInteger(iColumn)) {
896          if (value<lower+1.0e-6) {
897            nFix++;
898            newUpper=lower;
899          } else if (value>upper-1.0e-6) {
900            nFix++;
901            newLower=upper;
902          } else {
903            newLower = CoinMax(lower,floor(value)-2.0);
904            newUpper = CoinMin(upper,ceil(value)+2.0);
905          }
906        }
907        newSolver->setColLower(iColumn,newLower);
908        newSolver->setColUpper(iColumn,newUpper);
909      }
910      // try branch and bound
911      double * newSolution = new double [numberColumns];
912      printf("%d fixed after maximizing\n",nFix);
913      int returnCode = smallBranchAndBound(newSolver,
914                                           numberNodes_,newSolution,
915                                           solutionValue,
916                                           solutionValue,"CbcHeuristicNaive1");
917      if (returnCode<0)
918        returnCode=0; // returned on size
919      if ((returnCode&2)!=0) {
920        // could add cut
921        returnCode &= ~2;
922      }
923      if (returnCode==1) {
924        // solution
925        solutionFound=true;
926        memcpy(betterSolution,newSolution,
927               numberColumns*sizeof(double));
928        printf("Naive maximizing gave solution of %g\n",solutionValue);
929        cutoff = solutionValue - model_->getCutoffIncrement();
930      }
931      delete [] newSolution;
932    }
933  }
934#endif
935  delete newSolver;
936  return solutionFound ? 1 : 0;
937}
938// update model
939void CbcHeuristicNaive::setModel(CbcModel * model)
940{
941  model_ = model;
942}
943// Default Constructor
944CbcHeuristicCrossover::CbcHeuristicCrossover() 
945  :CbcHeuristic(),
946   numberSolutions_(0),
947   useNumber_(3)
948{
949  setWhen(1);
950}
951
952// Constructor with model - assumed before cuts
953
954CbcHeuristicCrossover::CbcHeuristicCrossover(CbcModel & model)
955  :CbcHeuristic(model),
956   numberSolutions_(0),
957   useNumber_(3)
958{
959  setWhen(1);
960  for (int i=0;i<10;i++)
961    random_[i]=model.randomNumberGenerator()->randomDouble();
962}
963
964// Destructor
965CbcHeuristicCrossover::~CbcHeuristicCrossover ()
966{
967}
968
969// Clone
970CbcHeuristic *
971CbcHeuristicCrossover::clone() const
972{
973  return new CbcHeuristicCrossover(*this);
974}
975// Create C++ lines to get to current state
976void 
977CbcHeuristicCrossover::generateCpp( FILE * fp) 
978{
979  CbcHeuristicCrossover other;
980  fprintf(fp,"0#include \"CbcHeuristicLocal.hpp\"\n");
981  fprintf(fp,"3  CbcHeuristicCrossover crossover(*cbcModel);\n");
982  CbcHeuristic::generateCpp(fp,"crossover");
983  if (useNumber_!=other.useNumber_)
984    fprintf(fp,"3  crossover.setNumberSolutions(%d);\n",useNumber_);
985  else
986    fprintf(fp,"4  crossover.setNumberSolutions(%d);\n",useNumber_);
987  fprintf(fp,"3  cbcModel->addHeuristic(&crossover);\n");
988}
989
990// Copy constructor
991CbcHeuristicCrossover::CbcHeuristicCrossover(const CbcHeuristicCrossover & rhs)
992:
993  CbcHeuristic(rhs),
994  attempts_(rhs.attempts_),
995  numberSolutions_(rhs.numberSolutions_),
996  useNumber_(rhs.useNumber_)
997{
998  memcpy(random_,rhs.random_,10*sizeof(double));
999}
1000
1001// Assignment operator
1002CbcHeuristicCrossover & 
1003CbcHeuristicCrossover::operator=( const CbcHeuristicCrossover& rhs)
1004{
1005  if (this!=&rhs) {
1006    CbcHeuristic::operator=(rhs);
1007    useNumber_ = rhs.useNumber_;
1008    attempts_ = rhs.attempts_;
1009    numberSolutions_ = rhs.numberSolutions_;
1010    memcpy(random_,rhs.random_,10*sizeof(double));
1011  }
1012  return *this;
1013}
1014
1015// Resets stuff if model changes
1016void 
1017CbcHeuristicCrossover::resetModel(CbcModel * model)
1018{
1019  CbcHeuristic::resetModel(model);
1020}
1021int
1022CbcHeuristicCrossover::solution(double & solutionValue,
1023                         double * betterSolution)
1024{
1025  if (when_==0)
1026    return 0;
1027  numCouldRun_++;
1028  bool useBest=(numberSolutions_!=model_->getSolutionCount());
1029  if (!useBest&&(when_%10)==1)
1030    return 0;
1031  numberSolutions_=model_->getSolutionCount();
1032  OsiSolverInterface * continuousSolver = model_->continuousSolver();
1033  int useNumber =CoinMin(model_->numberSavedSolutions(),useNumber_);
1034  if (useNumber<2||!continuousSolver)
1035    return 0;
1036  // Fix later
1037  if (!useBest)
1038    abort();
1039  numRuns_++;
1040  double cutoff;
1041  model_->solver()->getDblParam(OsiDualObjectiveLimit,cutoff);
1042  double direction = model_->solver()->getObjSense();
1043  cutoff *= direction;
1044  cutoff = CoinMin(cutoff,solutionValue);
1045  OsiSolverInterface * solver = cloneBut(2);
1046  // But reset bounds
1047  solver->setColLower(continuousSolver->getColLower());
1048  solver->setColUpper(continuousSolver->getColUpper());
1049  int numberColumns = solver->getNumCols();
1050  // Fixed
1051  double * fixed =new double [numberColumns];
1052  for (int i=0;i<numberColumns;i++)
1053    fixed[i]=-COIN_DBL_MAX;
1054  int whichSolution[10];
1055  for (int i=0;i<useNumber;i++)
1056    whichSolution[i]=i;
1057  for (int i=0;i<useNumber;i++) {
1058    int k = whichSolution[i];
1059    const double * solution = model_->savedSolution(k);
1060    for (int j=0;j<numberColumns;j++) {
1061      if (solver->isInteger(j)) {
1062        if (fixed[j]==-COIN_DBL_MAX) 
1063          fixed[j]=floor(solution[j]+0.5);
1064        else if (fabs(fixed[j]-solution[j])>1.0e-7)
1065          fixed[j]=COIN_DBL_MAX;
1066      }
1067    }
1068  }
1069  const double * colLower = solver->getColLower();
1070  for (int i=0;i<numberColumns;i++) {
1071    if (solver->isInteger(i)) {
1072      double value=fixed[i];
1073      if (value!=COIN_DBL_MAX) {
1074        if (when_<10) {
1075          solver->setColLower(i,value);
1076          solver->setColUpper(i,value);
1077        } else if (value==colLower[i]) {
1078          solver->setColUpper(i,value);
1079        }
1080      }
1081    }
1082  }
1083  int returnCode = smallBranchAndBound(solver,numberNodes_,betterSolution,
1084                                       solutionValue,
1085                                       solutionValue,"CbcHeuristicCrossover");
1086  if (returnCode<0)
1087    returnCode=0; // returned on size
1088  if ((returnCode&2)!=0) {
1089    // could add cut
1090    returnCode &= ~2;
1091  }
1092
1093  delete solver;
1094  return returnCode;
1095}
1096// update model
1097void CbcHeuristicCrossover::setModel(CbcModel * model)
1098{
1099  model_ = model;
1100  if (model) {
1101    for (int i=0;i<10;i++)
1102      random_[i]=model->randomNumberGenerator()->randomDouble();
1103  }
1104}
1105
1106 
Note: See TracBrowser for help on using the repository browser.