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

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

changes for heuristic clone

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