Changeset 139


Ignore:
Timestamp:
May 20, 2005 11:48:21 AM (16 years ago)
Author:
forrest
Message:

add equality

Location:
trunk/Samples
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • trunk/Samples/CbcHeuristicGreedy.cpp

    r83 r139  
    1111#include "OsiSolverInterface.hpp"
    1212#include "CbcModel.hpp"
     13#include "CbcStrategy.hpp"
    1314#include "CbcHeuristicGreedy.hpp"
    1415#include "CoinSort.hpp"
     16#include "CglPreProcess.hpp"
    1517// Default Constructor
    1618CbcGreedyCover::CbcGreedyCover()
     
    2022  originalNumberRows_=0;
    2123  algorithm_=0;
     24  numberTimes_=100;
    2225}
    2326
     
    3134  originalNumberRows_=model.solver()->getNumRows();
    3235  algorithm_=0;
     36  numberTimes_=100;
    3337}
    3438
     
    5155  matrix_(rhs.matrix_),
    5256  originalNumberRows_(rhs.originalNumberRows_),
    53   algorithm_(rhs.algorithm_)
     57  algorithm_(rhs.algorithm_),
     58  numberTimes_(rhs.numberTimes_)
    5459{
    5560  this->setWhen(rhs.when());
     
    6570    originalNumberRows_=rhs.originalNumberRows_;
    6671    algorithm_=rhs.algorithm_;
     72    numberTimes_=rhs.numberTimes_;
    6773  }
    6874  return *this;
     
    7884  if (!when()||(when()==1&&model_->phase()!=1))
    7985    return 0; // switched off
    80   if (model_->getNodeCount()>500)
     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)
    8192    return 0;
    8293  OsiSolverInterface * solver = model_->solver();
     
    219230            which[newNumber++]=iColumn;
    220231            double ratio = (cost/sum)*(1.0+perturb*CoinDrand48());
     232            // If at root choose first
     233            if (atRoot)
     234              ratio = iColumn;
    221235            if (ratio<bestRatio) {
    222236              bestRatio=ratio;
     
    387401  }
    388402}
     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}
    389861
    390862 
  • trunk/Samples/CbcHeuristicGreedy.hpp

    r83 r139  
    1 // Copyright (C) 2002, International Business Machines
     1// Copyright (C) 2005, International Business Machines
    22// Corporation and others.  All Rights Reserved.
    33#ifndef CbcHeuristicGreedy_H
     
    55
    66#include "CbcHeuristic.hpp"
    7 /** Greedy heuristic class
     7/** Greedy heuristic classes
    88 */
    99
     
    5757  inline void setAlgorithm(int value)
    5858  { algorithm_=value; };
     59  // Only do this many times
     60  inline int numberTimes() const
     61  { return numberTimes_; };
     62  inline void setNumberTimes(int value)
     63  { numberTimes_=value; };
    5964
    6065protected:
     
    7176  */
    7277  int algorithm_;
     78  /// Do this many times
     79  int numberTimes_;
    7380
    7481private:
     
    7885
    7986
     87class CbcGreedyEquality : public CbcHeuristic {
     88public:
     89
     90  // Default Constructor
     91  CbcGreedyEquality ();
     92
     93  /* Constructor with model - assumed before cuts
     94     Initial version does not do Lps
     95  */
     96  CbcGreedyEquality (CbcModel & model);
     97 
     98  // Copy constructor
     99  CbcGreedyEquality ( const CbcGreedyEquality &);
     100   
     101  // Destructor
     102  ~CbcGreedyEquality ();
     103 
     104  /// Clone
     105  virtual CbcHeuristic * clone() const;
     106
     107  /// update model (This is needed if cliques update matrix etc)
     108  virtual void setModel(CbcModel * model);
     109 
     110  /** returns 0 if no solution, 1 if valid solution.
     111      Sets solution values if good, sets objective value (only if good)
     112      We leave all variables which are at one at this node of the
     113      tree to that value and will
     114      initially set all others to zero.  We then sort all variables in order of their cost
     115      divided by the number of entries in rows which are not yet covered.  We randomize that
     116      value a bit so that ties will be broken in different ways on different runs of the heuristic.
     117      We then choose the best one and set it to one and repeat the exercise. 
     118
     119  */
     120  virtual int solution(double & objectiveValue,
     121                       double * newSolution);
     122  /// Validate model i.e. sets when_ to 0 if necessary (may be NULL)
     123  virtual void validate() ;
     124  /// Resets stuff if model changes
     125  virtual void resetModel(CbcModel * model);
     126  /* Algorithm
     127     0 - use current upper bounds
     128     1 - use original upper bounds
     129     If 10 added perturb ratios more
     130     if 100 added round up all >=0.5
     131  */
     132  inline int algorithm() const
     133  { return algorithm_; };
     134  inline void setAlgorithm(int value)
     135  { algorithm_=value; };
     136  // Fraction of rhs to cover before branch and cut
     137  inline void setFraction(double value)
     138  { fraction_ = value;};
     139  inline double fraction() const
     140  { return fraction_;};
     141  // Only do this many times
     142  inline int numberTimes() const
     143  { return numberTimes_; };
     144  inline void setNumberTimes(int value)
     145  { numberTimes_=value; };
     146protected:
     147  // Data
     148
     149  // Original matrix by column
     150  CoinPackedMatrix matrix_;
     151  // Fraction of rhs to cover before branch and cut
     152  double fraction_;
     153  // original number of rows
     154  int originalNumberRows_;
     155  /* Algorithm
     156     0 - use current upper bounds
     157     1 - use original upper bounds
     158     If 10 added perturb ratios more
     159  */
     160  int algorithm_;
     161  /// Do this many times
     162  int numberTimes_;
     163
     164private:
     165  /// Illegal Assignment operator
     166  CbcGreedyEquality & operator=(const CbcGreedyEquality& rhs);
     167};
     168
     169
    80170#endif
Note: See TracChangeset for help on using the changeset viewer.