source: branches/devel/Cbc/src/CbcHeuristic.cpp @ 648

Last change on this file since 648 was 642, checked in by forrest, 12 years ago

update branches/devel for threads

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 31.2 KB
Line 
1// Copyright (C) 2002, International Business Machines
2// Corporation and others.  All Rights Reserved.
3#if defined(_MSC_VER)
4// Turn off compiler warning about long names
5#  pragma warning(disable:4786)
6#endif
7
8#include "CbcConfig.h"
9
10#include <cassert>
11#include <cmath>
12#include <cfloat>
13
14#ifdef COIN_HAS_CLP
15#include "OsiClpSolverInterface.hpp"
16#endif
17#include "CbcModel.hpp"
18#include "CbcMessage.hpp"
19#include "CbcHeuristic.hpp"
20#include "CbcStrategy.hpp"
21#include "CglPreProcess.hpp"
22#include "OsiAuxInfo.hpp"
23
24// Default Constructor
25CbcHeuristic::CbcHeuristic() 
26  :model_(NULL),
27   when_(2),
28   numberNodes_(200),
29   fractionSmall_(1.0),
30   heuristicName_("Unknown")
31{
32  // As CbcHeuristic virtual need to modify .cpp if above change
33}
34
35// Constructor from model
36CbcHeuristic::CbcHeuristic(CbcModel & model)
37:
38  model_(&model),
39  when_(2),
40  numberNodes_(200),
41  fractionSmall_(1.0),
42  heuristicName_("Unknown")
43{
44  // As CbcHeuristic virtual need to modify .cpp if above change
45}
46// Copy constructor
47CbcHeuristic::CbcHeuristic(const CbcHeuristic & rhs)
48:
49  model_(rhs.model_),
50  when_(rhs.when_),
51  numberNodes_(rhs.numberNodes_),
52  fractionSmall_(rhs.fractionSmall_),
53  heuristicName_(rhs.heuristicName_)
54{
55}
56// Assignment operator
57CbcHeuristic & 
58CbcHeuristic::operator=( const CbcHeuristic& rhs)
59{
60  if (this!=&rhs) {
61    model_ = rhs.model_;
62    when_ = rhs.when_;
63    numberNodes_ = rhs.numberNodes_;
64    fractionSmall_ = rhs.fractionSmall_;
65    heuristicName_ = rhs.heuristicName_ ;
66  }
67  return *this;
68}
69
70// Resets stuff if model changes
71void 
72CbcHeuristic::resetModel(CbcModel * model)
73{
74  model_=model;
75}
76
77// Create C++ lines to get to current state
78void 
79CbcHeuristic::generateCpp( FILE * fp, const char * heuristic) 
80{
81  // hard coded as CbcHeuristic virtual
82  if (when_!=2)
83    fprintf(fp,"3  %s.setWhen(%d);\n",heuristic,when_);
84  else
85    fprintf(fp,"4  %s.setWhen(%d);\n",heuristic,when_);
86  if (numberNodes_!=200)
87    fprintf(fp,"3  %s.setNumberNodes(%d);\n",heuristic,numberNodes_);
88  else
89    fprintf(fp,"4  %s.setNumberNodes(%d);\n",heuristic,numberNodes_);
90  if (fractionSmall_!=1.0)
91    fprintf(fp,"3  %s.setFractionSmall(%g);\n",heuristic,fractionSmall_);
92  else
93    fprintf(fp,"4  %s.setFractionSmall(%g);\n",heuristic,fractionSmall_);
94  if (heuristicName_ != "Unknown")
95    fprintf(fp,"3  %s.setHeuristicName(\"%s\");\n",
96            heuristic,heuristicName_.c_str()) ;
97  else
98    fprintf(fp,"4  %s.setHeuristicName(\"%s\");\n",
99            heuristic,heuristicName_.c_str()) ;
100}
101// Destructor
102CbcHeuristic::~CbcHeuristic ()
103{
104}
105
106// update model
107void CbcHeuristic::setModel(CbcModel * model)
108{
109  model_ = model;
110}
111// Do mini branch and bound (return 1 if solution)
112int 
113CbcHeuristic::smallBranchAndBound(OsiSolverInterface * solver,int numberNodes,
114                                  double * newSolution, double & newSolutionValue,
115                                  double cutoff, std::string name) const
116{
117#ifdef COIN_HAS_CLP
118  OsiClpSolverInterface * osiclp = dynamic_cast< OsiClpSolverInterface*> (solver);
119  if (osiclp&&(osiclp->specialOptions()&65536)==0) {
120    // go faster stripes
121    if (osiclp->getNumRows()<300&&osiclp->getNumCols()<500) {
122      osiclp->setupForRepeatedUse(2,0);
123    } else {
124      osiclp->setupForRepeatedUse(0,0);
125    }
126    // Turn this off if you get problems
127    // Used to be automatically set
128    osiclp->setSpecialOptions(osiclp->specialOptions()|(128+64));
129  }
130#endif
131  // Reduce printout
132  solver->setHintParam(OsiDoReducePrint,true,OsiHintTry);
133  solver->setHintParam(OsiDoPresolveInInitial,false,OsiHintTry);
134  solver->setDblParam(OsiDualObjectiveLimit,cutoff*solver->getObjSense());
135  solver->initialSolve();
136  int returnCode=1;
137  int logLevel = model_->logLevel();
138  if (solver->isProvenOptimal()) {
139    CglPreProcess process;
140    /* Do not try and produce equality cliques and
141       do up to 5 passes */
142    if (logLevel<=1)
143      process.messageHandler()->setLogLevel(0);
144    OsiSolverInterface * solver2= process.preProcess(*solver);
145    if (!solver2) {
146      if (logLevel>1)
147        printf("Pre-processing says infeasible\n");
148      returnCode=2; // so will be infeasible
149    } else {
150      // see if too big
151      double before = solver->getNumRows()+solver->getNumCols();
152      double after = solver2->getNumRows()+solver2->getNumCols();
153      if (logLevel>1)
154        printf("before %d rows %d columns, after %d rows %d columns\n",
155               solver->getNumRows(),solver->getNumCols(),
156               solver2->getNumRows(),solver2->getNumCols());
157      if (after>fractionSmall_*before)
158        return 0;
159      solver2->resolve();
160      CbcModel model(*solver2);
161      if (logLevel<=1)
162        model.setLogLevel(0);
163      else
164        model.setLogLevel(logLevel);
165      model.setCutoff(cutoff);
166      model.setMaximumNodes(numberNodes);
167      model.solver()->setHintParam(OsiDoReducePrint,true,OsiHintTry);
168      // Lightweight
169      CbcStrategyDefaultSubTree strategy(model_,true,5,1,0);
170      model.setStrategy(strategy);
171      model.solver()->setIntParam(OsiMaxNumIterationHotStart,10);
172      // Do search
173      if (logLevel>1)
174        model_->messageHandler()->message(CBC_START_SUB,model_->messages())
175          << name
176          << model.getMaximumNodes()
177          <<CoinMessageEol;
178      // probably faster to use a basis to get integer solutions
179      model.setSpecialOptions(2);
180#ifdef CBC_THREAD
181      if (model_->getNumberThreads()>0&&(model_->getThreadMode()&1)!=0) {
182        // See if at root node
183        bool atRoot = model_->getNodeCount()==0;
184        int passNumber = model_->getCurrentPassNumber();
185        if (atRoot&&passNumber==1)
186          model.setNumberThreads(model_->getNumberThreads());
187      }
188#endif
189      model.branchAndBound();
190      if (logLevel>1)
191        model_->messageHandler()->message(CBC_END_SUB,model_->messages())
192          << name
193          <<CoinMessageEol;
194      if (model.getMinimizationObjValue()<CoinMin(cutoff,1.0e30)) {
195        // solution
196        returnCode=model.isProvenOptimal() ? 3 : 1;
197        // post process
198        process.postProcess(*model.solver());
199        if (solver->isProvenOptimal()) {
200          // Solution now back in solver
201          int numberColumns = solver->getNumCols();
202          memcpy(newSolution,solver->getColSolution(),
203                 numberColumns*sizeof(double));
204          newSolutionValue = model.getMinimizationObjValue();
205        } else {
206          // odd - but no good
207          returnCode=0; // so will be infeasible
208        }
209      } else {
210        // no good
211        returnCode=model.isProvenInfeasible() ? 2 : 0; // so will be infeasible
212      }
213    }
214  } else {
215    returnCode=2; // infeasible finished
216  }
217  return returnCode;
218}
219
220// Default Constructor
221CbcRounding::CbcRounding() 
222  :CbcHeuristic()
223{
224  // matrix and row copy will automatically be empty
225  seed_=1;
226}
227
228// Constructor from model
229CbcRounding::CbcRounding(CbcModel & model)
230  :CbcHeuristic(model)
231{
232  // Get a copy of original matrix (and by row for rounding);
233  assert(model.solver());
234  matrix_ = *model.solver()->getMatrixByCol();
235  matrixByRow_ = *model.solver()->getMatrixByRow();
236  seed_=1;
237}
238
239// Destructor
240CbcRounding::~CbcRounding ()
241{
242}
243
244// Clone
245CbcHeuristic *
246CbcRounding::clone() const
247{
248  return new CbcRounding(*this);
249}
250// Create C++ lines to get to current state
251void 
252CbcRounding::generateCpp( FILE * fp) 
253{
254  CbcRounding other;
255  fprintf(fp,"0#include \"CbcHeuristic.hpp\"\n");
256  fprintf(fp,"3  CbcRounding rounding(*cbcModel);\n");
257  CbcHeuristic::generateCpp(fp,"rounding");
258  if (seed_!=other.seed_)
259    fprintf(fp,"3  rounding.setSeed(%d);\n",seed_);
260  else
261    fprintf(fp,"4  rounding.setSeed(%d);\n",seed_);
262  fprintf(fp,"3  cbcModel->addHeuristic(&rounding);\n");
263}
264
265// Copy constructor
266CbcRounding::CbcRounding(const CbcRounding & rhs)
267:
268  CbcHeuristic(rhs),
269  matrix_(rhs.matrix_),
270  matrixByRow_(rhs.matrixByRow_),
271  seed_(rhs.seed_)
272{
273}
274
275// Assignment operator
276CbcRounding & 
277CbcRounding::operator=( const CbcRounding& rhs)
278{
279  if (this!=&rhs) {
280    CbcHeuristic::operator=(rhs);
281    matrix_ = rhs.matrix_;
282    matrixByRow_ = rhs.matrixByRow_;
283    seed_ = rhs.seed_;
284  }
285  return *this;
286}
287
288// Resets stuff if model changes
289void 
290CbcRounding::resetModel(CbcModel * model)
291{
292  model_=model;
293  // Get a copy of original matrix (and by row for rounding);
294  assert(model_->solver());
295  matrix_ = *model_->solver()->getMatrixByCol();
296  matrixByRow_ = *model_->solver()->getMatrixByRow();
297}
298// See if rounding will give solution
299// Sets value of solution
300// Assumes rhs for original matrix still okay
301// At present only works with integers
302// Fix values if asked for
303// Returns 1 if solution, 0 if not
304int
305CbcRounding::solution(double & solutionValue,
306                         double * betterSolution)
307{
308
309  // See if to do
310  if (!when()||(when()%10==1&&model_->phase()!=1)||
311      (when()%10==2&&(model_->phase()!=2&&model_->phase()!=3)))
312    return 0; // switched off
313  OsiSolverInterface * solver = model_->solver();
314  const double * lower = solver->getColLower();
315  const double * upper = solver->getColUpper();
316  const double * rowLower = solver->getRowLower();
317  const double * rowUpper = solver->getRowUpper();
318  const double * solution = solver->getColSolution();
319  const double * objective = solver->getObjCoefficients();
320  double integerTolerance = model_->getDblParam(CbcModel::CbcIntegerTolerance);
321  double primalTolerance;
322  solver->getDblParam(OsiPrimalTolerance,primalTolerance);
323
324  int numberRows = matrix_.getNumRows();
325  assert (numberRows<=solver->getNumRows());
326  int numberIntegers = model_->numberIntegers();
327  const int * integerVariable = model_->integerVariable();
328  int i;
329  double direction = solver->getObjSense();
330  double newSolutionValue = direction*solver->getObjValue();
331  int returnCode = 0;
332  // Column copy
333  const double * element = matrix_.getElements();
334  const int * row = matrix_.getIndices();
335  const CoinBigIndex * columnStart = matrix_.getVectorStarts();
336  const int * columnLength = matrix_.getVectorLengths();
337  // Row copy
338  const double * elementByRow = matrixByRow_.getElements();
339  const int * column = matrixByRow_.getIndices();
340  const CoinBigIndex * rowStart = matrixByRow_.getVectorStarts();
341  const int * rowLength = matrixByRow_.getVectorLengths();
342
343  // Get solution array for heuristic solution
344  int numberColumns = solver->getNumCols();
345  double * newSolution = new double [numberColumns];
346  memcpy(newSolution,solution,numberColumns*sizeof(double));
347
348  double * rowActivity = new double[numberRows];
349  memset(rowActivity,0,numberRows*sizeof(double));
350  for (i=0;i<numberColumns;i++) {
351    int j;
352    double value = newSolution[i];
353    if (value<lower[i]) {
354      value=lower[i];
355      newSolution[i]=value;
356    } else if (value>upper[i]) {
357      value=upper[i];
358      newSolution[i]=value;
359    }
360    if (value) {
361      for (j=columnStart[i];
362           j<columnStart[i]+columnLength[i];j++) {
363        int iRow=row[j];
364        rowActivity[iRow] += value*element[j];
365      }
366    }
367  }
368  // check was feasible - if not adjust (cleaning may move)
369  for (i=0;i<numberRows;i++) {
370    if(rowActivity[i]<rowLower[i]) {
371      //assert (rowActivity[i]>rowLower[i]-1000.0*primalTolerance);
372      rowActivity[i]=rowLower[i];
373    } else if(rowActivity[i]>rowUpper[i]) {
374      //assert (rowActivity[i]<rowUpper[i]+1000.0*primalTolerance);
375      rowActivity[i]=rowUpper[i];
376    }
377  }
378  for (i=0;i<numberIntegers;i++) {
379    int iColumn = integerVariable[i];
380    double value=newSolution[iColumn];
381    if (fabs(floor(value+0.5)-value)>integerTolerance) {
382      double below = floor(value);
383      double newValue=newSolution[iColumn];
384      double cost = direction * objective[iColumn];
385      double move;
386      if (cost>0.0) {
387        // try up
388        move = 1.0 -(value-below);
389      } else if (cost<0.0) {
390        // try down
391        move = below-value;
392      } else {
393        // won't be able to move unless we can grab another variable
394        double randomNumber = CoinDrand48();
395        // which way?
396        if (randomNumber<0.5) 
397          move = below-value;
398        else
399          move = 1.0 -(value-below);
400      }
401      newValue += move;
402      newSolution[iColumn] = newValue;
403      newSolutionValue += move*cost;
404      int j;
405      for (j=columnStart[iColumn];
406           j<columnStart[iColumn]+columnLength[iColumn];j++) {
407        int iRow = row[j];
408        rowActivity[iRow] += move*element[j];
409      }
410    }
411  }
412
413  double penalty=0.0;
414  const char * integerType = model_->integerType();
415  // see if feasible - just using singletons
416  for (i=0;i<numberRows;i++) {
417    double value = rowActivity[i];
418    double thisInfeasibility=0.0;
419    if (value<rowLower[i]-primalTolerance)
420      thisInfeasibility = value-rowLower[i];
421    else if (value>rowUpper[i]+primalTolerance)
422      thisInfeasibility = value-rowUpper[i];
423    if (thisInfeasibility) {
424      // See if there are any slacks I can use to fix up
425      // maybe put in coding for multiple slacks?
426      double bestCost = 1.0e50;
427      int k;
428      int iBest=-1;
429      double addCost=0.0;
430      double newValue=0.0;
431      double changeRowActivity=0.0;
432      double absInfeasibility = fabs(thisInfeasibility);
433      for (k=rowStart[i];k<rowStart[i]+rowLength[i];k++) {
434        int iColumn = column[k];
435        // See if all elements help
436        if (columnLength[iColumn]==1) {
437          double currentValue = newSolution[iColumn];
438          double elementValue = elementByRow[k];
439          double lowerValue = lower[iColumn];
440          double upperValue = upper[iColumn];
441          double gap = rowUpper[i]-rowLower[i];
442          double absElement=fabs(elementValue);
443          if (thisInfeasibility*elementValue>0.0) {
444            // we want to reduce
445            if ((currentValue-lowerValue)*absElement>=absInfeasibility) {
446              // possible - check if integer
447              double distance = absInfeasibility/absElement;
448              double thisCost = -direction*objective[iColumn]*distance;
449              if (integerType[iColumn]) {
450                distance = ceil(distance-primalTolerance);
451                if (currentValue-distance>=lowerValue-primalTolerance) {
452                  if (absInfeasibility-distance*absElement< -gap-primalTolerance)
453                    thisCost=1.0e100; // no good
454                  else
455                    thisCost = -direction*objective[iColumn]*distance;
456                } else {
457                  thisCost=1.0e100; // no good
458                }
459              }
460              if (thisCost<bestCost) {
461                bestCost=thisCost;
462                iBest=iColumn;
463                addCost = thisCost;
464                newValue = currentValue-distance;
465                changeRowActivity = -distance*elementValue;
466              }
467            }
468          } else {
469            // we want to increase
470            if ((upperValue-currentValue)*absElement>=absInfeasibility) {
471              // possible - check if integer
472              double distance = absInfeasibility/absElement;
473              double thisCost = direction*objective[iColumn]*distance;
474              if (integerType[iColumn]) {
475                distance = ceil(distance-1.0e-7);
476                assert (currentValue-distance<=upperValue+primalTolerance);
477                if (absInfeasibility-distance*absElement< -gap-primalTolerance)
478                  thisCost=1.0e100; // no good
479                else
480                  thisCost = direction*objective[iColumn]*distance;
481              }
482              if (thisCost<bestCost) {
483                bestCost=thisCost;
484                iBest=iColumn;
485                addCost = thisCost;
486                newValue = currentValue+distance;
487                changeRowActivity = distance*elementValue;
488              }
489            }
490          }
491        }
492      }
493      if (iBest>=0) {
494        /*printf("Infeasibility of %g on row %d cost %g\n",
495          thisInfeasibility,i,addCost);*/
496        newSolution[iBest]=newValue;
497        thisInfeasibility=0.0;
498        newSolutionValue += addCost;
499        rowActivity[i] += changeRowActivity;
500      }
501      penalty += fabs(thisInfeasibility);
502    }
503  }
504  if (penalty) {
505    // see if feasible using any
506    // first continuous
507    double penaltyChange=0.0;
508    int iColumn;
509    for (iColumn=0;iColumn<numberColumns;iColumn++) {
510      if (integerType[iColumn])
511        continue;
512      double currentValue = newSolution[iColumn];
513      double lowerValue = lower[iColumn];
514      double upperValue = upper[iColumn];
515      int j;
516      int anyBadDown=0;
517      int anyBadUp=0;
518      double upImprovement=0.0;
519      double downImprovement=0.0;
520      for (j=columnStart[iColumn];
521           j<columnStart[iColumn]+columnLength[iColumn];j++) {
522        int iRow = row[j];
523        if (rowUpper[iRow]>rowLower[iRow]) {
524          double value = element[j];
525          if (rowActivity[iRow]>rowUpper[iRow]+primalTolerance) {
526            // infeasible above
527            downImprovement += value;
528            upImprovement -= value;
529            if (value>0.0) 
530              anyBadUp++;
531            else 
532              anyBadDown++;
533          } else if (rowActivity[iRow]>rowUpper[iRow]-primalTolerance) {
534            // feasible at ub
535            if (value>0.0) {
536              upImprovement -= value;
537              anyBadUp++;
538            } else {
539              downImprovement += value;
540              anyBadDown++;
541            }
542          } else if (rowActivity[iRow]>rowLower[iRow]+primalTolerance) {
543            // feasible in interior
544          } else if (rowActivity[iRow]>rowLower[iRow]-primalTolerance) {
545            // feasible at lb
546            if (value<0.0) {
547              upImprovement += value;
548              anyBadUp++;
549            } else {
550              downImprovement -= value;
551              anyBadDown++;
552            }
553          } else {
554            // infeasible below
555            downImprovement -= value;
556            upImprovement += value;
557            if (value<0.0) 
558              anyBadUp++;
559            else 
560              anyBadDown++;
561          }
562        } else {
563          // equality row
564          double value = element[j];
565          if (rowActivity[iRow]>rowUpper[iRow]+primalTolerance) {
566            // infeasible above
567            downImprovement += value;
568            upImprovement -= value;
569            if (value>0.0) 
570              anyBadUp++;
571            else 
572              anyBadDown++;
573          } else if (rowActivity[iRow]<rowLower[iRow]-primalTolerance) {
574            // infeasible below
575            downImprovement -= value;
576            upImprovement += value;
577            if (value<0.0) 
578              anyBadUp++;
579            else 
580              anyBadDown++;
581          } else {
582            // feasible - no good
583            anyBadUp=-1;
584            anyBadDown=-1;
585            break;
586          }
587        }
588      }
589      // could change tests for anyBad
590      if (anyBadUp)
591        upImprovement=0.0;
592      if (anyBadDown)
593        downImprovement=0.0;
594      double way=0.0;
595      double improvement=0.0;
596      if (downImprovement>0.0&&currentValue>lowerValue) {
597        way=-1.0;
598        improvement = downImprovement;
599      } else if (upImprovement>0.0&&currentValue<upperValue) {
600        way=1.0;
601        improvement = upImprovement;
602      }
603      if (way) {
604        // can improve
605        double distance;
606        if (way>0.0)
607          distance = upperValue-currentValue;
608        else
609          distance = currentValue-lowerValue;
610        for (j=columnStart[iColumn];
611             j<columnStart[iColumn]+columnLength[iColumn];j++) {
612          int iRow = row[j];
613          double value = element[j]*way;
614          if (rowActivity[iRow]>rowUpper[iRow]+primalTolerance) {
615            // infeasible above
616            assert (value<0.0);
617            double gap = rowActivity[iRow]-rowUpper[iRow];
618            if (gap+value*distance<0.0) 
619              distance = -gap/value;
620          } else if (rowActivity[iRow]<rowLower[iRow]-primalTolerance) {
621            // infeasible below
622            assert (value>0.0);
623            double gap = rowActivity[iRow]-rowLower[iRow];
624            if (gap+value*distance>0.0) 
625              distance = -gap/value;
626          } else {
627            // feasible
628            if (value>0) {
629              double gap = rowActivity[iRow]-rowUpper[iRow];
630              if (gap+value*distance>0.0) 
631              distance = -gap/value;
632            } else {
633              double gap = rowActivity[iRow]-rowLower[iRow];
634              if (gap+value*distance<0.0) 
635                distance = -gap/value;
636            }
637          }
638        }
639        //move
640        penaltyChange += improvement*distance;
641        distance *= way;
642        newSolution[iColumn] += distance;
643        newSolutionValue += direction*objective[iColumn]*distance;
644        for (j=columnStart[iColumn];
645             j<columnStart[iColumn]+columnLength[iColumn];j++) {
646          int iRow = row[j];
647          double value = element[j];
648          rowActivity[iRow] += distance*value;
649        }
650      }
651    }
652    // and now all if improving
653    double lastChange= penaltyChange ? 1.0 : 0.0;
654    while (lastChange>1.0e-2) {
655      lastChange=0;
656      for (iColumn=0;iColumn<numberColumns;iColumn++) {
657        bool isInteger = integerType[iColumn];
658        double currentValue = newSolution[iColumn];
659        double lowerValue = lower[iColumn];
660        double upperValue = upper[iColumn];
661        int j;
662        int anyBadDown=0;
663        int anyBadUp=0;
664        double upImprovement=0.0;
665        double downImprovement=0.0;
666        for (j=columnStart[iColumn];
667             j<columnStart[iColumn]+columnLength[iColumn];j++) {
668          int iRow = row[j];
669          double value = element[j];
670          if (isInteger) {
671            if (value>0.0) {
672              if (rowActivity[iRow]+value>rowUpper[iRow]+primalTolerance)
673                anyBadUp++;
674              if (rowActivity[iRow]-value<rowLower[iRow]-primalTolerance)
675                anyBadDown++;
676            } else {
677              if (rowActivity[iRow]-value>rowUpper[iRow]+primalTolerance)
678                anyBadDown++;
679              if (rowActivity[iRow]+value<rowLower[iRow]-primalTolerance)
680                anyBadUp++;
681            }
682          }
683          if (rowUpper[iRow]>rowLower[iRow]) {
684            if (rowActivity[iRow]>rowUpper[iRow]+primalTolerance) {
685              // infeasible above
686              downImprovement += value;
687              upImprovement -= value;
688              if (value>0.0) 
689                anyBadUp++;
690              else 
691                anyBadDown++;
692            } else if (rowActivity[iRow]>rowUpper[iRow]-primalTolerance) {
693              // feasible at ub
694              if (value>0.0) {
695                upImprovement -= value;
696                anyBadUp++;
697              } else {
698                downImprovement += value;
699                anyBadDown++;
700              }
701            } else if (rowActivity[iRow]>rowLower[iRow]+primalTolerance) {
702              // feasible in interior
703            } else if (rowActivity[iRow]>rowLower[iRow]-primalTolerance) {
704              // feasible at lb
705              if (value<0.0) {
706                upImprovement += value;
707                anyBadUp++;
708              } else {
709                downImprovement -= value;
710                anyBadDown++;
711              }
712            } else {
713              // infeasible below
714              downImprovement -= value;
715              upImprovement += value;
716              if (value<0.0) 
717                anyBadUp++;
718              else 
719                anyBadDown++;
720            }
721          } else {
722            // equality row
723            if (rowActivity[iRow]>rowUpper[iRow]+primalTolerance) {
724              // infeasible above
725              downImprovement += value;
726              upImprovement -= value;
727              if (value>0.0) 
728                anyBadUp++;
729              else 
730                anyBadDown++;
731            } else if (rowActivity[iRow]<rowLower[iRow]-primalTolerance) {
732              // infeasible below
733              downImprovement -= value;
734              upImprovement += value;
735              if (value<0.0) 
736                anyBadUp++;
737              else 
738                anyBadDown++;
739            } else {
740              // feasible - no good
741              anyBadUp=-1;
742              anyBadDown=-1;
743              break;
744            }
745          }
746        }
747        // could change tests for anyBad
748        if (anyBadUp)
749          upImprovement=0.0;
750        if (anyBadDown)
751          downImprovement=0.0;
752        double way=0.0;
753        double improvement=0.0;
754        if (downImprovement>0.0&&currentValue>lowerValue) {
755          way=-1.0;
756          improvement = downImprovement;
757        } else if (upImprovement>0.0&&currentValue<upperValue) {
758          way=1.0;
759          improvement = upImprovement;
760        }
761        if (way) {
762          // can improve
763          double distance=COIN_DBL_MAX;
764          for (j=columnStart[iColumn];
765               j<columnStart[iColumn]+columnLength[iColumn];j++) {
766            int iRow = row[j];
767            double value = element[j]*way;
768            if (rowActivity[iRow]>rowUpper[iRow]+primalTolerance) {
769              // infeasible above
770              assert (value<0.0);
771              double gap = rowActivity[iRow]-rowUpper[iRow];
772              if (gap+value*distance<0.0) {
773                // If integer then has to move by 1
774                if (!isInteger)
775                  distance = -gap/value;
776                else
777                  distance = CoinMax(-gap/value,1.0);
778              }
779            } else if (rowActivity[iRow]<rowLower[iRow]-primalTolerance) {
780              // infeasible below
781              assert (value>0.0);
782              double gap = rowActivity[iRow]-rowLower[iRow];
783              if (gap+value*distance>0.0) {
784                // If integer then has to move by 1
785                if (!isInteger)
786                  distance = -gap/value;
787                else
788                  distance = CoinMax(-gap/value,1.0);
789              }
790            } else {
791              // feasible
792              if (value>0) {
793                double gap = rowActivity[iRow]-rowUpper[iRow];
794                if (gap+value*distance>0.0) 
795                  distance = -gap/value;
796              } else {
797                double gap = rowActivity[iRow]-rowLower[iRow];
798                if (gap+value*distance<0.0) 
799                  distance = -gap/value;
800              }
801            }
802          }
803          if (isInteger)
804            distance = floor(distance+1.05e-8);
805          if (!distance) {
806            // should never happen
807            //printf("zero distance in CbcRounding - debug\n");
808          }
809          //move
810          lastChange += improvement*distance;
811          distance *= way;
812          newSolution[iColumn] += distance;
813          newSolutionValue += direction*objective[iColumn]*distance;
814          for (j=columnStart[iColumn];
815               j<columnStart[iColumn]+columnLength[iColumn];j++) {
816            int iRow = row[j];
817            double value = element[j];
818            rowActivity[iRow] += distance*value;
819          }
820        }
821      }
822      penaltyChange += lastChange;
823    }
824    penalty -= penaltyChange;
825    if (penalty<1.0e-5*fabs(penaltyChange)) {
826      // recompute
827      penalty=0.0;
828      for (i=0;i<numberRows;i++) {
829        double value = rowActivity[i];
830        if (value<rowLower[i]-primalTolerance)
831          penalty += rowLower[i]-value;
832        else if (value>rowUpper[i]+primalTolerance)
833          penalty += value-rowUpper[i];
834      }
835    }
836  }
837
838  // Could also set SOS (using random) and repeat
839  if (!penalty) {
840    // See if we can do better
841    //seed_++;
842    //CoinSeedRandom(seed_);
843    // Random number between 0 and 1.
844    double randomNumber = CoinDrand48();
845    int iPass;
846    int start[2];
847    int end[2];
848    int iRandom = (int) (randomNumber*((double) numberIntegers));
849    start[0]=iRandom;
850    end[0]=numberIntegers;
851    start[1]=0;
852    end[1]=iRandom;
853    for (iPass=0;iPass<2;iPass++) {
854      int i;
855      for (i=start[iPass];i<end[iPass];i++) {
856        int iColumn = integerVariable[i];
857#ifndef NDEBUG
858        double value=newSolution[iColumn];
859        assert (fabs(floor(value+0.5)-value)<integerTolerance);
860#endif
861        double cost = direction * objective[iColumn];
862        double move=0.0;
863        if (cost>0.0)
864          move = -1.0;
865        else if (cost<0.0)
866          move=1.0;
867        while (move) {
868          bool good=true;
869          double newValue=newSolution[iColumn]+move;
870          if (newValue<lower[iColumn]-primalTolerance||
871              newValue>upper[iColumn]+primalTolerance) {
872            move=0.0;
873          } else {
874            // see if we can move
875            int j;
876            for (j=columnStart[iColumn];
877                 j<columnStart[iColumn]+columnLength[iColumn];j++) {
878              int iRow = row[j];
879              double newActivity = rowActivity[iRow] + move*element[j];
880              if (newActivity<rowLower[iRow]-primalTolerance||
881                  newActivity>rowUpper[iRow]+primalTolerance) {
882                good=false;
883                break;
884              }
885            }
886            if (good) {
887              newSolution[iColumn] = newValue;
888              newSolutionValue += move*cost;
889              int j;
890              for (j=columnStart[iColumn];
891                   j<columnStart[iColumn]+columnLength[iColumn];j++) {
892                int iRow = row[j];
893                rowActivity[iRow] += move*element[j];
894              }
895            } else {
896              move=0.0;
897            }
898          }
899        }
900      }
901    }
902    // Just in case of some stupidity
903    double objOffset=0.0;
904    solver->getDblParam(OsiObjOffset,objOffset);
905    newSolutionValue = -objOffset;
906    for ( i=0 ; i<numberColumns ; i++ )
907      newSolutionValue += objective[i]*newSolution[i];
908    newSolutionValue *= direction;
909    //printf("new solution value %g %g\n",newSolutionValue,solutionValue);
910    if (newSolutionValue<solutionValue) {
911      // paranoid check
912      memset(rowActivity,0,numberRows*sizeof(double));
913      for (i=0;i<numberColumns;i++) {
914        int j;
915        double value = newSolution[i];
916        if (value) {
917          for (j=columnStart[i];
918               j<columnStart[i]+columnLength[i];j++) {
919            int iRow=row[j];
920            rowActivity[iRow] += value*element[j];
921          }
922        }
923      }
924      // check was approximately feasible
925      bool feasible=true;
926      for (i=0;i<numberRows;i++) {
927        if(rowActivity[i]<rowLower[i]) {
928          if (rowActivity[i]<rowLower[i]-1000.0*primalTolerance)
929            feasible = false;
930        } else if(rowActivity[i]>rowUpper[i]) {
931          if (rowActivity[i]>rowUpper[i]+1000.0*primalTolerance)
932            feasible = false;
933        }
934      }
935      if (feasible) {
936        // new solution
937        memcpy(betterSolution,newSolution,numberColumns*sizeof(double));
938        solutionValue = newSolutionValue;
939        //printf("** Solution of %g found by rounding\n",newSolutionValue);
940        returnCode=1;
941      } else {
942        // Can easily happen
943        //printf("Debug CbcRounding giving bad solution\n");
944      }
945    }
946  }
947  delete [] newSolution;
948  delete [] rowActivity;
949  return returnCode;
950}
951// update model
952void CbcRounding::setModel(CbcModel * model)
953{
954  model_ = model;
955  // Get a copy of original matrix (and by row for rounding);
956  assert(model_->solver());
957  matrix_ = *model_->solver()->getMatrixByCol();
958  matrixByRow_ = *model_->solver()->getMatrixByRow();
959  // make sure model okay for heuristic
960  validate();
961}
962// Validate model i.e. sets when_ to 0 if necessary (may be NULL)
963void 
964CbcRounding::validate() 
965{
966  if (model_&&when()<10) {
967    if (model_->numberIntegers()!=
968        model_->numberObjects())
969      setWhen(0);
970  }
971}
972
973// Default Constructor
974CbcSerendipity::CbcSerendipity() 
975  :CbcHeuristic()
976{
977}
978
979// Constructor from model
980CbcSerendipity::CbcSerendipity(CbcModel & model)
981  :CbcHeuristic(model)
982{
983}
984
985// Destructor
986CbcSerendipity::~CbcSerendipity ()
987{
988}
989
990// Clone
991CbcHeuristic *
992CbcSerendipity::clone() const
993{
994  return new CbcSerendipity(*this);
995}
996// Create C++ lines to get to current state
997void 
998CbcSerendipity::generateCpp( FILE * fp) 
999{
1000  fprintf(fp,"0#include \"CbcHeuristic.hpp\"\n");
1001  fprintf(fp,"3  CbcSerendipity serendipity(*cbcModel);\n");
1002  CbcHeuristic::generateCpp(fp,"serendipity");
1003  fprintf(fp,"3  cbcModel->addHeuristic(&serendipity);\n");
1004}
1005
1006// Copy constructor
1007CbcSerendipity::CbcSerendipity(const CbcSerendipity & rhs)
1008:
1009  CbcHeuristic(rhs)
1010{
1011}
1012
1013// Assignment operator
1014CbcSerendipity & 
1015CbcSerendipity::operator=( const CbcSerendipity& rhs)
1016{
1017  if (this!=&rhs) {
1018    CbcHeuristic::operator=(rhs);
1019  }
1020  return *this;
1021}
1022
1023// Returns 1 if solution, 0 if not
1024int
1025CbcSerendipity::solution(double & solutionValue,
1026                         double * betterSolution)
1027{
1028  if (!model_)
1029    return 0;
1030  // get information on solver type
1031  OsiAuxInfo * auxInfo = model_->solver()->getAuxiliaryInfo();
1032  OsiBabSolver * auxiliaryInfo = dynamic_cast< OsiBabSolver *> (auxInfo);
1033  if (auxiliaryInfo)
1034    return auxiliaryInfo->solution(solutionValue,betterSolution,model_->solver()->getNumCols());
1035  else
1036    return 0;
1037}
1038// update model
1039void CbcSerendipity::setModel(CbcModel * model)
1040{
1041  model_ = model;
1042}
1043// Resets stuff if model changes
1044void 
1045CbcSerendipity::resetModel(CbcModel * model)
1046{
1047  model_ = model;
1048}
1049 
Note: See TracBrowser for help on using the repository browser.