source: trunk/CbcHeuristic.cpp @ 216

Last change on this file since 216 was 208, checked in by forrest, 15 years ago

large number of changes

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