source: trunk/Cbc/src/CbcHeuristicGreedy.cpp @ 1212

Last change on this file since 1212 was 1212, checked in by forrest, 10 years ago

fixes

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 26.8 KB
Line 
1/* $Id: CbcHeuristicGreedy.cpp 1212 2009-08-21 16:19:13Z forrest $ */
2// Copyright (C) 2005, International Business Machines
3// Corporation and others.  All Rights Reserved.
4#if defined(_MSC_VER)
5// Turn off compiler warning about long names
6#  pragma warning(disable:4786)
7#endif
8#include <cassert>
9#include <cstdlib>
10#include <cmath>
11#include <cfloat>
12
13#include "OsiSolverInterface.hpp"
14#include "CbcModel.hpp"
15#include "CbcStrategy.hpp"
16#include "CbcHeuristicGreedy.hpp"
17#include "CoinSort.hpp"
18#include "CglPreProcess.hpp"
19// Default Constructor
20CbcHeuristicGreedyCover::CbcHeuristicGreedyCover() 
21  :CbcHeuristic()
22{
23  // matrix  will automatically be empty
24  originalNumberRows_=0;
25  algorithm_=0;
26  numberTimes_=100;
27}
28
29// Constructor from model
30CbcHeuristicGreedyCover::CbcHeuristicGreedyCover(CbcModel & model)
31  :CbcHeuristic(model)
32{
33  gutsOfConstructor(&model);
34  algorithm_=0;
35  numberTimes_=100;
36  whereFrom_=1;
37}
38
39// Destructor
40CbcHeuristicGreedyCover::~CbcHeuristicGreedyCover ()
41{
42}
43
44// Clone
45CbcHeuristic *
46CbcHeuristicGreedyCover::clone() const
47{
48  return new CbcHeuristicGreedyCover(*this);
49}
50// Guts of constructor from a CbcModel
51void 
52CbcHeuristicGreedyCover::gutsOfConstructor(CbcModel * model)
53{
54  model_=model;
55  // Get a copy of original matrix
56  assert(model->solver());
57  if (model->solver()->getNumRows()) {
58    matrix_ = *model->solver()->getMatrixByCol();
59  }
60  originalNumberRows_=model->solver()->getNumRows();
61}
62// Create C++ lines to get to current state
63void 
64CbcHeuristicGreedyCover::generateCpp( FILE * fp) 
65{
66  CbcHeuristicGreedyCover other;
67  fprintf(fp,"0#include \"CbcHeuristicGreedy.hpp\"\n");
68  fprintf(fp,"3  CbcHeuristicGreedyCover heuristicGreedyCover(*cbcModel);\n");
69  CbcHeuristic::generateCpp(fp,"heuristicGreedyCover");
70  if (algorithm_!=other.algorithm_)
71    fprintf(fp,"3  heuristicGreedyCover.setAlgorithm(%d);\n",algorithm_);
72  else
73    fprintf(fp,"4  heuristicGreedyCover.setAlgorithm(%d);\n",algorithm_);
74  if (numberTimes_!=other.numberTimes_)
75    fprintf(fp,"3  heuristicGreedyCover.setNumberTimes(%d);\n",numberTimes_);
76  else
77    fprintf(fp,"4  heuristicGreedyCover.setNumberTimes(%d);\n",numberTimes_);
78  fprintf(fp,"3  cbcModel->addHeuristic(&heuristicGreedyCover);\n");
79}
80
81// Copy constructor
82CbcHeuristicGreedyCover::CbcHeuristicGreedyCover(const CbcHeuristicGreedyCover & rhs)
83:
84  CbcHeuristic(rhs),
85  matrix_(rhs.matrix_),
86  originalNumberRows_(rhs.originalNumberRows_),
87  algorithm_(rhs.algorithm_),
88  numberTimes_(rhs.numberTimes_)
89{
90}
91
92// Assignment operator
93CbcHeuristicGreedyCover & 
94CbcHeuristicGreedyCover::operator=( const CbcHeuristicGreedyCover& rhs)
95{
96  if (this != &rhs) {
97    CbcHeuristic::operator=(rhs);
98    matrix_=rhs.matrix_;
99    originalNumberRows_=rhs.originalNumberRows_;
100    algorithm_=rhs.algorithm_;
101    numberTimes_=rhs.numberTimes_;
102  }
103  return *this;
104}
105// Returns 1 if solution, 0 if not
106int
107CbcHeuristicGreedyCover::solution(double & solutionValue,
108                         double * betterSolution)
109{
110  numCouldRun_++;
111  if (!model_)
112    return 0;
113  // See if to do
114  if (!when()||(when()==1&&model_->phase()!=1))
115    return 0; // switched off
116  if (model_->getNodeCount()>numberTimes_)
117    return 0;
118  // See if at root node
119  bool atRoot = model_->getNodeCount()==0;
120  int passNumber = model_->getCurrentPassNumber();
121  if (atRoot&&passNumber!=1)
122    return 0;
123  OsiSolverInterface * solver = model_->solver();
124  const double * columnLower = solver->getColLower();
125  const double * columnUpper = solver->getColUpper();
126  // And original upper bounds in case we want to use them
127  const double * originalUpper = model_->continuousSolver()->getColUpper();
128  // But not if algorithm says so
129  if ((algorithm_%10)==0)
130    originalUpper = columnUpper;
131  const double * rowLower = solver->getRowLower();
132  const double * solution = solver->getColSolution();
133  const double * objective = solver->getObjCoefficients();
134  double integerTolerance = model_->getDblParam(CbcModel::CbcIntegerTolerance);
135  double primalTolerance;
136  solver->getDblParam(OsiPrimalTolerance,primalTolerance);
137
138  // This is number of rows when matrix was passed in
139  int numberRows = originalNumberRows_;
140  if (!numberRows)
141    return 0; // switched off
142
143  numRuns_++;
144  assert (numberRows==matrix_.getNumRows());
145  int iRow, iColumn;
146  double direction = solver->getObjSense();
147  double offset;
148  solver->getDblParam(OsiObjOffset,offset);
149  double newSolutionValue = -offset;
150  int returnCode = 0;
151
152  // Column copy
153  const double * element = matrix_.getElements();
154  const int * row = matrix_.getIndices();
155  const CoinBigIndex * columnStart = matrix_.getVectorStarts();
156  const int * columnLength = matrix_.getVectorLengths();
157
158  // Get solution array for heuristic solution
159  int numberColumns = solver->getNumCols();
160  double * newSolution = new double [numberColumns];
161  double * rowActivity = new double[numberRows];
162  memset(rowActivity,0,numberRows*sizeof(double));
163  bool allOnes=true;
164  // Get rounded down solution
165  for (iColumn=0;iColumn<numberColumns;iColumn++) {
166    CoinBigIndex j;
167    double value = solution[iColumn];
168    if (solver->isInteger(iColumn)) {
169      // Round down integer
170      if (fabs(floor(value+0.5)-value)<integerTolerance) {
171        value=floor(CoinMax(value+1.0e-3,columnLower[iColumn]));
172      } else {
173        value=CoinMax(floor(value),columnLower[iColumn]);
174      }
175    }
176    // make sure clean
177    value = CoinMin(value,columnUpper[iColumn]);
178    value = CoinMax(value,columnLower[iColumn]);
179    newSolution[iColumn]=value;
180    double cost = direction * objective[iColumn];
181    newSolutionValue += value*cost;
182    for (j=columnStart[iColumn];
183         j<columnStart[iColumn]+columnLength[iColumn];j++) {
184      int iRow=row[j];
185      rowActivity[iRow] += value*element[j];
186      if (element[j]!=1.0)
187        allOnes=false;
188    }
189  }
190  // See if we round up
191  bool roundup = ((algorithm_%100)!=0);
192  if (roundup&&allOnes) {
193    // Get rounded up solution
194    for (iColumn=0;iColumn<numberColumns;iColumn++) {
195      CoinBigIndex j;
196      double value = solution[iColumn];
197      if (solver->isInteger(iColumn)) {
198        // but round up if no activity
199        if (roundup&&value>=0.499999&&!newSolution[iColumn]) {
200          bool choose=true;
201          for (j=columnStart[iColumn];
202               j<columnStart[iColumn]+columnLength[iColumn];j++) {
203            int iRow=row[j];
204            if (rowActivity[iRow]) {
205              choose=false;
206              break;
207            }
208          }
209          if (choose) {
210            newSolution[iColumn]=1.0;
211            double cost = direction * objective[iColumn];
212            newSolutionValue += cost;
213            for (j=columnStart[iColumn];
214                 j<columnStart[iColumn]+columnLength[iColumn];j++) {
215              int iRow=row[j];
216              rowActivity[iRow] += 1.0;
217            }
218          }
219        }
220      }
221    }
222  }
223  // Get initial list
224  int * which = new int [numberColumns];
225  for (iColumn=0;iColumn<numberColumns;iColumn++) 
226    which[iColumn]=iColumn;
227  int numberLook=numberColumns;
228  // See if we want to perturb more
229  double perturb = ((algorithm_%10)==0) ? 0.1 : 0.25;
230  // Keep going round until a solution
231  while (true) {
232    // Get column with best ratio
233    int bestColumn=-1;
234    double bestRatio=COIN_DBL_MAX;
235    double bestStepSize = 0.0;
236    int newNumber=0;
237    for (int jColumn=0;jColumn<numberLook;jColumn++) {
238      int iColumn = which[jColumn];
239      CoinBigIndex j;
240      double value = newSolution[iColumn];
241      double cost = direction * objective[iColumn];
242      if (solver->isInteger(iColumn)) {
243        // use current upper or original upper
244        if (value+0.99<originalUpper[iColumn]) {
245          double sum=0.0;
246          int numberExact=0;
247          for (j=columnStart[iColumn];
248               j<columnStart[iColumn]+columnLength[iColumn];j++) {
249            int iRow=row[j];
250            double gap = rowLower[iRow]-rowActivity[iRow];
251            double elementValue = allOnes ? 1.0 : element[j];
252            if (gap>1.0e-7) {
253              sum += CoinMin(elementValue,gap);
254              if (fabs(elementValue-gap)<1.0e-7) 
255                numberExact++;
256            }
257          }
258          // could bias if exact
259          if (sum>0.0) {
260            // add to next time
261            which[newNumber++]=iColumn;
262            double ratio = (cost/sum)*(1.0+perturb*randomNumberGenerator_.randomDouble());
263            // If at root choose first
264            if (atRoot)
265              ratio = iColumn;
266            if (ratio<bestRatio) {
267              bestRatio=ratio;
268              bestColumn=iColumn;
269              bestStepSize=1.0;
270            }
271          }
272        }
273      } else {
274        // continuous
275        if (value<columnUpper[iColumn]) {
276          // Go through twice - first to get step length
277          double step=1.0e50;
278          for (j=columnStart[iColumn];
279               j<columnStart[iColumn]+columnLength[iColumn];j++) {
280            int iRow=row[j];
281            if (rowActivity[iRow]<rowLower[iRow]-1.0e-10&&
282                element[j]*step+rowActivity[iRow]>=rowLower[iRow]) {
283              step = (rowLower[iRow]-rowActivity[iRow])/element[j];;
284            }
285          }
286          // now ratio
287          if (step<1.0e50) {
288            // add to next time
289            which[newNumber++]=iColumn;
290            assert (step>0.0);
291            double sum=0.0;
292            for (j=columnStart[iColumn];
293                 j<columnStart[iColumn]+columnLength[iColumn];j++) {
294              int iRow=row[j];
295              double newActivity = element[j]*step+rowActivity[iRow];
296              if (rowActivity[iRow]<rowLower[iRow]-1.0e-10&&
297                newActivity>=rowLower[iRow]-1.0e-12) {
298                sum += element[j];
299              }
300            }
301            assert (sum>0.0);
302            double ratio = (cost/sum)*(1.0+perturb*randomNumberGenerator_.randomDouble());
303            if (ratio<bestRatio) {
304              bestRatio=ratio;
305              bestColumn=iColumn;
306              bestStepSize=step;
307            }
308          }
309        }
310      }
311    }
312    if (bestColumn<0)
313      break; // we have finished
314    // Increase chosen column
315    newSolution[bestColumn] += bestStepSize;
316    double cost = direction * objective[bestColumn];
317    newSolutionValue += bestStepSize*cost;
318    for (CoinBigIndex j=columnStart[bestColumn];
319         j<columnStart[bestColumn]+columnLength[bestColumn];j++) {
320      int iRow = row[j];
321      rowActivity[iRow] += bestStepSize*element[j];
322    }
323  }
324  delete [] which;
325  if (newSolutionValue<solutionValue) {
326    // check feasible
327    memset(rowActivity,0,numberRows*sizeof(double));
328    for (iColumn=0;iColumn<numberColumns;iColumn++) {
329      CoinBigIndex j;
330      double value = newSolution[iColumn];
331      if (value) {
332        for (j=columnStart[iColumn];
333             j<columnStart[iColumn]+columnLength[iColumn];j++) {
334          int iRow=row[j];
335          rowActivity[iRow] += value*element[j];
336        }
337      }
338    }
339    // check was approximately feasible
340    bool feasible=true;
341    for (iRow=0;iRow<numberRows;iRow++) {
342      if(rowActivity[iRow]<rowLower[iRow]) {
343        if (rowActivity[iRow]<rowLower[iRow]-10.0*primalTolerance)
344          feasible = false;
345      }
346    }
347    if (feasible) {
348      // new solution
349      memcpy(betterSolution,newSolution,numberColumns*sizeof(double));
350      solutionValue = newSolutionValue;
351      //printf("** Solution of %g found by rounding\n",newSolutionValue);
352      returnCode=1;
353    } else {
354      // Can easily happen
355      //printf("Debug CbcHeuristicGreedyCover giving bad solution\n");
356    }
357  }
358  delete [] newSolution;
359  delete [] rowActivity;
360  return returnCode;
361}
362// update model
363void CbcHeuristicGreedyCover::setModel(CbcModel * model)
364{
365  gutsOfConstructor(model);
366  validate();
367}
368// Resets stuff if model changes
369void 
370CbcHeuristicGreedyCover::resetModel(CbcModel * model)
371{
372  gutsOfConstructor(model);
373}
374// Validate model i.e. sets when_ to 0 if necessary (may be NULL)
375void 
376CbcHeuristicGreedyCover::validate() 
377{
378  if (model_&&when()<10) {
379    if (model_->numberIntegers()!=
380        model_->numberObjects()&&(model_->numberObjects()||
381                                  (model_->specialOptions()&1024)==0)) {
382      int numberOdd=0;
383      for (int i=0;i<model_->numberObjects();i++) {
384        if (!model_->object(i)->canDoHeuristics()) 
385          numberOdd++;
386      }
387      if (numberOdd)
388        setWhen(0);
389    }
390    // Only works if costs positive, coefficients positive and all rows G
391    OsiSolverInterface * solver = model_->solver();
392    const double * columnLower = solver->getColLower();
393    const double * rowUpper = solver->getRowUpper();
394    const double * objective = solver->getObjCoefficients();
395    double direction = solver->getObjSense();
396
397    int numberRows = solver->getNumRows();
398    // Column copy
399    const double * element = matrix_.getElements();
400    const CoinBigIndex * columnStart = matrix_.getVectorStarts();
401    const int * columnLength = matrix_.getVectorLengths();
402    bool good = true;
403    for (int iRow=0;iRow<numberRows;iRow++) {
404      if (rowUpper[iRow]<1.0e30)
405        good = false;
406    }
407    int numberColumns = solver->getNumCols();
408    for (int iColumn=0;iColumn<numberColumns;iColumn++) {
409      if (objective[iColumn]*direction<0.0)
410        good=false;
411      if (columnLower[iColumn]<0.0)
412        good=false;
413      CoinBigIndex j;
414      for (j=columnStart[iColumn];
415           j<columnStart[iColumn]+columnLength[iColumn];j++) {
416        if (element[j]<0.0)
417          good=false;
418      }
419    }
420    if (!good)
421      setWhen(0); // switch off
422  }
423}
424// Default Constructor
425CbcHeuristicGreedyEquality::CbcHeuristicGreedyEquality() 
426  :CbcHeuristic()
427{
428  // matrix  will automatically be empty
429  fraction_=1.0; // no branch and bound
430  originalNumberRows_=0;
431  algorithm_=0;
432  numberTimes_=100;
433  whereFrom_=1;
434}
435
436// Constructor from model
437CbcHeuristicGreedyEquality::CbcHeuristicGreedyEquality(CbcModel & model)
438  :CbcHeuristic(model)
439{
440  // Get a copy of original matrix
441  gutsOfConstructor(&model);
442  fraction_=1.0; // no branch and bound
443  algorithm_=0;
444  numberTimes_=100;
445  whereFrom_=1;
446}
447
448// Destructor
449CbcHeuristicGreedyEquality::~CbcHeuristicGreedyEquality ()
450{
451}
452
453// Clone
454CbcHeuristic *
455CbcHeuristicGreedyEquality::clone() const
456{
457  return new CbcHeuristicGreedyEquality(*this);
458}
459// Guts of constructor from a CbcModel
460void 
461CbcHeuristicGreedyEquality::gutsOfConstructor(CbcModel * model)
462{
463  model_=model;
464  // Get a copy of original matrix
465  assert(model->solver());
466  if (model->solver()->getNumRows()) {
467    matrix_ = *model->solver()->getMatrixByCol();
468  }
469  originalNumberRows_=model->solver()->getNumRows();
470}
471// Create C++ lines to get to current state
472void 
473CbcHeuristicGreedyEquality::generateCpp( FILE * fp) 
474{
475  CbcHeuristicGreedyEquality other;
476  fprintf(fp,"0#include \"CbcHeuristicGreedy.hpp\"\n");
477  fprintf(fp,"3  CbcHeuristicGreedyEquality heuristicGreedyEquality(*cbcModel);\n");
478  CbcHeuristic::generateCpp(fp,"heuristicGreedyEquality");
479  if (algorithm_!=other.algorithm_)
480    fprintf(fp,"3  heuristicGreedyEquality.setAlgorithm(%d);\n",algorithm_);
481  else
482    fprintf(fp,"4  heuristicGreedyEquality.setAlgorithm(%d);\n",algorithm_);
483  if (fraction_!=other.fraction_)
484    fprintf(fp,"3  heuristicGreedyEquality.setFraction(%g);\n",fraction_);
485  else
486    fprintf(fp,"4  heuristicGreedyEquality.setFraction(%g);\n",fraction_);
487  if (numberTimes_!=other.numberTimes_)
488    fprintf(fp,"3  heuristicGreedyEquality.setNumberTimes(%d);\n",numberTimes_);
489  else
490    fprintf(fp,"4  heuristicGreedyEquality.setNumberTimes(%d);\n",numberTimes_);
491  fprintf(fp,"3  cbcModel->addHeuristic(&heuristicGreedyEquality);\n");
492}
493
494// Copy constructor
495CbcHeuristicGreedyEquality::CbcHeuristicGreedyEquality(const CbcHeuristicGreedyEquality & rhs)
496:
497  CbcHeuristic(rhs),
498  matrix_(rhs.matrix_),
499  fraction_(rhs.fraction_),
500  originalNumberRows_(rhs.originalNumberRows_),
501  algorithm_(rhs.algorithm_),
502  numberTimes_(rhs.numberTimes_)
503{
504}
505
506// Assignment operator
507CbcHeuristicGreedyEquality & 
508CbcHeuristicGreedyEquality::operator=( const CbcHeuristicGreedyEquality& rhs)
509{
510  if (this != &rhs) {
511    CbcHeuristic::operator=(rhs);
512    matrix_=rhs.matrix_;
513    fraction_ = rhs.fraction_;
514    originalNumberRows_=rhs.originalNumberRows_;
515    algorithm_=rhs.algorithm_;
516    numberTimes_=rhs.numberTimes_;
517  }
518  return *this;
519}
520// Returns 1 if solution, 0 if not
521int
522CbcHeuristicGreedyEquality::solution(double & solutionValue,
523                         double * betterSolution)
524{
525  numCouldRun_++;
526  if (!model_)
527    return 0;
528  // See if to do
529  if (!when()||(when()==1&&model_->phase()!=1))
530    return 0; // switched off
531  if (model_->getNodeCount()>numberTimes_)
532    return 0;
533  // See if at root node
534  bool atRoot = model_->getNodeCount()==0;
535  int passNumber = model_->getCurrentPassNumber();
536  if (atRoot&&passNumber!=1)
537    return 0;
538  OsiSolverInterface * solver = model_->solver();
539  const double * columnLower = solver->getColLower();
540  const double * columnUpper = solver->getColUpper();
541  // And original upper bounds in case we want to use them
542  const double * originalUpper = model_->continuousSolver()->getColUpper();
543  // But not if algorithm says so
544  if ((algorithm_%10)==0)
545    originalUpper = columnUpper;
546  const double * rowLower = solver->getRowLower();
547  const double * rowUpper = solver->getRowUpper();
548  const double * solution = solver->getColSolution();
549  const double * objective = solver->getObjCoefficients();
550  double integerTolerance = model_->getDblParam(CbcModel::CbcIntegerTolerance);
551  double primalTolerance;
552  solver->getDblParam(OsiPrimalTolerance,primalTolerance);
553
554  // This is number of rows when matrix was passed in
555  int numberRows = originalNumberRows_;
556  if (!numberRows)
557    return 0; // switched off
558  numRuns_++;
559
560  assert (numberRows==matrix_.getNumRows());
561  int iRow, iColumn;
562  double direction = solver->getObjSense();
563  double offset;
564  solver->getDblParam(OsiObjOffset,offset);
565  double newSolutionValue = -offset;
566  int returnCode = 0;
567
568  // Column copy
569  const double * element = matrix_.getElements();
570  const int * row = matrix_.getIndices();
571  const CoinBigIndex * columnStart = matrix_.getVectorStarts();
572  const int * columnLength = matrix_.getVectorLengths();
573
574  // Get solution array for heuristic solution
575  int numberColumns = solver->getNumCols();
576  double * newSolution = new double [numberColumns];
577  double * rowActivity = new double[numberRows];
578  memset(rowActivity,0,numberRows*sizeof(double));
579  double rhsNeeded=0;
580  for (iRow=0;iRow<numberRows;iRow++) 
581    rhsNeeded += rowUpper[iRow];
582  rhsNeeded *= fraction_;
583  bool allOnes=true;
584  // Get rounded down solution
585  for (iColumn=0;iColumn<numberColumns;iColumn++) {
586    CoinBigIndex j;
587    double value = solution[iColumn];
588    if (solver->isInteger(iColumn)) {
589      // Round down integer
590      if (fabs(floor(value+0.5)-value)<integerTolerance) {
591        value=floor(CoinMax(value+1.0e-3,columnLower[iColumn]));
592      } else {
593        value=CoinMax(floor(value),columnLower[iColumn]);
594      }
595    }
596    // make sure clean
597    value = CoinMin(value,columnUpper[iColumn]);
598    value = CoinMax(value,columnLower[iColumn]);
599    newSolution[iColumn]=value;
600    double cost = direction * objective[iColumn];
601    newSolutionValue += value*cost;
602    for (j=columnStart[iColumn];
603         j<columnStart[iColumn]+columnLength[iColumn];j++) {
604      int iRow=row[j];
605      rowActivity[iRow] += value*element[j];
606      rhsNeeded -= value*element[j];
607      if (element[j]!=1.0)
608        allOnes=false;
609    }
610  }
611  // See if we round up
612  bool roundup = ((algorithm_%100)!=0);
613  if (roundup&&allOnes) {
614    // Get rounded up solution
615    for (iColumn=0;iColumn<numberColumns;iColumn++) {
616      CoinBigIndex j;
617      double value = solution[iColumn];
618      if (solver->isInteger(iColumn)) {
619        // but round up if no activity
620        if (roundup&&value>=0.6&&!newSolution[iColumn]) {
621          bool choose=true;
622          for (j=columnStart[iColumn];
623               j<columnStart[iColumn]+columnLength[iColumn];j++) {
624            int iRow=row[j];
625            if (rowActivity[iRow]) {
626              choose=false;
627              break;
628            }
629          }
630          if (choose) {
631            newSolution[iColumn]=1.0;
632            double cost = direction * objective[iColumn];
633            newSolutionValue += cost;
634            for (j=columnStart[iColumn];
635                 j<columnStart[iColumn]+columnLength[iColumn];j++) {
636              int iRow=row[j];
637              rowActivity[iRow] += 1.0;
638              rhsNeeded -= 1.0;
639            }
640          }
641        }
642      }
643    }
644  }
645  // Get initial list
646  int * which = new int [numberColumns];
647  for (iColumn=0;iColumn<numberColumns;iColumn++) 
648    which[iColumn]=iColumn;
649  int numberLook=numberColumns;
650  // See if we want to perturb more
651  double perturb = ((algorithm_%10)==0) ? 0.1 : 0.25;
652  // Keep going round until a solution
653  while (true) {
654    // Get column with best ratio
655    int bestColumn=-1;
656    double bestRatio=COIN_DBL_MAX;
657    double bestStepSize = 0.0;
658    int newNumber=0;
659    for (int jColumn=0;jColumn<numberLook;jColumn++) {
660      int iColumn = which[jColumn];
661      CoinBigIndex j;
662      double value = newSolution[iColumn];
663      double cost = direction * objective[iColumn];
664      if (solver->isInteger(iColumn)) {
665        // use current upper or original upper
666        if (value+0.9999<originalUpper[iColumn]) {
667          double movement = 1.0;
668          double sum=0.0;
669          for (j=columnStart[iColumn];
670               j<columnStart[iColumn]+columnLength[iColumn];j++) {
671            int iRow=row[j];
672            double gap = rowUpper[iRow]-rowActivity[iRow];
673            double elementValue = allOnes ? 1.0 : element[j];
674            sum += elementValue;
675            if (movement*elementValue>gap) {
676              movement = gap/elementValue;
677            }
678          }
679          if (movement>0.999999) {
680            // add to next time
681            which[newNumber++]=iColumn;
682            double ratio = (cost/sum)*(1.0+perturb*randomNumberGenerator_.randomDouble());
683            // If at root
684            if (atRoot) {
685              if (fraction_==1.0)
686                ratio = iColumn; // choose first
687              else
688                ratio = - solution[iColumn]; // choose largest
689            }
690            if (ratio<bestRatio) {
691              bestRatio=ratio;
692              bestColumn=iColumn;
693              bestStepSize=1.0;
694            }
695          }
696        }
697      } else {
698        // continuous
699        if (value<columnUpper[iColumn]) {
700          double movement=1.0e50;
701          double sum=0.0;
702          for (j=columnStart[iColumn];
703               j<columnStart[iColumn]+columnLength[iColumn];j++) {
704            int iRow=row[j];
705            if (element[j]*movement+rowActivity[iRow]>rowUpper[iRow]) {
706              movement = (rowUpper[iRow]-rowActivity[iRow])/element[j];;
707            }
708            sum += element[j];
709          }
710          // now ratio
711          if (movement>1.0e-7) {
712            // add to next time
713            which[newNumber++]=iColumn;
714            double ratio = (cost/sum)*(1.0+perturb*randomNumberGenerator_.randomDouble());
715            if (ratio<bestRatio) {
716              bestRatio=ratio;
717              bestColumn=iColumn;
718              bestStepSize=movement;
719            }
720          }
721        }
722      }
723    }
724    if (bestColumn<0)
725      break; // we have finished
726    // Increase chosen column
727    newSolution[bestColumn] += bestStepSize;
728    double cost = direction * objective[bestColumn];
729    newSolutionValue += bestStepSize*cost;
730    for (CoinBigIndex j=columnStart[bestColumn];
731         j<columnStart[bestColumn]+columnLength[bestColumn];j++) {
732      int iRow = row[j];
733      rowActivity[iRow] += bestStepSize*element[j];
734      rhsNeeded -= bestStepSize*element[j];
735    }
736    if (rhsNeeded<1.0e-8)
737      break;
738  }
739  delete [] which;
740  if (fraction_<1.0&&rhsNeeded<1.0e-8&&newSolutionValue<solutionValue) {
741    // do branch and cut
742    // fix all nonzero
743    OsiSolverInterface * newSolver = model_->continuousSolver()->clone();
744    for (iColumn=0;iColumn<numberColumns;iColumn++) {
745      if (newSolver->isInteger(iColumn))
746        newSolver->setColLower(iColumn,newSolution[iColumn]);
747    }
748    int returnCode = smallBranchAndBound(newSolver,200,newSolution,newSolutionValue,
749                                         solutionValue,"CbcHeuristicGreedy");
750    if (returnCode<0)
751      returnCode=0; // returned on size
752    if ((returnCode&2)!=0) {
753      // could add cut
754      returnCode &= ~2;
755    }
756    rhsNeeded = 1.0-returnCode;
757    delete newSolver;
758  }
759  if (newSolutionValue<solutionValue&&rhsNeeded<1.0e-8) {
760    // check feasible
761    memset(rowActivity,0,numberRows*sizeof(double));
762    for (iColumn=0;iColumn<numberColumns;iColumn++) {
763      CoinBigIndex j;
764      double value = newSolution[iColumn];
765      if (value) {
766        for (j=columnStart[iColumn];
767             j<columnStart[iColumn]+columnLength[iColumn];j++) {
768          int iRow=row[j];
769          rowActivity[iRow] += value*element[j];
770        }
771      }
772    }
773    // check was approximately feasible
774    bool feasible=true;
775    for (iRow=0;iRow<numberRows;iRow++) {
776      if(rowActivity[iRow]<rowLower[iRow]) {
777        if (rowActivity[iRow]<rowLower[iRow]-10.0*primalTolerance)
778          feasible = false;
779      }
780    }
781    if (feasible) {
782      // new solution
783      memcpy(betterSolution,newSolution,numberColumns*sizeof(double));
784      solutionValue = newSolutionValue;
785      returnCode=1;
786    }
787  }
788  delete [] newSolution;
789  delete [] rowActivity;
790  if (atRoot&&fraction_==1.0) {
791    // try quick search
792    fraction_=0.4;
793    int newCode=this->solution(solutionValue,betterSolution);
794    if (newCode)
795      returnCode=1;
796    fraction_=1.0;
797  }
798  return returnCode;
799}
800// update model
801void CbcHeuristicGreedyEquality::setModel(CbcModel * model)
802{
803  gutsOfConstructor(model);
804  validate();
805}
806// Resets stuff if model changes
807void 
808CbcHeuristicGreedyEquality::resetModel(CbcModel * model)
809{
810  gutsOfConstructor(model);
811}
812// Validate model i.e. sets when_ to 0 if necessary (may be NULL)
813void 
814CbcHeuristicGreedyEquality::validate() 
815{
816  if (model_&&when()<10) {
817    if (model_->numberIntegers()!=
818        model_->numberObjects())
819      setWhen(0);
820    // Only works if costs positive, coefficients positive and all rows E or L
821    // And if values are integer
822    OsiSolverInterface * solver = model_->solver();
823    const double * columnLower = solver->getColLower();
824    const double * rowUpper = solver->getRowUpper();
825    const double * rowLower = solver->getRowLower();
826    const double * objective = solver->getObjCoefficients();
827    double direction = solver->getObjSense();
828
829    int numberRows = solver->getNumRows();
830    // Column copy
831    const double * element = matrix_.getElements();
832    const CoinBigIndex * columnStart = matrix_.getVectorStarts();
833    const int * columnLength = matrix_.getVectorLengths();
834    bool good = true;
835    for (int iRow=0;iRow<numberRows;iRow++) {
836      if (rowUpper[iRow]>1.0e30)
837        good = false;
838      if (rowLower[iRow]>0.0&&rowLower[iRow]!=rowUpper[iRow])
839        good = false;
840      if (floor(rowUpper[iRow]+0.5)!=rowUpper[iRow])
841        good = false;
842    }
843    int numberColumns = solver->getNumCols();
844    for (int iColumn=0;iColumn<numberColumns;iColumn++) {
845      if (objective[iColumn]*direction<0.0)
846        good=false;
847      if (columnLower[iColumn]<0.0)
848        good=false;
849      CoinBigIndex j;
850      for (j=columnStart[iColumn];
851           j<columnStart[iColumn]+columnLength[iColumn];j++) {
852        if (element[j]<0.0)
853          good=false;
854        if (floor(element[j]+0.5)!=element[j])
855          good = false;
856      }
857    }
858    if (!good)
859      setWhen(0); // switch off
860  }
861}
862
863 
Note: See TracBrowser for help on using the repository browser.