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

Last change on this file since 502 was 502, checked in by forrest, 14 years ago

mostly for rins

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