source: trunk/Cbc/src/CbcHeuristicLocal.cpp @ 1212

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

fixes

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 28.5 KB
Line 
1/* $Id: CbcHeuristicLocal.cpp 1212 2009-08-21 16:19:13Z 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 char[numberColumns];
44  memset(used_,0,numberColumns);
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_ = new char[numberColumns];
85    memcpy(used_,rhs.used_,numberColumns);
86  } else {
87    used_=NULL;
88  }
89}
90
91// Assignment operator
92CbcHeuristicLocal & 
93CbcHeuristicLocal::operator=( const CbcHeuristicLocal& rhs)
94{
95  if (this!=&rhs) {
96    CbcHeuristic::operator=(rhs);
97    matrix_ = rhs.matrix_;
98    numberSolutions_ = rhs.numberSolutions_;
99    swap_ = rhs.swap_;
100    delete [] used_;
101    if (model_&&rhs.used_) {
102      int numberColumns = model_->solver()->getNumCols();
103      used_ = new char[numberColumns];
104      memcpy(used_,rhs.used_,numberColumns);
105    } else {
106      used_=NULL;
107    }
108  }
109  return *this;
110}
111
112// Resets stuff if model changes
113void 
114CbcHeuristicLocal::resetModel(CbcModel * /*model*/)
115{
116  //CbcHeuristic::resetModel(model);
117  delete [] used_;
118  if (model_&&used_) {
119    int numberColumns = model_->solver()->getNumCols();
120    used_ = new char[numberColumns];
121    memset(used_,0,numberColumns);
122  } else {
123    used_=NULL;
124  }
125}
126// This version fixes stuff and does IP
127int 
128CbcHeuristicLocal::solutionFix(double & objectiveValue,
129                            double * newSolution,
130                               const int * /*keep*/)
131{
132  numCouldRun_++;
133  // See if to do
134  if (!when()||(when()==1&&model_->phase()!=1))
135    return 0; // switched off
136  // Don't do if it was this heuristic which found solution!
137  if (this==model_->lastHeuristic())
138    return 0;
139  OsiSolverInterface * newSolver = model_->continuousSolver()->clone();
140  const double * colLower = newSolver->getColLower();
141  //const double * colUpper = newSolver->getColUpper();
142
143  int numberIntegers = model_->numberIntegers();
144  const int * integerVariable = model_->integerVariable();
145 
146  int i;
147  int nFix=0;
148  for (i=0;i<numberIntegers;i++) {
149    int iColumn=integerVariable[i];
150    const OsiObject * object = model_->object(i);
151    // get original bounds
152    double originalLower;
153    double originalUpper;
154    getIntegerInformation( object,originalLower, originalUpper); 
155    newSolver->setColLower(iColumn,CoinMax(colLower[iColumn],originalLower));
156    if (!used_[iColumn]) {
157      newSolver->setColUpper(iColumn,colLower[iColumn]);
158      nFix++;
159    }
160  }
161  int returnCode = smallBranchAndBound(newSolver,numberNodes_,newSolution,objectiveValue,
162                                         objectiveValue,"CbcHeuristicLocal");
163  if (returnCode<0)
164    returnCode=0; // returned on size
165  if ((returnCode&2)!=0) {
166    // could add cut
167    returnCode &= ~2;
168  }
169
170  delete newSolver;
171  return returnCode;
172}
173/*
174  First tries setting a variable to better value.  If feasible then
175  tries setting others.  If not feasible then tries swaps
176  Returns 1 if solution, 0 if not */
177int
178CbcHeuristicLocal::solution(double & solutionValue,
179                         double * betterSolution)
180{
181
182  numCouldRun_++;
183  if (numberSolutions_==model_->getSolutionCount())
184    return 0;
185  numberSolutions_=model_->getSolutionCount();
186  if ((model_->getNumCols()>100000&&model_->getNumCols()>
187       10*model_->getNumRows())||numberSolutions_<=1)
188    return 0; // probably not worth it
189  // worth trying
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  if (!solution)
196    return 0; // No solution found yet
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  numRuns_++;
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 OsiObject * object = model_->object(i);
236    // get original bounds
237    double originalLower;
238    double originalUpper;
239    getIntegerInformation( object,originalLower, originalUpper); 
240    double value=newSolution[iColumn];
241    if (value<originalLower) {
242      value=originalLower;
243      newSolution[iColumn]=value;
244    } else if (value>originalUpper) {
245      value=originalUpper;
246      newSolution[iColumn]=value;
247    }
248    double nearest=floor(value+0.5);
249    //assert(fabs(value-nearest)<10.0*primalTolerance);
250    value=nearest;
251    newSolution[iColumn]=nearest;
252    // if away from lower bound mark that fact
253    if (nearest>originalLower) {
254      used_[iColumn]=1;
255    }
256    cost[i] = direction*objective[iColumn];
257    int iway=0;
258   
259    if (value>originalLower+0.5) 
260      iway = 1;
261    if (value<originalUpper-0.5) 
262      iway |= 2;
263    way[i]=static_cast<char>(iway);
264  }
265  // get row activities
266  double * rowActivity = new double[numberRows];
267  memset(rowActivity,0,numberRows*sizeof(double));
268
269  for (i=0;i<numberColumns;i++) {
270    int j;
271    double value = newSolution[i];
272    if (value) {
273      for (j=columnStart[i];
274           j<columnStart[i]+columnLength[i];j++) {
275        int iRow=row[j];
276        rowActivity[iRow] += value*element[j];
277      }
278    }
279  }
280  // check was feasible - if not adjust (cleaning may move)
281  // if very infeasible then give up
282  bool tryHeuristic=true;
283  for (i=0;i<numberRows;i++) {
284    if(rowActivity[i]<rowLower[i]) {
285      if (rowActivity[i]<rowLower[i]-10.0*primalTolerance)
286        tryHeuristic=false;
287      rowActivity[i]=rowLower[i];
288    } else if(rowActivity[i]>rowUpper[i]) {
289      if (rowActivity[i]<rowUpper[i]+10.0*primalTolerance)
290        tryHeuristic=false;
291      rowActivity[i]=rowUpper[i];
292    }
293  }
294  // Switch off if may take too long
295  if (model_->getNumCols()>10000&&model_->getNumCols()>
296      10*model_->getNumRows())
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 OsiObject * object = model_->object(goodK);
496        // get original bounds
497        double originalLower;
498        double originalUpper;
499        getIntegerInformation( object,originalLower, originalUpper); 
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]=static_cast<char>(iway);
509      }
510    }
511    if (bestChange+newSolutionValue<solutionValue) {
512      // paranoid check
513      memset(rowActivity,0,numberRows*sizeof(double));
514     
515      for (i=0;i<numberColumns;i++) {
516        int j;
517        double value = newSolution[i];
518        if (value) {
519          for (j=columnStart[i];
520               j<columnStart[i]+columnLength[i];j++) {
521            int iRow=row[j];
522            rowActivity[iRow] += value*element[j];
523          }
524        }
525      }
526      int numberBad=0;
527      double sumBad=0.0;
528      // check was approximately feasible
529      for (i=0;i<numberRows;i++) {
530        if(rowActivity[i]<rowLower[i]) {
531          sumBad += rowLower[i]-rowActivity[i];
532          if (rowActivity[i]<rowLower[i]-10.0*primalTolerance)
533            numberBad++;
534        } else if(rowActivity[i]>rowUpper[i]) {
535          sumBad += rowUpper[i]-rowActivity[i];
536          if (rowActivity[i]>rowUpper[i]+10.0*primalTolerance)
537            numberBad++;
538        }
539      }
540      if (!numberBad) {
541        for (i=0;i<numberIntegers;i++) {
542          int iColumn = integerVariable[i];
543          const OsiObject * object = model_->object(i);
544          // get original bounds
545          double originalLower;
546          double originalUpper;
547          getIntegerInformation( object,originalLower, originalUpper); 
548         
549          double value=newSolution[iColumn];
550          // if away from lower bound mark that fact
551          if (value>originalLower) {
552            used_[iColumn]=1;
553          }
554        }
555        // new solution
556        memcpy(betterSolution,newSolution,numberColumns*sizeof(double));
557        CoinWarmStartBasis * basis =
558          dynamic_cast<CoinWarmStartBasis *>(solver->getWarmStart()) ;
559        if (basis) {
560          model_->setBestSolutionBasis(* basis);
561          delete basis;
562        }
563        returnCode=1;
564        solutionValue = newSolutionValue + bestChange;
565      } else {
566        // bad solution - should not happen so debug if see message
567        printf("Local search got bad solution with %d infeasibilities summing to %g\n",
568               numberBad,sumBad);
569      }
570    }
571  }
572  delete [] newSolution;
573  delete [] rowActivity;
574  delete [] way;
575  delete [] cost;
576  delete [] save;
577  delete [] mark;
578  if (numberSolutions_>1&&swap_==1) {
579    int i;
580    for ( i=0;i<numberColumns;i++) {
581      if (used_[i]>1)
582        break;
583    }
584    if (i==numberColumns) {
585      // modify used_ if just one
586      const int * used = model_->usedInSolution();
587      for (int i=0;i<numberColumns;i++)
588        used_[i]= static_cast<char>(CoinMin(used[i],255));
589    }
590    // try merge
591    int returnCode2=solutionFix( solutionValue, betterSolution,NULL);
592    if (returnCode2)
593      returnCode=1;
594  }
595  return returnCode;
596}
597// update model
598void CbcHeuristicLocal::setModel(CbcModel * model)
599{
600  model_ = model;
601  // Get a copy of original matrix
602  assert(model_->solver());
603  if (model_->solver()->getNumRows()) {
604    matrix_ = *model_->solver()->getMatrixByCol();
605  }
606  delete [] used_;
607  int numberColumns = model->solver()->getNumCols();
608  used_ = new char[numberColumns];
609  memset(used_,0,numberColumns);
610}
611// Default Constructor
612CbcHeuristicNaive::CbcHeuristicNaive() 
613  :CbcHeuristic()
614{
615  large_=1.0e6;
616}
617
618// Constructor with model - assumed before cuts
619
620CbcHeuristicNaive::CbcHeuristicNaive(CbcModel & model)
621  :CbcHeuristic(model)
622{
623  large_=1.0e6;
624}
625
626// Destructor
627CbcHeuristicNaive::~CbcHeuristicNaive ()
628{
629}
630
631// Clone
632CbcHeuristic *
633CbcHeuristicNaive::clone() const
634{
635  return new CbcHeuristicNaive(*this);
636}
637// Create C++ lines to get to current state
638void 
639CbcHeuristicNaive::generateCpp( FILE * fp) 
640{
641  CbcHeuristicNaive other;
642  fprintf(fp,"0#include \"CbcHeuristicLocal.hpp\"\n");
643  fprintf(fp,"3  CbcHeuristicNaive naive(*cbcModel);\n");
644  CbcHeuristic::generateCpp(fp,"naive");
645  if (large_!=other.large_)
646    fprintf(fp,"3  naive.setLarge(%g);\n",large_);
647  else
648    fprintf(fp,"4  naive.setLarge(%g);\n",large_);
649  fprintf(fp,"3  cbcModel->addHeuristic(&naive);\n");
650}
651
652// Copy constructor
653CbcHeuristicNaive::CbcHeuristicNaive(const CbcHeuristicNaive & rhs)
654:
655  CbcHeuristic(rhs),
656  large_(rhs.large_)
657{
658}
659
660// Assignment operator
661CbcHeuristicNaive & 
662CbcHeuristicNaive::operator=( const CbcHeuristicNaive& rhs)
663{
664  if (this!=&rhs) {
665    CbcHeuristic::operator=(rhs);
666    large_ = rhs.large_;
667  }
668  return *this;
669}
670
671// Resets stuff if model changes
672void 
673CbcHeuristicNaive::resetModel(CbcModel * model)
674{
675  CbcHeuristic::resetModel(model);
676}
677int
678CbcHeuristicNaive::solution(double & solutionValue,
679                         double * betterSolution)
680{
681  numCouldRun_++;
682  // See if to do
683  bool atRoot = model_->getNodeCount()==0;
684  int passNumber = model_->getCurrentPassNumber();
685  if (!when()||(when()==1&&model_->phase()!=1)||!atRoot||passNumber!=1)
686    return 0; // switched off
687  // Don't do if it was this heuristic which found solution!
688  if (this==model_->lastHeuristic())
689    return 0;
690  numRuns_++;
691  double cutoff;
692  model_->solver()->getDblParam(OsiDualObjectiveLimit,cutoff);
693  double direction = model_->solver()->getObjSense();
694  cutoff *= direction;
695  cutoff = CoinMin(cutoff,solutionValue);
696  OsiSolverInterface * solver = model_->continuousSolver();
697  if (!solver)
698    solver = model_->solver();
699  const double * colLower = solver->getColLower();
700  const double * colUpper = solver->getColUpper();
701  const double * objective = solver->getObjCoefficients();
702
703  int numberColumns = model_->getNumCols();
704  int numberIntegers = model_->numberIntegers();
705  const int * integerVariable = model_->integerVariable();
706 
707  int i;
708  bool solutionFound=false;
709  CoinWarmStartBasis saveBasis;
710  CoinWarmStartBasis * basis =
711    dynamic_cast<CoinWarmStartBasis *>(solver->getWarmStart()) ;
712  if (basis) {
713    saveBasis = * basis;
714    delete basis;
715  }
716  // First just fix all integers as close to zero as possible
717  OsiSolverInterface * newSolver = cloneBut(7); // wassolver->clone();
718  for (i=0;i<numberIntegers;i++) {
719    int iColumn=integerVariable[i];
720    double lower = colLower[iColumn];
721    double upper = colUpper[iColumn];
722    double value;
723    if (lower>0.0)
724      value=lower;
725    else if (upper<0.0)
726      value=upper;
727    else
728      value=0.0;
729    newSolver->setColLower(iColumn,value);
730    newSolver->setColUpper(iColumn,value);
731  }
732  newSolver->initialSolve();
733  if (newSolver->isProvenOptimal()) {
734    double solValue = newSolver->getObjValue()*direction ;
735    if (solValue<cutoff) {
736      // we have a solution
737      solutionFound=true;
738      solutionValue=solValue;
739      memcpy(betterSolution,newSolver->getColSolution(),
740             numberColumns*sizeof(double));
741      printf("Naive fixing close to zero gave solution of %g\n",solutionValue);
742      cutoff = solValue - model_->getCutoffIncrement();
743    }
744  }
745  // Now fix all integers as close to zero if zero or large cost
746  int nFix=0;
747  for (i=0;i<numberIntegers;i++) {
748    int iColumn=integerVariable[i];
749    double lower = colLower[iColumn];
750    double upper = colUpper[iColumn];
751    double value;
752    if (fabs(objective[i])>0.0&&fabs(objective[i])<large_) {
753      nFix++;
754      if (lower>0.0)
755        value=lower;
756      else if (upper<0.0)
757        value=upper;
758      else
759        value=0.0;
760      newSolver->setColLower(iColumn,value);
761      newSolver->setColUpper(iColumn,value);
762    } else {
763      // set back to original
764      newSolver->setColLower(iColumn,lower);
765      newSolver->setColUpper(iColumn,upper);
766    }
767  }
768  const double * solution = solver->getColSolution();
769  if (nFix) {
770    newSolver->setWarmStart(&saveBasis);
771    newSolver->setColSolution(solution);
772    newSolver->initialSolve();
773    if (newSolver->isProvenOptimal()) {
774      double solValue = newSolver->getObjValue()*direction ;
775      if (solValue<cutoff) {
776        // try branch and bound
777        double * newSolution = new double [numberColumns];
778        printf("%d fixed after fixing costs\n",nFix);
779        int returnCode = smallBranchAndBound(newSolver,
780                                             numberNodes_,newSolution,
781                                             solutionValue,
782                                             solutionValue,"CbcHeuristicNaive1");
783        if (returnCode<0)
784          returnCode=0; // returned on size
785        if ((returnCode&2)!=0) {
786          // could add cut
787          returnCode &= ~2;
788        }
789        if (returnCode==1) {
790          // solution
791          solutionFound=true;
792          memcpy(betterSolution,newSolution,
793                 numberColumns*sizeof(double));
794          printf("Naive fixing zeros gave solution of %g\n",solutionValue);
795          cutoff = solutionValue - model_->getCutoffIncrement();
796        }
797        delete [] newSolution;
798      }
799    }
800  }
801#if 1
802  newSolver->setObjSense(-direction); // maximize
803  newSolver->setWarmStart(&saveBasis);
804  newSolver->setColSolution(solution);
805  for (int iColumn=0;iColumn<numberColumns;iColumn++) {
806    double value = solution[iColumn];
807    double lower = colLower[iColumn];
808    double upper = colUpper[iColumn];
809    double newLower;
810    double newUpper;
811    if (newSolver->isInteger(iColumn)) {
812      newLower = CoinMax(lower,floor(value)-2.0);
813      newUpper = CoinMin(upper,ceil(value)+2.0);
814    } else {
815      newLower = CoinMax(lower,value-1.0e5);
816      newUpper = CoinMin(upper,value+1.0e-5);
817    }
818    newSolver->setColLower(iColumn,newLower);
819    newSolver->setColUpper(iColumn,newUpper);
820  }
821  newSolver->initialSolve();
822  if (newSolver->isProvenOptimal()) {
823    double solValue = newSolver->getObjValue()*direction ;
824    if (solValue<cutoff) {
825      nFix=0;
826      newSolver->setObjSense(direction); // correct direction
827      //const double * thisSolution = newSolver->getColSolution();
828      for (int iColumn=0;iColumn<numberColumns;iColumn++) {
829        double value = solution[iColumn];
830        double lower = colLower[iColumn];
831        double upper = colUpper[iColumn];
832        double newLower=lower;
833        double newUpper=upper;
834        if (newSolver->isInteger(iColumn)) {
835          if (value<lower+1.0e-6) {
836            nFix++;
837            newUpper=lower;
838          } else if (value>upper-1.0e-6) {
839            nFix++;
840            newLower=upper;
841          } else {
842            newLower = CoinMax(lower,floor(value)-2.0);
843            newUpper = CoinMin(upper,ceil(value)+2.0);
844          }
845        }
846        newSolver->setColLower(iColumn,newLower);
847        newSolver->setColUpper(iColumn,newUpper);
848      }
849      // try branch and bound
850      double * newSolution = new double [numberColumns];
851      printf("%d fixed after maximizing\n",nFix);
852      int returnCode = smallBranchAndBound(newSolver,
853                                           numberNodes_,newSolution,
854                                           solutionValue,
855                                           solutionValue,"CbcHeuristicNaive1");
856      if (returnCode<0)
857        returnCode=0; // returned on size
858      if ((returnCode&2)!=0) {
859        // could add cut
860        returnCode &= ~2;
861      }
862      if (returnCode==1) {
863        // solution
864        solutionFound=true;
865        memcpy(betterSolution,newSolution,
866               numberColumns*sizeof(double));
867        printf("Naive maximizing gave solution of %g\n",solutionValue);
868        cutoff = solutionValue - model_->getCutoffIncrement();
869      }
870      delete [] newSolution;
871    }
872  }
873#endif
874  delete newSolver;
875  return solutionFound ? 1 : 0;
876}
877// update model
878void CbcHeuristicNaive::setModel(CbcModel * model)
879{
880  model_ = model;
881}
882// Default Constructor
883CbcHeuristicCrossover::CbcHeuristicCrossover() 
884  :CbcHeuristic(),
885   numberSolutions_(0),
886   useNumber_(3)
887{
888  setWhen(1);
889}
890
891// Constructor with model - assumed before cuts
892
893CbcHeuristicCrossover::CbcHeuristicCrossover(CbcModel & model)
894  :CbcHeuristic(model),
895   numberSolutions_(0),
896   useNumber_(3)
897{
898  setWhen(1);
899  for (int i=0;i<10;i++)
900    random_[i]=model.randomNumberGenerator()->randomDouble();
901}
902
903// Destructor
904CbcHeuristicCrossover::~CbcHeuristicCrossover ()
905{
906}
907
908// Clone
909CbcHeuristic *
910CbcHeuristicCrossover::clone() const
911{
912  return new CbcHeuristicCrossover(*this);
913}
914// Create C++ lines to get to current state
915void 
916CbcHeuristicCrossover::generateCpp( FILE * fp) 
917{
918  CbcHeuristicCrossover other;
919  fprintf(fp,"0#include \"CbcHeuristicLocal.hpp\"\n");
920  fprintf(fp,"3  CbcHeuristicCrossover crossover(*cbcModel);\n");
921  CbcHeuristic::generateCpp(fp,"crossover");
922  if (useNumber_!=other.useNumber_)
923    fprintf(fp,"3  crossover.setNumberSolutions(%d);\n",useNumber_);
924  else
925    fprintf(fp,"4  crossover.setNumberSolutions(%d);\n",useNumber_);
926  fprintf(fp,"3  cbcModel->addHeuristic(&crossover);\n");
927}
928
929// Copy constructor
930CbcHeuristicCrossover::CbcHeuristicCrossover(const CbcHeuristicCrossover & rhs)
931:
932  CbcHeuristic(rhs),
933  attempts_(rhs.attempts_),
934  numberSolutions_(rhs.numberSolutions_),
935  useNumber_(rhs.useNumber_)
936{
937  memcpy(random_,rhs.random_,10*sizeof(double));
938}
939
940// Assignment operator
941CbcHeuristicCrossover & 
942CbcHeuristicCrossover::operator=( const CbcHeuristicCrossover& rhs)
943{
944  if (this!=&rhs) {
945    CbcHeuristic::operator=(rhs);
946    useNumber_ = rhs.useNumber_;
947    attempts_ = rhs.attempts_;
948    numberSolutions_ = rhs.numberSolutions_;
949    memcpy(random_,rhs.random_,10*sizeof(double));
950  }
951  return *this;
952}
953
954// Resets stuff if model changes
955void 
956CbcHeuristicCrossover::resetModel(CbcModel * model)
957{
958  CbcHeuristic::resetModel(model);
959}
960int
961CbcHeuristicCrossover::solution(double & solutionValue,
962                         double * betterSolution)
963{
964  if (when_==0)
965    return 0;
966  numCouldRun_++;
967  bool useBest=(numberSolutions_!=model_->getSolutionCount());
968  if (!useBest&&(when_%10)==1)
969    return 0;
970  numberSolutions_=model_->getSolutionCount();
971  OsiSolverInterface * continuousSolver = model_->continuousSolver();
972  int useNumber =CoinMin(model_->numberSavedSolutions(),useNumber_);
973  if (useNumber<2||!continuousSolver)
974    return 0;
975  // Fix later
976  if (!useBest)
977    abort();
978  numRuns_++;
979  double cutoff;
980  model_->solver()->getDblParam(OsiDualObjectiveLimit,cutoff);
981  double direction = model_->solver()->getObjSense();
982  cutoff *= direction;
983  cutoff = CoinMin(cutoff,solutionValue);
984  OsiSolverInterface * solver = cloneBut(2);
985  // But reset bounds
986  solver->setColLower(continuousSolver->getColLower());
987  solver->setColUpper(continuousSolver->getColUpper());
988  int numberColumns = solver->getNumCols();
989  // Fixed
990  double * fixed =new double [numberColumns];
991  for (int i=0;i<numberColumns;i++)
992    fixed[i]=-COIN_DBL_MAX;
993  int whichSolution[10];
994  for (int i=0;i<useNumber;i++)
995    whichSolution[i]=i;
996  for (int i=0;i<useNumber;i++) {
997    int k = whichSolution[i];
998    const double * solution = model_->savedSolution(k);
999    for (int j=0;j<numberColumns;j++) {
1000      if (solver->isInteger(j)) {
1001        if (fixed[j]==-COIN_DBL_MAX) 
1002          fixed[j]=floor(solution[j]+0.5);
1003        else if (fabs(fixed[j]-solution[j])>1.0e-7)
1004          fixed[j]=COIN_DBL_MAX;
1005      }
1006    }
1007  }
1008  const double * colLower = solver->getColLower();
1009  for (int i=0;i<numberColumns;i++) {
1010    if (solver->isInteger(i)) {
1011      double value=fixed[i];
1012      if (value!=COIN_DBL_MAX) {
1013        if (when_<10) {
1014          solver->setColLower(i,value);
1015          solver->setColUpper(i,value);
1016        } else if (value==colLower[i]) {
1017          solver->setColUpper(i,value);
1018        }
1019      }
1020    }
1021  }
1022  int returnCode = smallBranchAndBound(solver,numberNodes_,betterSolution,
1023                                       solutionValue,
1024                                       solutionValue,"CbcHeuristicCrossover");
1025  if (returnCode<0)
1026    returnCode=0; // returned on size
1027  if ((returnCode&2)!=0) {
1028    // could add cut
1029    returnCode &= ~2;
1030  }
1031
1032  delete solver;
1033  return returnCode;
1034}
1035// update model
1036void CbcHeuristicCrossover::setModel(CbcModel * model)
1037{
1038  model_ = model;
1039  if (model) {
1040    for (int i=0;i<10;i++)
1041      random_[i]=model->randomNumberGenerator()->randomDouble();
1042  }
1043}
1044
1045 
Note: See TracBrowser for help on using the repository browser.