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

Last change on this file since 616 was 607, checked in by lou, 13 years ago

Add names for heuristic, analogous to names for cut generators.

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 30.6 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      model.branchAndBound();
181      if (logLevel>1)
182        model_->messageHandler()->message(CBC_END_SUB,model_->messages())
183          << name
184          <<CoinMessageEol;
185      if (model.getMinimizationObjValue()<CoinMin(cutoff,1.0e30)) {
186        // solution
187        returnCode=model.isProvenOptimal() ? 3 : 1;
188        // post process
189        process.postProcess(*model.solver());
190        if (solver->isProvenOptimal()) {
191          // Solution now back in solver
192          int numberColumns = solver->getNumCols();
193          memcpy(newSolution,solver->getColSolution(),
194                 numberColumns*sizeof(double));
195          newSolutionValue = model.getMinimizationObjValue();
196        } else {
197          // odd - but no good
198          returnCode=0; // so will be infeasible
199        }
200      } else {
201        // no good
202        returnCode=model.isProvenInfeasible() ? 2 : 0; // so will be infeasible
203      }
204    }
205  } else {
206    returnCode=2; // infeasible finished
207  }
208  return returnCode;
209}
210
211// Default Constructor
212CbcRounding::CbcRounding() 
213  :CbcHeuristic()
214{
215  // matrix and row copy will automatically be empty
216  seed_=1;
217}
218
219// Constructor from model
220CbcRounding::CbcRounding(CbcModel & model)
221  :CbcHeuristic(model)
222{
223  // Get a copy of original matrix (and by row for rounding);
224  assert(model.solver());
225  matrix_ = *model.solver()->getMatrixByCol();
226  matrixByRow_ = *model.solver()->getMatrixByRow();
227  seed_=1;
228}
229
230// Destructor
231CbcRounding::~CbcRounding ()
232{
233}
234
235// Clone
236CbcHeuristic *
237CbcRounding::clone() const
238{
239  return new CbcRounding(*this);
240}
241// Create C++ lines to get to current state
242void 
243CbcRounding::generateCpp( FILE * fp) 
244{
245  CbcRounding other;
246  fprintf(fp,"0#include \"CbcHeuristic.hpp\"\n");
247  fprintf(fp,"3  CbcRounding rounding(*cbcModel);\n");
248  CbcHeuristic::generateCpp(fp,"rounding");
249  if (seed_!=other.seed_)
250    fprintf(fp,"3  rounding.setSeed(%d);\n",seed_);
251  else
252    fprintf(fp,"4  rounding.setSeed(%d);\n",seed_);
253  fprintf(fp,"3  cbcModel->addHeuristic(&rounding);\n");
254}
255
256// Copy constructor
257CbcRounding::CbcRounding(const CbcRounding & rhs)
258:
259  CbcHeuristic(rhs),
260  matrix_(rhs.matrix_),
261  matrixByRow_(rhs.matrixByRow_),
262  seed_(rhs.seed_)
263{
264}
265
266// Assignment operator
267CbcRounding & 
268CbcRounding::operator=( const CbcRounding& rhs)
269{
270  if (this!=&rhs) {
271    CbcHeuristic::operator=(rhs);
272    matrix_ = rhs.matrix_;
273    matrixByRow_ = rhs.matrixByRow_;
274    seed_ = rhs.seed_;
275  }
276  return *this;
277}
278
279// Resets stuff if model changes
280void 
281CbcRounding::resetModel(CbcModel * model)
282{
283  model_=model;
284  // Get a copy of original matrix (and by row for rounding);
285  assert(model_->solver());
286  matrix_ = *model_->solver()->getMatrixByCol();
287  matrixByRow_ = *model_->solver()->getMatrixByRow();
288}
289// See if rounding will give solution
290// Sets value of solution
291// Assumes rhs for original matrix still okay
292// At present only works with integers
293// Fix values if asked for
294// Returns 1 if solution, 0 if not
295int
296CbcRounding::solution(double & solutionValue,
297                         double * betterSolution)
298{
299
300  // See if to do
301  if (!when()||(when()%10==1&&model_->phase()!=1)||
302      (when()%10==2&&(model_->phase()!=2&&model_->phase()!=3)))
303    return 0; // switched off
304
305  OsiSolverInterface * solver = model_->solver();
306  const double * lower = solver->getColLower();
307  const double * upper = solver->getColUpper();
308  const double * rowLower = solver->getRowLower();
309  const double * rowUpper = solver->getRowUpper();
310  const double * solution = solver->getColSolution();
311  const double * objective = solver->getObjCoefficients();
312  double integerTolerance = model_->getDblParam(CbcModel::CbcIntegerTolerance);
313  double primalTolerance;
314  solver->getDblParam(OsiPrimalTolerance,primalTolerance);
315
316  int numberRows = matrix_.getNumRows();
317
318  int numberIntegers = model_->numberIntegers();
319  const int * integerVariable = model_->integerVariable();
320  int i;
321  double direction = solver->getObjSense();
322  double newSolutionValue = direction*solver->getObjValue();
323  int returnCode = 0;
324
325  // Column copy
326  const double * element = matrix_.getElements();
327  const int * row = matrix_.getIndices();
328  const CoinBigIndex * columnStart = matrix_.getVectorStarts();
329  const int * columnLength = matrix_.getVectorLengths();
330  // Row copy
331  const double * elementByRow = matrixByRow_.getElements();
332  const int * column = matrixByRow_.getIndices();
333  const CoinBigIndex * rowStart = matrixByRow_.getVectorStarts();
334  const int * rowLength = matrixByRow_.getVectorLengths();
335
336  // Get solution array for heuristic solution
337  int numberColumns = solver->getNumCols();
338  double * newSolution = new double [numberColumns];
339  memcpy(newSolution,solution,numberColumns*sizeof(double));
340
341  double * rowActivity = new double[numberRows];
342  memset(rowActivity,0,numberRows*sizeof(double));
343  for (i=0;i<numberColumns;i++) {
344    int j;
345    double value = newSolution[i];
346    if (value<lower[i]) {
347      value=lower[i];
348      newSolution[i]=value;
349    } else if (value>upper[i]) {
350      value=upper[i];
351      newSolution[i]=value;
352    }
353    if (value) {
354      for (j=columnStart[i];
355           j<columnStart[i]+columnLength[i];j++) {
356        int iRow=row[j];
357        rowActivity[iRow] += value*element[j];
358      }
359    }
360  }
361  // check was feasible - if not adjust (cleaning may move)
362  for (i=0;i<numberRows;i++) {
363    if(rowActivity[i]<rowLower[i]) {
364      //assert (rowActivity[i]>rowLower[i]-1000.0*primalTolerance);
365      rowActivity[i]=rowLower[i];
366    } else if(rowActivity[i]>rowUpper[i]) {
367      //assert (rowActivity[i]<rowUpper[i]+1000.0*primalTolerance);
368      rowActivity[i]=rowUpper[i];
369    }
370  }
371  for (i=0;i<numberIntegers;i++) {
372    int iColumn = integerVariable[i];
373    double value=newSolution[iColumn];
374    if (fabs(floor(value+0.5)-value)>integerTolerance) {
375      double below = floor(value);
376      double newValue=newSolution[iColumn];
377      double cost = direction * objective[iColumn];
378      double move;
379      if (cost>0.0) {
380        // try up
381        move = 1.0 -(value-below);
382      } else if (cost<0.0) {
383        // try down
384        move = below-value;
385      } else {
386        // won't be able to move unless we can grab another variable
387        double randomNumber = CoinDrand48();
388        // which way?
389        if (randomNumber<0.5) 
390          move = below-value;
391        else
392          move = 1.0 -(value-below);
393      }
394      newValue += move;
395      newSolution[iColumn] = newValue;
396      newSolutionValue += move*cost;
397      int j;
398      for (j=columnStart[iColumn];
399           j<columnStart[iColumn]+columnLength[iColumn];j++) {
400        int iRow = row[j];
401        rowActivity[iRow] += move*element[j];
402      }
403    }
404  }
405
406  double penalty=0.0;
407  const char * integerType = model_->integerType();
408  // see if feasible - just using singletons
409  for (i=0;i<numberRows;i++) {
410    double value = rowActivity[i];
411    double thisInfeasibility=0.0;
412    if (value<rowLower[i]-primalTolerance)
413      thisInfeasibility = value-rowLower[i];
414    else if (value>rowUpper[i]+primalTolerance)
415      thisInfeasibility = value-rowUpper[i];
416    if (thisInfeasibility) {
417      // See if there are any slacks I can use to fix up
418      // maybe put in coding for multiple slacks?
419      double bestCost = 1.0e50;
420      int k;
421      int iBest=-1;
422      double addCost=0.0;
423      double newValue=0.0;
424      double changeRowActivity=0.0;
425      double absInfeasibility = fabs(thisInfeasibility);
426      for (k=rowStart[i];k<rowStart[i]+rowLength[i];k++) {
427        int iColumn = column[k];
428        // See if all elements help
429        if (columnLength[iColumn]==1) {
430          double currentValue = newSolution[iColumn];
431          double elementValue = elementByRow[k];
432          double lowerValue = lower[iColumn];
433          double upperValue = upper[iColumn];
434          double gap = rowUpper[i]-rowLower[i];
435          double absElement=fabs(elementValue);
436          if (thisInfeasibility*elementValue>0.0) {
437            // we want to reduce
438            if ((currentValue-lowerValue)*absElement>=absInfeasibility) {
439              // possible - check if integer
440              double distance = absInfeasibility/absElement;
441              double thisCost = -direction*objective[iColumn]*distance;
442              if (integerType[iColumn]) {
443                distance = ceil(distance-primalTolerance);
444                if (currentValue-distance>=lowerValue-primalTolerance) {
445                  if (absInfeasibility-distance*absElement< -gap-primalTolerance)
446                    thisCost=1.0e100; // no good
447                  else
448                    thisCost = -direction*objective[iColumn]*distance;
449                } else {
450                  thisCost=1.0e100; // no good
451                }
452              }
453              if (thisCost<bestCost) {
454                bestCost=thisCost;
455                iBest=iColumn;
456                addCost = thisCost;
457                newValue = currentValue-distance;
458                changeRowActivity = -distance*elementValue;
459              }
460            }
461          } else {
462            // we want to increase
463            if ((upperValue-currentValue)*absElement>=absInfeasibility) {
464              // possible - check if integer
465              double distance = absInfeasibility/absElement;
466              double thisCost = direction*objective[iColumn]*distance;
467              if (integerType[iColumn]) {
468                distance = ceil(distance-1.0e-7);
469                assert (currentValue-distance<=upperValue+primalTolerance);
470                if (absInfeasibility-distance*absElement< -gap-primalTolerance)
471                  thisCost=1.0e100; // no good
472                else
473                  thisCost = direction*objective[iColumn]*distance;
474              }
475              if (thisCost<bestCost) {
476                bestCost=thisCost;
477                iBest=iColumn;
478                addCost = thisCost;
479                newValue = currentValue+distance;
480                changeRowActivity = distance*elementValue;
481              }
482            }
483          }
484        }
485      }
486      if (iBest>=0) {
487        /*printf("Infeasibility of %g on row %d cost %g\n",
488          thisInfeasibility,i,addCost);*/
489        newSolution[iBest]=newValue;
490        thisInfeasibility=0.0;
491        newSolutionValue += addCost;
492        rowActivity[i] += changeRowActivity;
493      }
494      penalty += fabs(thisInfeasibility);
495    }
496  }
497  if (penalty) {
498    // see if feasible using any
499    // first continuous
500    double penaltyChange=0.0;
501    int iColumn;
502    for (iColumn=0;iColumn<numberColumns;iColumn++) {
503      if (integerType[iColumn])
504        continue;
505      double currentValue = newSolution[iColumn];
506      double lowerValue = lower[iColumn];
507      double upperValue = upper[iColumn];
508      int j;
509      int anyBadDown=0;
510      int anyBadUp=0;
511      double upImprovement=0.0;
512      double downImprovement=0.0;
513      for (j=columnStart[iColumn];
514           j<columnStart[iColumn]+columnLength[iColumn];j++) {
515        int iRow = row[j];
516        if (rowUpper[iRow]>rowLower[iRow]) {
517          double value = element[j];
518          if (rowActivity[iRow]>rowUpper[iRow]+primalTolerance) {
519            // infeasible above
520            downImprovement += value;
521            upImprovement -= value;
522            if (value>0.0) 
523              anyBadUp++;
524            else 
525              anyBadDown++;
526          } else if (rowActivity[iRow]>rowUpper[iRow]-primalTolerance) {
527            // feasible at ub
528            if (value>0.0) {
529              upImprovement -= value;
530              anyBadUp++;
531            } else {
532              downImprovement += value;
533              anyBadDown++;
534            }
535          } else if (rowActivity[iRow]>rowLower[iRow]+primalTolerance) {
536            // feasible in interior
537          } else if (rowActivity[iRow]>rowLower[iRow]-primalTolerance) {
538            // feasible at lb
539            if (value<0.0) {
540              upImprovement += value;
541              anyBadUp++;
542            } else {
543              downImprovement -= value;
544              anyBadDown++;
545            }
546          } else {
547            // infeasible below
548            downImprovement -= value;
549            upImprovement += value;
550            if (value<0.0) 
551              anyBadUp++;
552            else 
553              anyBadDown++;
554          }
555        } else {
556          // equality row
557          double value = element[j];
558          if (rowActivity[iRow]>rowUpper[iRow]+primalTolerance) {
559            // infeasible above
560            downImprovement += value;
561            upImprovement -= value;
562            if (value>0.0) 
563              anyBadUp++;
564            else 
565              anyBadDown++;
566          } else if (rowActivity[iRow]<rowLower[iRow]-primalTolerance) {
567            // infeasible below
568            downImprovement -= value;
569            upImprovement += value;
570            if (value<0.0) 
571              anyBadUp++;
572            else 
573              anyBadDown++;
574          } else {
575            // feasible - no good
576            anyBadUp=-1;
577            anyBadDown=-1;
578            break;
579          }
580        }
581      }
582      // could change tests for anyBad
583      if (anyBadUp)
584        upImprovement=0.0;
585      if (anyBadDown)
586        downImprovement=0.0;
587      double way=0.0;
588      double improvement=0.0;
589      if (downImprovement>0.0&&currentValue>lowerValue) {
590        way=-1.0;
591        improvement = downImprovement;
592      } else if (upImprovement>0.0&&currentValue<upperValue) {
593        way=1.0;
594        improvement = upImprovement;
595      }
596      if (way) {
597        // can improve
598        double distance;
599        if (way>0.0)
600          distance = upperValue-currentValue;
601        else
602          distance = currentValue-lowerValue;
603        for (j=columnStart[iColumn];
604             j<columnStart[iColumn]+columnLength[iColumn];j++) {
605          int iRow = row[j];
606          double value = element[j]*way;
607          if (rowActivity[iRow]>rowUpper[iRow]+primalTolerance) {
608            // infeasible above
609            assert (value<0.0);
610            double gap = rowActivity[iRow]-rowUpper[iRow];
611            if (gap+value*distance<0.0) 
612              distance = -gap/value;
613          } else if (rowActivity[iRow]<rowLower[iRow]-primalTolerance) {
614            // infeasible below
615            assert (value>0.0);
616            double gap = rowActivity[iRow]-rowLower[iRow];
617            if (gap+value*distance>0.0) 
618              distance = -gap/value;
619          } else {
620            // feasible
621            if (value>0) {
622              double gap = rowActivity[iRow]-rowUpper[iRow];
623              if (gap+value*distance>0.0) 
624              distance = -gap/value;
625            } else {
626              double gap = rowActivity[iRow]-rowLower[iRow];
627              if (gap+value*distance<0.0) 
628                distance = -gap/value;
629            }
630          }
631        }
632        //move
633        penaltyChange += improvement*distance;
634        distance *= way;
635        newSolution[iColumn] += distance;
636        newSolutionValue += direction*objective[iColumn]*distance;
637        for (j=columnStart[iColumn];
638             j<columnStart[iColumn]+columnLength[iColumn];j++) {
639          int iRow = row[j];
640          double value = element[j];
641          rowActivity[iRow] += distance*value;
642        }
643      }
644    }
645    // and now all if improving
646    double lastChange= penaltyChange ? 1.0 : 0.0;
647    while (lastChange>1.0e-2) {
648      lastChange=0;
649      for (iColumn=0;iColumn<numberColumns;iColumn++) {
650        bool isInteger = integerType[iColumn];
651        double currentValue = newSolution[iColumn];
652        double lowerValue = lower[iColumn];
653        double upperValue = upper[iColumn];
654        int j;
655        int anyBadDown=0;
656        int anyBadUp=0;
657        double upImprovement=0.0;
658        double downImprovement=0.0;
659        for (j=columnStart[iColumn];
660             j<columnStart[iColumn]+columnLength[iColumn];j++) {
661          int iRow = row[j];
662          double value = element[j];
663          if (isInteger) {
664            if (value>0.0) {
665              if (rowActivity[iRow]+value>rowUpper[iRow]+primalTolerance)
666                anyBadUp++;
667              if (rowActivity[iRow]-value<rowLower[iRow]-primalTolerance)
668                anyBadDown++;
669            } else {
670              if (rowActivity[iRow]-value>rowUpper[iRow]+primalTolerance)
671                anyBadDown++;
672              if (rowActivity[iRow]+value<rowLower[iRow]-primalTolerance)
673                anyBadUp++;
674            }
675          }
676          if (rowUpper[iRow]>rowLower[iRow]) {
677            if (rowActivity[iRow]>rowUpper[iRow]+primalTolerance) {
678              // infeasible above
679              downImprovement += value;
680              upImprovement -= value;
681              if (value>0.0) 
682                anyBadUp++;
683              else 
684                anyBadDown++;
685            } else if (rowActivity[iRow]>rowUpper[iRow]-primalTolerance) {
686              // feasible at ub
687              if (value>0.0) {
688                upImprovement -= value;
689                anyBadUp++;
690              } else {
691                downImprovement += value;
692                anyBadDown++;
693              }
694            } else if (rowActivity[iRow]>rowLower[iRow]+primalTolerance) {
695              // feasible in interior
696            } else if (rowActivity[iRow]>rowLower[iRow]-primalTolerance) {
697              // feasible at lb
698              if (value<0.0) {
699                upImprovement += value;
700                anyBadUp++;
701              } else {
702                downImprovement -= value;
703                anyBadDown++;
704              }
705            } else {
706              // infeasible below
707              downImprovement -= value;
708              upImprovement += value;
709              if (value<0.0) 
710                anyBadUp++;
711              else 
712                anyBadDown++;
713            }
714          } else {
715            // equality row
716            if (rowActivity[iRow]>rowUpper[iRow]+primalTolerance) {
717              // infeasible above
718              downImprovement += value;
719              upImprovement -= value;
720              if (value>0.0) 
721                anyBadUp++;
722              else 
723                anyBadDown++;
724            } else if (rowActivity[iRow]<rowLower[iRow]-primalTolerance) {
725              // infeasible below
726              downImprovement -= value;
727              upImprovement += value;
728              if (value<0.0) 
729                anyBadUp++;
730              else 
731                anyBadDown++;
732            } else {
733              // feasible - no good
734              anyBadUp=-1;
735              anyBadDown=-1;
736              break;
737            }
738          }
739        }
740        // could change tests for anyBad
741        if (anyBadUp)
742          upImprovement=0.0;
743        if (anyBadDown)
744          downImprovement=0.0;
745        double way=0.0;
746        double improvement=0.0;
747        if (downImprovement>0.0&&currentValue>lowerValue) {
748          way=-1.0;
749          improvement = downImprovement;
750        } else if (upImprovement>0.0&&currentValue<upperValue) {
751          way=1.0;
752          improvement = upImprovement;
753        }
754        if (way) {
755          // can improve
756          double distance=COIN_DBL_MAX;
757          for (j=columnStart[iColumn];
758               j<columnStart[iColumn]+columnLength[iColumn];j++) {
759            int iRow = row[j];
760            double value = element[j]*way;
761            if (rowActivity[iRow]>rowUpper[iRow]+primalTolerance) {
762              // infeasible above
763              assert (value<0.0);
764              double gap = rowActivity[iRow]-rowUpper[iRow];
765              if (gap+value*distance<0.0) {
766                // If integer then has to move by 1
767                if (!isInteger)
768                  distance = -gap/value;
769                else
770                  distance = CoinMax(-gap/value,1.0);
771              }
772            } else if (rowActivity[iRow]<rowLower[iRow]-primalTolerance) {
773              // infeasible below
774              assert (value>0.0);
775              double gap = rowActivity[iRow]-rowLower[iRow];
776              if (gap+value*distance>0.0) {
777                // If integer then has to move by 1
778                if (!isInteger)
779                  distance = -gap/value;
780                else
781                  distance = CoinMax(-gap/value,1.0);
782              }
783            } else {
784              // feasible
785              if (value>0) {
786                double gap = rowActivity[iRow]-rowUpper[iRow];
787                if (gap+value*distance>0.0) 
788                  distance = -gap/value;
789              } else {
790                double gap = rowActivity[iRow]-rowLower[iRow];
791                if (gap+value*distance<0.0) 
792                  distance = -gap/value;
793              }
794            }
795          }
796          if (isInteger)
797            distance = floor(distance+1.05e-8);
798          if (!distance) {
799            // should never happen
800            //printf("zero distance in CbcRounding - debug\n");
801          }
802          //move
803          lastChange += improvement*distance;
804          distance *= way;
805          newSolution[iColumn] += distance;
806          newSolutionValue += direction*objective[iColumn]*distance;
807          for (j=columnStart[iColumn];
808               j<columnStart[iColumn]+columnLength[iColumn];j++) {
809            int iRow = row[j];
810            double value = element[j];
811            rowActivity[iRow] += distance*value;
812          }
813        }
814      }
815      penaltyChange += lastChange;
816    }
817    penalty -= penaltyChange;
818    if (penalty<1.0e-5*fabs(penaltyChange)) {
819      // recompute
820      penalty=0.0;
821      for (i=0;i<numberRows;i++) {
822        double value = rowActivity[i];
823        if (value<rowLower[i]-primalTolerance)
824          penalty += rowLower[i]-value;
825        else if (value>rowUpper[i]+primalTolerance)
826          penalty += value-rowUpper[i];
827      }
828    }
829  }
830
831  // Could also set SOS (using random) and repeat
832  if (!penalty) {
833    // See if we can do better
834    //seed_++;
835    //CoinSeedRandom(seed_);
836    // Random number between 0 and 1.
837    double randomNumber = CoinDrand48();
838    int iPass;
839    int start[2];
840    int end[2];
841    int iRandom = (int) (randomNumber*((double) numberIntegers));
842    start[0]=iRandom;
843    end[0]=numberIntegers;
844    start[1]=0;
845    end[1]=iRandom;
846    for (iPass=0;iPass<2;iPass++) {
847      int i;
848      for (i=start[iPass];i<end[iPass];i++) {
849        int iColumn = integerVariable[i];
850#ifndef NDEBUG
851        double value=newSolution[iColumn];
852        assert (fabs(floor(value+0.5)-value)<integerTolerance);
853#endif
854        double cost = direction * objective[iColumn];
855        double move=0.0;
856        if (cost>0.0)
857          move = -1.0;
858        else if (cost<0.0)
859          move=1.0;
860        while (move) {
861          bool good=true;
862          double newValue=newSolution[iColumn]+move;
863          if (newValue<lower[iColumn]-primalTolerance||
864              newValue>upper[iColumn]+primalTolerance) {
865            move=0.0;
866          } else {
867            // see if we can move
868            int j;
869            for (j=columnStart[iColumn];
870                 j<columnStart[iColumn]+columnLength[iColumn];j++) {
871              int iRow = row[j];
872              double newActivity = rowActivity[iRow] + move*element[j];
873              if (newActivity<rowLower[iRow]-primalTolerance||
874                  newActivity>rowUpper[iRow]+primalTolerance) {
875                good=false;
876                break;
877              }
878            }
879            if (good) {
880              newSolution[iColumn] = newValue;
881              newSolutionValue += move*cost;
882              int j;
883              for (j=columnStart[iColumn];
884                   j<columnStart[iColumn]+columnLength[iColumn];j++) {
885                int iRow = row[j];
886                rowActivity[iRow] += move*element[j];
887              }
888            } else {
889              move=0.0;
890            }
891          }
892        }
893      }
894    }
895    if (newSolutionValue<solutionValue) {
896      // paranoid check
897      memset(rowActivity,0,numberRows*sizeof(double));
898      for (i=0;i<numberColumns;i++) {
899        int j;
900        double value = newSolution[i];
901        if (value) {
902          for (j=columnStart[i];
903               j<columnStart[i]+columnLength[i];j++) {
904            int iRow=row[j];
905            rowActivity[iRow] += value*element[j];
906          }
907        }
908      }
909      // check was approximately feasible
910      bool feasible=true;
911      for (i=0;i<numberRows;i++) {
912        if(rowActivity[i]<rowLower[i]) {
913          if (rowActivity[i]<rowLower[i]-1000.0*primalTolerance)
914            feasible = false;
915        } else if(rowActivity[i]>rowUpper[i]) {
916          if (rowActivity[i]>rowUpper[i]+1000.0*primalTolerance)
917            feasible = false;
918        }
919      }
920      if (feasible) {
921        // new solution
922        memcpy(betterSolution,newSolution,numberColumns*sizeof(double));
923        solutionValue = newSolutionValue;
924        //printf("** Solution of %g found by rounding\n",newSolutionValue);
925        returnCode=1;
926      } else {
927        // Can easily happen
928        //printf("Debug CbcRounding giving bad solution\n");
929      }
930    }
931  }
932  delete [] newSolution;
933  delete [] rowActivity;
934  return returnCode;
935}
936// update model
937void CbcRounding::setModel(CbcModel * model)
938{
939  model_ = model;
940  // Get a copy of original matrix (and by row for rounding);
941  assert(model_->solver());
942  matrix_ = *model_->solver()->getMatrixByCol();
943  matrixByRow_ = *model_->solver()->getMatrixByRow();
944  // make sure model okay for heuristic
945  validate();
946}
947// Validate model i.e. sets when_ to 0 if necessary (may be NULL)
948void 
949CbcRounding::validate() 
950{
951  if (model_&&when()<10) {
952    if (model_->numberIntegers()!=
953        model_->numberObjects())
954      setWhen(0);
955  }
956}
957
958// Default Constructor
959CbcSerendipity::CbcSerendipity() 
960  :CbcHeuristic()
961{
962}
963
964// Constructor from model
965CbcSerendipity::CbcSerendipity(CbcModel & model)
966  :CbcHeuristic(model)
967{
968}
969
970// Destructor
971CbcSerendipity::~CbcSerendipity ()
972{
973}
974
975// Clone
976CbcHeuristic *
977CbcSerendipity::clone() const
978{
979  return new CbcSerendipity(*this);
980}
981// Create C++ lines to get to current state
982void 
983CbcSerendipity::generateCpp( FILE * fp) 
984{
985  fprintf(fp,"0#include \"CbcHeuristic.hpp\"\n");
986  fprintf(fp,"3  CbcSerendipity serendipity(*cbcModel);\n");
987  CbcHeuristic::generateCpp(fp,"serendipity");
988  fprintf(fp,"3  cbcModel->addHeuristic(&serendipity);\n");
989}
990
991// Copy constructor
992CbcSerendipity::CbcSerendipity(const CbcSerendipity & rhs)
993:
994  CbcHeuristic(rhs)
995{
996}
997
998// Assignment operator
999CbcSerendipity & 
1000CbcSerendipity::operator=( const CbcSerendipity& rhs)
1001{
1002  if (this!=&rhs) {
1003    CbcHeuristic::operator=(rhs);
1004  }
1005  return *this;
1006}
1007
1008// Returns 1 if solution, 0 if not
1009int
1010CbcSerendipity::solution(double & solutionValue,
1011                         double * betterSolution)
1012{
1013  if (!model_)
1014    return 0;
1015  // get information on solver type
1016  OsiAuxInfo * auxInfo = model_->solver()->getAuxiliaryInfo();
1017  OsiBabSolver * auxiliaryInfo = dynamic_cast< OsiBabSolver *> (auxInfo);
1018  if (auxiliaryInfo)
1019    return auxiliaryInfo->solution(solutionValue,betterSolution,model_->solver()->getNumCols());
1020  else
1021    return 0;
1022}
1023// update model
1024void CbcSerendipity::setModel(CbcModel * model)
1025{
1026  model_ = model;
1027}
1028// Resets stuff if model changes
1029void 
1030CbcSerendipity::resetModel(CbcModel * model)
1031{
1032  model_ = model;
1033}
1034 
Note: See TracBrowser for help on using the repository browser.