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 ClpSimplexPrimal_H |
---|
11 | #define ClpSimplexPrimal_H |
---|
12 | |
---|
13 | #include "ClpSimplex.hpp" |
---|
14 | |
---|
15 | /** This solves LPs using the primal simplex method |
---|
16 | |
---|
17 | It inherits from ClpSimplex. It has no data of its own and |
---|
18 | is never created - only cast from a ClpSimplex object at algorithm time. |
---|
19 | |
---|
20 | */ |
---|
21 | |
---|
22 | class ClpSimplexPrimal : public ClpSimplex { |
---|
23 | |
---|
24 | public: |
---|
25 | |
---|
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 | |
---|
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 | /** |
---|
176 | Row array has pivot column |
---|
177 | This chooses pivot row. |
---|
178 | Rhs array is used for distance to next bound (for speed) |
---|
179 | For speed, we may need to go to a bucket approach when many |
---|
180 | variables go through bounds |
---|
181 | If valuesPass non-zero then compute dj for direction |
---|
182 | */ |
---|
183 | void primalRow(CoinIndexedVector * rowArray, |
---|
184 | CoinIndexedVector * rhsArray, |
---|
185 | CoinIndexedVector * spareArray, |
---|
186 | CoinIndexedVector * spareArray2, |
---|
187 | int valuesPass); |
---|
188 | /** |
---|
189 | Chooses primal pivot column |
---|
190 | updateArray has cost updates (also use pivotRow_ from last iteration) |
---|
191 | Would be faster with separate region to scan |
---|
192 | and will have this (with square of infeasibility) when steepest |
---|
193 | For easy problems we can just choose one of the first columns we look at |
---|
194 | */ |
---|
195 | void primalColumn(CoinIndexedVector * updateArray, |
---|
196 | CoinIndexedVector * spareRow1, |
---|
197 | CoinIndexedVector * spareRow2, |
---|
198 | CoinIndexedVector * spareColumn1, |
---|
199 | CoinIndexedVector * spareColumn2); |
---|
200 | |
---|
201 | /** Checks if tentative optimal actually means unbounded in primal |
---|
202 | Returns -3 if not, 2 if is unbounded */ |
---|
203 | int checkUnbounded(CoinIndexedVector * ray,CoinIndexedVector * spare, |
---|
204 | double changeCost); |
---|
205 | /** Refactorizes if necessary |
---|
206 | Checks if finished. Updates status. |
---|
207 | lastCleaned refers to iteration at which some objective/feasibility |
---|
208 | cleaning too place. |
---|
209 | |
---|
210 | type - 0 initial so set up save arrays etc |
---|
211 | - 1 normal -if good update save |
---|
212 | - 2 restoring from saved |
---|
213 | saveModel is normally NULL but may not be if doing Sprint |
---|
214 | */ |
---|
215 | void statusOfProblemInPrimal(int & lastCleaned, int type, |
---|
216 | ClpSimplexProgress * progress, |
---|
217 | bool doFactorization, |
---|
218 | int ifValuesPass, |
---|
219 | ClpSimplex * saveModel=NULL); |
---|
220 | /// Perturbs problem (method depends on perturbation()) |
---|
221 | void perturb(int type); |
---|
222 | /// Take off effect of perturbation and say whether to try dual |
---|
223 | bool unPerturb(); |
---|
224 | /// Unflag all variables and return number unflagged |
---|
225 | int unflag(); |
---|
226 | /** Get next superbasic -1 if none, |
---|
227 | Normal type is 1 |
---|
228 | If type is 3 then initializes sorted list |
---|
229 | if 2 uses list. |
---|
230 | */ |
---|
231 | int nextSuperBasic(int superBasicType,CoinIndexedVector * columnArray); |
---|
232 | |
---|
233 | /// Create primal ray |
---|
234 | void primalRay(CoinIndexedVector * rowArray); |
---|
235 | /// Clears all bits and clears rowArray[1] etc |
---|
236 | void clearAll(); |
---|
237 | |
---|
238 | /// Sort of lexicographic resolve |
---|
239 | int lexSolve(); |
---|
240 | |
---|
241 | //@} |
---|
242 | }; |
---|
243 | #endif |
---|
244 | |
---|