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

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

try and make a bit faster
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 47.7 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        if (model.getNumCols()) {
391          setCutAndHeuristicOptions(model);
392          model.branchAndBound();
393        } else {
394          // empty model
395          model.setMinimizationObjValue(model.solver()->getObjSense()*model.solver()->getObjValue());
396        }
397        if (logLevel>1)
398          model_->messageHandler()->message(CBC_END_SUB,model_->messages())
399            << name
400            <<CoinMessageEol;
401        if (model.getMinimizationObjValue()<CoinMin(cutoff,1.0e30)) {
402          // solution
403          if (model.getNumCols())
404            returnCode=model.isProvenOptimal() ? 3 : 1;
405          else
406            returnCode=3;
407          // post process
408#ifdef COIN_HAS_CLP
409          OsiClpSolverInterface * clpSolver = dynamic_cast< OsiClpSolverInterface*> (model.solver());
410          if (clpSolver) {
411            ClpSimplex * lpSolver = clpSolver->getModelPtr();
412            lpSolver->setSpecialOptions(lpSolver->specialOptions()|0x01000000); // say is Cbc (and in branch and bound)
413          }
414#endif
415          process.postProcess(*model.solver());
416          if (solver->isProvenOptimal()) {
417            // Solution now back in solver
418            int numberColumns = solver->getNumCols();
419            memcpy(newSolution,solver->getColSolution(),
420                   numberColumns*sizeof(double));
421            newSolutionValue = model.getMinimizationObjValue();
422          } else {
423            // odd - but no good
424            returnCode=0; // so will be infeasible
425          }
426        } else {
427        // no good
428          returnCode=model.isProvenInfeasible() ? 2 : 0; // so will be infeasible
429        }
430        if (model.status()==5)
431          returnCode=-2; // stop
432        if (model.isProvenInfeasible())
433          status=1;
434        else if (model.isProvenOptimal())
435          status=2;
436      }
437    }
438  } else {
439    returnCode=2; // infeasible finished
440  }
441  model_->setLogLevel(logLevel);
442  if (reset) {
443    for (int iColumn=0;iColumn<numberColumns;iColumn++) {
444      if (reset[iColumn])
445        solver->setColLower(iColumn,0.0);
446    }
447    delete [] reset;
448  }
449#ifdef COIN_DEVELOP
450  getHistoryStatistics_=true;
451  if (returnCode==1||returnCode==2) {
452    if (status==1)
453      printf("heuristic could add cut because infeasible (%s)\n",heuristicName_.c_str()); 
454    else if (status==2)
455      printf("heuristic could add cut because optimal (%s)\n",heuristicName_.c_str());
456  } 
457#endif
458  return returnCode;
459}
460
461// Default Constructor
462CbcRounding::CbcRounding() 
463  :CbcHeuristic()
464{
465  // matrix and row copy will automatically be empty
466  seed_=1;
467  down_ = NULL;
468  up_ = NULL;
469  equal_ = NULL;
470}
471
472// Constructor from model
473CbcRounding::CbcRounding(CbcModel & model)
474  :CbcHeuristic(model)
475{
476  // Get a copy of original matrix (and by row for rounding);
477  assert(model.solver());
478  matrix_ = *model.solver()->getMatrixByCol();
479  matrixByRow_ = *model.solver()->getMatrixByRow();
480  validate();
481  seed_=1;
482}
483
484// Destructor
485CbcRounding::~CbcRounding ()
486{
487  delete [] down_;
488  delete [] up_;
489  delete [] equal_;
490}
491
492// Clone
493CbcHeuristic *
494CbcRounding::clone() const
495{
496  return new CbcRounding(*this);
497}
498// Create C++ lines to get to current state
499void 
500CbcRounding::generateCpp( FILE * fp) 
501{
502  CbcRounding other;
503  fprintf(fp,"0#include \"CbcHeuristic.hpp\"\n");
504  fprintf(fp,"3  CbcRounding rounding(*cbcModel);\n");
505  CbcHeuristic::generateCpp(fp,"rounding");
506  if (seed_!=other.seed_)
507    fprintf(fp,"3  rounding.setSeed(%d);\n",seed_);
508  else
509    fprintf(fp,"4  rounding.setSeed(%d);\n",seed_);
510  fprintf(fp,"3  cbcModel->addHeuristic(&rounding);\n");
511}
512//#define NEW_ROUNDING
513// Copy constructor
514CbcRounding::CbcRounding(const CbcRounding & rhs)
515:
516  CbcHeuristic(rhs),
517  matrix_(rhs.matrix_),
518  matrixByRow_(rhs.matrixByRow_),
519  seed_(rhs.seed_)
520{
521#ifdef NEW_ROUNDING
522  int numberColumns = matrix_.getNumCols();
523  down_ = CoinCopyOfArray(rhs.down_,numberColumns);
524  up_ = CoinCopyOfArray(rhs.up_,numberColumns);
525  equal_ = CoinCopyOfArray(rhs.equal_,numberColumns);
526#else
527  down_ = NULL;
528  up_ = NULL;
529  equal_ = NULL;
530#endif 
531}
532
533// Assignment operator
534CbcRounding & 
535CbcRounding::operator=( const CbcRounding& rhs)
536{
537  if (this!=&rhs) {
538    CbcHeuristic::operator=(rhs);
539    matrix_ = rhs.matrix_;
540    matrixByRow_ = rhs.matrixByRow_;
541#ifdef NEW_ROUNDING
542    delete [] down_;
543    delete [] up_;
544    delete [] equal_;
545    int numberColumns = matrix_.getNumCols();
546    down_ = CoinCopyOfArray(rhs.down_,numberColumns);
547    up_ = CoinCopyOfArray(rhs.up_,numberColumns);
548    equal_ = CoinCopyOfArray(rhs.equal_,numberColumns);
549#else
550    down_ = NULL;
551    up_ = NULL;
552    equal_ = NULL;
553#endif 
554    seed_ = rhs.seed_;
555  }
556  return *this;
557}
558
559// Resets stuff if model changes
560void 
561CbcRounding::resetModel(CbcModel * model)
562{
563  model_=model;
564  // Get a copy of original matrix (and by row for rounding);
565  assert(model_->solver());
566  matrix_ = *model_->solver()->getMatrixByCol();
567  matrixByRow_ = *model_->solver()->getMatrixByRow();
568  validate();
569}
570// See if rounding will give solution
571// Sets value of solution
572// Assumes rhs for original matrix still okay
573// At present only works with integers
574// Fix values if asked for
575// Returns 1 if solution, 0 if not
576int
577CbcRounding::solution(double & solutionValue,
578                      double * betterSolution)
579{
580
581  // See if to do
582  if (!when()||(when()%10==1&&model_->phase()!=1)||
583      (when()%10==2&&(model_->phase()!=2&&model_->phase()!=3)))
584    return 0; // switched off
585  OsiSolverInterface * solver = model_->solver();
586  double direction = solver->getObjSense();
587  double newSolutionValue = direction*solver->getObjValue();
588  return solution(solutionValue,betterSolution,newSolutionValue);
589}
590// See if rounding will give solution
591// Sets value of solution
592// Assumes rhs for original matrix still okay
593// At present only works with integers
594// Fix values if asked for
595// Returns 1 if solution, 0 if not
596int
597CbcRounding::solution(double & solutionValue,
598                      double * betterSolution,
599                      double newSolutionValue)
600{
601
602  // See if to do
603  if (!when()||(when()%10==1&&model_->phase()!=1)||
604      (when()%10==2&&(model_->phase()!=2&&model_->phase()!=3)))
605    return 0; // switched off
606  OsiSolverInterface * solver = model_->solver();
607  const double * lower = solver->getColLower();
608  const double * upper = solver->getColUpper();
609  const double * rowLower = solver->getRowLower();
610  const double * rowUpper = solver->getRowUpper();
611  const double * solution = solver->getColSolution();
612  const double * objective = solver->getObjCoefficients();
613  double integerTolerance = model_->getDblParam(CbcModel::CbcIntegerTolerance);
614  double primalTolerance;
615  solver->getDblParam(OsiPrimalTolerance,primalTolerance);
616
617  int numberRows = matrix_.getNumRows();
618  assert (numberRows<=solver->getNumRows());
619  int numberIntegers = model_->numberIntegers();
620  const int * integerVariable = model_->integerVariable();
621  int i;
622  double direction = solver->getObjSense();
623  //double newSolutionValue = direction*solver->getObjValue();
624  int returnCode = 0;
625  // Column copy
626  const double * element = matrix_.getElements();
627  const int * row = matrix_.getIndices();
628  const CoinBigIndex * columnStart = matrix_.getVectorStarts();
629  const int * columnLength = matrix_.getVectorLengths();
630  // Row copy
631  const double * elementByRow = matrixByRow_.getElements();
632  const int * column = matrixByRow_.getIndices();
633  const CoinBigIndex * rowStart = matrixByRow_.getVectorStarts();
634  const int * rowLength = matrixByRow_.getVectorLengths();
635
636  // Get solution array for heuristic solution
637  int numberColumns = solver->getNumCols();
638  double * newSolution = new double [numberColumns];
639  memcpy(newSolution,solution,numberColumns*sizeof(double));
640
641  double * rowActivity = new double[numberRows];
642  memset(rowActivity,0,numberRows*sizeof(double));
643  for (i=0;i<numberColumns;i++) {
644    int j;
645    double value = newSolution[i];
646    if (value<lower[i]) {
647      value=lower[i];
648      newSolution[i]=value;
649    } else if (value>upper[i]) {
650      value=upper[i];
651      newSolution[i]=value;
652    }
653    if (value) {
654      for (j=columnStart[i];
655           j<columnStart[i]+columnLength[i];j++) {
656        int iRow=row[j];
657        rowActivity[iRow] += value*element[j];
658      }
659    }
660  }
661  // check was feasible - if not adjust (cleaning may move)
662  for (i=0;i<numberRows;i++) {
663    if(rowActivity[i]<rowLower[i]) {
664      //assert (rowActivity[i]>rowLower[i]-1000.0*primalTolerance);
665      rowActivity[i]=rowLower[i];
666    } else if(rowActivity[i]>rowUpper[i]) {
667      //assert (rowActivity[i]<rowUpper[i]+1000.0*primalTolerance);
668      rowActivity[i]=rowUpper[i];
669    }
670  }
671  for (i=0;i<numberIntegers;i++) {
672    int iColumn = integerVariable[i];
673    double value=newSolution[iColumn];
674    if (fabs(floor(value+0.5)-value)>integerTolerance) {
675      double below = floor(value);
676      double newValue=newSolution[iColumn];
677      double cost = direction * objective[iColumn];
678      double move;
679      if (cost>0.0) {
680        // try up
681        move = 1.0 -(value-below);
682      } else if (cost<0.0) {
683        // try down
684        move = below-value;
685      } else {
686        // won't be able to move unless we can grab another variable
687        double randomNumber = randomNumberGenerator_.randomDouble();
688        // which way?
689        if (randomNumber<0.5) 
690          move = below-value;
691        else
692          move = 1.0 -(value-below);
693      }
694      newValue += move;
695      newSolution[iColumn] = newValue;
696      newSolutionValue += move*cost;
697      int j;
698      for (j=columnStart[iColumn];
699           j<columnStart[iColumn]+columnLength[iColumn];j++) {
700        int iRow = row[j];
701        rowActivity[iRow] += move*element[j];
702      }
703    }
704  }
705
706  double penalty=0.0;
707  const char * integerType = model_->integerType();
708  // see if feasible - just using singletons
709  for (i=0;i<numberRows;i++) {
710    double value = rowActivity[i];
711    double thisInfeasibility=0.0;
712    if (value<rowLower[i]-primalTolerance)
713      thisInfeasibility = value-rowLower[i];
714    else if (value>rowUpper[i]+primalTolerance)
715      thisInfeasibility = value-rowUpper[i];
716    if (thisInfeasibility) {
717      // See if there are any slacks I can use to fix up
718      // maybe put in coding for multiple slacks?
719      double bestCost = 1.0e50;
720      int k;
721      int iBest=-1;
722      double addCost=0.0;
723      double newValue=0.0;
724      double changeRowActivity=0.0;
725      double absInfeasibility = fabs(thisInfeasibility);
726      for (k=rowStart[i];k<rowStart[i]+rowLength[i];k++) {
727        int iColumn = column[k];
728        // See if all elements help
729        if (columnLength[iColumn]==1) {
730          double currentValue = newSolution[iColumn];
731          double elementValue = elementByRow[k];
732          double lowerValue = lower[iColumn];
733          double upperValue = upper[iColumn];
734          double gap = rowUpper[i]-rowLower[i];
735          double absElement=fabs(elementValue);
736          if (thisInfeasibility*elementValue>0.0) {
737            // we want to reduce
738            if ((currentValue-lowerValue)*absElement>=absInfeasibility) {
739              // possible - check if integer
740              double distance = absInfeasibility/absElement;
741              double thisCost = -direction*objective[iColumn]*distance;
742              if (integerType[iColumn]) {
743                distance = ceil(distance-primalTolerance);
744                if (currentValue-distance>=lowerValue-primalTolerance) {
745                  if (absInfeasibility-distance*absElement< -gap-primalTolerance)
746                    thisCost=1.0e100; // no good
747                  else
748                    thisCost = -direction*objective[iColumn]*distance;
749                } else {
750                  thisCost=1.0e100; // no good
751                }
752              }
753              if (thisCost<bestCost) {
754                bestCost=thisCost;
755                iBest=iColumn;
756                addCost = thisCost;
757                newValue = currentValue-distance;
758                changeRowActivity = -distance*elementValue;
759              }
760            }
761          } else {
762            // we want to increase
763            if ((upperValue-currentValue)*absElement>=absInfeasibility) {
764              // possible - check if integer
765              double distance = absInfeasibility/absElement;
766              double thisCost = direction*objective[iColumn]*distance;
767              if (integerType[iColumn]) {
768                distance = ceil(distance-1.0e-7);
769                assert (currentValue-distance<=upperValue+primalTolerance);
770                if (absInfeasibility-distance*absElement< -gap-primalTolerance)
771                  thisCost=1.0e100; // no good
772                else
773                  thisCost = direction*objective[iColumn]*distance;
774              }
775              if (thisCost<bestCost) {
776                bestCost=thisCost;
777                iBest=iColumn;
778                addCost = thisCost;
779                newValue = currentValue+distance;
780                changeRowActivity = distance*elementValue;
781              }
782            }
783          }
784        }
785      }
786      if (iBest>=0) {
787        /*printf("Infeasibility of %g on row %d cost %g\n",
788          thisInfeasibility,i,addCost);*/
789        newSolution[iBest]=newValue;
790        thisInfeasibility=0.0;
791        newSolutionValue += addCost;
792        rowActivity[i] += changeRowActivity;
793      }
794      penalty += fabs(thisInfeasibility);
795    }
796  }
797  if (penalty) {
798    // see if feasible using any
799    // first continuous
800    double penaltyChange=0.0;
801    int iColumn;
802    for (iColumn=0;iColumn<numberColumns;iColumn++) {
803      if (integerType[iColumn])
804        continue;
805      double currentValue = newSolution[iColumn];
806      double lowerValue = lower[iColumn];
807      double upperValue = upper[iColumn];
808      int j;
809      int anyBadDown=0;
810      int anyBadUp=0;
811      double upImprovement=0.0;
812      double downImprovement=0.0;
813      for (j=columnStart[iColumn];
814           j<columnStart[iColumn]+columnLength[iColumn];j++) {
815        int iRow = row[j];
816        if (rowUpper[iRow]>rowLower[iRow]) {
817          double value = element[j];
818          if (rowActivity[iRow]>rowUpper[iRow]+primalTolerance) {
819            // infeasible above
820            downImprovement += value;
821            upImprovement -= value;
822            if (value>0.0) 
823              anyBadUp++;
824            else 
825              anyBadDown++;
826          } else if (rowActivity[iRow]>rowUpper[iRow]-primalTolerance) {
827            // feasible at ub
828            if (value>0.0) {
829              upImprovement -= value;
830              anyBadUp++;
831            } else {
832              downImprovement += value;
833              anyBadDown++;
834            }
835          } else if (rowActivity[iRow]>rowLower[iRow]+primalTolerance) {
836            // feasible in interior
837          } else if (rowActivity[iRow]>rowLower[iRow]-primalTolerance) {
838            // feasible at lb
839            if (value<0.0) {
840              upImprovement += value;
841              anyBadUp++;
842            } else {
843              downImprovement -= value;
844              anyBadDown++;
845            }
846          } else {
847            // infeasible below
848            downImprovement -= value;
849            upImprovement += value;
850            if (value<0.0) 
851              anyBadUp++;
852            else 
853              anyBadDown++;
854          }
855        } else {
856          // equality row
857          double value = element[j];
858          if (rowActivity[iRow]>rowUpper[iRow]+primalTolerance) {
859            // infeasible above
860            downImprovement += value;
861            upImprovement -= value;
862            if (value>0.0) 
863              anyBadUp++;
864            else 
865              anyBadDown++;
866          } else if (rowActivity[iRow]<rowLower[iRow]-primalTolerance) {
867            // infeasible below
868            downImprovement -= value;
869            upImprovement += value;
870            if (value<0.0) 
871              anyBadUp++;
872            else 
873              anyBadDown++;
874          } else {
875            // feasible - no good
876            anyBadUp=-1;
877            anyBadDown=-1;
878            break;
879          }
880        }
881      }
882      // could change tests for anyBad
883      if (anyBadUp)
884        upImprovement=0.0;
885      if (anyBadDown)
886        downImprovement=0.0;
887      double way=0.0;
888      double improvement=0.0;
889      if (downImprovement>0.0&&currentValue>lowerValue) {
890        way=-1.0;
891        improvement = downImprovement;
892      } else if (upImprovement>0.0&&currentValue<upperValue) {
893        way=1.0;
894        improvement = upImprovement;
895      }
896      if (way) {
897        // can improve
898        double distance;
899        if (way>0.0)
900          distance = upperValue-currentValue;
901        else
902          distance = currentValue-lowerValue;
903        for (j=columnStart[iColumn];
904             j<columnStart[iColumn]+columnLength[iColumn];j++) {
905          int iRow = row[j];
906          double value = element[j]*way;
907          if (rowActivity[iRow]>rowUpper[iRow]+primalTolerance) {
908            // infeasible above
909            assert (value<0.0);
910            double gap = rowActivity[iRow]-rowUpper[iRow];
911            if (gap+value*distance<0.0) 
912              distance = -gap/value;
913          } else if (rowActivity[iRow]<rowLower[iRow]-primalTolerance) {
914            // infeasible below
915            assert (value>0.0);
916            double gap = rowActivity[iRow]-rowLower[iRow];
917            if (gap+value*distance>0.0) 
918              distance = -gap/value;
919          } else {
920            // feasible
921            if (value>0) {
922              double gap = rowActivity[iRow]-rowUpper[iRow];
923              if (gap+value*distance>0.0) 
924              distance = -gap/value;
925            } else {
926              double gap = rowActivity[iRow]-rowLower[iRow];
927              if (gap+value*distance<0.0) 
928                distance = -gap/value;
929            }
930          }
931        }
932        //move
933        penaltyChange += improvement*distance;
934        distance *= way;
935        newSolution[iColumn] += distance;
936        newSolutionValue += direction*objective[iColumn]*distance;
937        for (j=columnStart[iColumn];
938             j<columnStart[iColumn]+columnLength[iColumn];j++) {
939          int iRow = row[j];
940          double value = element[j];
941          rowActivity[iRow] += distance*value;
942        }
943      }
944    }
945    // and now all if improving
946    double lastChange= penaltyChange ? 1.0 : 0.0;
947    while (lastChange>1.0e-2) {
948      lastChange=0;
949      for (iColumn=0;iColumn<numberColumns;iColumn++) {
950        bool isInteger = (integerType[iColumn]!=0);
951        double currentValue = newSolution[iColumn];
952        double lowerValue = lower[iColumn];
953        double upperValue = upper[iColumn];
954        int j;
955        int anyBadDown=0;
956        int anyBadUp=0;
957        double upImprovement=0.0;
958        double downImprovement=0.0;
959        for (j=columnStart[iColumn];
960             j<columnStart[iColumn]+columnLength[iColumn];j++) {
961          int iRow = row[j];
962          double value = element[j];
963          if (isInteger) {
964            if (value>0.0) {
965              if (rowActivity[iRow]+value>rowUpper[iRow]+primalTolerance)
966                anyBadUp++;
967              if (rowActivity[iRow]-value<rowLower[iRow]-primalTolerance)
968                anyBadDown++;
969            } else {
970              if (rowActivity[iRow]-value>rowUpper[iRow]+primalTolerance)
971                anyBadDown++;
972              if (rowActivity[iRow]+value<rowLower[iRow]-primalTolerance)
973                anyBadUp++;
974            }
975          }
976          if (rowUpper[iRow]>rowLower[iRow]) {
977            if (rowActivity[iRow]>rowUpper[iRow]+primalTolerance) {
978              // infeasible above
979              downImprovement += value;
980              upImprovement -= value;
981              if (value>0.0) 
982                anyBadUp++;
983              else 
984                anyBadDown++;
985            } else if (rowActivity[iRow]>rowUpper[iRow]-primalTolerance) {
986              // feasible at ub
987              if (value>0.0) {
988                upImprovement -= value;
989                anyBadUp++;
990              } else {
991                downImprovement += value;
992                anyBadDown++;
993              }
994            } else if (rowActivity[iRow]>rowLower[iRow]+primalTolerance) {
995              // feasible in interior
996            } else if (rowActivity[iRow]>rowLower[iRow]-primalTolerance) {
997              // feasible at lb
998              if (value<0.0) {
999                upImprovement += value;
1000                anyBadUp++;
1001              } else {
1002                downImprovement -= value;
1003                anyBadDown++;
1004              }
1005            } else {
1006              // infeasible below
1007              downImprovement -= value;
1008              upImprovement += value;
1009              if (value<0.0) 
1010                anyBadUp++;
1011              else 
1012                anyBadDown++;
1013            }
1014          } else {
1015            // equality row
1016            if (rowActivity[iRow]>rowUpper[iRow]+primalTolerance) {
1017              // infeasible above
1018              downImprovement += value;
1019              upImprovement -= value;
1020              if (value>0.0) 
1021                anyBadUp++;
1022              else 
1023                anyBadDown++;
1024            } else if (rowActivity[iRow]<rowLower[iRow]-primalTolerance) {
1025              // infeasible below
1026              downImprovement -= value;
1027              upImprovement += value;
1028              if (value<0.0) 
1029                anyBadUp++;
1030              else 
1031                anyBadDown++;
1032            } else {
1033              // feasible - no good
1034              anyBadUp=-1;
1035              anyBadDown=-1;
1036              break;
1037            }
1038          }
1039        }
1040        // could change tests for anyBad
1041        if (anyBadUp)
1042          upImprovement=0.0;
1043        if (anyBadDown)
1044          downImprovement=0.0;
1045        double way=0.0;
1046        double improvement=0.0;
1047        if (downImprovement>0.0&&currentValue>lowerValue) {
1048          way=-1.0;
1049          improvement = downImprovement;
1050        } else if (upImprovement>0.0&&currentValue<upperValue) {
1051          way=1.0;
1052          improvement = upImprovement;
1053        }
1054        if (way) {
1055          // can improve
1056          double distance=COIN_DBL_MAX;
1057          for (j=columnStart[iColumn];
1058               j<columnStart[iColumn]+columnLength[iColumn];j++) {
1059            int iRow = row[j];
1060            double value = element[j]*way;
1061            if (rowActivity[iRow]>rowUpper[iRow]+primalTolerance) {
1062              // infeasible above
1063              assert (value<0.0);
1064              double gap = rowActivity[iRow]-rowUpper[iRow];
1065              if (gap+value*distance<0.0) {
1066                // If integer then has to move by 1
1067                if (!isInteger)
1068                  distance = -gap/value;
1069                else
1070                  distance = CoinMax(-gap/value,1.0);
1071              }
1072            } else if (rowActivity[iRow]<rowLower[iRow]-primalTolerance) {
1073              // infeasible below
1074              assert (value>0.0);
1075              double gap = rowActivity[iRow]-rowLower[iRow];
1076              if (gap+value*distance>0.0) {
1077                // If integer then has to move by 1
1078                if (!isInteger)
1079                  distance = -gap/value;
1080                else
1081                  distance = CoinMax(-gap/value,1.0);
1082              }
1083            } else {
1084              // feasible
1085              if (value>0) {
1086                double gap = rowActivity[iRow]-rowUpper[iRow];
1087                if (gap+value*distance>0.0) 
1088                  distance = -gap/value;
1089              } else {
1090                double gap = rowActivity[iRow]-rowLower[iRow];
1091                if (gap+value*distance<0.0) 
1092                  distance = -gap/value;
1093              }
1094            }
1095          }
1096          if (isInteger)
1097            distance = floor(distance+1.05e-8);
1098          if (!distance) {
1099            // should never happen
1100            //printf("zero distance in CbcRounding - debug\n");
1101          }
1102          //move
1103          lastChange += improvement*distance;
1104          distance *= way;
1105          newSolution[iColumn] += distance;
1106          newSolutionValue += direction*objective[iColumn]*distance;
1107          for (j=columnStart[iColumn];
1108               j<columnStart[iColumn]+columnLength[iColumn];j++) {
1109            int iRow = row[j];
1110            double value = element[j];
1111            rowActivity[iRow] += distance*value;
1112          }
1113        }
1114      }
1115      penaltyChange += lastChange;
1116    }
1117    penalty -= penaltyChange;
1118    if (penalty<1.0e-5*fabs(penaltyChange)) {
1119      // recompute
1120      penalty=0.0;
1121      for (i=0;i<numberRows;i++) {
1122        double value = rowActivity[i];
1123        if (value<rowLower[i]-primalTolerance)
1124          penalty += rowLower[i]-value;
1125        else if (value>rowUpper[i]+primalTolerance)
1126          penalty += value-rowUpper[i];
1127      }
1128    }
1129  }
1130
1131  // Could also set SOS (using random) and repeat
1132  if (!penalty) {
1133    // See if we can do better
1134    //seed_++;
1135    //CoinSeedRandom(seed_);
1136    // Random number between 0 and 1.
1137    double randomNumber = randomNumberGenerator_.randomDouble();
1138    int iPass;
1139    int start[2];
1140    int end[2];
1141    int iRandom = (int) (randomNumber*((double) numberIntegers));
1142    start[0]=iRandom;
1143    end[0]=numberIntegers;
1144    start[1]=0;
1145    end[1]=iRandom;
1146    for (iPass=0;iPass<2;iPass++) {
1147      int i;
1148      for (i=start[iPass];i<end[iPass];i++) {
1149        int iColumn = integerVariable[i];
1150#ifndef NDEBUG
1151        double value=newSolution[iColumn];
1152        assert (fabs(floor(value+0.5)-value)<integerTolerance);
1153#endif
1154        double cost = direction * objective[iColumn];
1155        double move=0.0;
1156        if (cost>0.0)
1157          move = -1.0;
1158        else if (cost<0.0)
1159          move=1.0;
1160        while (move) {
1161          bool good=true;
1162          double newValue=newSolution[iColumn]+move;
1163          if (newValue<lower[iColumn]-primalTolerance||
1164              newValue>upper[iColumn]+primalTolerance) {
1165            move=0.0;
1166          } else {
1167            // see if we can move
1168            int j;
1169            for (j=columnStart[iColumn];
1170                 j<columnStart[iColumn]+columnLength[iColumn];j++) {
1171              int iRow = row[j];
1172              double newActivity = rowActivity[iRow] + move*element[j];
1173              if (newActivity<rowLower[iRow]-primalTolerance||
1174                  newActivity>rowUpper[iRow]+primalTolerance) {
1175                good=false;
1176                break;
1177              }
1178            }
1179            if (good) {
1180              newSolution[iColumn] = newValue;
1181              newSolutionValue += move*cost;
1182              int j;
1183              for (j=columnStart[iColumn];
1184                   j<columnStart[iColumn]+columnLength[iColumn];j++) {
1185                int iRow = row[j];
1186                rowActivity[iRow] += move*element[j];
1187              }
1188            } else {
1189              move=0.0;
1190            }
1191          }
1192        }
1193      }
1194    }
1195    // Just in case of some stupidity
1196    double objOffset=0.0;
1197    solver->getDblParam(OsiObjOffset,objOffset);
1198    newSolutionValue = -objOffset;
1199    for ( i=0 ; i<numberColumns ; i++ )
1200      newSolutionValue += objective[i]*newSolution[i];
1201    newSolutionValue *= direction;
1202    //printf("new solution value %g %g\n",newSolutionValue,solutionValue);
1203    if (newSolutionValue<solutionValue) {
1204      // paranoid check
1205      memset(rowActivity,0,numberRows*sizeof(double));
1206      for (i=0;i<numberColumns;i++) {
1207        int j;
1208        double value = newSolution[i];
1209        if (value) {
1210          for (j=columnStart[i];
1211               j<columnStart[i]+columnLength[i];j++) {
1212            int iRow=row[j];
1213            rowActivity[iRow] += value*element[j];
1214          }
1215        }
1216      }
1217      // check was approximately feasible
1218      bool feasible=true;
1219      for (i=0;i<numberRows;i++) {
1220        if(rowActivity[i]<rowLower[i]) {
1221          if (rowActivity[i]<rowLower[i]-1000.0*primalTolerance)
1222            feasible = false;
1223        } else if(rowActivity[i]>rowUpper[i]) {
1224          if (rowActivity[i]>rowUpper[i]+1000.0*primalTolerance)
1225            feasible = false;
1226        }
1227      }
1228      if (feasible) {
1229        // new solution
1230        memcpy(betterSolution,newSolution,numberColumns*sizeof(double));
1231        solutionValue = newSolutionValue;
1232        //printf("** Solution of %g found by rounding\n",newSolutionValue);
1233        returnCode=1;
1234      } else {
1235        // Can easily happen
1236        //printf("Debug CbcRounding giving bad solution\n");
1237      }
1238    }
1239  }
1240#ifdef NEW_ROUNDING
1241  if (!returnCode) {
1242#if 0
1243    // back to starting point
1244    memcpy(newSolution,solution,numberColumns*sizeof(double));
1245    memset(rowActivity,0,numberRows*sizeof(double));
1246    for (i=0;i<numberColumns;i++) {
1247      int j;
1248      double value = newSolution[i];
1249      if (value<lower[i]) {
1250        value=lower[i];
1251        newSolution[i]=value;
1252      } else if (value>upper[i]) {
1253        value=upper[i];
1254        newSolution[i]=value;
1255      }
1256      if (value) {
1257        for (j=columnStart[i];
1258             j<columnStart[i]+columnLength[i];j++) {
1259          int iRow=row[j];
1260          rowActivity[iRow] += value*element[j];
1261        }
1262      }
1263    }
1264    // check was feasible - if not adjust (cleaning may move)
1265    for (i=0;i<numberRows;i++) {
1266      if(rowActivity[i]<rowLower[i]) {
1267        //assert (rowActivity[i]>rowLower[i]-1000.0*primalTolerance);
1268        rowActivity[i]=rowLower[i];
1269      } else if(rowActivity[i]>rowUpper[i]) {
1270        //assert (rowActivity[i]<rowUpper[i]+1000.0*primalTolerance);
1271        rowActivity[i]=rowUpper[i];
1272      }
1273    }
1274#endif
1275    int * candidate = new int [numberColumns];
1276    int nCandidate=0;
1277    for (iColumn=0;iColumn<numberColumns;iColumn++) {
1278      bool isInteger = (integerType[iColumn]!=0);
1279      if (isInteger) {
1280        double currentValue = newSolution[iColumn];
1281        if (fabs(currentValue-floor(currentValue+0.5))>1.0e-8)
1282          candidate[nCandidate++]=iColumn;
1283      }
1284    }
1285    if (true) {
1286      // Rounding as in Berthold
1287      while (nCandidate) {
1288        double infeasibility =1.0e-7;
1289        int iRow=-1;
1290        for (i=0;i<numberRows;i++) {
1291          double value=0.0;
1292          if(rowActivity[i]<rowLower[i]) {
1293            value = rowLower[i]-rowActivity[i];
1294          } else if(rowActivity[i]>rowUpper[i]) {
1295            value = rowActivity[i]-rowUpper[i];
1296          }
1297          if (value>infeasibility) {
1298            infeasibility = value;
1299            iRow=i;
1300          }
1301        }
1302        if (iRow>=0) {
1303          // infeasible
1304        } else {
1305          // feasible
1306        }
1307      }
1308    } else {
1309      // Shifting as in Berthold
1310    }
1311    delete [] candidate;
1312  }
1313#endif
1314  delete [] newSolution;
1315  delete [] rowActivity;
1316  return returnCode;
1317}
1318// update model
1319void CbcRounding::setModel(CbcModel * model)
1320{
1321  model_ = model;
1322  // Get a copy of original matrix (and by row for rounding);
1323  assert(model_->solver());
1324  matrix_ = *model_->solver()->getMatrixByCol();
1325  matrixByRow_ = *model_->solver()->getMatrixByRow();
1326  // make sure model okay for heuristic
1327  validate();
1328}
1329// Validate model i.e. sets when_ to 0 if necessary (may be NULL)
1330void 
1331CbcRounding::validate() 
1332{
1333  if (model_&&when()<10) {
1334    if (model_->numberIntegers()!=
1335        model_->numberObjects())
1336      setWhen(0);
1337  }
1338#ifdef NEW_ROUNDING
1339  int numberColumns = matrix_.getNumCols();
1340  down_ = new unsigned short [numberColumns];
1341  up_ = new unsigned short [numberColumns];
1342  equal_ = new unsigned short [numberColumns];
1343  // Column copy
1344  const double * element = matrix_.getElements();
1345  const int * row = matrix_.getIndices();
1346  const CoinBigIndex * columnStart = matrix_.getVectorStarts();
1347  const int * columnLength = matrix_.getVectorLengths();
1348  const double * rowLower = model.solver()->getRowLower();
1349  const double * rowUpper = model.solver()->getRowUpper();
1350  for (int i=0;i<numberColumns;i++) {
1351    int down=0;
1352    int up=0;
1353    int equal=0;
1354    if (columnLength[i]>65535) {
1355      equal[0]=65535; 
1356      break; // unlikely to work
1357    }
1358    for (CoinBigIndex j=columnStart[i];
1359         j<columnStart[i]+columnLength[i];j++) {
1360      int iRow=row[j];
1361      if (rowLower[iRow]>-1.0e20&&rowUpper[iRow]<1.0e20) {
1362        equal++;
1363      } else if (element[j]>0.0) {
1364        if (rowUpper[iRow]<1.0e20)
1365          up++;
1366        else
1367          down--;
1368      } else {
1369        if (rowLower[iRow]>-1.0e20)
1370          up++;
1371        else
1372          down--;
1373      }
1374    }
1375    down_[i] = (unsigned short) down;
1376    up_[i] = (unsigned short) up;
1377    equal_[i] = (unsigned short) equal;
1378  }
1379#else
1380  down_ = NULL;
1381  up_ = NULL;
1382  equal_ = NULL;
1383#endif 
1384}
1385
1386// Default Constructor
1387CbcHeuristicPartial::CbcHeuristicPartial() 
1388  :CbcHeuristic()
1389{
1390  fixPriority_ = 10000;
1391}
1392
1393// Constructor from model
1394CbcHeuristicPartial::CbcHeuristicPartial(CbcModel & model, int fixPriority, int numberNodes)
1395  :CbcHeuristic(model)
1396{
1397  fixPriority_ = fixPriority;
1398  setNumberNodes(numberNodes);
1399  validate();
1400}
1401
1402// Destructor
1403CbcHeuristicPartial::~CbcHeuristicPartial ()
1404{
1405}
1406
1407// Clone
1408CbcHeuristic *
1409CbcHeuristicPartial::clone() const
1410{
1411  return new CbcHeuristicPartial(*this);
1412}
1413// Create C++ lines to get to current state
1414void 
1415CbcHeuristicPartial::generateCpp( FILE * fp) 
1416{
1417  CbcHeuristicPartial other;
1418  fprintf(fp,"0#include \"CbcHeuristic.hpp\"\n");
1419  fprintf(fp,"3  CbcHeuristicPartial partial(*cbcModel);\n");
1420  CbcHeuristic::generateCpp(fp,"partial");
1421  if (fixPriority_!=other.fixPriority_)
1422    fprintf(fp,"3  partial.setFixPriority(%d);\n",fixPriority_);
1423  else
1424    fprintf(fp,"4  partial.setFixPriority(%d);\n",fixPriority_);
1425  fprintf(fp,"3  cbcModel->addHeuristic(&partial);\n");
1426}
1427//#define NEW_PARTIAL
1428// Copy constructor
1429CbcHeuristicPartial::CbcHeuristicPartial(const CbcHeuristicPartial & rhs)
1430:
1431  CbcHeuristic(rhs),
1432  fixPriority_(rhs.fixPriority_)
1433{
1434}
1435
1436// Assignment operator
1437CbcHeuristicPartial & 
1438CbcHeuristicPartial::operator=( const CbcHeuristicPartial& rhs)
1439{
1440  if (this!=&rhs) {
1441    CbcHeuristic::operator=(rhs);
1442    fixPriority_ = rhs.fixPriority_;
1443  }
1444  return *this;
1445}
1446
1447// Resets stuff if model changes
1448void 
1449CbcHeuristicPartial::resetModel(CbcModel * model)
1450{
1451  model_=model;
1452  // Get a copy of original matrix (and by row for partial);
1453  assert(model_->solver());
1454  validate();
1455}
1456// See if partial will give solution
1457// Sets value of solution
1458// Assumes rhs for original matrix still okay
1459// At present only works with integers
1460// Fix values if asked for
1461// Returns 1 if solution, 0 if not
1462int
1463CbcHeuristicPartial::solution(double & solutionValue,
1464                      double * betterSolution)
1465{
1466  // Return if already done
1467  if (fixPriority_<0)
1468    return 0; // switched off
1469  const double * hotstartSolution = model_->hotstartSolution();
1470  const int * hotstartPriorities = model_->hotstartPriorities();
1471  if (!hotstartSolution)
1472    return 0;
1473  OsiSolverInterface * solver = model_->solver();
1474 
1475  int numberIntegers = model_->numberIntegers();
1476  const int * integerVariable = model_->integerVariable();
1477 
1478  OsiSolverInterface * newSolver = model_->continuousSolver()->clone();
1479  const double * colLower = newSolver->getColLower();
1480  const double * colUpper = newSolver->getColUpper();
1481
1482  double primalTolerance;
1483  solver->getDblParam(OsiPrimalTolerance,primalTolerance);
1484   
1485  int i;
1486  int numberFixed=0;
1487  int returnCode=0;
1488
1489  for (i=0;i<numberIntegers;i++) {
1490    int iColumn=integerVariable[i];
1491    if (abs(hotstartPriorities[iColumn])<=fixPriority_) {
1492      double value = hotstartSolution[iColumn];
1493      double lower = colLower[iColumn];
1494      double upper = colUpper[iColumn];
1495      value = CoinMax(value,lower);
1496      value = CoinMin(value,upper);
1497      if (fabs(value-floor(value+0.5))<1.0e-8) {
1498        value = floor(value+0.5);
1499        newSolver->setColLower(iColumn,value);
1500        newSolver->setColUpper(iColumn,value);
1501        numberFixed++;
1502      }
1503    }
1504  }
1505  if (numberFixed>numberIntegers/5-100000000) {
1506#ifdef COIN_DEVELOP
1507    printf("%d integers fixed\n",numberFixed);
1508#endif
1509    returnCode = smallBranchAndBound(newSolver,numberNodes_,betterSolution,solutionValue,
1510                                     model_->getCutoff(),"CbcHeuristicPartial");
1511    if (returnCode<0)
1512      returnCode=0; // returned on size
1513    //printf("return code %d",returnCode);
1514    if ((returnCode&2)!=0) {
1515      // could add cut
1516      returnCode &= ~2;
1517      //printf("could add cut with %d elements (if all 0-1)\n",nFix);
1518    } else {
1519      //printf("\n");
1520    }
1521  }
1522  fixPriority_=-1; // switch off
1523 
1524  delete newSolver;
1525  return returnCode;
1526}
1527// update model
1528void CbcHeuristicPartial::setModel(CbcModel * model)
1529{
1530  model_ = model;
1531  assert(model_->solver());
1532  // make sure model okay for heuristic
1533  validate();
1534}
1535// Validate model i.e. sets when_ to 0 if necessary (may be NULL)
1536void 
1537CbcHeuristicPartial::validate() 
1538{
1539  if (model_&&when()<10) {
1540    if (model_->numberIntegers()!=
1541        model_->numberObjects())
1542      setWhen(0);
1543  }
1544}
1545
1546// Default Constructor
1547CbcSerendipity::CbcSerendipity() 
1548  :CbcHeuristic()
1549{
1550}
1551
1552// Constructor from model
1553CbcSerendipity::CbcSerendipity(CbcModel & model)
1554  :CbcHeuristic(model)
1555{
1556}
1557
1558// Destructor
1559CbcSerendipity::~CbcSerendipity ()
1560{
1561}
1562
1563// Clone
1564CbcHeuristic *
1565CbcSerendipity::clone() const
1566{
1567  return new CbcSerendipity(*this);
1568}
1569// Create C++ lines to get to current state
1570void 
1571CbcSerendipity::generateCpp( FILE * fp) 
1572{
1573  fprintf(fp,"0#include \"CbcHeuristic.hpp\"\n");
1574  fprintf(fp,"3  CbcSerendipity serendipity(*cbcModel);\n");
1575  CbcHeuristic::generateCpp(fp,"serendipity");
1576  fprintf(fp,"3  cbcModel->addHeuristic(&serendipity);\n");
1577}
1578
1579// Copy constructor
1580CbcSerendipity::CbcSerendipity(const CbcSerendipity & rhs)
1581:
1582  CbcHeuristic(rhs)
1583{
1584}
1585
1586// Assignment operator
1587CbcSerendipity & 
1588CbcSerendipity::operator=( const CbcSerendipity& rhs)
1589{
1590  if (this!=&rhs) {
1591    CbcHeuristic::operator=(rhs);
1592  }
1593  return *this;
1594}
1595
1596// Returns 1 if solution, 0 if not
1597int
1598CbcSerendipity::solution(double & solutionValue,
1599                         double * betterSolution)
1600{
1601  if (!model_)
1602    return 0;
1603  // get information on solver type
1604  OsiAuxInfo * auxInfo = model_->solver()->getAuxiliaryInfo();
1605  OsiBabSolver * auxiliaryInfo = dynamic_cast< OsiBabSolver *> (auxInfo);
1606  if (auxiliaryInfo)
1607    return auxiliaryInfo->solution(solutionValue,betterSolution,model_->solver()->getNumCols());
1608  else
1609    return 0;
1610}
1611// update model
1612void CbcSerendipity::setModel(CbcModel * model)
1613{
1614  model_ = model;
1615}
1616// Resets stuff if model changes
1617void 
1618CbcSerendipity::resetModel(CbcModel * model)
1619{
1620  model_ = model;
1621}
1622 
Note: See TracBrowser for help on using the repository browser.