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

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

max times and fix bug in fpump

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