source: stable/2.0/Cbc/src/CbcHeuristicGreedy.cpp @ 905

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

include cstdlib before cmath to get things to compile on AIX with xlC (same as changeset 904 in trunk)

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