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

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

changes to heuristics and allow collection of more statistics

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