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

Last change on this file since 1148 was 1148, checked in by forrest, 11 years ago

changes to use heuristics with SOS etc

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 26.7 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      int numberOdd=0;
379      for (int i=0;i<model_->numberObjects();i++) {
380        if (!model_->object(i)->canDoHeuristics()) 
381          numberOdd++;
382      }
383      if (numberOdd)
384        setWhen(0);
385    }
386    // Only works if costs positive, coefficients positive and all rows G
387    OsiSolverInterface * solver = model_->solver();
388    const double * columnLower = solver->getColLower();
389    const double * rowUpper = solver->getRowUpper();
390    const double * objective = solver->getObjCoefficients();
391    double direction = solver->getObjSense();
392
393    int numberRows = solver->getNumRows();
394    // Column copy
395    const double * element = matrix_.getElements();
396    const CoinBigIndex * columnStart = matrix_.getVectorStarts();
397    const int * columnLength = matrix_.getVectorLengths();
398    bool good = true;
399    for (int iRow=0;iRow<numberRows;iRow++) {
400      if (rowUpper[iRow]<1.0e30)
401        good = false;
402    }
403    int numberColumns = solver->getNumCols();
404    for (int iColumn=0;iColumn<numberColumns;iColumn++) {
405      if (objective[iColumn]*direction<0.0)
406        good=false;
407      if (columnLower[iColumn]<0.0)
408        good=false;
409      CoinBigIndex j;
410      for (j=columnStart[iColumn];
411           j<columnStart[iColumn]+columnLength[iColumn];j++) {
412        if (element[j]<0.0)
413          good=false;
414      }
415    }
416    if (!good)
417      setWhen(0); // switch off
418  }
419}
420// Default Constructor
421CbcHeuristicGreedyEquality::CbcHeuristicGreedyEquality() 
422  :CbcHeuristic()
423{
424  // matrix  will automatically be empty
425  fraction_=1.0; // no branch and bound
426  originalNumberRows_=0;
427  algorithm_=0;
428  numberTimes_=100;
429}
430
431// Constructor from model
432CbcHeuristicGreedyEquality::CbcHeuristicGreedyEquality(CbcModel & model)
433  :CbcHeuristic(model)
434{
435  // Get a copy of original matrix
436  gutsOfConstructor(&model);
437  fraction_=1.0; // no branch and bound
438  algorithm_=0;
439  numberTimes_=100;
440}
441
442// Destructor
443CbcHeuristicGreedyEquality::~CbcHeuristicGreedyEquality ()
444{
445}
446
447// Clone
448CbcHeuristic *
449CbcHeuristicGreedyEquality::clone() const
450{
451  return new CbcHeuristicGreedyEquality(*this);
452}
453// Guts of constructor from a CbcModel
454void 
455CbcHeuristicGreedyEquality::gutsOfConstructor(CbcModel * model)
456{
457  model_=model;
458  // Get a copy of original matrix
459  assert(model->solver());
460  matrix_ = *model->solver()->getMatrixByCol();
461  originalNumberRows_=model->solver()->getNumRows();
462}
463// Create C++ lines to get to current state
464void 
465CbcHeuristicGreedyEquality::generateCpp( FILE * fp) 
466{
467  CbcHeuristicGreedyEquality other;
468  fprintf(fp,"0#include \"CbcHeuristicGreedy.hpp\"\n");
469  fprintf(fp,"3  CbcHeuristicGreedyEquality heuristicGreedyEquality(*cbcModel);\n");
470  CbcHeuristic::generateCpp(fp,"heuristicGreedyEquality");
471  if (algorithm_!=other.algorithm_)
472    fprintf(fp,"3  heuristicGreedyEquality.setAlgorithm(%d);\n",algorithm_);
473  else
474    fprintf(fp,"4  heuristicGreedyEquality.setAlgorithm(%d);\n",algorithm_);
475  if (fraction_!=other.fraction_)
476    fprintf(fp,"3  heuristicGreedyEquality.setFraction(%g);\n",fraction_);
477  else
478    fprintf(fp,"4  heuristicGreedyEquality.setFraction(%g);\n",fraction_);
479  if (numberTimes_!=other.numberTimes_)
480    fprintf(fp,"3  heuristicGreedyEquality.setNumberTimes(%d);\n",numberTimes_);
481  else
482    fprintf(fp,"4  heuristicGreedyEquality.setNumberTimes(%d);\n",numberTimes_);
483  fprintf(fp,"3  cbcModel->addHeuristic(&heuristicGreedyEquality);\n");
484}
485
486// Copy constructor
487CbcHeuristicGreedyEquality::CbcHeuristicGreedyEquality(const CbcHeuristicGreedyEquality & rhs)
488:
489  CbcHeuristic(rhs),
490  matrix_(rhs.matrix_),
491  fraction_(rhs.fraction_),
492  originalNumberRows_(rhs.originalNumberRows_),
493  algorithm_(rhs.algorithm_),
494  numberTimes_(rhs.numberTimes_)
495{
496}
497
498// Assignment operator
499CbcHeuristicGreedyEquality & 
500CbcHeuristicGreedyEquality::operator=( const CbcHeuristicGreedyEquality& rhs)
501{
502  if (this != &rhs) {
503    CbcHeuristic::operator=(rhs);
504    matrix_=rhs.matrix_;
505    fraction_ = rhs.fraction_;
506    originalNumberRows_=rhs.originalNumberRows_;
507    algorithm_=rhs.algorithm_;
508    numberTimes_=rhs.numberTimes_;
509  }
510  return *this;
511}
512// Returns 1 if solution, 0 if not
513int
514CbcHeuristicGreedyEquality::solution(double & solutionValue,
515                         double * betterSolution)
516{
517  numCouldRun_++;
518  if (!model_)
519    return 0;
520  // See if to do
521  if (!when()||(when()==1&&model_->phase()!=1))
522    return 0; // switched off
523  if (model_->getNodeCount()>numberTimes_)
524    return 0;
525  // See if at root node
526  bool atRoot = model_->getNodeCount()==0;
527  int passNumber = model_->getCurrentPassNumber();
528  if (atRoot&&passNumber!=1)
529    return 0;
530  OsiSolverInterface * solver = model_->solver();
531  const double * columnLower = solver->getColLower();
532  const double * columnUpper = solver->getColUpper();
533  // And original upper bounds in case we want to use them
534  const double * originalUpper = model_->continuousSolver()->getColUpper();
535  // But not if algorithm says so
536  if ((algorithm_%10)==0)
537    originalUpper = columnUpper;
538  const double * rowLower = solver->getRowLower();
539  const double * rowUpper = solver->getRowUpper();
540  const double * solution = solver->getColSolution();
541  const double * objective = solver->getObjCoefficients();
542  double integerTolerance = model_->getDblParam(CbcModel::CbcIntegerTolerance);
543  double primalTolerance;
544  solver->getDblParam(OsiPrimalTolerance,primalTolerance);
545
546  // This is number of rows when matrix was passed in
547  int numberRows = originalNumberRows_;
548  if (!numberRows)
549    return 0; // switched off
550  numRuns_++;
551
552  assert (numberRows==matrix_.getNumRows());
553  int iRow, iColumn;
554  double direction = solver->getObjSense();
555  double offset;
556  solver->getDblParam(OsiObjOffset,offset);
557  double newSolutionValue = -offset;
558  int returnCode = 0;
559
560  // Column copy
561  const double * element = matrix_.getElements();
562  const int * row = matrix_.getIndices();
563  const CoinBigIndex * columnStart = matrix_.getVectorStarts();
564  const int * columnLength = matrix_.getVectorLengths();
565
566  // Get solution array for heuristic solution
567  int numberColumns = solver->getNumCols();
568  double * newSolution = new double [numberColumns];
569  double * rowActivity = new double[numberRows];
570  memset(rowActivity,0,numberRows*sizeof(double));
571  double rhsNeeded=0;
572  for (iRow=0;iRow<numberRows;iRow++) 
573    rhsNeeded += rowUpper[iRow];
574  rhsNeeded *= fraction_;
575  bool allOnes=true;
576  // Get rounded down solution
577  for (iColumn=0;iColumn<numberColumns;iColumn++) {
578    CoinBigIndex j;
579    double value = solution[iColumn];
580    if (solver->isInteger(iColumn)) {
581      // Round down integer
582      if (fabs(floor(value+0.5)-value)<integerTolerance) {
583        value=floor(CoinMax(value+1.0e-3,columnLower[iColumn]));
584      } else {
585        value=CoinMax(floor(value),columnLower[iColumn]);
586      }
587    }
588    // make sure clean
589    value = CoinMin(value,columnUpper[iColumn]);
590    value = CoinMax(value,columnLower[iColumn]);
591    newSolution[iColumn]=value;
592    double cost = direction * objective[iColumn];
593    newSolutionValue += value*cost;
594    for (j=columnStart[iColumn];
595         j<columnStart[iColumn]+columnLength[iColumn];j++) {
596      int iRow=row[j];
597      rowActivity[iRow] += value*element[j];
598      rhsNeeded -= value*element[j];
599      if (element[j]!=1.0)
600        allOnes=false;
601    }
602  }
603  // See if we round up
604  bool roundup = ((algorithm_%100)!=0);
605  if (roundup&&allOnes) {
606    // Get rounded up solution
607    for (iColumn=0;iColumn<numberColumns;iColumn++) {
608      CoinBigIndex j;
609      double value = solution[iColumn];
610      if (solver->isInteger(iColumn)) {
611        // but round up if no activity
612        if (roundup&&value>=0.6&&!newSolution[iColumn]) {
613          bool choose=true;
614          for (j=columnStart[iColumn];
615               j<columnStart[iColumn]+columnLength[iColumn];j++) {
616            int iRow=row[j];
617            if (rowActivity[iRow]) {
618              choose=false;
619              break;
620            }
621          }
622          if (choose) {
623            newSolution[iColumn]=1.0;
624            double cost = direction * objective[iColumn];
625            newSolutionValue += cost;
626            for (j=columnStart[iColumn];
627                 j<columnStart[iColumn]+columnLength[iColumn];j++) {
628              int iRow=row[j];
629              rowActivity[iRow] += 1.0;
630              rhsNeeded -= 1.0;
631            }
632          }
633        }
634      }
635    }
636  }
637  // Get initial list
638  int * which = new int [numberColumns];
639  for (iColumn=0;iColumn<numberColumns;iColumn++) 
640    which[iColumn]=iColumn;
641  int numberLook=numberColumns;
642  // See if we want to perturb more
643  double perturb = ((algorithm_%10)==0) ? 0.1 : 0.25;
644  // Keep going round until a solution
645  while (true) {
646    // Get column with best ratio
647    int bestColumn=-1;
648    double bestRatio=COIN_DBL_MAX;
649    double bestStepSize = 0.0;
650    int newNumber=0;
651    for (int jColumn=0;jColumn<numberLook;jColumn++) {
652      int iColumn = which[jColumn];
653      CoinBigIndex j;
654      double value = newSolution[iColumn];
655      double cost = direction * objective[iColumn];
656      if (solver->isInteger(iColumn)) {
657        // use current upper or original upper
658        if (value+0.9999<originalUpper[iColumn]) {
659          double movement = 1.0;
660          double sum=0.0;
661          for (j=columnStart[iColumn];
662               j<columnStart[iColumn]+columnLength[iColumn];j++) {
663            int iRow=row[j];
664            double gap = rowUpper[iRow]-rowActivity[iRow];
665            double elementValue = allOnes ? 1.0 : element[j];
666            sum += elementValue;
667            if (movement*elementValue>gap) {
668              movement = gap/elementValue;
669            }
670          }
671          if (movement>0.999999) {
672            // add to next time
673            which[newNumber++]=iColumn;
674            double ratio = (cost/sum)*(1.0+perturb*randomNumberGenerator_.randomDouble());
675            // If at root
676            if (atRoot) {
677              if (fraction_==1.0)
678                ratio = iColumn; // choose first
679              else
680                ratio = - solution[iColumn]; // choose largest
681            }
682            if (ratio<bestRatio) {
683              bestRatio=ratio;
684              bestColumn=iColumn;
685              bestStepSize=1.0;
686            }
687          }
688        }
689      } else {
690        // continuous
691        if (value<columnUpper[iColumn]) {
692          double movement=1.0e50;
693          double sum=0.0;
694          for (j=columnStart[iColumn];
695               j<columnStart[iColumn]+columnLength[iColumn];j++) {
696            int iRow=row[j];
697            if (element[j]*movement+rowActivity[iRow]>rowUpper[iRow]) {
698              movement = (rowUpper[iRow]-rowActivity[iRow])/element[j];;
699            }
700            sum += element[j];
701          }
702          // now ratio
703          if (movement>1.0e-7) {
704            // add to next time
705            which[newNumber++]=iColumn;
706            double ratio = (cost/sum)*(1.0+perturb*randomNumberGenerator_.randomDouble());
707            if (ratio<bestRatio) {
708              bestRatio=ratio;
709              bestColumn=iColumn;
710              bestStepSize=movement;
711            }
712          }
713        }
714      }
715    }
716    if (bestColumn<0)
717      break; // we have finished
718    // Increase chosen column
719    newSolution[bestColumn] += bestStepSize;
720    double cost = direction * objective[bestColumn];
721    newSolutionValue += bestStepSize*cost;
722    for (CoinBigIndex j=columnStart[bestColumn];
723         j<columnStart[bestColumn]+columnLength[bestColumn];j++) {
724      int iRow = row[j];
725      rowActivity[iRow] += bestStepSize*element[j];
726      rhsNeeded -= bestStepSize*element[j];
727    }
728    if (rhsNeeded<1.0e-8)
729      break;
730  }
731  delete [] which;
732  if (fraction_<1.0&&rhsNeeded<1.0e-8&&newSolutionValue<solutionValue) {
733    // do branch and cut
734    // fix all nonzero
735    OsiSolverInterface * newSolver = model_->continuousSolver()->clone();
736    for (iColumn=0;iColumn<numberColumns;iColumn++) {
737      if (newSolver->isInteger(iColumn))
738        newSolver->setColLower(iColumn,newSolution[iColumn]);
739    }
740    int returnCode = smallBranchAndBound(newSolver,200,newSolution,newSolutionValue,
741                                         solutionValue,"CbcHeuristicGreedy");
742    if (returnCode<0)
743      returnCode=0; // returned on size
744    if ((returnCode&2)!=0) {
745      // could add cut
746      returnCode &= ~2;
747    }
748    rhsNeeded = 1.0-returnCode;
749    delete newSolver;
750  }
751  if (newSolutionValue<solutionValue&&rhsNeeded<1.0e-8) {
752    // check feasible
753    memset(rowActivity,0,numberRows*sizeof(double));
754    for (iColumn=0;iColumn<numberColumns;iColumn++) {
755      CoinBigIndex j;
756      double value = newSolution[iColumn];
757      if (value) {
758        for (j=columnStart[iColumn];
759             j<columnStart[iColumn]+columnLength[iColumn];j++) {
760          int iRow=row[j];
761          rowActivity[iRow] += value*element[j];
762        }
763      }
764    }
765    // check was approximately feasible
766    bool feasible=true;
767    for (iRow=0;iRow<numberRows;iRow++) {
768      if(rowActivity[iRow]<rowLower[iRow]) {
769        if (rowActivity[iRow]<rowLower[iRow]-10.0*primalTolerance)
770          feasible = false;
771      }
772    }
773    if (feasible) {
774      // new solution
775      memcpy(betterSolution,newSolution,numberColumns*sizeof(double));
776      solutionValue = newSolutionValue;
777      returnCode=1;
778    }
779  }
780  delete [] newSolution;
781  delete [] rowActivity;
782  if (atRoot&&fraction_==1.0) {
783    // try quick search
784    fraction_=0.4;
785    int newCode=this->solution(solutionValue,betterSolution);
786    if (newCode)
787      returnCode=1;
788    fraction_=1.0;
789  }
790  return returnCode;
791}
792// update model
793void CbcHeuristicGreedyEquality::setModel(CbcModel * model)
794{
795  gutsOfConstructor(model);
796  validate();
797}
798// Resets stuff if model changes
799void 
800CbcHeuristicGreedyEquality::resetModel(CbcModel * model)
801{
802  gutsOfConstructor(model);
803}
804// Validate model i.e. sets when_ to 0 if necessary (may be NULL)
805void 
806CbcHeuristicGreedyEquality::validate() 
807{
808  if (model_&&when()<10) {
809    if (model_->numberIntegers()!=
810        model_->numberObjects())
811      setWhen(0);
812    // Only works if costs positive, coefficients positive and all rows E or L
813    // And if values are integer
814    OsiSolverInterface * solver = model_->solver();
815    const double * columnLower = solver->getColLower();
816    const double * rowUpper = solver->getRowUpper();
817    const double * rowLower = solver->getRowLower();
818    const double * objective = solver->getObjCoefficients();
819    double direction = solver->getObjSense();
820
821    int numberRows = solver->getNumRows();
822    // Column copy
823    const double * element = matrix_.getElements();
824    const CoinBigIndex * columnStart = matrix_.getVectorStarts();
825    const int * columnLength = matrix_.getVectorLengths();
826    bool good = true;
827    for (int iRow=0;iRow<numberRows;iRow++) {
828      if (rowUpper[iRow]>1.0e30)
829        good = false;
830      if (rowLower[iRow]>0.0&&rowLower[iRow]!=rowUpper[iRow])
831        good = false;
832      if (floor(rowUpper[iRow]+0.5)!=rowUpper[iRow])
833        good = false;
834    }
835    int numberColumns = solver->getNumCols();
836    for (int iColumn=0;iColumn<numberColumns;iColumn++) {
837      if (objective[iColumn]*direction<0.0)
838        good=false;
839      if (columnLower[iColumn]<0.0)
840        good=false;
841      CoinBigIndex j;
842      for (j=columnStart[iColumn];
843           j<columnStart[iColumn]+columnLength[iColumn];j++) {
844        if (element[j]<0.0)
845          good=false;
846        if (floor(element[j]+0.5)!=element[j])
847          good = false;
848      }
849    }
850    if (!good)
851      setWhen(0); // switch off
852  }
853}
854
855 
Note: See TracBrowser for help on using the repository browser.