source: trunk/CbcHeuristic.cpp @ 277

Last change on this file since 277 was 277, checked in by lou, 14 years ago

Revise cbc build process for better control of solvers included in build
(COIN_USE_XXX -> CBC_USE_XXX). Add CbcEventHandler? for independence from
clp.

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