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

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

for deterministic parallel
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 43.3 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 "CbcHeuristicFPump.hpp"
21#include "CbcStrategy.hpp"
22#include "CglPreProcess.hpp"
23#include "OsiAuxInfo.hpp"
24#include "OsiPresolve.hpp"
25
26// Default Constructor
27CbcHeuristic::CbcHeuristic() 
28  :model_(NULL),
29   when_(2),
30   numberNodes_(200),
31   feasibilityPumpOptions_(-1),
32   fractionSmall_(1.0),
33   heuristicName_("Unknown")
34{
35  // As CbcHeuristic virtual need to modify .cpp if above change
36}
37
38// Constructor from model
39CbcHeuristic::CbcHeuristic(CbcModel & model)
40:
41  model_(&model),
42  when_(2),
43  numberNodes_(200),
44  feasibilityPumpOptions_(-1),
45  fractionSmall_(1.0),
46  heuristicName_("Unknown")
47{
48  // As CbcHeuristic virtual need to modify .cpp if above change
49}
50// Copy constructor
51CbcHeuristic::CbcHeuristic(const CbcHeuristic & rhs)
52:
53  model_(rhs.model_),
54  when_(rhs.when_),
55  numberNodes_(rhs.numberNodes_),
56  feasibilityPumpOptions_(rhs.feasibilityPumpOptions_),
57  fractionSmall_(rhs.fractionSmall_),
58  randomNumberGenerator_(rhs.randomNumberGenerator_),
59  heuristicName_(rhs.heuristicName_)
60{
61}
62// Assignment operator
63CbcHeuristic & 
64CbcHeuristic::operator=( const CbcHeuristic& rhs)
65{
66  if (this!=&rhs) {
67    model_ = rhs.model_;
68    when_ = rhs.when_;
69    numberNodes_ = rhs.numberNodes_;
70    feasibilityPumpOptions_ = rhs.feasibilityPumpOptions_;
71    fractionSmall_ = rhs.fractionSmall_;
72    randomNumberGenerator_ = rhs.randomNumberGenerator_;
73    heuristicName_ = rhs.heuristicName_ ;
74  }
75  return *this;
76}
77
78// Resets stuff if model changes
79void 
80CbcHeuristic::resetModel(CbcModel * model)
81{
82  model_=model;
83}
84// Set seed
85void
86CbcHeuristic::setSeed(int value)
87{
88  randomNumberGenerator_.setSeed(value);
89}
90
91// Create C++ lines to get to current state
92void 
93CbcHeuristic::generateCpp( FILE * fp, const char * heuristic) 
94{
95  // hard coded as CbcHeuristic virtual
96  if (when_!=2)
97    fprintf(fp,"3  %s.setWhen(%d);\n",heuristic,when_);
98  else
99    fprintf(fp,"4  %s.setWhen(%d);\n",heuristic,when_);
100  if (numberNodes_!=200)
101    fprintf(fp,"3  %s.setNumberNodes(%d);\n",heuristic,numberNodes_);
102  else
103    fprintf(fp,"4  %s.setNumberNodes(%d);\n",heuristic,numberNodes_);
104  if (fractionSmall_!=1.0)
105    fprintf(fp,"3  %s.setFractionSmall(%g);\n",heuristic,fractionSmall_);
106  else
107    fprintf(fp,"4  %s.setFractionSmall(%g);\n",heuristic,fractionSmall_);
108  if (heuristicName_ != "Unknown")
109    fprintf(fp,"3  %s.setHeuristicName(\"%s\");\n",
110            heuristic,heuristicName_.c_str()) ;
111  else
112    fprintf(fp,"4  %s.setHeuristicName(\"%s\");\n",
113            heuristic,heuristicName_.c_str()) ;
114}
115// Destructor
116CbcHeuristic::~CbcHeuristic ()
117{
118}
119
120// update model
121void CbcHeuristic::setModel(CbcModel * model)
122{
123  model_ = model;
124}
125#ifdef COIN_DEVELOP
126extern bool getHistoryStatistics_;
127#endif
128// Do mini branch and bound (return 1 if solution)
129int 
130CbcHeuristic::smallBranchAndBound(OsiSolverInterface * solver,int numberNodes,
131                                  double * newSolution, double & newSolutionValue,
132                                  double cutoff, std::string name) const
133{
134#ifdef COIN_HAS_CLP
135  OsiClpSolverInterface * osiclp = dynamic_cast< OsiClpSolverInterface*> (solver);
136  if (osiclp&&(osiclp->specialOptions()&65536)==0) {
137    // go faster stripes
138    if (osiclp->getNumRows()<300&&osiclp->getNumCols()<500) {
139      osiclp->setupForRepeatedUse(2,0);
140    } else {
141      osiclp->setupForRepeatedUse(0,0);
142    }
143    // Turn this off if you get problems
144    // Used to be automatically set
145    osiclp->setSpecialOptions(osiclp->specialOptions()|(128+64));
146    ClpSimplex * lpSolver = osiclp->getModelPtr();
147    lpSolver->setSpecialOptions(lpSolver->specialOptions()|0x01000000); // say is Cbc (and in branch and bound)
148  }
149#endif
150#ifdef COIN_DEVELOP
151  getHistoryStatistics_=false;
152#endif
153  int status=0;
154  int logLevel = model_->logLevel();
155  char generalPrint[100];
156  // Do presolve to see if possible
157  int numberColumns = solver->getNumCols();
158  char * reset = NULL;
159  int returnCode=1;
160  {
161    int saveLogLevel = solver->messageHandler()->logLevel();
162    if (saveLogLevel==1)
163      solver->messageHandler()->setLogLevel(0);
164    OsiPresolve * pinfo = new OsiPresolve();
165    int presolveActions=0;
166    // Allow dual stuff on integers
167    presolveActions=1;
168    // Do not allow all +1 to be tampered with
169    //if (allPlusOnes)
170    //presolveActions |= 2;
171    // allow transfer of costs
172    // presolveActions |= 4;
173    pinfo->setPresolveActions(presolveActions);
174    OsiSolverInterface * presolvedModel = pinfo->presolvedModel(*solver,1.0e-8,true,2);
175    delete pinfo;
176    // see if too big
177    double before = 2*solver->getNumRows()+solver->getNumCols();
178    if (presolvedModel) {
179      int afterRows = presolvedModel->getNumRows();
180      int afterCols = presolvedModel->getNumCols();
181      delete presolvedModel;
182      double after = 2*afterRows+afterCols;
183      if (after>fractionSmall_*before&&after>300) {
184        // Need code to try again to compress further using used
185        const int * used =  model_->usedInSolution();
186        int maxUsed=0;
187        int iColumn;
188        const double * lower = solver->getColLower();
189        const double * upper = solver->getColUpper();
190        for (iColumn=0;iColumn<numberColumns;iColumn++) {
191          if (upper[iColumn]>lower[iColumn]) {
192            if (solver->isBinary(iColumn))
193              maxUsed = CoinMax(maxUsed,used[iColumn]);
194          }
195        }
196        if (maxUsed) {
197          reset = new char [numberColumns];
198          int nFix=0;
199          for (iColumn=0;iColumn<numberColumns;iColumn++) {
200            reset[iColumn]=0;
201            if (upper[iColumn]>lower[iColumn]) {
202              if (solver->isBinary(iColumn)&&used[iColumn]==maxUsed) {
203                bool setValue=true;
204                if (maxUsed==1) {
205                  double randomNumber = randomNumberGenerator_.randomDouble();
206                  if (randomNumber>0.3)
207                    setValue=false;
208                }
209                if (setValue) {
210                  reset[iColumn]=1;
211                  solver->setColLower(iColumn,1.0);
212                  nFix++;
213                }
214              }
215            }
216          }
217          pinfo = new OsiPresolve();
218          presolveActions=0;
219          // Allow dual stuff on integers
220          presolveActions=1;
221          // Do not allow all +1 to be tampered with
222          //if (allPlusOnes)
223          //presolveActions |= 2;
224          // allow transfer of costs
225          // presolveActions |= 4;
226          pinfo->setPresolveActions(presolveActions);
227          presolvedModel = pinfo->presolvedModel(*solver,1.0e-8,true,2);
228          delete pinfo;
229          if(presolvedModel) {
230            // see if too big
231            int afterRows2 = presolvedModel->getNumRows();
232            int afterCols2 = presolvedModel->getNumCols();
233            delete presolvedModel;
234            double after = 2*afterRows2+afterCols2;
235            if (after>fractionSmall_*before&&after>300) {
236              sprintf(generalPrint,"Full problem %d rows %d columns, reduced to %d rows %d columns - %d fixed gives %d, %d - still too large",
237                      solver->getNumRows(),solver->getNumCols(),
238                      afterRows,afterCols,nFix,afterRows2,afterCols2);
239            } else {
240              sprintf(generalPrint,"Full problem %d rows %d columns, reduced to %d rows %d columns - %d fixed gives %d, %d - ok now",
241                      solver->getNumRows(),solver->getNumCols(),
242                      afterRows,afterCols,nFix,afterRows2,afterCols2);
243            }
244            model_->messageHandler()->message(CBC_FPUMP1,model_->messages())
245              << generalPrint
246              <<CoinMessageEol;
247          } else {
248            returnCode=-1; // infeasible
249          }
250        }
251      }
252    } else {
253      returnCode=-1; // infeasible
254    }
255    solver->messageHandler()->setLogLevel(saveLogLevel);
256  }
257  if (returnCode==-1) {
258    delete [] reset;
259#ifdef COIN_DEVELOP
260    getHistoryStatistics_=true;
261#endif
262    return returnCode;
263  }
264  // Reduce printout
265  solver->setHintParam(OsiDoReducePrint,true,OsiHintTry);
266  solver->setHintParam(OsiDoPresolveInInitial,false,OsiHintTry);
267  solver->setDblParam(OsiDualObjectiveLimit,cutoff*solver->getObjSense());
268  solver->initialSolve();
269  if (solver->isProvenOptimal()) {
270    CglPreProcess process;
271    /* Do not try and produce equality cliques and
272       do up to 2 passes */
273    if (logLevel<=1)
274      process.messageHandler()->setLogLevel(0);
275    OsiSolverInterface * solver2= process.preProcessNonDefault(*solver,false,2);
276    if (!solver2) {
277      if (logLevel>1)
278        printf("Pre-processing says infeasible\n");
279      returnCode=2; // so will be infeasible
280    } else {
281      // see if too big
282      double before = 2*solver->getNumRows()+solver->getNumCols();
283      double after = 2*solver2->getNumRows()+solver2->getNumCols();
284      if (after>fractionSmall_*before&&after>300) {
285        sprintf(generalPrint,"Full problem %d rows %d columns, reduced to %d rows %d columns - too large",
286                solver->getNumRows(),solver->getNumCols(),
287                solver2->getNumRows(),solver2->getNumCols());
288        model_->messageHandler()->message(CBC_FPUMP1,model_->messages())
289          << generalPrint
290          <<CoinMessageEol;
291        returnCode = -1;
292      } else {
293        sprintf(generalPrint,"Full problem %d rows %d columns, reduced to %d rows %d columns",
294                solver->getNumRows(),solver->getNumCols(),
295                solver2->getNumRows(),solver2->getNumCols());
296        model_->messageHandler()->message(CBC_FPUMP1,model_->messages())
297          << generalPrint
298          <<CoinMessageEol;
299      }
300      if (returnCode==1) {
301        solver2->resolve();
302        CbcModel model(*solver2);
303        if (logLevel<=1)
304          model.setLogLevel(0);
305        else
306          model.setLogLevel(logLevel);
307        if (feasibilityPumpOptions_>=0) {
308          CbcHeuristicFPump heuristic4;
309          int pumpTune=feasibilityPumpOptions_;
310          if (pumpTune>0) {
311            /*
312              >=10000000 for using obj
313              >=1000000 use as accumulate switch
314              >=1000 use index+1 as number of large loops
315              >=100 use 0.05 objvalue as increment
316              >=10 use +0.1 objvalue for cutoff (add)
317              1 == fix ints at bounds, 2 fix all integral ints, 3 and continuous at bounds
318              4 and static continuous, 5 as 3 but no internal integers
319              6 as 3 but all slack basis!
320            */
321            double value = solver2->getObjSense()*solver2->getObjValue();
322            int w = pumpTune/10;
323            int c = w % 10;
324            w /= 10;
325            int i = w % 10;
326            w /= 10;
327            int r = w;
328            int accumulate = r/1000;
329            r -= 1000*accumulate;
330            if (accumulate>=10) {
331              int which = accumulate/10;
332              accumulate -= 10*which;
333              which--;
334              // weights and factors
335              double weight[]={0.1,0.1,0.5,0.5,1.0,1.0,5.0,5.0};
336              double factor[] = {0.1,0.5,0.1,0.5,0.1,0.5,0.1,0.5};
337              heuristic4.setInitialWeight(weight[which]);
338              heuristic4.setWeightFactor(factor[which]);
339            }
340            // fake cutoff
341            if (c) {
342              double cutoff;
343              solver2->getDblParam(OsiDualObjectiveLimit,cutoff);
344              cutoff = CoinMin(cutoff,value + 0.1*fabs(value)*c);
345              heuristic4.setFakeCutoff(cutoff);
346            }
347            if (i||r) {
348              // also set increment
349              //double increment = (0.01*i+0.005)*(fabs(value)+1.0e-12);
350              double increment = 0.0;
351              heuristic4.setAbsoluteIncrement(increment);
352              heuristic4.setAccumulate(accumulate);
353              heuristic4.setMaximumRetries(r+1);
354            }
355            pumpTune = pumpTune%100;
356            if (pumpTune==6)
357              pumpTune =13;
358            heuristic4.setWhen(pumpTune+10);
359          }
360          heuristic4.setHeuristicName("feasibility pump");
361          model.addHeuristic(&heuristic4);
362        }
363        model.setCutoff(cutoff);
364        model.setMaximumNodes(numberNodes);
365        model.solver()->setHintParam(OsiDoReducePrint,true,OsiHintTry);
366        // Lightweight
367        CbcStrategyDefaultSubTree strategy(model_,true,5,1,0);
368        model.setStrategy(strategy);
369        model.solver()->setIntParam(OsiMaxNumIterationHotStart,10);
370        // Do search
371        if (logLevel>1)
372          model_->messageHandler()->message(CBC_START_SUB,model_->messages())
373            << name
374            << model.getMaximumNodes()
375            <<CoinMessageEol;
376        // probably faster to use a basis to get integer solutions
377        model.setSpecialOptions(2);
378#ifdef CBC_THREAD
379        if (model_->getNumberThreads()>0&&(model_->getThreadMode()&1)!=0) {
380          // See if at root node
381          bool atRoot = model_->getNodeCount()==0;
382          int passNumber = model_->getCurrentPassNumber();
383          if (atRoot&&passNumber==1)
384            model.setNumberThreads(model_->getNumberThreads());
385        }
386#endif
387        model.setMaximumCutPassesAtRoot(CoinMin(20,model_->getMaximumCutPassesAtRoot()));
388        model.setParentModel(*model_);
389        model.setOriginalColumns(process.originalColumns());
390        model.branchAndBound();
391        if (logLevel>1)
392          model_->messageHandler()->message(CBC_END_SUB,model_->messages())
393            << name
394            <<CoinMessageEol;
395        if (model.getMinimizationObjValue()<CoinMin(cutoff,1.0e30)) {
396          // solution
397          returnCode=model.isProvenOptimal() ? 3 : 1;
398          // post process
399#ifdef COIN_HAS_CLP
400          OsiClpSolverInterface * clpSolver = dynamic_cast< OsiClpSolverInterface*> (model.solver());
401          if (clpSolver) {
402            ClpSimplex * lpSolver = clpSolver->getModelPtr();
403            lpSolver->setSpecialOptions(lpSolver->specialOptions()|0x01000000); // say is Cbc (and in branch and bound)
404          }
405#endif
406          process.postProcess(*model.solver());
407          if (solver->isProvenOptimal()) {
408            // Solution now back in solver
409            int numberColumns = solver->getNumCols();
410            memcpy(newSolution,solver->getColSolution(),
411                   numberColumns*sizeof(double));
412            newSolutionValue = model.getMinimizationObjValue();
413          } else {
414            // odd - but no good
415            returnCode=0; // so will be infeasible
416          }
417        } else {
418        // no good
419          returnCode=model.isProvenInfeasible() ? 2 : 0; // so will be infeasible
420        }
421        if (model.status()==5)
422          returnCode=-2; // stop
423        if (model.isProvenInfeasible())
424          status=1;
425        else if (model.isProvenOptimal())
426          status=2;
427      }
428    }
429  } else {
430    returnCode=2; // infeasible finished
431  }
432  model_->setLogLevel(logLevel);
433  if (reset) {
434    for (int iColumn=0;iColumn<numberColumns;iColumn++) {
435      if (reset[iColumn])
436        solver->setColLower(iColumn,0.0);
437    }
438    delete [] reset;
439  }
440#ifdef COIN_DEVELOP
441  getHistoryStatistics_=true;
442  if (returnCode==1||returnCode==2) {
443    if (status==1)
444      printf("heuristic could add cut because infeasible (%s)\n",heuristicName_.c_str()); 
445    else if (status==2)
446      printf("heuristic could add cut because optimal (%s)\n",heuristicName_.c_str());
447  } 
448#endif
449  return returnCode;
450}
451
452// Default Constructor
453CbcRounding::CbcRounding() 
454  :CbcHeuristic()
455{
456  // matrix and row copy will automatically be empty
457  seed_=1;
458  down_ = NULL;
459  up_ = NULL;
460  equal_ = NULL;
461}
462
463// Constructor from model
464CbcRounding::CbcRounding(CbcModel & model)
465  :CbcHeuristic(model)
466{
467  // Get a copy of original matrix (and by row for rounding);
468  assert(model.solver());
469  matrix_ = *model.solver()->getMatrixByCol();
470  matrixByRow_ = *model.solver()->getMatrixByRow();
471  validate();
472  seed_=1;
473}
474
475// Destructor
476CbcRounding::~CbcRounding ()
477{
478  delete [] down_;
479  delete [] up_;
480  delete [] equal_;
481}
482
483// Clone
484CbcHeuristic *
485CbcRounding::clone() const
486{
487  return new CbcRounding(*this);
488}
489// Create C++ lines to get to current state
490void 
491CbcRounding::generateCpp( FILE * fp) 
492{
493  CbcRounding other;
494  fprintf(fp,"0#include \"CbcHeuristic.hpp\"\n");
495  fprintf(fp,"3  CbcRounding rounding(*cbcModel);\n");
496  CbcHeuristic::generateCpp(fp,"rounding");
497  if (seed_!=other.seed_)
498    fprintf(fp,"3  rounding.setSeed(%d);\n",seed_);
499  else
500    fprintf(fp,"4  rounding.setSeed(%d);\n",seed_);
501  fprintf(fp,"3  cbcModel->addHeuristic(&rounding);\n");
502}
503//#define NEW_ROUNDING
504// Copy constructor
505CbcRounding::CbcRounding(const CbcRounding & rhs)
506:
507  CbcHeuristic(rhs),
508  matrix_(rhs.matrix_),
509  matrixByRow_(rhs.matrixByRow_),
510  seed_(rhs.seed_)
511{
512#ifdef NEW_ROUNDING
513  int numberColumns = matrix_.getNumCols();
514  down_ = CoinCopyOfArray(rhs.down_,numberColumns);
515  up_ = CoinCopyOfArray(rhs.up_,numberColumns);
516  equal_ = CoinCopyOfArray(rhs.equal_,numberColumns);
517#else
518  down_ = NULL;
519  up_ = NULL;
520  equal_ = NULL;
521#endif 
522}
523
524// Assignment operator
525CbcRounding & 
526CbcRounding::operator=( const CbcRounding& rhs)
527{
528  if (this!=&rhs) {
529    CbcHeuristic::operator=(rhs);
530    matrix_ = rhs.matrix_;
531    matrixByRow_ = rhs.matrixByRow_;
532#ifdef NEW_ROUNDING
533    delete [] down_;
534    delete [] up_;
535    delete [] equal_;
536    int numberColumns = matrix_.getNumCols();
537    down_ = CoinCopyOfArray(rhs.down_,numberColumns);
538    up_ = CoinCopyOfArray(rhs.up_,numberColumns);
539    equal_ = CoinCopyOfArray(rhs.equal_,numberColumns);
540#else
541    down_ = NULL;
542    up_ = NULL;
543    equal_ = NULL;
544#endif 
545    seed_ = rhs.seed_;
546  }
547  return *this;
548}
549
550// Resets stuff if model changes
551void 
552CbcRounding::resetModel(CbcModel * model)
553{
554  model_=model;
555  // Get a copy of original matrix (and by row for rounding);
556  assert(model_->solver());
557  matrix_ = *model_->solver()->getMatrixByCol();
558  matrixByRow_ = *model_->solver()->getMatrixByRow();
559  validate();
560}
561// See if rounding will give solution
562// Sets value of solution
563// Assumes rhs for original matrix still okay
564// At present only works with integers
565// Fix values if asked for
566// Returns 1 if solution, 0 if not
567int
568CbcRounding::solution(double & solutionValue,
569                      double * betterSolution)
570{
571
572  // See if to do
573  if (!when()||(when()%10==1&&model_->phase()!=1)||
574      (when()%10==2&&(model_->phase()!=2&&model_->phase()!=3)))
575    return 0; // switched off
576  OsiSolverInterface * solver = model_->solver();
577  double direction = solver->getObjSense();
578  double newSolutionValue = direction*solver->getObjValue();
579  return solution(solutionValue,betterSolution,newSolutionValue);
580}
581// See if rounding will give solution
582// Sets value of solution
583// Assumes rhs for original matrix still okay
584// At present only works with integers
585// Fix values if asked for
586// Returns 1 if solution, 0 if not
587int
588CbcRounding::solution(double & solutionValue,
589                      double * betterSolution,
590                      double newSolutionValue)
591{
592
593  // See if to do
594  if (!when()||(when()%10==1&&model_->phase()!=1)||
595      (when()%10==2&&(model_->phase()!=2&&model_->phase()!=3)))
596    return 0; // switched off
597  OsiSolverInterface * solver = model_->solver();
598  const double * lower = solver->getColLower();
599  const double * upper = solver->getColUpper();
600  const double * rowLower = solver->getRowLower();
601  const double * rowUpper = solver->getRowUpper();
602  const double * solution = solver->getColSolution();
603  const double * objective = solver->getObjCoefficients();
604  double integerTolerance = model_->getDblParam(CbcModel::CbcIntegerTolerance);
605  double primalTolerance;
606  solver->getDblParam(OsiPrimalTolerance,primalTolerance);
607
608  int numberRows = matrix_.getNumRows();
609  assert (numberRows<=solver->getNumRows());
610  int numberIntegers = model_->numberIntegers();
611  const int * integerVariable = model_->integerVariable();
612  int i;
613  double direction = solver->getObjSense();
614  //double newSolutionValue = direction*solver->getObjValue();
615  int returnCode = 0;
616  // Column copy
617  const double * element = matrix_.getElements();
618  const int * row = matrix_.getIndices();
619  const CoinBigIndex * columnStart = matrix_.getVectorStarts();
620  const int * columnLength = matrix_.getVectorLengths();
621  // Row copy
622  const double * elementByRow = matrixByRow_.getElements();
623  const int * column = matrixByRow_.getIndices();
624  const CoinBigIndex * rowStart = matrixByRow_.getVectorStarts();
625  const int * rowLength = matrixByRow_.getVectorLengths();
626
627  // Get solution array for heuristic solution
628  int numberColumns = solver->getNumCols();
629  double * newSolution = new double [numberColumns];
630  memcpy(newSolution,solution,numberColumns*sizeof(double));
631
632  double * rowActivity = new double[numberRows];
633  memset(rowActivity,0,numberRows*sizeof(double));
634  for (i=0;i<numberColumns;i++) {
635    int j;
636    double value = newSolution[i];
637    if (value<lower[i]) {
638      value=lower[i];
639      newSolution[i]=value;
640    } else if (value>upper[i]) {
641      value=upper[i];
642      newSolution[i]=value;
643    }
644    if (value) {
645      for (j=columnStart[i];
646           j<columnStart[i]+columnLength[i];j++) {
647        int iRow=row[j];
648        rowActivity[iRow] += value*element[j];
649      }
650    }
651  }
652  // check was feasible - if not adjust (cleaning may move)
653  for (i=0;i<numberRows;i++) {
654    if(rowActivity[i]<rowLower[i]) {
655      //assert (rowActivity[i]>rowLower[i]-1000.0*primalTolerance);
656      rowActivity[i]=rowLower[i];
657    } else if(rowActivity[i]>rowUpper[i]) {
658      //assert (rowActivity[i]<rowUpper[i]+1000.0*primalTolerance);
659      rowActivity[i]=rowUpper[i];
660    }
661  }
662  for (i=0;i<numberIntegers;i++) {
663    int iColumn = integerVariable[i];
664    double value=newSolution[iColumn];
665    if (fabs(floor(value+0.5)-value)>integerTolerance) {
666      double below = floor(value);
667      double newValue=newSolution[iColumn];
668      double cost = direction * objective[iColumn];
669      double move;
670      if (cost>0.0) {
671        // try up
672        move = 1.0 -(value-below);
673      } else if (cost<0.0) {
674        // try down
675        move = below-value;
676      } else {
677        // won't be able to move unless we can grab another variable
678        double randomNumber = randomNumberGenerator_.randomDouble();
679        // which way?
680        if (randomNumber<0.5) 
681          move = below-value;
682        else
683          move = 1.0 -(value-below);
684      }
685      newValue += move;
686      newSolution[iColumn] = newValue;
687      newSolutionValue += move*cost;
688      int j;
689      for (j=columnStart[iColumn];
690           j<columnStart[iColumn]+columnLength[iColumn];j++) {
691        int iRow = row[j];
692        rowActivity[iRow] += move*element[j];
693      }
694    }
695  }
696
697  double penalty=0.0;
698  const char * integerType = model_->integerType();
699  // see if feasible - just using singletons
700  for (i=0;i<numberRows;i++) {
701    double value = rowActivity[i];
702    double thisInfeasibility=0.0;
703    if (value<rowLower[i]-primalTolerance)
704      thisInfeasibility = value-rowLower[i];
705    else if (value>rowUpper[i]+primalTolerance)
706      thisInfeasibility = value-rowUpper[i];
707    if (thisInfeasibility) {
708      // See if there are any slacks I can use to fix up
709      // maybe put in coding for multiple slacks?
710      double bestCost = 1.0e50;
711      int k;
712      int iBest=-1;
713      double addCost=0.0;
714      double newValue=0.0;
715      double changeRowActivity=0.0;
716      double absInfeasibility = fabs(thisInfeasibility);
717      for (k=rowStart[i];k<rowStart[i]+rowLength[i];k++) {
718        int iColumn = column[k];
719        // See if all elements help
720        if (columnLength[iColumn]==1) {
721          double currentValue = newSolution[iColumn];
722          double elementValue = elementByRow[k];
723          double lowerValue = lower[iColumn];
724          double upperValue = upper[iColumn];
725          double gap = rowUpper[i]-rowLower[i];
726          double absElement=fabs(elementValue);
727          if (thisInfeasibility*elementValue>0.0) {
728            // we want to reduce
729            if ((currentValue-lowerValue)*absElement>=absInfeasibility) {
730              // possible - check if integer
731              double distance = absInfeasibility/absElement;
732              double thisCost = -direction*objective[iColumn]*distance;
733              if (integerType[iColumn]) {
734                distance = ceil(distance-primalTolerance);
735                if (currentValue-distance>=lowerValue-primalTolerance) {
736                  if (absInfeasibility-distance*absElement< -gap-primalTolerance)
737                    thisCost=1.0e100; // no good
738                  else
739                    thisCost = -direction*objective[iColumn]*distance;
740                } else {
741                  thisCost=1.0e100; // no good
742                }
743              }
744              if (thisCost<bestCost) {
745                bestCost=thisCost;
746                iBest=iColumn;
747                addCost = thisCost;
748                newValue = currentValue-distance;
749                changeRowActivity = -distance*elementValue;
750              }
751            }
752          } else {
753            // we want to increase
754            if ((upperValue-currentValue)*absElement>=absInfeasibility) {
755              // possible - check if integer
756              double distance = absInfeasibility/absElement;
757              double thisCost = direction*objective[iColumn]*distance;
758              if (integerType[iColumn]) {
759                distance = ceil(distance-1.0e-7);
760                assert (currentValue-distance<=upperValue+primalTolerance);
761                if (absInfeasibility-distance*absElement< -gap-primalTolerance)
762                  thisCost=1.0e100; // no good
763                else
764                  thisCost = direction*objective[iColumn]*distance;
765              }
766              if (thisCost<bestCost) {
767                bestCost=thisCost;
768                iBest=iColumn;
769                addCost = thisCost;
770                newValue = currentValue+distance;
771                changeRowActivity = distance*elementValue;
772              }
773            }
774          }
775        }
776      }
777      if (iBest>=0) {
778        /*printf("Infeasibility of %g on row %d cost %g\n",
779          thisInfeasibility,i,addCost);*/
780        newSolution[iBest]=newValue;
781        thisInfeasibility=0.0;
782        newSolutionValue += addCost;
783        rowActivity[i] += changeRowActivity;
784      }
785      penalty += fabs(thisInfeasibility);
786    }
787  }
788  if (penalty) {
789    // see if feasible using any
790    // first continuous
791    double penaltyChange=0.0;
792    int iColumn;
793    for (iColumn=0;iColumn<numberColumns;iColumn++) {
794      if (integerType[iColumn])
795        continue;
796      double currentValue = newSolution[iColumn];
797      double lowerValue = lower[iColumn];
798      double upperValue = upper[iColumn];
799      int j;
800      int anyBadDown=0;
801      int anyBadUp=0;
802      double upImprovement=0.0;
803      double downImprovement=0.0;
804      for (j=columnStart[iColumn];
805           j<columnStart[iColumn]+columnLength[iColumn];j++) {
806        int iRow = row[j];
807        if (rowUpper[iRow]>rowLower[iRow]) {
808          double value = element[j];
809          if (rowActivity[iRow]>rowUpper[iRow]+primalTolerance) {
810            // infeasible above
811            downImprovement += value;
812            upImprovement -= value;
813            if (value>0.0) 
814              anyBadUp++;
815            else 
816              anyBadDown++;
817          } else if (rowActivity[iRow]>rowUpper[iRow]-primalTolerance) {
818            // feasible at ub
819            if (value>0.0) {
820              upImprovement -= value;
821              anyBadUp++;
822            } else {
823              downImprovement += value;
824              anyBadDown++;
825            }
826          } else if (rowActivity[iRow]>rowLower[iRow]+primalTolerance) {
827            // feasible in interior
828          } else if (rowActivity[iRow]>rowLower[iRow]-primalTolerance) {
829            // feasible at lb
830            if (value<0.0) {
831              upImprovement += value;
832              anyBadUp++;
833            } else {
834              downImprovement -= value;
835              anyBadDown++;
836            }
837          } else {
838            // infeasible below
839            downImprovement -= value;
840            upImprovement += value;
841            if (value<0.0) 
842              anyBadUp++;
843            else 
844              anyBadDown++;
845          }
846        } else {
847          // equality row
848          double value = element[j];
849          if (rowActivity[iRow]>rowUpper[iRow]+primalTolerance) {
850            // infeasible above
851            downImprovement += value;
852            upImprovement -= value;
853            if (value>0.0) 
854              anyBadUp++;
855            else 
856              anyBadDown++;
857          } else if (rowActivity[iRow]<rowLower[iRow]-primalTolerance) {
858            // infeasible below
859            downImprovement -= value;
860            upImprovement += value;
861            if (value<0.0) 
862              anyBadUp++;
863            else 
864              anyBadDown++;
865          } else {
866            // feasible - no good
867            anyBadUp=-1;
868            anyBadDown=-1;
869            break;
870          }
871        }
872      }
873      // could change tests for anyBad
874      if (anyBadUp)
875        upImprovement=0.0;
876      if (anyBadDown)
877        downImprovement=0.0;
878      double way=0.0;
879      double improvement=0.0;
880      if (downImprovement>0.0&&currentValue>lowerValue) {
881        way=-1.0;
882        improvement = downImprovement;
883      } else if (upImprovement>0.0&&currentValue<upperValue) {
884        way=1.0;
885        improvement = upImprovement;
886      }
887      if (way) {
888        // can improve
889        double distance;
890        if (way>0.0)
891          distance = upperValue-currentValue;
892        else
893          distance = currentValue-lowerValue;
894        for (j=columnStart[iColumn];
895             j<columnStart[iColumn]+columnLength[iColumn];j++) {
896          int iRow = row[j];
897          double value = element[j]*way;
898          if (rowActivity[iRow]>rowUpper[iRow]+primalTolerance) {
899            // infeasible above
900            assert (value<0.0);
901            double gap = rowActivity[iRow]-rowUpper[iRow];
902            if (gap+value*distance<0.0) 
903              distance = -gap/value;
904          } else if (rowActivity[iRow]<rowLower[iRow]-primalTolerance) {
905            // infeasible below
906            assert (value>0.0);
907            double gap = rowActivity[iRow]-rowLower[iRow];
908            if (gap+value*distance>0.0) 
909              distance = -gap/value;
910          } else {
911            // feasible
912            if (value>0) {
913              double gap = rowActivity[iRow]-rowUpper[iRow];
914              if (gap+value*distance>0.0) 
915              distance = -gap/value;
916            } else {
917              double gap = rowActivity[iRow]-rowLower[iRow];
918              if (gap+value*distance<0.0) 
919                distance = -gap/value;
920            }
921          }
922        }
923        //move
924        penaltyChange += improvement*distance;
925        distance *= way;
926        newSolution[iColumn] += distance;
927        newSolutionValue += direction*objective[iColumn]*distance;
928        for (j=columnStart[iColumn];
929             j<columnStart[iColumn]+columnLength[iColumn];j++) {
930          int iRow = row[j];
931          double value = element[j];
932          rowActivity[iRow] += distance*value;
933        }
934      }
935    }
936    // and now all if improving
937    double lastChange= penaltyChange ? 1.0 : 0.0;
938    while (lastChange>1.0e-2) {
939      lastChange=0;
940      for (iColumn=0;iColumn<numberColumns;iColumn++) {
941        bool isInteger = (integerType[iColumn]!=0);
942        double currentValue = newSolution[iColumn];
943        double lowerValue = lower[iColumn];
944        double upperValue = upper[iColumn];
945        int j;
946        int anyBadDown=0;
947        int anyBadUp=0;
948        double upImprovement=0.0;
949        double downImprovement=0.0;
950        for (j=columnStart[iColumn];
951             j<columnStart[iColumn]+columnLength[iColumn];j++) {
952          int iRow = row[j];
953          double value = element[j];
954          if (isInteger) {
955            if (value>0.0) {
956              if (rowActivity[iRow]+value>rowUpper[iRow]+primalTolerance)
957                anyBadUp++;
958              if (rowActivity[iRow]-value<rowLower[iRow]-primalTolerance)
959                anyBadDown++;
960            } else {
961              if (rowActivity[iRow]-value>rowUpper[iRow]+primalTolerance)
962                anyBadDown++;
963              if (rowActivity[iRow]+value<rowLower[iRow]-primalTolerance)
964                anyBadUp++;
965            }
966          }
967          if (rowUpper[iRow]>rowLower[iRow]) {
968            if (rowActivity[iRow]>rowUpper[iRow]+primalTolerance) {
969              // infeasible above
970              downImprovement += value;
971              upImprovement -= value;
972              if (value>0.0) 
973                anyBadUp++;
974              else 
975                anyBadDown++;
976            } else if (rowActivity[iRow]>rowUpper[iRow]-primalTolerance) {
977              // feasible at ub
978              if (value>0.0) {
979                upImprovement -= value;
980                anyBadUp++;
981              } else {
982                downImprovement += value;
983                anyBadDown++;
984              }
985            } else if (rowActivity[iRow]>rowLower[iRow]+primalTolerance) {
986              // feasible in interior
987            } else if (rowActivity[iRow]>rowLower[iRow]-primalTolerance) {
988              // feasible at lb
989              if (value<0.0) {
990                upImprovement += value;
991                anyBadUp++;
992              } else {
993                downImprovement -= value;
994                anyBadDown++;
995              }
996            } else {
997              // infeasible below
998              downImprovement -= value;
999              upImprovement += value;
1000              if (value<0.0) 
1001                anyBadUp++;
1002              else 
1003                anyBadDown++;
1004            }
1005          } else {
1006            // equality row
1007            if (rowActivity[iRow]>rowUpper[iRow]+primalTolerance) {
1008              // infeasible above
1009              downImprovement += value;
1010              upImprovement -= value;
1011              if (value>0.0) 
1012                anyBadUp++;
1013              else 
1014                anyBadDown++;
1015            } else if (rowActivity[iRow]<rowLower[iRow]-primalTolerance) {
1016              // infeasible below
1017              downImprovement -= value;
1018              upImprovement += value;
1019              if (value<0.0) 
1020                anyBadUp++;
1021              else 
1022                anyBadDown++;
1023            } else {
1024              // feasible - no good
1025              anyBadUp=-1;
1026              anyBadDown=-1;
1027              break;
1028            }
1029          }
1030        }
1031        // could change tests for anyBad
1032        if (anyBadUp)
1033          upImprovement=0.0;
1034        if (anyBadDown)
1035          downImprovement=0.0;
1036        double way=0.0;
1037        double improvement=0.0;
1038        if (downImprovement>0.0&&currentValue>lowerValue) {
1039          way=-1.0;
1040          improvement = downImprovement;
1041        } else if (upImprovement>0.0&&currentValue<upperValue) {
1042          way=1.0;
1043          improvement = upImprovement;
1044        }
1045        if (way) {
1046          // can improve
1047          double distance=COIN_DBL_MAX;
1048          for (j=columnStart[iColumn];
1049               j<columnStart[iColumn]+columnLength[iColumn];j++) {
1050            int iRow = row[j];
1051            double value = element[j]*way;
1052            if (rowActivity[iRow]>rowUpper[iRow]+primalTolerance) {
1053              // infeasible above
1054              assert (value<0.0);
1055              double gap = rowActivity[iRow]-rowUpper[iRow];
1056              if (gap+value*distance<0.0) {
1057                // If integer then has to move by 1
1058                if (!isInteger)
1059                  distance = -gap/value;
1060                else
1061                  distance = CoinMax(-gap/value,1.0);
1062              }
1063            } else if (rowActivity[iRow]<rowLower[iRow]-primalTolerance) {
1064              // infeasible below
1065              assert (value>0.0);
1066              double gap = rowActivity[iRow]-rowLower[iRow];
1067              if (gap+value*distance>0.0) {
1068                // If integer then has to move by 1
1069                if (!isInteger)
1070                  distance = -gap/value;
1071                else
1072                  distance = CoinMax(-gap/value,1.0);
1073              }
1074            } else {
1075              // feasible
1076              if (value>0) {
1077                double gap = rowActivity[iRow]-rowUpper[iRow];
1078                if (gap+value*distance>0.0) 
1079                  distance = -gap/value;
1080              } else {
1081                double gap = rowActivity[iRow]-rowLower[iRow];
1082                if (gap+value*distance<0.0) 
1083                  distance = -gap/value;
1084              }
1085            }
1086          }
1087          if (isInteger)
1088            distance = floor(distance+1.05e-8);
1089          if (!distance) {
1090            // should never happen
1091            //printf("zero distance in CbcRounding - debug\n");
1092          }
1093          //move
1094          lastChange += improvement*distance;
1095          distance *= way;
1096          newSolution[iColumn] += distance;
1097          newSolutionValue += direction*objective[iColumn]*distance;
1098          for (j=columnStart[iColumn];
1099               j<columnStart[iColumn]+columnLength[iColumn];j++) {
1100            int iRow = row[j];
1101            double value = element[j];
1102            rowActivity[iRow] += distance*value;
1103          }
1104        }
1105      }
1106      penaltyChange += lastChange;
1107    }
1108    penalty -= penaltyChange;
1109    if (penalty<1.0e-5*fabs(penaltyChange)) {
1110      // recompute
1111      penalty=0.0;
1112      for (i=0;i<numberRows;i++) {
1113        double value = rowActivity[i];
1114        if (value<rowLower[i]-primalTolerance)
1115          penalty += rowLower[i]-value;
1116        else if (value>rowUpper[i]+primalTolerance)
1117          penalty += value-rowUpper[i];
1118      }
1119    }
1120  }
1121
1122  // Could also set SOS (using random) and repeat
1123  if (!penalty) {
1124    // See if we can do better
1125    //seed_++;
1126    //CoinSeedRandom(seed_);
1127    // Random number between 0 and 1.
1128    double randomNumber = randomNumberGenerator_.randomDouble();
1129    int iPass;
1130    int start[2];
1131    int end[2];
1132    int iRandom = (int) (randomNumber*((double) numberIntegers));
1133    start[0]=iRandom;
1134    end[0]=numberIntegers;
1135    start[1]=0;
1136    end[1]=iRandom;
1137    for (iPass=0;iPass<2;iPass++) {
1138      int i;
1139      for (i=start[iPass];i<end[iPass];i++) {
1140        int iColumn = integerVariable[i];
1141#ifndef NDEBUG
1142        double value=newSolution[iColumn];
1143        assert (fabs(floor(value+0.5)-value)<integerTolerance);
1144#endif
1145        double cost = direction * objective[iColumn];
1146        double move=0.0;
1147        if (cost>0.0)
1148          move = -1.0;
1149        else if (cost<0.0)
1150          move=1.0;
1151        while (move) {
1152          bool good=true;
1153          double newValue=newSolution[iColumn]+move;
1154          if (newValue<lower[iColumn]-primalTolerance||
1155              newValue>upper[iColumn]+primalTolerance) {
1156            move=0.0;
1157          } else {
1158            // see if we can move
1159            int j;
1160            for (j=columnStart[iColumn];
1161                 j<columnStart[iColumn]+columnLength[iColumn];j++) {
1162              int iRow = row[j];
1163              double newActivity = rowActivity[iRow] + move*element[j];
1164              if (newActivity<rowLower[iRow]-primalTolerance||
1165                  newActivity>rowUpper[iRow]+primalTolerance) {
1166                good=false;
1167                break;
1168              }
1169            }
1170            if (good) {
1171              newSolution[iColumn] = newValue;
1172              newSolutionValue += move*cost;
1173              int j;
1174              for (j=columnStart[iColumn];
1175                   j<columnStart[iColumn]+columnLength[iColumn];j++) {
1176                int iRow = row[j];
1177                rowActivity[iRow] += move*element[j];
1178              }
1179            } else {
1180              move=0.0;
1181            }
1182          }
1183        }
1184      }
1185    }
1186    // Just in case of some stupidity
1187    double objOffset=0.0;
1188    solver->getDblParam(OsiObjOffset,objOffset);
1189    newSolutionValue = -objOffset;
1190    for ( i=0 ; i<numberColumns ; i++ )
1191      newSolutionValue += objective[i]*newSolution[i];
1192    newSolutionValue *= direction;
1193    //printf("new solution value %g %g\n",newSolutionValue,solutionValue);
1194    if (newSolutionValue<solutionValue) {
1195      // paranoid check
1196      memset(rowActivity,0,numberRows*sizeof(double));
1197      for (i=0;i<numberColumns;i++) {
1198        int j;
1199        double value = newSolution[i];
1200        if (value) {
1201          for (j=columnStart[i];
1202               j<columnStart[i]+columnLength[i];j++) {
1203            int iRow=row[j];
1204            rowActivity[iRow] += value*element[j];
1205          }
1206        }
1207      }
1208      // check was approximately feasible
1209      bool feasible=true;
1210      for (i=0;i<numberRows;i++) {
1211        if(rowActivity[i]<rowLower[i]) {
1212          if (rowActivity[i]<rowLower[i]-1000.0*primalTolerance)
1213            feasible = false;
1214        } else if(rowActivity[i]>rowUpper[i]) {
1215          if (rowActivity[i]>rowUpper[i]+1000.0*primalTolerance)
1216            feasible = false;
1217        }
1218      }
1219      if (feasible) {
1220        // new solution
1221        memcpy(betterSolution,newSolution,numberColumns*sizeof(double));
1222        solutionValue = newSolutionValue;
1223        //printf("** Solution of %g found by rounding\n",newSolutionValue);
1224        returnCode=1;
1225      } else {
1226        // Can easily happen
1227        //printf("Debug CbcRounding giving bad solution\n");
1228      }
1229    }
1230  }
1231#ifdef NEW_ROUNDING
1232  if (!returnCode) {
1233#if 0
1234    // back to starting point
1235    memcpy(newSolution,solution,numberColumns*sizeof(double));
1236    memset(rowActivity,0,numberRows*sizeof(double));
1237    for (i=0;i<numberColumns;i++) {
1238      int j;
1239      double value = newSolution[i];
1240      if (value<lower[i]) {
1241        value=lower[i];
1242        newSolution[i]=value;
1243      } else if (value>upper[i]) {
1244        value=upper[i];
1245        newSolution[i]=value;
1246      }
1247      if (value) {
1248        for (j=columnStart[i];
1249             j<columnStart[i]+columnLength[i];j++) {
1250          int iRow=row[j];
1251          rowActivity[iRow] += value*element[j];
1252        }
1253      }
1254    }
1255    // check was feasible - if not adjust (cleaning may move)
1256    for (i=0;i<numberRows;i++) {
1257      if(rowActivity[i]<rowLower[i]) {
1258        //assert (rowActivity[i]>rowLower[i]-1000.0*primalTolerance);
1259        rowActivity[i]=rowLower[i];
1260      } else if(rowActivity[i]>rowUpper[i]) {
1261        //assert (rowActivity[i]<rowUpper[i]+1000.0*primalTolerance);
1262        rowActivity[i]=rowUpper[i];
1263      }
1264    }
1265#endif
1266    int * candidate = new int [numberColumns];
1267    int nCandidate=0;
1268    for (iColumn=0;iColumn<numberColumns;iColumn++) {
1269      bool isInteger = (integerType[iColumn]!=0);
1270      if (isInteger) {
1271        double currentValue = newSolution[iColumn];
1272        if (fabs(currentValue-floor(currentValue+0.5))>1.0e-8)
1273          candidate[nCandidate++]=iColumn;
1274      }
1275    }
1276    if (true) {
1277      // Rounding as in Berthold
1278      while (nCandidate) {
1279        double infeasibility =1.0e-7;
1280        int iRow=-1;
1281        for (i=0;i<numberRows;i++) {
1282          double value=0.0;
1283          if(rowActivity[i]<rowLower[i]) {
1284            value = rowLower[i]-rowActivity[i];
1285          } else if(rowActivity[i]>rowUpper[i]) {
1286            value = rowActivity[i]-rowUpper[i];
1287          }
1288          if (value>infeasibility) {
1289            infeasibility = value;
1290            iRow=i;
1291          }
1292        }
1293        if (iRow>=0) {
1294          // infeasible
1295        } else {
1296          // feasible
1297        }
1298      }
1299    } else {
1300      // Shifting as in Berthold
1301    }
1302    delete [] candidate;
1303  }
1304#endif
1305  delete [] newSolution;
1306  delete [] rowActivity;
1307  return returnCode;
1308}
1309// update model
1310void CbcRounding::setModel(CbcModel * model)
1311{
1312  model_ = model;
1313  // Get a copy of original matrix (and by row for rounding);
1314  assert(model_->solver());
1315  matrix_ = *model_->solver()->getMatrixByCol();
1316  matrixByRow_ = *model_->solver()->getMatrixByRow();
1317  // make sure model okay for heuristic
1318  validate();
1319}
1320// Validate model i.e. sets when_ to 0 if necessary (may be NULL)
1321void 
1322CbcRounding::validate() 
1323{
1324  if (model_&&when()<10) {
1325    if (model_->numberIntegers()!=
1326        model_->numberObjects())
1327      setWhen(0);
1328  }
1329#ifdef NEW_ROUNDING
1330  int numberColumns = matrix_.getNumCols();
1331  down_ = new unsigned short [numberColumns];
1332  up_ = new unsigned short [numberColumns];
1333  equal_ = new unsigned short [numberColumns];
1334  // Column copy
1335  const double * element = matrix_.getElements();
1336  const int * row = matrix_.getIndices();
1337  const CoinBigIndex * columnStart = matrix_.getVectorStarts();
1338  const int * columnLength = matrix_.getVectorLengths();
1339  const double * rowLower = model.solver()->getRowLower();
1340  const double * rowUpper = model.solver()->getRowUpper();
1341  for (int i=0;i<numberColumns;i++) {
1342    int down=0;
1343    int up=0;
1344    int equal=0;
1345    if (columnLength[i]>65535) {
1346      equal[0]=65535; 
1347      break; // unlikely to work
1348    }
1349    for (CoinBigIndex j=columnStart[i];
1350         j<columnStart[i]+columnLength[i];j++) {
1351      int iRow=row[j];
1352      if (rowLower[iRow]>-1.0e20&&rowUpper[iRow]<1.0e20) {
1353        equal++;
1354      } else if (element[j]>0.0) {
1355        if (rowUpper[iRow]<1.0e20)
1356          up++;
1357        else
1358          down--;
1359      } else {
1360        if (rowLower[iRow]>-1.0e20)
1361          up++;
1362        else
1363          down--;
1364      }
1365    }
1366    down_[i] = (unsigned short) down;
1367    up_[i] = (unsigned short) up;
1368    equal_[i] = (unsigned short) equal;
1369  }
1370#else
1371  down_ = NULL;
1372  up_ = NULL;
1373  equal_ = NULL;
1374#endif 
1375}
1376
1377// Default Constructor
1378CbcSerendipity::CbcSerendipity() 
1379  :CbcHeuristic()
1380{
1381}
1382
1383// Constructor from model
1384CbcSerendipity::CbcSerendipity(CbcModel & model)
1385  :CbcHeuristic(model)
1386{
1387}
1388
1389// Destructor
1390CbcSerendipity::~CbcSerendipity ()
1391{
1392}
1393
1394// Clone
1395CbcHeuristic *
1396CbcSerendipity::clone() const
1397{
1398  return new CbcSerendipity(*this);
1399}
1400// Create C++ lines to get to current state
1401void 
1402CbcSerendipity::generateCpp( FILE * fp) 
1403{
1404  fprintf(fp,"0#include \"CbcHeuristic.hpp\"\n");
1405  fprintf(fp,"3  CbcSerendipity serendipity(*cbcModel);\n");
1406  CbcHeuristic::generateCpp(fp,"serendipity");
1407  fprintf(fp,"3  cbcModel->addHeuristic(&serendipity);\n");
1408}
1409
1410// Copy constructor
1411CbcSerendipity::CbcSerendipity(const CbcSerendipity & rhs)
1412:
1413  CbcHeuristic(rhs)
1414{
1415}
1416
1417// Assignment operator
1418CbcSerendipity & 
1419CbcSerendipity::operator=( const CbcSerendipity& rhs)
1420{
1421  if (this!=&rhs) {
1422    CbcHeuristic::operator=(rhs);
1423  }
1424  return *this;
1425}
1426
1427// Returns 1 if solution, 0 if not
1428int
1429CbcSerendipity::solution(double & solutionValue,
1430                         double * betterSolution)
1431{
1432  if (!model_)
1433    return 0;
1434  // get information on solver type
1435  OsiAuxInfo * auxInfo = model_->solver()->getAuxiliaryInfo();
1436  OsiBabSolver * auxiliaryInfo = dynamic_cast< OsiBabSolver *> (auxInfo);
1437  if (auxiliaryInfo)
1438    return auxiliaryInfo->solution(solutionValue,betterSolution,model_->solver()->getNumCols());
1439  else
1440    return 0;
1441}
1442// update model
1443void CbcSerendipity::setModel(CbcModel * model)
1444{
1445  model_ = model;
1446}
1447// Resets stuff if model changes
1448void 
1449CbcSerendipity::resetModel(CbcModel * model)
1450{
1451  model_ = model;
1452}
1453 
Note: See TracBrowser for help on using the repository browser.