source: branches/devel/Cbc/src/CbcHeuristicGreedy.cpp @ 439

Last change on this file since 439 was 356, checked in by ladanyi, 13 years ago

finishing conversion to svn

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