source: trunk/Cbc/src/CbcHeuristic.cpp @ 741

Last change on this file since 741 was 741, checked in by forrest, 13 years ago

event handling

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