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

Last change on this file since 863 was 863, checked in by forrest, 12 years ago

make print line longer

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