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

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

changes for postprocess loop and LOS

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