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

Last change on this file since 355 was 355, checked in by andreasw, 13 years ago

undid last commit (patches incorrectly applied)

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