source: trunk/CbcHeuristicGreedy.cpp @ 279

Last change on this file since 279 was 279, checked in by forrest, 13 years ago

for timing

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