1 // Last edit: 02/05/2013 2 // 3 // Name: CglGMI.hpp 4 // Author: Giacomo Nannicini 5 // Singapore University of Technology and Design, Singapore 6 // email: nannicini@sutd.edu.sg 7 // Date: 11/17/09 8 //----------------------------------------------------------------------------- 9 // Copyright (C) 2009, Giacomo Nannicini. All Rights Reserved. 10 11 #ifndef CglGMI_H 12 #define CglGMI_H 13 14 #include "CglCutGenerator.hpp" 15 #include "CglGMIParam.hpp" 16 #include "CoinWarmStartBasis.hpp" 17 #include "CoinFactorization.hpp" 18 19 /* Enable tracking of rejection of cutting planes. If this is disabled, 20 the cut generator is slightly faster. If defined, it enables proper use 21 of setTrackRejection and related functions. */ 22 //#define TRACK_REJECT 23 24 /* Debug output */ 25 //#define GMI_TRACE 26 27 /* Debug output: print optimal tableau */ 28 //#define GMI_TRACETAB 29 30 /* Print reason for cut rejection, whenever a cut is discarded */ 31 //#define GMI_TRACE_CLEAN 32 33 /** Gomory cut generator with several cleaning procedures, used to test 34 * the numerical safety of the resulting cuts 35 */ 36 37 class CglGMI : public CglCutGenerator { 38 39 friend void CglGMIUnitTest(const OsiSolverInterface * siP, 40 const std::string mpdDir); 41 public: 42 43 /** Public enum: all possible reasons for cut rejection */ 44 enum RejectionType{ 45 failureFractionality, 46 failureDynamism, 47 failureViolation, 48 failureSupport, 49 failureScale 50 }; 51 52 /**@name generateCuts */ 53 //@{ 54 /** Generate Gomory Mixed-Integer cuts for the model of the solver 55 interface si. 56 57 Insert the generated cuts into OsiCuts cs. 58 59 Warning: This generator currently works only with the Lp solvers Clp or 60 Cplex9.0 or higher. It requires access to the optimal tableau and 61 optimal basis inverse and makes assumptions on the way slack variables 62 are added by the solver. The Osi implementations for Clp and Cplex 63 verify these assumptions. 64 65 When calling the generator, the solver interface si must contain 66 an optimized problem and information related to the optimal 67 basis must be available through the OsiSolverInterface methods 68 (si->optimalBasisIsAvailable() must return 'true'). It is also 69 essential that the integrality of structural variable i can be 70 obtained using si->isInteger(i). 71 72 */ 73 virtual void generateCuts(const OsiSolverInterface & si, OsiCuts & cs, 74 const CglTreeInfo info = CglTreeInfo()); 75 76 /// Return true if needs optimal basis to do cuts (will return true) needsOptimalBasis() const77 virtual bool needsOptimalBasis() const { return true; } 78 //@} 79 80 /**@name Common Methods */ 81 //@{ 82 // Function for checking equality with user tolerance areEqual(double x,double y,double epsAbs=1e-12,double epsRel=1e-12)83 inline bool areEqual(double x, double y, 84 double epsAbs = 1e-12, 85 double epsRel = 1e-12) { 86 return (fabs((x) - (y)) <= 87 std::max(epsAbs, epsRel * std::max(fabs(x), fabs(y)))); 88 } 89 90 // Function for checking is a number is zero isZero(double x,double epsZero=1e-20)91 inline bool isZero(double x, double epsZero = 1e-20) { 92 return (fabs(x) <= epsZero); 93 } 94 95 96 // Function for checking if a number is integer isIntegerValue(double x,double intEpsAbs=1e-9,double intEpsRel=1e-15)97 inline bool isIntegerValue(double x, 98 double intEpsAbs = 1e-9, 99 double intEpsRel = 1e-15) { 100 return (fabs((x) - floor((x)+0.5)) <= 101 std::max(intEpsAbs, intEpsRel * fabs(x))); 102 } 103 104 105 //@} 106 107 108 /**@name Public Methods */ 109 //@{ 110 111 // Set the parameters to the values of the given CglGMIParam object. 112 void setParam(const CglGMIParam &source); 113 // Return the CglGMIParam object of the generator. getParam() const114 inline CglGMIParam getParam() const {return param;} getParam()115 inline CglGMIParam & getParam() {return param;} 116 117 // Compute entries of is_integer. 118 void computeIsInteger(); 119 120 /// Print the current simplex tableau 121 void printOptTab(OsiSolverInterface *solver) const; 122 123 /// Set/get tracking of the rejection of cutting planes. 124 /// Note that all rejection related functions will not do anything 125 /// unless the generator is compiled with the define GMI_TRACK_REJECTION 126 void setTrackRejection(bool value); 127 bool getTrackRejection(); 128 129 /// Get number of cuts rejected for given reason; see above 130 int getNumberRejectedCuts(RejectionType reason); 131 132 /// Reset counters for cut rejection tracking; see above 133 void resetRejectionCounters(); 134 135 /// Get total number of generated cuts since last resetRejectionCounters() 136 int getNumberGeneratedCuts(); 137 138 //@} 139 140 /**@name Constructors and destructors */ 141 //@{ 142 /// Default constructor 143 CglGMI(); 144 145 /// Constructor with specified parameters 146 CglGMI(const CglGMIParam ¶m); 147 148 /// Copy constructor 149 CglGMI(const CglGMI &); 150 151 /// Clone 152 virtual CglCutGenerator * clone() const; 153 154 /// Assignment operator 155 CglGMI & operator=(const CglGMI& rhs); 156 157 /// Destructor 158 virtual ~CglGMI(); 159 /// Create C++ lines to get to current state 160 virtual std::string generateCpp( FILE * fp); 161 162 //@} 163 164 private: 165 166 // Private member methods 167 168 /**@name Private member methods */ 169 170 //@{ 171 172 // Method generating the cuts after all CglGMI members are properly set. 173 void generateCuts(OsiCuts & cs); 174 175 /// Compute the fractional part of value, allowing for small error. 176 inline double aboveInteger(double value) const; 177 178 /// Compute the fractionalities involved in the cut, and the cut rhs. 179 /// Returns true if cut is accepted, false if discarded 180 inline bool computeCutFractionality(double varRhs, double& cutRhs); 181 182 /// Compute the cut coefficient on a given variable 183 inline double computeCutCoefficient(double rowElem, int index); 184 185 /// Use multiples of the initial inequalities to cancel out the coefficient 186 /// on a slack variables. 187 inline void eliminateSlack(double cutElem, int cutIndex, double* cut, 188 double& cutRhs, const double *elements, 189 const CoinBigIndex *rowStart, const int *indices, 190 const int *rowLength, const double *rhs); 191 192 /// Change the sign of the coefficients of the non basic 193 /// variables at their upper bound. 194 inline void flip(double& rowElem, int rowIndex); 195 196 /// Change the sign of the coefficients of the non basic 197 /// variables at their upper bound and do the translations restoring 198 /// the original bounds. Modify the right hand side 199 /// accordingly. Two functions: one for original variables, one for slacks. 200 inline void unflipOrig(double& rowElem, int rowIndex, double& rowRhs); 201 inline void unflipSlack(double& rowElem, int rowIndex, double& rowRhs, 202 const double* slack_val); 203 204 /// Pack a row of ncol elements 205 inline void packRow(double* row, double* rowElem, int* rowIndex, 206 int& rowNz); 207 208 /// Clean the cutting plane; the cleaning procedure does several things 209 /// like removing small coefficients, scaling, and checks several 210 /// acceptance criteria. If this returns false, the cut should be discarded. 211 /// There are several cleaning procedures available, that can be selected 212 /// via the parameter param.setCLEANING_PROCEDURE(int value) 213 bool cleanCut(double* cutElem, int* cutIndex, int& cutNz, 214 double& cutRhs, const double* xbar); 215 216 /// Cut cleaning procedures: return true if successfull, false if 217 /// cut should be discarded by the caller of if problems encountered 218 219 /// Check the violation 220 bool checkViolation(const double* cutElem, const int* cutIndex, 221 int cutNz, double cutrhs, const double* xbar); 222 223 /// Check the dynamism 224 bool checkDynamism(const double* cutElem, const int* cutIndex, 225 int cutNz); 226 227 /// Check the support 228 bool checkSupport(int cutNz); 229 230 /// Remove small coefficients and adjust the rhs accordingly 231 bool removeSmallCoefficients(double* cutElem, int* cutIndex, 232 int& cutNz, double& cutRhs); 233 234 /// Adjust the rhs by relaxing by a small amount (relative or absolute) 235 void relaxRhs(double& rhs); 236 237 /// Scale the cutting plane in different ways; 238 /// scaling_type possible values: 239 /// 0 : scale to obtain integral cut 240 /// 1 : scale based on norm, to obtain cut norm equal to ncol 241 /// 2 : scale to obtain largest coefficient equal to 1 242 bool scaleCut(double* cutElem, int* cutIndex, int cutNz, 243 double& cutRhs, int scalingType); 244 245 /// Scale the cutting plane in order to generate integral coefficients 246 bool scaleCutIntegral(double* cutElem, int* cutIndex, int cutNz, 247 double& cutRhs); 248 249 /// Compute the nearest rational number; used by scale_row_integral 250 bool nearestRational(double val, double maxdelta, long maxdnom, 251 long& numerator, long& denominator); 252 253 /// Compute the greatest common divisor 254 long computeGcd(long a, long b); 255 256 /// print a vector of integers 257 void printvecINT(const char *vecstr, const int *x, int n) const; 258 /// print a vector of doubles: dense form 259 void printvecDBL(const char *vecstr, const double *x, int n) const; 260 /// print a vector of doubles: sparse form 261 void printvecDBL(const char *vecstr, const double *elem, const int * index, 262 int nz) const; 263 264 /// Recompute the simplex tableau for want of a better accuracy. 265 /// Requires an empty CoinFactorization object to do the computations, 266 /// and two empty (already allocated) arrays which will contain 267 /// the basis indices on exit. Returns 0 if successfull. 268 int factorize(CoinFactorization & factorization, 269 int* colBasisIndex, int* rowBasisIndex); 270 271 272 //@} 273 274 275 // Private member data 276 277 /**@name Private member data */ 278 279 //@{ 280 281 /// Object with CglGMIParam members. 282 CglGMIParam param; 283 284 /// Number of rows ( = number of slack variables) in the current LP. 285 int nrow; 286 287 /// Number of structural variables in the current LP. 288 int ncol; 289 290 /// Lower bounds for structural variables 291 const double *colLower; 292 293 /// Upper bounds for structural variables 294 const double *colUpper; 295 296 /// Lower bounds for constraints 297 const double *rowLower; 298 299 /// Upper bounds for constraints 300 const double *rowUpper; 301 302 /// Righ hand side for constraints (upper bound for ranged constraints). 303 const double *rowRhs; 304 305 /// Characteristic vectors of structural integer variables or continuous 306 /// variables currently fixed to integer values. 307 bool *isInteger; 308 309 /// Current basis status: columns 310 int *cstat; 311 312 /// Current basis status: rows 313 int *rstat; 314 315 /// Pointer on solver. Reset by each call to generateCuts(). 316 OsiSolverInterface *solver; 317 318 /// Pointer on point to separate. Reset by each call to generateCuts(). 319 const double *xlp; 320 321 /// Pointer on row activity. Reset by each call to generateCuts(). 322 const double *rowActivity; 323 324 /// Pointer on matrix of coefficient ordered by rows. 325 /// Reset by each call to generateCuts(). 326 const CoinPackedMatrix *byRow; 327 328 /// Pointer on matrix of coefficient ordered by columns. 329 /// Reset by each call to generateCuts(). 330 const CoinPackedMatrix *byCol; 331 332 /// Fractionality of the cut and related quantities. 333 double f0; 334 double f0compl; 335 double ratiof0compl; 336 337 #if defined(TRACK_REJECT) || defined (TRACK_REJECT_SIMPLE) 338 /// Should we track the reason of each cut rejection? 339 bool trackRejection; 340 /// Number of failures by type 341 int fracFail; 342 int dynFail; 343 int violFail; 344 int suppFail; 345 int smallCoeffFail; 346 int scaleFail; 347 /// Total number of generated cuts 348 int numGeneratedCuts; 349 #endif 350 351 //@} 352 }; 353 354 //############################################################################# 355 /** A function that tests the methods in the CglGMI class. The 356 only reason for it not to be a member method is that this way it doesn't 357 have to be compiled into the library. And that's a gain, because the 358 library should be compiled with optimization on, but this method should be 359 compiled with debugging. */ 360 void CglGMIUnitTest(const OsiSolverInterface * siP, 361 const std::string mpdDir ); 362 363 364 #endif 365