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