source: stable/2.0/Cbc/src/CbcHeuristic.cpp @ 905

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

include cstdlib before cmath to get things to compile on AIX with xlC (same as changeset 904 in trunk)

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