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

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

Changes to try and make faster.
ALSO createBranch is now createCbcBranch. I don't think anyone uses derived classes from CbcBranchBase? but if they do they will get annoyed at me and have to change a few names. I am just trying to get rid of all but all warnings

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