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

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

include cstdlib before cmath to get things to compile on AIX with xlC

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