1 // Copyright (C) 2002, International Business Machines 2 // Corporation and others. All Rights Reserved. 3 4 /* 5 Authors 6 7 John Forrest 8 9 */ 10 #ifndef IClpSimplexPrimal_H 11 #define IClpSimplexPrimal_H 12 13 #include "IClpSimplex.hpp" 14 //#define NPY_NO_DEPRECATED_API 15 16 17 /** This solves LPs using the primal simplex method 18 19 It inherits from IClpSimplex. It has no data of its own and 20 is never created - only cast from a IClpSimplex object at algorithm time. 21 22 */ 23 24 class IClpSimplexPrimal_Wolfe : public IClpSimplex{ 25 26 public: 27 28 /**@name Description of algorithm */ 29 //@{ 30 /** Primal algorithm 31 32 Method 33 34 It tries to be a single phase approach with a weight of 1.0 being 35 given to getting optimal and a weight of infeasibilityCost_ being 36 given to getting primal feasible. In this version I have tried to 37 be clever in a stupid way. The idea of fake bounds in dual 38 seems to work so the primal analogue would be that of getting 39 bounds on reduced costs (by a presolve approach) and using 40 these for being above or below feasible region. I decided to waste 41 memory and keep these explicitly. This allows for non-linear 42 costs! I have not tested non-linear costs but will be glad 43 to do something if a reasonable example is provided. 44 45 The code is designed to take advantage of sparsity so arrays are 46 seldom zeroed out from scratch or gone over in their entirety. 47 The only exception is a full scan to find incoming variable for 48 Dantzig row choice. For steepest edge we keep an updated list 49 of dual infeasibilities (actually squares). 50 On easy problems we don't need full scan - just 51 pick first reasonable. This method has not been coded. 52 53 One problem is how to tackle degeneracy and accuracy. At present 54 I am using the modification of costs which I put in OSL and which was 55 extended by Gill et al. I am still not sure whether we will also 56 need explicit perturbation. 57 58 The flow of primal is three while loops as follows: 59 60 while (not finished) { 61 62 while (not clean solution) { 63 64 Factorize and/or clean up solution by changing bounds so 65 primal feasible. If looks finished check fake primal bounds. 66 Repeat until status is iterating (-1) or finished (0,1,2) 67 68 } 69 70 while (status==-1) { 71 72 Iterate until no pivot in or out or time to re-factorize. 73 74 Flow is: 75 76 choose pivot column (incoming variable). if none then 77 we are primal feasible so looks as if done but we need to 78 break and check bounds etc. 79 80 Get pivot column in tableau 81 82 Choose outgoing row. If we don't find one then we look 83 primal unbounded so break and check bounds etc. (Also the 84 pivot tolerance is larger after any iterations so that may be 85 reason) 86 87 If we do find outgoing row, we may have to adjust costs to 88 keep going forwards (anti-degeneracy). Check pivot will be stable 89 and if unstable throw away iteration and break to re-factorize. 90 If minor error re-factorize after iteration. 91 92 Update everything (this may involve changing bounds on 93 variables to stay primal feasible. 94 95 } 96 97 } 98 99 TODO's (or maybe not) 100 101 At present we never check we are going forwards. I overdid that in 102 OSL so will try and make a last resort. 103 104 Needs partial scan pivot in option. 105 106 May need other anti-degeneracy measures, especially if we try and use 107 loose tolerances as a way to solve in fewer iterations. 108 109 I like idea of dynamic scaling. This gives opportunity to decouple 110 different implications of scaling for accuracy, iteration count and 111 feasibility tolerance. 112 113 for use of exotic parameter startFinishoptions see IClpSimplex.hpp 114 */ 115 116 int primal(int ifValuesPass=0, int startFinishOptions=0); 117 //@} 118 119 /**@name For advanced users */ 120 //@{ 121 /// Do not change infeasibility cost and always say optimal 122 void alwaysOptimal(bool onOff); 123 bool alwaysOptimal() const; 124 /** Normally outgoing variables can go out to slightly negative 125 values (but within tolerance) - this is to help stability and 126 and degeneracy. This can be switched off 127 */ 128 void exactOutgoing(bool onOff); 129 bool exactOutgoing() const; 130 //@} 131 132 /**@name Functions used in primal */ 133 //@{ 134 /** This has the flow between re-factorizations 135 136 Returns a code to say where decision to exit was made 137 Problem status set to: 138 139 -2 re-factorize 140 -4 Looks optimal/infeasible 141 -5 Looks unbounded 142 +3 max iterations 143 144 valuesOption has original value of valuesPass 145 */ 146 int whileIterating(int valuesOption); 147 148 /** Do last half of an iteration. This is split out so people can 149 force incoming variable. If solveType_ is 2 then this may 150 re-factorize while normally it would exit to re-factorize. 151 Return codes 152 Reasons to come out (normal mode/user mode): 153 -1 normal 154 -2 factorize now - good iteration/ NA 155 -3 slight inaccuracy - refactorize - iteration done/ same but factor done 156 -4 inaccuracy - refactorize - no iteration/ NA 157 -5 something flagged - go round again/ pivot not possible 158 +2 looks unbounded 159 +3 max iterations (iteration done) 160 161 With solveType_ ==2 this should 162 Pivot in a variable and choose an outgoing one. Assumes primal 163 feasible - will not go through a bound. Returns step length in theta 164 Returns ray in ray_ 165 */ 166 int pivotResult(int ifValuesPass=0); 167 168 169 /** The primals are updated by the given array. 170 Returns number of infeasibilities. 171 After rowArray will have cost changes for use next iteration 172 */ 173 int updatePrimalsInPrimal(CoinIndexedVector * rowArray, 174 double theta, 175 double & objectiveChange, 176 int valuesPass); 177 /** 178 Row array has pivot column 179 This chooses pivot row. 180 Rhs array is used for distance to next bound (for speed) 181 For speed, we may need to go to a bucket approach when many 182 variables go through bounds 183 If valuesPass non-zero then compute dj for direction 184 */ 185 void primalRow(CoinIndexedVector * rowArray, 186 CoinIndexedVector * rhsArray, 187 CoinIndexedVector * spareArray, 188 CoinIndexedVector * spareArray2, 189 int valuesPass); 190 /** 191 Chooses primal pivot column 192 updateArray has cost updates (also use pivotRow_ from last iteration) 193 Would be faster with separate region to scan 194 and will have this (with square of infeasibility) when steepest 195 For easy problems we can just choose one of the first columns we look at 196 */ 197 void primalColumn(CoinIndexedVector * updateArray, 198 CoinIndexedVector * spareRow1, 199 CoinIndexedVector * spareRow2, 200 CoinIndexedVector * spareColumn1, 201 CoinIndexedVector * spareColumn2); 202 203 /** Checks if tentative optimal actually means unbounded in primal 204 Returns -3 if not, 2 if is unbounded */ 205 int checkUnbounded(CoinIndexedVector * ray,CoinIndexedVector * spare, 206 double changeCost); 207 /** Refactorizes if necessary 208 Checks if finished. Updates status. 209 lastCleaned refers to iteration at which some objective/feasibility 210 cleaning too place. 211 212 type - 0 initial so set up save arrays etc 213 - 1 normal -if good update save 214 - 2 restoring from saved 215 saveModel is normally NULL but may not be if doing Sprint 216 */ 217 void statusOfProblemInPrimal(int & lastCleaned, int type, 218 ClpSimplexProgress * progress, 219 bool doFactorization, 220 int ifValuesPass, 221 IClpSimplex * saveModel=NULL); 222 /// Perturbs problem (method depends on perturbation()) 223 void perturb(int type); 224 /// Take off effect of perturbation and say whether to try dual 225 bool unPerturb(); 226 /// Unflag all variables and return number unflagged 227 int unflag(); 228 /** Get next superbasic -1 if none, 229 Normal type is 1 230 If type is 3 then initializes sorted list 231 if 2 uses list. 232 */ 233 int nextSuperBasic(int superBasicType,CoinIndexedVector * columnArray); 234 235 /// Create primal ray 236 void primalRay(CoinIndexedVector * rowArray); 237 /// Clears all bits and clears rowArray[1] etc 238 void clearAll(); 239 240 /// Sort of lexicographic resolve 241 int lexSolve(); 242 243 //@} 244 }; 245 #endif 246 247