source: trunk/Samples/CbcHeuristicGreedy.cpp @ 139

Last change on this file since 139 was 139, checked in by forrest, 16 years ago

add equality

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 25.9 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
18CbcGreedyCover::CbcGreedyCover() 
19  :CbcHeuristic()
20{
21  // matrix  will automatically be empty
22  originalNumberRows_=0;
23  algorithm_=0;
24  numberTimes_=100;
25}
26
27// Constructor from model
28CbcGreedyCover::CbcGreedyCover(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
40CbcGreedyCover::~CbcGreedyCover ()
41{
42}
43
44// Clone
45CbcHeuristic *
46CbcGreedyCover::clone() const
47{
48  return new CbcGreedyCover(*this);
49}
50
51// Copy constructor
52CbcGreedyCover::CbcGreedyCover(const CbcGreedyCover & 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
64CbcGreedyCover & 
65CbcGreedyCover::operator=( const CbcGreedyCover& 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
78CbcGreedyCover::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              if (rowActivity[iRow]<rowLower[iRow]-1.0e-10&&
265                element[j]*step+rowActivity[iRow]>=rowLower[iRow]) {
266                sum += element[j];
267              }
268              assert (sum>0.0);
269              double ratio = (cost/sum)*(1.0+perturb*CoinDrand48());
270              if (ratio<bestRatio) {
271                bestRatio=ratio;
272                bestColumn=iColumn;
273                bestStepSize=step;
274              }
275            }
276          }
277        }
278      }
279    }
280    if (bestColumn<0)
281      break; // we have finished
282    // Increase chosen column
283    newSolution[bestColumn] += bestStepSize;
284    double cost = direction * objective[bestColumn];
285    newSolutionValue += bestStepSize*cost;
286    for (CoinBigIndex j=columnStart[bestColumn];
287         j<columnStart[bestColumn]+columnLength[bestColumn];j++) {
288      int iRow = row[j];
289      rowActivity[iRow] += bestStepSize*element[j];
290    }
291  }
292  delete [] which;
293  if (newSolutionValue<solutionValue) {
294    // check feasible
295    memset(rowActivity,0,numberRows*sizeof(double));
296    for (iColumn=0;iColumn<numberColumns;iColumn++) {
297      CoinBigIndex j;
298      double value = newSolution[iColumn];
299      if (value) {
300        for (j=columnStart[iColumn];
301             j<columnStart[iColumn]+columnLength[iColumn];j++) {
302          int iRow=row[j];
303          rowActivity[iRow] += value*element[j];
304        }
305      }
306    }
307    // check was approximately feasible
308    bool feasible=true;
309    for (iRow=0;iRow<numberRows;iRow++) {
310      if(rowActivity[iRow]<rowLower[iRow]) {
311        if (rowActivity[iRow]<rowLower[iRow]-10.0*primalTolerance)
312          feasible = false;
313      }
314    }
315    if (feasible) {
316      // new solution
317      memcpy(betterSolution,newSolution,numberColumns*sizeof(double));
318      solutionValue = newSolutionValue;
319      //printf("** Solution of %g found by rounding\n",newSolutionValue);
320      returnCode=1;
321    } else {
322      // Can easily happen
323      //printf("Debug CbcGreedyCover giving bad solution\n");
324    }
325  }
326  delete [] newSolution;
327  delete [] rowActivity;
328  return returnCode;
329}
330// update model
331void CbcGreedyCover::setModel(CbcModel * model)
332{
333#define SLOPPY
334#ifndef SLOPPY
335  model_ = model;
336  assert(model_->solver());
337  *this = CbcGreedyCover(*model); 
338#else
339  if (model_&&model!=model_) {
340    model_ = model;
341    assert(model_->solver());
342    *this = CbcGreedyCover(*model); 
343  }
344#endif
345  validate();
346}
347// Resets stuff if model changes
348void 
349CbcGreedyCover::resetModel(CbcModel * model)
350{
351#ifndef SLOPPY
352  model_ = model;
353  assert(model_->solver());
354  *this = CbcGreedyCover(*model); 
355#else
356  // switch off
357  model_ = NULL;
358  matrix_ = CoinPackedMatrix();
359#endif
360}
361// Validate model i.e. sets when_ to 0 if necessary (may be NULL)
362void 
363CbcGreedyCover::validate() 
364{
365  if (model_&&when()<10) {
366    if (model_->numberIntegers()!=
367        model_->numberObjects())
368      setWhen(0);
369    // Only works if costs positive, coefficients positive and all rows G
370    OsiSolverInterface * solver = model_->solver();
371    const double * columnLower = solver->getColLower();
372    const double * rowUpper = solver->getRowUpper();
373    const double * objective = solver->getObjCoefficients();
374    double direction = solver->getObjSense();
375
376    int numberRows = solver->getNumRows();
377    // Column copy
378    const double * element = matrix_.getElements();
379    const CoinBigIndex * columnStart = matrix_.getVectorStarts();
380    const int * columnLength = matrix_.getVectorLengths();
381    bool good = true;
382    for (int iRow=0;iRow<numberRows;iRow++) {
383      if (rowUpper[iRow]<1.0e30)
384        good = false;
385    }
386    int numberColumns = solver->getNumCols();
387    for (int iColumn=0;iColumn<numberColumns;iColumn++) {
388      if (objective[iColumn]*direction<0.0)
389        good=false;
390      if (columnLower[iColumn]<0.0)
391        good=false;
392      CoinBigIndex j;
393      for (j=columnStart[iColumn];
394           j<columnStart[iColumn]+columnLength[iColumn];j++) {
395        if (element[j]<0.0)
396          good=false;
397      }
398    }
399    if (!good)
400      setWhen(0); // switch off
401  }
402}
403// Default Constructor
404CbcGreedyEquality::CbcGreedyEquality() 
405  :CbcHeuristic()
406{
407  // matrix  will automatically be empty
408  fraction_=1.0; // no branch and bound
409  originalNumberRows_=0;
410  algorithm_=0;
411  numberTimes_=100;
412}
413
414// Constructor from model
415CbcGreedyEquality::CbcGreedyEquality(CbcModel & model)
416  :CbcHeuristic(model)
417{
418  // Get a copy of original matrix
419  assert(model.solver());
420  matrix_ = *model.solver()->getMatrixByCol();
421  fraction_=1.0; // no branch and bound
422  originalNumberRows_=model.solver()->getNumRows();
423  algorithm_=0;
424  numberTimes_=100;
425}
426
427// Destructor
428CbcGreedyEquality::~CbcGreedyEquality ()
429{
430}
431
432// Clone
433CbcHeuristic *
434CbcGreedyEquality::clone() const
435{
436  return new CbcGreedyEquality(*this);
437}
438
439// Copy constructor
440CbcGreedyEquality::CbcGreedyEquality(const CbcGreedyEquality & rhs)
441:
442  CbcHeuristic(rhs),
443  matrix_(rhs.matrix_),
444  fraction_(rhs.fraction_),
445  originalNumberRows_(rhs.originalNumberRows_),
446  algorithm_(rhs.algorithm_),
447  numberTimes_(rhs.numberTimes_)
448{
449  this->setWhen(rhs.when());
450}
451
452// Assignment operator
453CbcGreedyEquality & 
454CbcGreedyEquality::operator=( const CbcGreedyEquality& rhs)
455{
456  if (this != &rhs) {
457    setWhen(rhs.when());
458    matrix_=rhs.matrix_;
459    fraction_ = rhs.fraction_;
460    originalNumberRows_=rhs.originalNumberRows_;
461    algorithm_=rhs.algorithm_;
462    numberTimes_=rhs.numberTimes_;
463  }
464  return *this;
465}
466// Returns 1 if solution, 0 if not
467int
468CbcGreedyEquality::solution(double & solutionValue,
469                         double * betterSolution)
470{
471  if (!model_)
472    return 0;
473  // See if to do
474  if (!when()||(when()==1&&model_->phase()!=1))
475    return 0; // switched off
476  if (model_->getNodeCount()>numberTimes_)
477    return 0;
478  // See if at root node
479  bool atRoot = model_->getNodeCount()==0;
480  int passNumber = model_->getCurrentPassNumber();
481  if (atRoot&&passNumber!=1)
482    return 0;
483  OsiSolverInterface * solver = model_->solver();
484  const double * columnLower = solver->getColLower();
485  const double * columnUpper = solver->getColUpper();
486  // And original upper bounds in case we want to use them
487  const double * originalUpper = model_->continuousSolver()->getColUpper();
488  // But not if algorithm says so
489  if ((algorithm_%10)==0)
490    originalUpper = columnUpper;
491  const double * rowLower = solver->getRowLower();
492  const double * rowUpper = solver->getRowUpper();
493  const double * solution = solver->getColSolution();
494  const double * objective = solver->getObjCoefficients();
495  double integerTolerance = model_->getDblParam(CbcModel::CbcIntegerTolerance);
496  double primalTolerance;
497  solver->getDblParam(OsiPrimalTolerance,primalTolerance);
498
499  // This is number of rows when matrix was passed in
500  int numberRows = originalNumberRows_;
501  if (!numberRows)
502    return 0; // switched off
503
504  assert (numberRows==matrix_.getNumRows());
505  int iRow, iColumn;
506  double direction = solver->getObjSense();
507  double offset;
508  solver->getDblParam(OsiObjOffset,offset);
509  double newSolutionValue = -offset;
510  int returnCode = 0;
511
512  // Column copy
513  const double * element = matrix_.getElements();
514  const int * row = matrix_.getIndices();
515  const CoinBigIndex * columnStart = matrix_.getVectorStarts();
516  const int * columnLength = matrix_.getVectorLengths();
517
518  // Get solution array for heuristic solution
519  int numberColumns = solver->getNumCols();
520  double * newSolution = new double [numberColumns];
521  double * rowActivity = new double[numberRows];
522  memset(rowActivity,0,numberRows*sizeof(double));
523  double rhsNeeded=0;
524  for (iRow=0;iRow<numberRows;iRow++) 
525    rhsNeeded += rowUpper[iRow];
526  rhsNeeded *= fraction_;
527  bool allOnes=true;
528  // Get rounded down solution
529  for (iColumn=0;iColumn<numberColumns;iColumn++) {
530    CoinBigIndex j;
531    double value = solution[iColumn];
532    if (solver->isInteger(iColumn)) {
533      // Round down integer
534      if (fabs(floor(value+0.5)-value)<integerTolerance) {
535        value=floor(CoinMax(value+1.0e-3,columnLower[iColumn]));
536      } else {
537        value=CoinMax(floor(value),columnLower[iColumn]);
538      }
539    }
540    // make sure clean
541    value = CoinMin(value,columnUpper[iColumn]);
542    value = CoinMax(value,columnLower[iColumn]);
543    newSolution[iColumn]=value;
544    double cost = direction * objective[iColumn];
545    newSolutionValue += value*cost;
546    for (j=columnStart[iColumn];
547         j<columnStart[iColumn]+columnLength[iColumn];j++) {
548      int iRow=row[j];
549      rowActivity[iRow] += value*element[j];
550      rhsNeeded -= value*element[j];
551      if (element[j]!=1.0)
552        allOnes=false;
553    }
554  }
555  // See if we round up
556  bool roundup = ((algorithm_%100)!=0);
557  if (roundup&&allOnes) {
558    // Get rounded up solution
559    for (iColumn=0;iColumn<numberColumns;iColumn++) {
560      CoinBigIndex j;
561      double value = solution[iColumn];
562      if (solver->isInteger(iColumn)) {
563        // but round up if no activity
564        if (roundup&&value>=0.6&&!newSolution[iColumn]) {
565          bool choose=true;
566          for (j=columnStart[iColumn];
567               j<columnStart[iColumn]+columnLength[iColumn];j++) {
568            int iRow=row[j];
569            if (rowActivity[iRow]) {
570              choose=false;
571              break;
572            }
573          }
574          if (choose) {
575            newSolution[iColumn]=1.0;
576            double cost = direction * objective[iColumn];
577            newSolutionValue += cost;
578            for (j=columnStart[iColumn];
579                 j<columnStart[iColumn]+columnLength[iColumn];j++) {
580              int iRow=row[j];
581              rowActivity[iRow] += 1.0;
582              rhsNeeded -= 1.0;
583            }
584          }
585        }
586      }
587    }
588  }
589  // Get initial list
590  int * which = new int [numberColumns];
591  for (iColumn=0;iColumn<numberColumns;iColumn++) 
592    which[iColumn]=iColumn;
593  int numberLook=numberColumns;
594  // See if we want to perturb more
595  double perturb = ((algorithm_%10)==0) ? 0.1 : 0.25;
596  // Keep going round until a solution
597  while (true) {
598    // Get column with best ratio
599    int bestColumn=-1;
600    double bestRatio=COIN_DBL_MAX;
601    double bestStepSize = 0.0;
602    int newNumber=0;
603    for (int jColumn=0;jColumn<numberLook;jColumn++) {
604      int iColumn = which[jColumn];
605      CoinBigIndex j;
606      double value = newSolution[iColumn];
607      double cost = direction * objective[iColumn];
608      if (solver->isInteger(iColumn)) {
609        // use current upper or original upper
610        if (value+0.9999<originalUpper[iColumn]) {
611          double movement = 1.0;
612          double sum=0.0;
613          for (j=columnStart[iColumn];
614               j<columnStart[iColumn]+columnLength[iColumn];j++) {
615            int iRow=row[j];
616            double gap = rowUpper[iRow]-rowActivity[iRow];
617            double elementValue = allOnes ? 1.0 : element[j];
618            sum += elementValue;
619            if (movement*elementValue>gap) {
620              movement = gap/elementValue;
621            }
622          }
623          if (movement>0.999999) {
624            // add to next time
625            which[newNumber++]=iColumn;
626            double ratio = (cost/sum)*(1.0+perturb*CoinDrand48());
627            // If at root
628            if (atRoot) {
629              if (fraction_==1.0)
630                ratio = iColumn; // choose first
631              else
632                ratio = - solution[iColumn]; // choose largest
633            }
634            if (ratio<bestRatio) {
635              bestRatio=ratio;
636              bestColumn=iColumn;
637              bestStepSize=1.0;
638            }
639          }
640        }
641      } else {
642        // continuous
643        if (value<columnUpper[iColumn]) {
644          double movement=1.0e50;
645          double sum=0.0;
646          for (j=columnStart[iColumn];
647               j<columnStart[iColumn]+columnLength[iColumn];j++) {
648            int iRow=row[j];
649            if (element[j]*movement+rowActivity[iRow]>rowUpper[iRow]) {
650              movement = (rowUpper[iRow]-rowActivity[iRow])/element[j];;
651            }
652            sum += element[j];
653          }
654          // now ratio
655          if (movement>1.0e-7) {
656            // add to next time
657            which[newNumber++]=iColumn;
658            double ratio = (cost/sum)*(1.0+perturb*CoinDrand48());
659            if (ratio<bestRatio) {
660              bestRatio=ratio;
661              bestColumn=iColumn;
662              bestStepSize=movement;
663            }
664          }
665        }
666      }
667    }
668    if (bestColumn<0)
669      break; // we have finished
670    // Increase chosen column
671    newSolution[bestColumn] += bestStepSize;
672    double cost = direction * objective[bestColumn];
673    newSolutionValue += bestStepSize*cost;
674    for (CoinBigIndex j=columnStart[bestColumn];
675         j<columnStart[bestColumn]+columnLength[bestColumn];j++) {
676      int iRow = row[j];
677      rowActivity[iRow] += bestStepSize*element[j];
678      rhsNeeded -= bestStepSize*element[j];
679    }
680    if (rhsNeeded<1.0e-8)
681      break;
682  }
683  delete [] which;
684  if (fraction_<1.0&&rhsNeeded<1.0e-8&&newSolutionValue<solutionValue) {
685    // do branch and cut
686    // fix all nonzero
687    OsiSolverInterface * newSolver = model_->continuousSolver()->clone();
688    for (iColumn=0;iColumn<numberColumns;iColumn++) {
689      if (newSolver->isInteger(iColumn))
690        newSolver->setColLower(iColumn,newSolution[iColumn]);
691    }
692    // Reduce printout
693    newSolver->setHintParam(OsiDoReducePrint,true,OsiHintTry);
694    newSolver->setHintParam(OsiDoPresolveInInitial,true,OsiHintTry);
695    newSolver->setDblParam(OsiDualObjectiveLimit,solutionValue);
696    newSolver->initialSolve();
697    if (newSolver->isProvenOptimal()) {
698      CglPreProcess process;
699      /* Do not try and produce equality cliques and
700         do up to 5 passes */
701      OsiSolverInterface * solver2= process.preProcess(*newSolver);
702      if (!solver2) {
703        printf("Pre-processing says infeasible\n");
704        rhsNeeded=1.0; // so will be infeasible
705      } else {
706        solver2->resolve();
707        CbcModel model(*solver2);
708        model.setCutoff(solutionValue);
709        model.setMaximumNodes(200);
710        model.solver()->setHintParam(OsiDoReducePrint,true,OsiHintTry);
711        CbcStrategyDefault strategy(true,5,5,0);
712        model.setStrategy(strategy);
713        // Lightweight
714        model.setNumberStrong(5);
715        model.setNumberBeforeTrust(1);
716        model.solver()->setIntParam(OsiMaxNumIterationHotStart,10);
717        // Do complete search
718        model.branchAndBound();
719        if (model.getMinimizationObjValue()<solutionValue) {
720          // solution
721          rhsNeeded=0.0;
722          // post process
723          process.postProcess(*model.solver());
724          // Solution now back in newSolver
725          memcpy(newSolution,newSolver->getColSolution(),
726                 numberColumns*sizeof(double));
727          newSolutionValue = model.getMinimizationObjValue();
728        } else {
729          // no good
730          rhsNeeded=1.0; // so will be infeasible
731        }
732      }
733    } else {
734      rhsNeeded=1.0;
735    }
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      printf("** Solution of %g found by greedy\n",newSolutionValue);
765      returnCode=1;
766    }
767  }
768  delete [] newSolution;
769  delete [] rowActivity;
770  if (atRoot&&fraction_==1.0) {
771    // try quick search
772    fraction_=0.3;
773    int newCode=this->solution(solutionValue,betterSolution);
774    if (newCode)
775      returnCode=1;
776    fraction_=1.0;
777  }
778  return returnCode;
779}
780// update model
781void CbcGreedyEquality::setModel(CbcModel * model)
782{
783#define SLOPPY
784#ifndef SLOPPY
785  model_ = model;
786  assert(model_->solver());
787  *this = CbcGreedyEquality(*model); 
788#else
789  if (model_&&model!=model_) {
790    model_ = model;
791    assert(model_->solver());
792    *this = CbcGreedyEquality(*model); 
793  }
794#endif
795  validate();
796}
797// Resets stuff if model changes
798void 
799CbcGreedyEquality::resetModel(CbcModel * model)
800{
801#ifndef SLOPPY
802  model_ = model;
803  assert(model_->solver());
804  *this = CbcGreedyEquality(*model); 
805#else
806  // switch off
807  model_ = NULL;
808  matrix_ = CoinPackedMatrix();
809#endif
810}
811// Validate model i.e. sets when_ to 0 if necessary (may be NULL)
812void 
813CbcGreedyEquality::validate() 
814{
815  if (model_&&when()<10) {
816    if (model_->numberIntegers()!=
817        model_->numberObjects())
818      setWhen(0);
819    // Only works if costs positive, coefficients positive and all rows E or L
820    // And if values are integer
821    OsiSolverInterface * solver = model_->solver();
822    const double * columnLower = solver->getColLower();
823    const double * rowUpper = solver->getRowUpper();
824    const double * rowLower = solver->getRowLower();
825    const double * objective = solver->getObjCoefficients();
826    double direction = solver->getObjSense();
827
828    int numberRows = solver->getNumRows();
829    // Column copy
830    const double * element = matrix_.getElements();
831    const CoinBigIndex * columnStart = matrix_.getVectorStarts();
832    const int * columnLength = matrix_.getVectorLengths();
833    bool good = true;
834    for (int iRow=0;iRow<numberRows;iRow++) {
835      if (rowUpper[iRow]>1.0e30)
836        good = false;
837      if (rowLower[iRow]>0.0&&rowLower[iRow]!=rowUpper[iRow])
838        good = false;
839      if (floor(rowUpper[iRow]+0.5)!=rowUpper[iRow])
840        good = false;
841    }
842    int numberColumns = solver->getNumCols();
843    for (int iColumn=0;iColumn<numberColumns;iColumn++) {
844      if (objective[iColumn]*direction<0.0)
845        good=false;
846      if (columnLower[iColumn]<0.0)
847        good=false;
848      CoinBigIndex j;
849      for (j=columnStart[iColumn];
850           j<columnStart[iColumn]+columnLength[iColumn];j++) {
851        if (element[j]<0.0)
852          good=false;
853        if (floor(element[j]+0.5)!=element[j])
854          good = false;
855      }
856    }
857    if (!good)
858      setWhen(0); // switch off
859  }
860}
861
862 
Note: See TracBrowser for help on using the repository browser.