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

Last change on this file since 1103 was 1103, checked in by forrest, 11 years ago

compiler warnings

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 23.7 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 <cstdlib>
9#include <cmath>
10#include <cfloat>
11
12#include "OsiSolverInterface.hpp"
13#include "CbcModel.hpp"
14#include "CbcMessage.hpp"
15#include "CbcHeuristicLocal.hpp"
16#include "CbcBranchActual.hpp"
17#include "CbcStrategy.hpp"
18#include "CglPreProcess.hpp"
19
20// Default Constructor
21CbcHeuristicLocal::CbcHeuristicLocal() 
22  :CbcHeuristic()
23{
24  numberSolutions_=0;
25  swap_=0;
26  used_=NULL;
27}
28
29// Constructor with model - assumed before cuts
30
31CbcHeuristicLocal::CbcHeuristicLocal(CbcModel & model)
32  :CbcHeuristic(model)
33{
34  numberSolutions_=0;
35  swap_=0;
36  // Get a copy of original matrix
37  assert(model.solver());
38  matrix_ = *model.solver()->getMatrixByCol();
39  int numberColumns = model.solver()->getNumCols();
40  used_ = new char[numberColumns];
41  memset(used_,0,numberColumns);
42}
43
44// Destructor
45CbcHeuristicLocal::~CbcHeuristicLocal ()
46{
47  delete [] used_;
48}
49
50// Clone
51CbcHeuristic *
52CbcHeuristicLocal::clone() const
53{
54  return new CbcHeuristicLocal(*this);
55}
56// Create C++ lines to get to current state
57void 
58CbcHeuristicLocal::generateCpp( FILE * fp) 
59{
60  CbcHeuristicLocal other;
61  fprintf(fp,"0#include \"CbcHeuristicLocal.hpp\"\n");
62  fprintf(fp,"3  CbcHeuristicLocal heuristicLocal(*cbcModel);\n");
63  CbcHeuristic::generateCpp(fp,"heuristicLocal");
64  if (swap_!=other.swap_)
65    fprintf(fp,"3  heuristicLocal.setSearchType(%d);\n",swap_);
66  else
67    fprintf(fp,"4  heuristicLocal.setSearchType(%d);\n",swap_);
68  fprintf(fp,"3  cbcModel->addHeuristic(&heuristicLocal);\n");
69}
70
71// Copy constructor
72CbcHeuristicLocal::CbcHeuristicLocal(const CbcHeuristicLocal & rhs)
73:
74  CbcHeuristic(rhs),
75  matrix_(rhs.matrix_),
76  numberSolutions_(rhs.numberSolutions_),
77  swap_(rhs.swap_)
78{
79  if (model_&&rhs.used_) {
80    int numberColumns = model_->solver()->getNumCols();
81    used_ = new char[numberColumns];
82    memcpy(used_,rhs.used_,numberColumns);
83  } else {
84    used_=NULL;
85  }
86}
87
88// Assignment operator
89CbcHeuristicLocal & 
90CbcHeuristicLocal::operator=( const CbcHeuristicLocal& rhs)
91{
92  if (this!=&rhs) {
93    CbcHeuristic::operator=(rhs);
94    matrix_ = rhs.matrix_;
95    numberSolutions_ = rhs.numberSolutions_;
96    swap_ = rhs.swap_;
97    delete [] used_;
98    if (model_&&rhs.used_) {
99      int numberColumns = model_->solver()->getNumCols();
100      used_ = new char[numberColumns];
101      memcpy(used_,rhs.used_,numberColumns);
102    } else {
103      used_=NULL;
104    }
105  }
106  return *this;
107}
108
109// Resets stuff if model changes
110void 
111CbcHeuristicLocal::resetModel(CbcModel * model)
112{
113  //CbcHeuristic::resetModel(model);
114  delete [] used_;
115  if (model_&&used_) {
116    int numberColumns = model_->solver()->getNumCols();
117    used_ = new char[numberColumns];
118    memset(used_,0,numberColumns);
119  } else {
120    used_=NULL;
121  }
122}
123// This version fixes stuff and does IP
124int 
125CbcHeuristicLocal::solutionFix(double & objectiveValue,
126                            double * newSolution,
127                            const int * keep)
128{
129  numCouldRun_++;
130  // See if to do
131  if (!when()||(when()==1&&model_->phase()!=1))
132    return 0; // switched off
133  // Don't do if it was this heuristic which found solution!
134  if (this==model_->lastHeuristic())
135    return 0;
136  OsiSolverInterface * newSolver = model_->continuousSolver()->clone();
137  const double * colLower = newSolver->getColLower();
138  //const double * colUpper = newSolver->getColUpper();
139
140  int numberIntegers = model_->numberIntegers();
141  const int * integerVariable = model_->integerVariable();
142 
143  int i;
144  int nFix=0;
145  for (i=0;i<numberIntegers;i++) {
146    int iColumn=integerVariable[i];
147    const OsiObject * object = model_->object(i);
148    // get original bounds
149    double originalLower;
150    double originalUpper;
151    getIntegerInformation( object,originalLower, originalUpper); 
152    newSolver->setColLower(iColumn,CoinMax(colLower[iColumn],originalLower));
153    if (!used_[iColumn]) {
154      newSolver->setColUpper(iColumn,colLower[iColumn]);
155      nFix++;
156    }
157  }
158  int returnCode = smallBranchAndBound(newSolver,numberNodes_,newSolution,objectiveValue,
159                                         objectiveValue,"CbcHeuristicLocal");
160  if (returnCode<0)
161    returnCode=0; // returned on size
162  if ((returnCode&2)!=0) {
163    // could add cut
164    returnCode &= ~2;
165  }
166
167  delete newSolver;
168  return returnCode;
169}
170/*
171  First tries setting a variable to better value.  If feasible then
172  tries setting others.  If not feasible then tries swaps
173  Returns 1 if solution, 0 if not */
174int
175CbcHeuristicLocal::solution(double & solutionValue,
176                         double * betterSolution)
177{
178
179  numCouldRun_++;
180  if (numberSolutions_==model_->getSolutionCount())
181    return 0;
182  numberSolutions_=model_->getSolutionCount();
183  if ((model_->getNumCols()>100000&&model_->getNumCols()>
184       10*model_->getNumRows())||numberSolutions_<=1)
185    return 0; // probably not worth it
186  // worth trying
187
188  OsiSolverInterface * solver = model_->solver();
189  const double * rowLower = solver->getRowLower();
190  const double * rowUpper = solver->getRowUpper();
191  const double * solution = model_->bestSolution();
192  if (!solution)
193    return 0; // No solution found yet
194  const double * objective = solver->getObjCoefficients();
195  double primalTolerance;
196  solver->getDblParam(OsiPrimalTolerance,primalTolerance);
197
198  int numberRows = matrix_.getNumRows();
199
200  int numberIntegers = model_->numberIntegers();
201  const int * integerVariable = model_->integerVariable();
202 
203  int i;
204  double direction = solver->getObjSense();
205  double newSolutionValue = model_->getObjValue()*direction;
206  int returnCode = 0;
207  numRuns_++;
208  // Column copy
209  const double * element = matrix_.getElements();
210  const int * row = matrix_.getIndices();
211  const CoinBigIndex * columnStart = matrix_.getVectorStarts();
212  const int * columnLength = matrix_.getVectorLengths();
213
214  // Get solution array for heuristic solution
215  int numberColumns = solver->getNumCols();
216  double * newSolution = new double [numberColumns];
217  memcpy(newSolution,solution,numberColumns*sizeof(double));
218
219  // way is 1 if down possible, 2 if up possible, 3 if both possible
220  char * way = new char[numberIntegers];
221  // corrected costs
222  double * cost = new double[numberIntegers];
223  // for array to mark infeasible rows after iColumn branch
224  char * mark = new char[numberRows];
225  memset(mark,0,numberRows);
226  // space to save values so we don't introduce rounding errors
227  double * save = new double[numberRows];
228
229  // clean solution
230  for (i=0;i<numberIntegers;i++) {
231    int iColumn = integerVariable[i];
232    const OsiObject * object = model_->object(i);
233    // get original bounds
234    double originalLower;
235    double originalUpper;
236    getIntegerInformation( object,originalLower, originalUpper); 
237    double value=newSolution[iColumn];
238    if (value<originalLower) {
239      value=originalLower;
240      newSolution[iColumn]=value;
241    } else if (value>originalUpper) {
242      value=originalUpper;
243      newSolution[iColumn]=value;
244    }
245    double nearest=floor(value+0.5);
246    //assert(fabs(value-nearest)<10.0*primalTolerance);
247    value=nearest;
248    newSolution[iColumn]=nearest;
249    // if away from lower bound mark that fact
250    if (nearest>originalLower) {
251      used_[iColumn]=1;
252    }
253    cost[i] = direction*objective[iColumn];
254    int iway=0;
255   
256    if (value>originalLower+0.5) 
257      iway = 1;
258    if (value<originalUpper-0.5) 
259      iway |= 2;
260    way[i]=iway;
261  }
262  // get row activities
263  double * rowActivity = new double[numberRows];
264  memset(rowActivity,0,numberRows*sizeof(double));
265
266  for (i=0;i<numberColumns;i++) {
267    int j;
268    double value = newSolution[i];
269    if (value) {
270      for (j=columnStart[i];
271           j<columnStart[i]+columnLength[i];j++) {
272        int iRow=row[j];
273        rowActivity[iRow] += value*element[j];
274      }
275    }
276  }
277  // check was feasible - if not adjust (cleaning may move)
278  // if very infeasible then give up
279  bool tryHeuristic=true;
280  for (i=0;i<numberRows;i++) {
281    if(rowActivity[i]<rowLower[i]) {
282      if (rowActivity[i]<rowLower[i]-10.0*primalTolerance)
283        tryHeuristic=false;
284      rowActivity[i]=rowLower[i];
285    } else if(rowActivity[i]>rowUpper[i]) {
286      if (rowActivity[i]<rowUpper[i]+10.0*primalTolerance)
287        tryHeuristic=false;
288      rowActivity[i]=rowUpper[i];
289    }
290  }
291  // Switch off if may take too long
292  if (model_->getNumCols()>10000&&model_->getNumCols()>
293      10*model_->getNumRows())
294    tryHeuristic=false;
295  if (tryHeuristic) {
296   
297    // best change in objective
298    double bestChange=0.0;
299   
300    for (i=0;i<numberIntegers;i++) {
301      int iColumn = integerVariable[i];
302     
303      double objectiveCoefficient = cost[i];
304      int k;
305      int j;
306      int goodK=-1;
307      int wayK=-1,wayI=-1;
308      if ((way[i]&1)!=0) {
309        int numberInfeasible=0;
310        // save row activities and adjust
311        for (j=columnStart[iColumn];
312             j<columnStart[iColumn]+columnLength[iColumn];j++) {
313          int iRow = row[j];
314          save[iRow]=rowActivity[iRow];
315          rowActivity[iRow] -= element[j];
316          if(rowActivity[iRow]<rowLower[iRow]-primalTolerance||
317             rowActivity[iRow]>rowUpper[iRow]+primalTolerance) {
318            // mark row
319            mark[iRow]=1;
320            numberInfeasible++;
321          }
322        }
323        // try down
324        for (k=i+1;k<numberIntegers;k++) {
325          if ((way[k]&1)!=0) {
326            // try down
327            if (-objectiveCoefficient-cost[k]<bestChange) {
328              // see if feasible down
329              bool good=true;
330              int numberMarked=0;
331              int kColumn = integerVariable[k];
332              for (j=columnStart[kColumn];
333                   j<columnStart[kColumn]+columnLength[kColumn];j++) {
334                int iRow = row[j];
335                double newValue = rowActivity[iRow] - element[j];
336                if(newValue<rowLower[iRow]-primalTolerance||
337                   newValue>rowUpper[iRow]+primalTolerance) {
338                  good=false;
339                  break;
340                } else if (mark[iRow]) {
341                  // made feasible
342                  numberMarked++;
343                }
344              }
345              if (good&&numberMarked==numberInfeasible) {
346                // better solution
347                goodK=k;
348                wayK=-1;
349                wayI=-1;
350                bestChange = -objectiveCoefficient-cost[k];
351              }
352            }
353          }
354          if ((way[k]&2)!=0) {
355            // try up
356            if (-objectiveCoefficient+cost[k]<bestChange) {
357              // see if feasible up
358              bool good=true;
359              int numberMarked=0;
360              int kColumn = integerVariable[k];
361              for (j=columnStart[kColumn];
362                   j<columnStart[kColumn]+columnLength[kColumn];j++) {
363                int iRow = row[j];
364                double newValue = rowActivity[iRow] + element[j];
365                if(newValue<rowLower[iRow]-primalTolerance||
366                   newValue>rowUpper[iRow]+primalTolerance) {
367                  good=false;
368                  break;
369                } else if (mark[iRow]) {
370                  // made feasible
371                  numberMarked++;
372                }
373              }
374              if (good&&numberMarked==numberInfeasible) {
375                // better solution
376                goodK=k;
377                wayK=1;
378                wayI=-1;
379                bestChange = -objectiveCoefficient+cost[k];
380              }
381            }
382          }
383        }
384        // restore row activities
385        for (j=columnStart[iColumn];
386             j<columnStart[iColumn]+columnLength[iColumn];j++) {
387          int iRow = row[j];
388          rowActivity[iRow] = save[iRow];
389          mark[iRow]=0;
390        }
391      }
392      if ((way[i]&2)!=0) {
393        int numberInfeasible=0;
394        // save row activities and adjust
395        for (j=columnStart[iColumn];
396             j<columnStart[iColumn]+columnLength[iColumn];j++) {
397          int iRow = row[j];
398          save[iRow]=rowActivity[iRow];
399          rowActivity[iRow] += element[j];
400          if(rowActivity[iRow]<rowLower[iRow]-primalTolerance||
401             rowActivity[iRow]>rowUpper[iRow]+primalTolerance) {
402            // mark row
403            mark[iRow]=1;
404            numberInfeasible++;
405          }
406        }
407        // try up
408        for (k=i+1;k<numberIntegers;k++) {
409          if ((way[k]&1)!=0) {
410            // try down
411            if (objectiveCoefficient-cost[k]<bestChange) {
412              // see if feasible down
413              bool good=true;
414              int numberMarked=0;
415              int kColumn = integerVariable[k];
416              for (j=columnStart[kColumn];
417                   j<columnStart[kColumn]+columnLength[kColumn];j++) {
418                int iRow = row[j];
419                double newValue = rowActivity[iRow] - element[j];
420                if(newValue<rowLower[iRow]-primalTolerance||
421                   newValue>rowUpper[iRow]+primalTolerance) {
422                  good=false;
423                  break;
424                } else if (mark[iRow]) {
425                  // made feasible
426                  numberMarked++;
427                }
428              }
429              if (good&&numberMarked==numberInfeasible) {
430                // better solution
431                goodK=k;
432                wayK=-1;
433                wayI=1;
434                bestChange = objectiveCoefficient-cost[k];
435              }
436            }
437          }
438          if ((way[k]&2)!=0) {
439            // try up
440            if (objectiveCoefficient+cost[k]<bestChange) {
441              // see if feasible up
442              bool good=true;
443              int numberMarked=0;
444              int kColumn = integerVariable[k];
445              for (j=columnStart[kColumn];
446                   j<columnStart[kColumn]+columnLength[kColumn];j++) {
447                int iRow = row[j];
448                double newValue = rowActivity[iRow] + element[j];
449                if(newValue<rowLower[iRow]-primalTolerance||
450                   newValue>rowUpper[iRow]+primalTolerance) {
451                  good=false;
452                  break;
453                } else if (mark[iRow]) {
454                  // made feasible
455                  numberMarked++;
456                }
457              }
458              if (good&&numberMarked==numberInfeasible) {
459                // better solution
460                goodK=k;
461                wayK=1;
462                wayI=1;
463                bestChange = objectiveCoefficient+cost[k];
464              }
465            }
466          }
467        }
468        // restore row activities
469        for (j=columnStart[iColumn];
470             j<columnStart[iColumn]+columnLength[iColumn];j++) {
471          int iRow = row[j];
472          rowActivity[iRow] = save[iRow];
473          mark[iRow]=0;
474        }
475      }
476      if (goodK>=0) {
477        // we found something - update solution
478        for (j=columnStart[iColumn];
479             j<columnStart[iColumn]+columnLength[iColumn];j++) {
480          int iRow = row[j];
481          rowActivity[iRow]  += wayI * element[j];
482        }
483        newSolution[iColumn] += wayI;
484        int kColumn = integerVariable[goodK];
485        for (j=columnStart[kColumn];
486             j<columnStart[kColumn]+columnLength[kColumn];j++) {
487          int iRow = row[j];
488          rowActivity[iRow]  += wayK * element[j];
489        }
490        newSolution[kColumn] += wayK;
491        // See if k can go further ?
492        const OsiObject * object = model_->object(goodK);
493        // get original bounds
494        double originalLower;
495        double originalUpper;
496        getIntegerInformation( object,originalLower, originalUpper); 
497       
498        double value=newSolution[kColumn];
499        int iway=0;
500       
501        if (value>originalLower+0.5) 
502          iway = 1;
503        if (value<originalUpper-0.5) 
504          iway |= 2;
505        way[goodK]=iway;
506      }
507    }
508    if (bestChange+newSolutionValue<solutionValue) {
509      // paranoid check
510      memset(rowActivity,0,numberRows*sizeof(double));
511     
512      for (i=0;i<numberColumns;i++) {
513        int j;
514        double value = newSolution[i];
515        if (value) {
516          for (j=columnStart[i];
517               j<columnStart[i]+columnLength[i];j++) {
518            int iRow=row[j];
519            rowActivity[iRow] += value*element[j];
520          }
521        }
522      }
523      int numberBad=0;
524      double sumBad=0.0;
525      // check was approximately feasible
526      for (i=0;i<numberRows;i++) {
527        if(rowActivity[i]<rowLower[i]) {
528          sumBad += rowLower[i]-rowActivity[i];
529          if (rowActivity[i]<rowLower[i]-10.0*primalTolerance)
530            numberBad++;
531        } else if(rowActivity[i]>rowUpper[i]) {
532          sumBad += rowUpper[i]-rowActivity[i];
533          if (rowActivity[i]>rowUpper[i]+10.0*primalTolerance)
534            numberBad++;
535        }
536      }
537      if (!numberBad) {
538        for (i=0;i<numberIntegers;i++) {
539          int iColumn = integerVariable[i];
540          const OsiObject * object = model_->object(i);
541          // get original bounds
542          double originalLower;
543          double originalUpper;
544          getIntegerInformation( object,originalLower, originalUpper); 
545         
546          double value=newSolution[iColumn];
547          // if away from lower bound mark that fact
548          if (value>originalLower) {
549            used_[iColumn]=1;
550          }
551        }
552        // new solution
553        memcpy(betterSolution,newSolution,numberColumns*sizeof(double));
554        CoinWarmStartBasis * basis =
555          dynamic_cast<CoinWarmStartBasis *>(solver->getWarmStart()) ;
556        if (basis) {
557          model_->setBestSolutionBasis(* basis);
558          delete basis;
559        }
560        returnCode=1;
561        solutionValue = newSolutionValue + bestChange;
562      } else {
563        // bad solution - should not happen so debug if see message
564        printf("Local search got bad solution with %d infeasibilities summing to %g\n",
565               numberBad,sumBad);
566      }
567    }
568  }
569  delete [] newSolution;
570  delete [] rowActivity;
571  delete [] way;
572  delete [] cost;
573  delete [] save;
574  delete [] mark;
575  if (numberSolutions_>1&&swap_==1) {
576    // try merge
577    int returnCode2=solutionFix( solutionValue, betterSolution,NULL);
578    if (returnCode2)
579      returnCode=1;
580  }
581  return returnCode;
582}
583// update model
584void CbcHeuristicLocal::setModel(CbcModel * model)
585{
586  model_ = model;
587  // Get a copy of original matrix
588  assert(model_->solver());
589  matrix_ = *model_->solver()->getMatrixByCol();
590  delete [] used_;
591  int numberColumns = model->solver()->getNumCols();
592  used_ = new char[numberColumns];
593  memset(used_,0,numberColumns);
594}
595// Default Constructor
596CbcHeuristicNaive::CbcHeuristicNaive() 
597  :CbcHeuristic()
598{
599  large_=1.0e6;
600}
601
602// Constructor with model - assumed before cuts
603
604CbcHeuristicNaive::CbcHeuristicNaive(CbcModel & model)
605  :CbcHeuristic(model)
606{
607  large_=1.0e6;
608}
609
610// Destructor
611CbcHeuristicNaive::~CbcHeuristicNaive ()
612{
613}
614
615// Clone
616CbcHeuristic *
617CbcHeuristicNaive::clone() const
618{
619  return new CbcHeuristicNaive(*this);
620}
621// Create C++ lines to get to current state
622void 
623CbcHeuristicNaive::generateCpp( FILE * fp) 
624{
625  CbcHeuristicNaive other;
626  fprintf(fp,"0#include \"CbcHeuristicNaive.hpp\"\n");
627  fprintf(fp,"3  CbcHeuristicNaive naive(*cbcModel);\n");
628  CbcHeuristic::generateCpp(fp,"naive");
629  if (large_!=other.large_)
630    fprintf(fp,"3  naive.setLarge(%g);\n",large_);
631  else
632    fprintf(fp,"4  naive.setLarge(%g);\n",large_);
633  fprintf(fp,"3  cbcModel->addHeuristic(&naive);\n");
634}
635
636// Copy constructor
637CbcHeuristicNaive::CbcHeuristicNaive(const CbcHeuristicNaive & rhs)
638:
639  CbcHeuristic(rhs),
640  large_(rhs.large_)
641{
642}
643
644// Assignment operator
645CbcHeuristicNaive & 
646CbcHeuristicNaive::operator=( const CbcHeuristicNaive& rhs)
647{
648  if (this!=&rhs) {
649    CbcHeuristic::operator=(rhs);
650    large_ = rhs.large_;
651  }
652  return *this;
653}
654
655// Resets stuff if model changes
656void 
657CbcHeuristicNaive::resetModel(CbcModel * model)
658{
659  CbcHeuristic::resetModel(model);
660}
661int
662CbcHeuristicNaive::solution(double & solutionValue,
663                         double * betterSolution)
664{
665  numCouldRun_++;
666  // See if to do
667  bool atRoot = model_->getNodeCount()==0;
668  int passNumber = model_->getCurrentPassNumber();
669  if (!when()||(when()==1&&model_->phase()!=1)||!atRoot||passNumber!=1)
670    return 0; // switched off
671  // Don't do if it was this heuristic which found solution!
672  if (this==model_->lastHeuristic())
673    return 0;
674  numRuns_++;
675  double cutoff;
676  model_->solver()->getDblParam(OsiDualObjectiveLimit,cutoff);
677  double direction = model_->solver()->getObjSense();
678  cutoff *= direction;
679  cutoff = CoinMin(cutoff,solutionValue);
680  OsiSolverInterface * solver = model_->continuousSolver();
681  if (!solver)
682    solver = model_->solver();
683  const double * colLower = solver->getColLower();
684  const double * colUpper = solver->getColUpper();
685  const double * objective = solver->getObjCoefficients();
686
687  int numberColumns = model_->getNumCols();
688  int numberIntegers = model_->numberIntegers();
689  const int * integerVariable = model_->integerVariable();
690 
691  int i;
692  bool solutionFound=false;
693  CoinWarmStartBasis saveBasis;
694  CoinWarmStartBasis * basis =
695    dynamic_cast<CoinWarmStartBasis *>(solver->getWarmStart()) ;
696  if (basis) {
697    saveBasis = * basis;
698    delete basis;
699  }
700  // First just fix all integers as close to zero as possible
701  OsiSolverInterface * newSolver = solver->clone();
702  for (i=0;i<numberIntegers;i++) {
703    int iColumn=integerVariable[i];
704    double lower = colLower[iColumn];
705    double upper = colUpper[iColumn];
706    double value;
707    if (lower>0.0)
708      value=lower;
709    else if (upper<0.0)
710      value=upper;
711    else
712      value=0.0;
713    newSolver->setColLower(iColumn,value);
714    newSolver->setColUpper(iColumn,value);
715  }
716  newSolver->initialSolve();
717  if (newSolver->isProvenOptimal()) {
718    double solValue = newSolver->getObjValue()*direction ;
719    if (solValue<cutoff) {
720      // we have a solution
721      solutionFound=true;
722      solutionValue=solValue;
723      memcpy(betterSolution,newSolver->getColSolution(),
724             numberColumns*sizeof(double));
725      printf("Naive fixing close to zero gave solution of %g\n",solutionValue);
726      cutoff = solValue - model_->getCutoffIncrement();
727    }
728  }
729  // Now fix all integers as close to zero if zero or large cost
730  int nFix=0;
731  for (i=0;i<numberIntegers;i++) {
732    int iColumn=integerVariable[i];
733    double lower = colLower[iColumn];
734    double upper = colUpper[iColumn];
735    double value;
736    if (fabs(objective[i])>0.0&&fabs(objective[i])<large_) {
737      nFix++;
738      if (lower>0.0)
739        value=lower;
740      else if (upper<0.0)
741        value=upper;
742      else
743        value=0.0;
744      newSolver->setColLower(iColumn,value);
745      newSolver->setColUpper(iColumn,value);
746    }
747  }
748  const double * solution = solver->getColSolution();
749  if (nFix) {
750    newSolver->setWarmStart(&saveBasis);
751    newSolver->setColSolution(solution);
752    newSolver->initialSolve();
753    if (newSolver->isProvenOptimal()) {
754      double solValue = newSolver->getObjValue()*direction ;
755      if (solValue<cutoff) {
756        // try branch and bound
757        double * newSolution = new double [numberColumns];
758        printf("%d fixed after fixing costs\n",nFix);
759        int returnCode = smallBranchAndBound(newSolver,
760                                             numberNodes_,newSolution,
761                                             solutionValue,
762                                             solutionValue,"CbcHeuristicNaive1");
763        if (returnCode<0)
764          returnCode=0; // returned on size
765        if ((returnCode&2)!=0) {
766          // could add cut
767          returnCode &= ~2;
768        }
769        if (returnCode==1) {
770          // solution
771          solutionFound=true;
772          memcpy(betterSolution,newSolution,
773                 numberColumns*sizeof(double));
774          printf("Naive fixing zeros gave solution of %g\n",solutionValue);
775          cutoff = solutionValue - model_->getCutoffIncrement();
776        }
777        delete [] newSolution;
778      }
779    }
780  }
781#if 1
782  newSolver->setObjSense(-direction); // maximize
783  newSolver->setWarmStart(&saveBasis);
784  newSolver->setColSolution(solution);
785  for (int iColumn=0;iColumn<numberColumns;iColumn++) {
786    double value = solution[iColumn];
787    double lower = colLower[iColumn];
788    double upper = colUpper[iColumn];
789    double newLower;
790    double newUpper;
791    if (newSolver->isInteger(iColumn)) {
792      newLower = CoinMax(lower,floor(value)-2.0);
793      newUpper = CoinMin(upper,ceil(value)+2.0);
794    } else {
795      newLower = CoinMax(lower,value-1.0e5);
796      newUpper = CoinMin(upper,value+1.0e-5);
797    }
798    newSolver->setColLower(iColumn,newLower);
799    newSolver->setColUpper(iColumn,newUpper);
800  }
801  newSolver->initialSolve();
802  if (newSolver->isProvenOptimal()) {
803    double solValue = newSolver->getObjValue()*direction ;
804    if (solValue<cutoff) {
805      nFix=0;
806      newSolver->setObjSense(direction); // correct direction
807      //const double * thisSolution = newSolver->getColSolution();
808      for (int iColumn=0;iColumn<numberColumns;iColumn++) {
809        double value = solution[iColumn];
810        double lower = colLower[iColumn];
811        double upper = colUpper[iColumn];
812        double newLower=lower;
813        double newUpper=upper;
814        if (newSolver->isInteger(iColumn)) {
815          if (value<lower+1.0e-6) {
816            nFix++;
817            newUpper=lower;
818          } else if (value>upper-1.0e-6) {
819            nFix++;
820            newLower=upper;
821          } else {
822            newLower = CoinMax(lower,floor(value)-2.0);
823            newUpper = CoinMin(upper,ceil(value)+2.0);
824          }
825        }
826        newSolver->setColLower(iColumn,newLower);
827        newSolver->setColUpper(iColumn,newUpper);
828      }
829      // try branch and bound
830      double * newSolution = new double [numberColumns];
831      printf("%d fixed after maximizing\n",nFix);
832      int returnCode = smallBranchAndBound(newSolver,
833                                           numberNodes_,newSolution,
834                                           solutionValue,
835                                           solutionValue,"CbcHeuristicNaive1");
836      if (returnCode<0)
837        returnCode=0; // returned on size
838      if ((returnCode&2)!=0) {
839        // could add cut
840        returnCode &= ~2;
841      }
842      if (returnCode==1) {
843        // solution
844        solutionFound=true;
845        memcpy(betterSolution,newSolution,
846               numberColumns*sizeof(double));
847        printf("Naive maximizing gave solution of %g\n",solutionValue);
848        cutoff = solutionValue - model_->getCutoffIncrement();
849      }
850      delete [] newSolution;
851    }
852  }
853#endif
854  delete newSolver;
855  return solutionFound ? 1 : 0;
856}
857// update model
858void CbcHeuristicNaive::setModel(CbcModel * model)
859{
860  model_ = model;
861}
862
863 
Note: See TracBrowser for help on using the repository browser.