source: branches/devel/Cbc/src/CbcHeuristic.cpp @ 539

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

for nonlinear

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