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

Last change on this file since 356 was 356, checked in by ladanyi, 13 years ago

finishing conversion to svn

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