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

Last change on this file since 310 was 310, checked in by andreasw, 13 years ago

first commit for autotools conversion to be able to move more files

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 27.2 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#include <cassert>
8#include <cmath>
9#include <cfloat>
10
11#ifdef CBC_USE_CLP
12#include "OsiClpSolverInterface.hpp"
13#endif
14#include "CbcModel.hpp"
15#include "CbcMessage.hpp"
16#include "CbcHeuristic.hpp"
17#include "CbcStrategy.hpp"
18#include "CglPreProcess.hpp"
19#include "OsiAuxInfo.hpp"
20
21// Default Constructor
22CbcHeuristic::CbcHeuristic() 
23  :model_(NULL),
24   when_(2)
25{
26}
27
28// Constructor from model
29CbcHeuristic::CbcHeuristic(CbcModel & model)
30:
31  model_(&model),
32  when_(2)
33{
34}
35// Resets stuff if model changes
36void 
37CbcHeuristic::resetModel(CbcModel * model)
38{
39  model_=model;
40}
41
42// Destructor
43CbcHeuristic::~CbcHeuristic ()
44{
45}
46
47// update model
48void CbcHeuristic::setModel(CbcModel * model)
49{
50  model_ = model;
51}
52// Do mini branch and bound (return 1 if solution)
53int 
54CbcHeuristic::smallBranchAndBound(OsiSolverInterface * solver,int numberNodes,
55                                  double * newSolution, double & newSolutionValue,
56                                  double cutoff, std::string name) const
57{
58#ifdef CBC_USE_CLP
59  OsiClpSolverInterface * osiclp = dynamic_cast< OsiClpSolverInterface*> (solver);
60  if (osiclp) {
61    // go faster stripes
62    if (osiclp->getNumRows()<300&&osiclp->getNumCols()<500) {
63      osiclp->setupForRepeatedUse(2,0);
64    } else {
65      osiclp->setupForRepeatedUse(0,0);
66    }
67    // Turn this off if you get problems
68    // Used to be automatically set
69    osiclp->setSpecialOptions(osiclp->specialOptions()|(128+64));
70  }
71#endif
72  // Reduce printout
73  solver->setHintParam(OsiDoReducePrint,true,OsiHintTry);
74  solver->setHintParam(OsiDoPresolveInInitial,false,OsiHintTry);
75  solver->setDblParam(OsiDualObjectiveLimit,cutoff);
76  solver->initialSolve();
77  int returnCode=1;
78  int logLevel = model_->logLevel();
79  if (solver->isProvenOptimal()) {
80    CglPreProcess process;
81    /* Do not try and produce equality cliques and
82       do up to 5 passes */
83    if (logLevel<=1)
84      process.messageHandler()->setLogLevel(0);
85    OsiSolverInterface * solver2= process.preProcess(*solver);
86    if (!solver2) {
87      if (logLevel>1)
88        printf("Pre-processing says infeasible\n");
89      returnCode=0; // so will be infeasible
90    } else {
91      solver2->resolve();
92      CbcModel model(*solver2);
93      if (logLevel<=1)
94        model.setLogLevel(0);
95      else
96        model.setLogLevel(logLevel);
97      model.setCutoff(cutoff);
98      model.setMaximumNodes(numberNodes);
99      model.solver()->setHintParam(OsiDoReducePrint,true,OsiHintTry);
100      // Lightweight
101      CbcStrategyDefaultSubTree strategy(model_,true,5,1,0);
102      model.setStrategy(strategy);
103      model.solver()->setIntParam(OsiMaxNumIterationHotStart,10);
104      // Do search
105      if (logLevel>1)
106        model_->messageHandler()->message(CBC_START_SUB,model_->messages())
107          << name
108          << model.getMaximumNodes()
109          <<CoinMessageEol;
110      // probably faster to use a basis to get integer solutions
111      model.setSpecialOptions(2);
112      model.branchAndBound();
113      if (logLevel>1)
114        model_->messageHandler()->message(CBC_END_SUB,model_->messages())
115          << name
116          <<CoinMessageEol;
117      if (model.getMinimizationObjValue()<cutoff) {
118        // solution
119        returnCode=1;
120        // post process
121        process.postProcess(*model.solver());
122        if (solver->isProvenOptimal()) {
123          // Solution now back in solver
124          int numberColumns = solver->getNumCols();
125          memcpy(newSolution,solver->getColSolution(),
126                 numberColumns*sizeof(double));
127          newSolutionValue = model.getMinimizationObjValue();
128        } else {
129          // odd - but no good
130          returnCode=0; // so will be infeasible
131        }
132      } else {
133        // no good
134        returnCode=0; // so will be infeasible
135      }
136    }
137  } else {
138    returnCode=0;
139  }
140  return returnCode;
141}
142
143// Default Constructor
144CbcRounding::CbcRounding() 
145  :CbcHeuristic()
146{
147  // matrix and row copy will automatically be empty
148}
149
150// Constructor from model
151CbcRounding::CbcRounding(CbcModel & model)
152  :CbcHeuristic(model)
153{
154  // Get a copy of original matrix (and by row for rounding);
155  assert(model.solver());
156  matrix_ = *model.solver()->getMatrixByCol();
157  matrixByRow_ = *model.solver()->getMatrixByRow();
158  seed_=1;
159}
160
161// Destructor
162CbcRounding::~CbcRounding ()
163{
164}
165
166// Clone
167CbcHeuristic *
168CbcRounding::clone() const
169{
170  return new CbcRounding(*this);
171}
172
173// Copy constructor
174CbcRounding::CbcRounding(const CbcRounding & rhs)
175:
176  CbcHeuristic(rhs),
177  matrix_(rhs.matrix_),
178  matrixByRow_(rhs.matrixByRow_),
179  seed_(rhs.seed_)
180{
181  setWhen(rhs.when());
182}
183
184// Resets stuff if model changes
185void 
186CbcRounding::resetModel(CbcModel * model)
187{
188  model_=model;
189  // Get a copy of original matrix (and by row for rounding);
190  assert(model_->solver());
191  matrix_ = *model_->solver()->getMatrixByCol();
192  matrixByRow_ = *model_->solver()->getMatrixByRow();
193}
194// See if rounding will give solution
195// Sets value of solution
196// Assumes rhs for original matrix still okay
197// At present only works with integers
198// Fix values if asked for
199// Returns 1 if solution, 0 if not
200int
201CbcRounding::solution(double & solutionValue,
202                         double * betterSolution)
203{
204
205  // See if to do
206  if (!when()||(when()%10==1&&model_->phase()!=1)||
207      (when()%10==2&&(model_->phase()!=2&&model_->phase()!=3)))
208    return 0; // switched off
209
210  OsiSolverInterface * solver = model_->solver();
211  const double * lower = solver->getColLower();
212  const double * upper = solver->getColUpper();
213  const double * rowLower = solver->getRowLower();
214  const double * rowUpper = solver->getRowUpper();
215  const double * solution = solver->getColSolution();
216  const double * objective = solver->getObjCoefficients();
217  double integerTolerance = model_->getDblParam(CbcModel::CbcIntegerTolerance);
218  double primalTolerance;
219  solver->getDblParam(OsiPrimalTolerance,primalTolerance);
220
221  int numberRows = matrix_.getNumRows();
222
223  int numberIntegers = model_->numberIntegers();
224  const int * integerVariable = model_->integerVariable();
225  int i;
226  double direction = solver->getObjSense();
227  double newSolutionValue = direction*solver->getObjValue();
228  int returnCode = 0;
229
230  // Column copy
231  const double * element = matrix_.getElements();
232  const int * row = matrix_.getIndices();
233  const CoinBigIndex * columnStart = matrix_.getVectorStarts();
234  const int * columnLength = matrix_.getVectorLengths();
235  // Row copy
236  const double * elementByRow = matrixByRow_.getElements();
237  const int * column = matrixByRow_.getIndices();
238  const CoinBigIndex * rowStart = matrixByRow_.getVectorStarts();
239  const int * rowLength = matrixByRow_.getVectorLengths();
240
241  // Get solution array for heuristic solution
242  int numberColumns = solver->getNumCols();
243  double * newSolution = new double [numberColumns];
244  memcpy(newSolution,solution,numberColumns*sizeof(double));
245
246  double * rowActivity = new double[numberRows];
247  memset(rowActivity,0,numberRows*sizeof(double));
248  for (i=0;i<numberColumns;i++) {
249    int j;
250    double value = newSolution[i];
251    if (value<lower[i]) {
252      value=lower[i];
253      newSolution[i]=value;
254    } else if (value>upper[i]) {
255      value=upper[i];
256      newSolution[i]=value;
257    }
258    if (value) {
259      for (j=columnStart[i];
260           j<columnStart[i]+columnLength[i];j++) {
261        int iRow=row[j];
262        rowActivity[iRow] += value*element[j];
263      }
264    }
265  }
266  // check was feasible - if not adjust (cleaning may move)
267  for (i=0;i<numberRows;i++) {
268    if(rowActivity[i]<rowLower[i]) {
269      //assert (rowActivity[i]>rowLower[i]-1000.0*primalTolerance);
270      rowActivity[i]=rowLower[i];
271    } else if(rowActivity[i]>rowUpper[i]) {
272      //assert (rowActivity[i]<rowUpper[i]+1000.0*primalTolerance);
273      rowActivity[i]=rowUpper[i];
274    }
275  }
276  for (i=0;i<numberIntegers;i++) {
277    int iColumn = integerVariable[i];
278    double value=newSolution[iColumn];
279    if (fabs(floor(value+0.5)-value)>integerTolerance) {
280      double below = floor(value);
281      double newValue=newSolution[iColumn];
282      double cost = direction * objective[iColumn];
283      double move;
284      if (cost>0.0) {
285        // try up
286        move = 1.0 -(value-below);
287      } else if (cost<0.0) {
288        // try down
289        move = below-value;
290      } else {
291        // won't be able to move unless we can grab another variable
292        double randomNumber = CoinDrand48();
293        // which way?
294        if (randomNumber<0.5) 
295          move = below-value;
296        else
297          move = 1.0 -(value-below);
298      }
299      newValue += move;
300      newSolution[iColumn] = newValue;
301      newSolutionValue += move*cost;
302      int j;
303      for (j=columnStart[iColumn];
304           j<columnStart[iColumn]+columnLength[iColumn];j++) {
305        int iRow = row[j];
306        rowActivity[iRow] += move*element[j];
307      }
308    }
309  }
310
311  double penalty=0.0;
312  const char * integerType = model_->integerType();
313  // see if feasible - just using singletons
314  for (i=0;i<numberRows;i++) {
315    double value = rowActivity[i];
316    double thisInfeasibility=0.0;
317    if (value<rowLower[i]-primalTolerance)
318      thisInfeasibility = value-rowLower[i];
319    else if (value>rowUpper[i]+primalTolerance)
320      thisInfeasibility = value-rowUpper[i];
321    if (thisInfeasibility) {
322      // See if there are any slacks I can use to fix up
323      // maybe put in coding for multiple slacks?
324      double bestCost = 1.0e50;
325      int k;
326      int iBest=-1;
327      double addCost=0.0;
328      double newValue=0.0;
329      double changeRowActivity=0.0;
330      double absInfeasibility = fabs(thisInfeasibility);
331      for (k=rowStart[i];k<rowStart[i]+rowLength[i];k++) {
332        int iColumn = column[k];
333        // See if all elements help
334        if (columnLength[iColumn]==1) {
335          double currentValue = newSolution[iColumn];
336          double elementValue = elementByRow[k];
337          double lowerValue = lower[iColumn];
338          double upperValue = upper[iColumn];
339          double gap = rowUpper[i]-rowLower[i];
340          double absElement=fabs(elementValue);
341          if (thisInfeasibility*elementValue>0.0) {
342            // we want to reduce
343            if ((currentValue-lowerValue)*absElement>=absInfeasibility) {
344              // possible - check if integer
345              double distance = absInfeasibility/absElement;
346              double thisCost = -direction*objective[iColumn]*distance;
347              if (integerType[iColumn]) {
348                distance = ceil(distance-primalTolerance);
349                if (currentValue-distance>=lowerValue-primalTolerance) {
350                  if (absInfeasibility-distance*absElement< -gap-primalTolerance)
351                    thisCost=1.0e100; // no good
352                  else
353                    thisCost = -direction*objective[iColumn]*distance;
354                } else {
355                  thisCost=1.0e100; // no good
356                }
357              }
358              if (thisCost<bestCost) {
359                bestCost=thisCost;
360                iBest=iColumn;
361                addCost = thisCost;
362                newValue = currentValue-distance;
363                changeRowActivity = -distance*elementValue;
364              }
365            }
366          } else {
367            // we want to increase
368            if ((upperValue-currentValue)*absElement>=absInfeasibility) {
369              // possible - check if integer
370              double distance = absInfeasibility/absElement;
371              double thisCost = direction*objective[iColumn]*distance;
372              if (integerType[iColumn]) {
373                distance = ceil(distance-1.0e-7);
374                assert (currentValue-distance<=upperValue+primalTolerance);
375                if (absInfeasibility-distance*absElement< -gap-primalTolerance)
376                  thisCost=1.0e100; // no good
377                else
378                  thisCost = direction*objective[iColumn]*distance;
379              }
380              if (thisCost<bestCost) {
381                bestCost=thisCost;
382                iBest=iColumn;
383                addCost = thisCost;
384                newValue = currentValue+distance;
385                changeRowActivity = distance*elementValue;
386              }
387            }
388          }
389        }
390      }
391      if (iBest>=0) {
392        /*printf("Infeasibility of %g on row %d cost %g\n",
393          thisInfeasibility,i,addCost);*/
394        newSolution[iBest]=newValue;
395        thisInfeasibility=0.0;
396        newSolutionValue += addCost;
397        rowActivity[i] += changeRowActivity;
398      }
399      penalty += fabs(thisInfeasibility);
400    }
401  }
402  if (penalty) {
403    // see if feasible using any
404    // first continuous
405    double penaltyChange=0.0;
406    int iColumn;
407    for (iColumn=0;iColumn<numberColumns;iColumn++) {
408      if (integerType[iColumn])
409        continue;
410      double currentValue = newSolution[iColumn];
411      double lowerValue = lower[iColumn];
412      double upperValue = upper[iColumn];
413      int j;
414      int anyBadDown=0;
415      int anyBadUp=0;
416      double upImprovement=0.0;
417      double downImprovement=0.0;
418      for (j=columnStart[iColumn];
419           j<columnStart[iColumn]+columnLength[iColumn];j++) {
420        int iRow = row[j];
421        if (rowUpper[iRow]>rowLower[iRow]) {
422          double value = element[j];
423          if (rowActivity[iRow]>rowUpper[iRow]+primalTolerance) {
424            // infeasible above
425            downImprovement += value;
426            upImprovement -= value;
427            if (value>0.0) 
428              anyBadUp++;
429            else 
430              anyBadDown++;
431          } else if (rowActivity[iRow]>rowUpper[iRow]-primalTolerance) {
432            // feasible at ub
433            if (value>0.0) {
434              upImprovement -= value;
435              anyBadUp++;
436            } else {
437              downImprovement += value;
438              anyBadDown++;
439            }
440          } else if (rowActivity[iRow]>rowLower[iRow]+primalTolerance) {
441            // feasible in interior
442          } else if (rowActivity[iRow]>rowLower[iRow]-primalTolerance) {
443            // feasible at lb
444            if (value<0.0) {
445              upImprovement += value;
446              anyBadUp++;
447            } else {
448              downImprovement -= value;
449              anyBadDown++;
450            }
451          } else {
452            // infeasible below
453            downImprovement -= value;
454            upImprovement += value;
455            if (value<0.0) 
456              anyBadUp++;
457            else 
458              anyBadDown++;
459          }
460        } else {
461          // equality row
462          double value = element[j];
463          if (rowActivity[iRow]>rowUpper[iRow]+primalTolerance) {
464            // infeasible above
465            downImprovement += value;
466            upImprovement -= value;
467            if (value>0.0) 
468              anyBadUp++;
469            else 
470              anyBadDown++;
471          } else if (rowActivity[iRow]<rowLower[iRow]-primalTolerance) {
472            // infeasible below
473            downImprovement -= value;
474            upImprovement += value;
475            if (value<0.0) 
476              anyBadUp++;
477            else 
478              anyBadDown++;
479          } else {
480            // feasible - no good
481            anyBadUp=-1;
482            anyBadDown=-1;
483            break;
484          }
485        }
486      }
487      // could change tests for anyBad
488      if (anyBadUp)
489        upImprovement=0.0;
490      if (anyBadDown)
491        downImprovement=0.0;
492      double way=0.0;
493      double improvement=0.0;
494      if (downImprovement>0.0&&currentValue>lowerValue) {
495        way=-1.0;
496        improvement = downImprovement;
497      } else if (upImprovement>0.0&&currentValue<upperValue) {
498        way=1.0;
499        improvement = upImprovement;
500      }
501      if (way) {
502        // can improve
503        double distance;
504        if (way>0.0)
505          distance = upperValue-currentValue;
506        else
507          distance = currentValue-lowerValue;
508        for (j=columnStart[iColumn];
509             j<columnStart[iColumn]+columnLength[iColumn];j++) {
510          int iRow = row[j];
511          double value = element[j]*way;
512          if (rowActivity[iRow]>rowUpper[iRow]+primalTolerance) {
513            // infeasible above
514            assert (value<0.0);
515            double gap = rowActivity[iRow]-rowUpper[iRow];
516            if (gap+value*distance<0.0) 
517              distance = -gap/value;
518          } else if (rowActivity[iRow]<rowLower[iRow]-primalTolerance) {
519            // infeasible below
520            assert (value>0.0);
521            double gap = rowActivity[iRow]-rowLower[iRow];
522            if (gap+value*distance>0.0) 
523              distance = -gap/value;
524          } else {
525            // feasible
526            if (value>0) {
527              double gap = rowActivity[iRow]-rowUpper[iRow];
528              if (gap+value*distance>0.0) 
529              distance = -gap/value;
530            } else {
531              double gap = rowActivity[iRow]-rowLower[iRow];
532              if (gap+value*distance<0.0) 
533                distance = -gap/value;
534            }
535          }
536        }
537        //move
538        penaltyChange += improvement*distance;
539        distance *= way;
540        newSolution[iColumn] += distance;
541        newSolutionValue += direction*objective[iColumn]*distance;
542        for (j=columnStart[iColumn];
543             j<columnStart[iColumn]+columnLength[iColumn];j++) {
544          int iRow = row[j];
545          double value = element[j];
546          rowActivity[iRow] += distance*value;
547        }
548      }
549    }
550    // and now all if improving
551    double lastChange= penaltyChange ? 1.0 : 0.0;
552    while (lastChange>1.0e-2) {
553      lastChange=0;
554      for (iColumn=0;iColumn<numberColumns;iColumn++) {
555        bool isInteger = integerType[iColumn];
556        double currentValue = newSolution[iColumn];
557        double lowerValue = lower[iColumn];
558        double upperValue = upper[iColumn];
559        int j;
560        int anyBadDown=0;
561        int anyBadUp=0;
562        double upImprovement=0.0;
563        double downImprovement=0.0;
564        for (j=columnStart[iColumn];
565             j<columnStart[iColumn]+columnLength[iColumn];j++) {
566          int iRow = row[j];
567          double value = element[j];
568          if (isInteger) {
569            if (value>0.0) {
570              if (rowActivity[iRow]+value>rowUpper[iRow]+primalTolerance)
571                anyBadUp++;
572              if (rowActivity[iRow]-value<rowLower[iRow]-primalTolerance)
573                anyBadDown++;
574            } else {
575              if (rowActivity[iRow]-value>rowUpper[iRow]+primalTolerance)
576                anyBadDown++;
577              if (rowActivity[iRow]+value<rowLower[iRow]-primalTolerance)
578                anyBadUp++;
579            }
580          }
581          if (rowUpper[iRow]>rowLower[iRow]) {
582            if (rowActivity[iRow]>rowUpper[iRow]+primalTolerance) {
583              // infeasible above
584              downImprovement += value;
585              upImprovement -= value;
586              if (value>0.0) 
587                anyBadUp++;
588              else 
589                anyBadDown++;
590            } else if (rowActivity[iRow]>rowUpper[iRow]-primalTolerance) {
591              // feasible at ub
592              if (value>0.0) {
593                upImprovement -= value;
594                anyBadUp++;
595              } else {
596                downImprovement += value;
597                anyBadDown++;
598              }
599            } else if (rowActivity[iRow]>rowLower[iRow]+primalTolerance) {
600              // feasible in interior
601            } else if (rowActivity[iRow]>rowLower[iRow]-primalTolerance) {
602              // feasible at lb
603              if (value<0.0) {
604                upImprovement += value;
605                anyBadUp++;
606              } else {
607                downImprovement -= value;
608                anyBadDown++;
609              }
610            } else {
611              // infeasible below
612              downImprovement -= value;
613              upImprovement += value;
614              if (value<0.0) 
615                anyBadUp++;
616              else 
617                anyBadDown++;
618            }
619          } else {
620            // equality row
621            if (rowActivity[iRow]>rowUpper[iRow]+primalTolerance) {
622              // infeasible above
623              downImprovement += value;
624              upImprovement -= value;
625              if (value>0.0) 
626                anyBadUp++;
627              else 
628                anyBadDown++;
629            } else if (rowActivity[iRow]<rowLower[iRow]-primalTolerance) {
630              // infeasible below
631              downImprovement -= value;
632              upImprovement += value;
633              if (value<0.0) 
634                anyBadUp++;
635              else 
636                anyBadDown++;
637            } else {
638              // feasible - no good
639              anyBadUp=-1;
640              anyBadDown=-1;
641              break;
642            }
643          }
644        }
645        // could change tests for anyBad
646        if (anyBadUp)
647          upImprovement=0.0;
648        if (anyBadDown)
649          downImprovement=0.0;
650        double way=0.0;
651        double improvement=0.0;
652        if (downImprovement>0.0&&currentValue>lowerValue) {
653          way=-1.0;
654          improvement = downImprovement;
655        } else if (upImprovement>0.0&&currentValue<upperValue) {
656          way=1.0;
657          improvement = upImprovement;
658        }
659        if (way) {
660          // can improve
661          double distance=COIN_DBL_MAX;
662          for (j=columnStart[iColumn];
663               j<columnStart[iColumn]+columnLength[iColumn];j++) {
664            int iRow = row[j];
665            double value = element[j]*way;
666            if (rowActivity[iRow]>rowUpper[iRow]+primalTolerance) {
667              // infeasible above
668              assert (value<0.0);
669              double gap = rowActivity[iRow]-rowUpper[iRow];
670              if (gap+value*distance<0.0) {
671                // If integer then has to move by 1
672                if (!isInteger)
673                  distance = -gap/value;
674                else
675                  distance = CoinMax(-gap/value,1.0);
676              }
677            } else if (rowActivity[iRow]<rowLower[iRow]-primalTolerance) {
678              // infeasible below
679              assert (value>0.0);
680              double gap = rowActivity[iRow]-rowLower[iRow];
681              if (gap+value*distance>0.0) {
682                // If integer then has to move by 1
683                if (!isInteger)
684                  distance = -gap/value;
685                else
686                  distance = CoinMax(-gap/value,1.0);
687              }
688            } else {
689              // feasible
690              if (value>0) {
691                double gap = rowActivity[iRow]-rowUpper[iRow];
692                if (gap+value*distance>0.0) 
693                  distance = -gap/value;
694              } else {
695                double gap = rowActivity[iRow]-rowLower[iRow];
696                if (gap+value*distance<0.0) 
697                  distance = -gap/value;
698              }
699            }
700          }
701          if (isInteger)
702            distance = floor(distance+1.05e-8);
703          if (!distance) {
704            // should never happen
705            //printf("zero distance in CbcRounding - debug\n");
706          }
707          //move
708          lastChange += improvement*distance;
709          distance *= way;
710          newSolution[iColumn] += distance;
711          newSolutionValue += direction*objective[iColumn]*distance;
712          for (j=columnStart[iColumn];
713               j<columnStart[iColumn]+columnLength[iColumn];j++) {
714            int iRow = row[j];
715            double value = element[j];
716            rowActivity[iRow] += distance*value;
717          }
718        }
719      }
720      penaltyChange += lastChange;
721    }
722    penalty -= penaltyChange;
723    if (penalty<1.0e-5*fabs(penaltyChange)) {
724      // recompute
725      penalty=0.0;
726      for (i=0;i<numberRows;i++) {
727        double value = rowActivity[i];
728        if (value<rowLower[i]-primalTolerance)
729          penalty += rowLower[i]-value;
730        else if (value>rowUpper[i]+primalTolerance)
731          penalty += value-rowUpper[i];
732      }
733    }
734  }
735
736  // Could also set SOS (using random) and repeat
737  if (!penalty) {
738    // See if we can do better
739    //seed_++;
740    //CoinSeedRandom(seed_);
741    // Random number between 0 and 1.
742    double randomNumber = CoinDrand48();
743    int iPass;
744    int start[2];
745    int end[2];
746    int iRandom = (int) (randomNumber*((double) numberIntegers));
747    start[0]=iRandom;
748    end[0]=numberIntegers;
749    start[1]=0;
750    end[1]=iRandom;
751    for (iPass=0;iPass<2;iPass++) {
752      int i;
753      for (i=start[iPass];i<end[iPass];i++) {
754        int iColumn = integerVariable[i];
755#ifndef NDEBUG
756        double value=newSolution[iColumn];
757        assert (fabs(floor(value+0.5)-value)<integerTolerance);
758#endif
759        double cost = direction * objective[iColumn];
760        double move=0.0;
761        if (cost>0.0)
762          move = -1.0;
763        else if (cost<0.0)
764          move=1.0;
765        while (move) {
766          bool good=true;
767          double newValue=newSolution[iColumn]+move;
768          if (newValue<lower[iColumn]-primalTolerance||
769              newValue>upper[iColumn]+primalTolerance) {
770            move=0.0;
771          } else {
772            // see if we can move
773            int j;
774            for (j=columnStart[iColumn];
775                 j<columnStart[iColumn]+columnLength[iColumn];j++) {
776              int iRow = row[j];
777              double newActivity = rowActivity[iRow] + move*element[j];
778              if (newActivity<rowLower[iRow]-primalTolerance||
779                  newActivity>rowUpper[iRow]+primalTolerance) {
780                good=false;
781                break;
782              }
783            }
784            if (good) {
785              newSolution[iColumn] = newValue;
786              newSolutionValue += move*cost;
787              int j;
788              for (j=columnStart[iColumn];
789                   j<columnStart[iColumn]+columnLength[iColumn];j++) {
790                int iRow = row[j];
791                rowActivity[iRow] += move*element[j];
792              }
793            } else {
794              move=0.0;
795            }
796          }
797        }
798      }
799    }
800    if (newSolutionValue<solutionValue) {
801      // paranoid check
802      memset(rowActivity,0,numberRows*sizeof(double));
803      for (i=0;i<numberColumns;i++) {
804        int j;
805        double value = newSolution[i];
806        if (value) {
807          for (j=columnStart[i];
808               j<columnStart[i]+columnLength[i];j++) {
809            int iRow=row[j];
810            rowActivity[iRow] += value*element[j];
811          }
812        }
813      }
814      // check was approximately feasible
815      bool feasible=true;
816      for (i=0;i<numberRows;i++) {
817        if(rowActivity[i]<rowLower[i]) {
818          if (rowActivity[i]<rowLower[i]-1000.0*primalTolerance)
819            feasible = false;
820        } else if(rowActivity[i]>rowUpper[i]) {
821          if (rowActivity[i]>rowUpper[i]+1000.0*primalTolerance)
822            feasible = false;
823        }
824      }
825      if (feasible) {
826        // new solution
827        memcpy(betterSolution,newSolution,numberColumns*sizeof(double));
828        solutionValue = newSolutionValue;
829        //printf("** Solution of %g found by rounding\n",newSolutionValue);
830        returnCode=1;
831      } else {
832        // Can easily happen
833        //printf("Debug CbcRounding giving bad solution\n");
834      }
835    }
836  }
837  delete [] newSolution;
838  delete [] rowActivity;
839  return returnCode;
840}
841// update model
842void CbcRounding::setModel(CbcModel * model)
843{
844  model_ = model;
845  // Get a copy of original matrix (and by row for rounding);
846  assert(model_->solver());
847  matrix_ = *model_->solver()->getMatrixByCol();
848  matrixByRow_ = *model_->solver()->getMatrixByRow();
849  // make sure model okay for heuristic
850  validate();
851}
852// Validate model i.e. sets when_ to 0 if necessary (may be NULL)
853void 
854CbcRounding::validate() 
855{
856  if (model_&&when()<10) {
857    if (model_->numberIntegers()!=
858        model_->numberObjects())
859      setWhen(0);
860  }
861}
862
863// Default Constructor
864CbcSerendipity::CbcSerendipity() 
865  :CbcHeuristic()
866{
867}
868
869// Constructor from model
870CbcSerendipity::CbcSerendipity(CbcModel & model)
871  :CbcHeuristic(model)
872{
873}
874
875// Destructor
876CbcSerendipity::~CbcSerendipity ()
877{
878}
879
880// Clone
881CbcHeuristic *
882CbcSerendipity::clone() const
883{
884  return new CbcSerendipity(*this);
885}
886
887// Copy constructor
888CbcSerendipity::CbcSerendipity(const CbcSerendipity & rhs)
889:
890  CbcHeuristic(rhs)
891{
892}
893
894// Returns 1 if solution, 0 if not
895int
896CbcSerendipity::solution(double & solutionValue,
897                         double * betterSolution)
898{
899  if (!model_)
900    return 0;
901  // get information on solver type
902  OsiAuxInfo * auxInfo = model_->solver()->getAuxiliaryInfo();
903  OsiBabSolver * auxiliaryInfo = dynamic_cast< OsiBabSolver *> (auxInfo);
904  if (auxiliaryInfo)
905    return auxiliaryInfo->solution(solutionValue,betterSolution,model_->solver()->getNumCols());
906  else
907    return 0;
908}
909// update model
910void CbcSerendipity::setModel(CbcModel * model)
911{
912  model_ = model;
913}
914// Resets stuff if model changes
915void 
916CbcSerendipity::resetModel(CbcModel * model)
917{
918  model_ = model;
919}
920 
Note: See TracBrowser for help on using the repository browser.