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

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

fix seg fault on preprocessed problem

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