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

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

to stop heuristics faster

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