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

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

add multi heuristic

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