1/* $Id: ClpSimplexDual.hpp 1665 2011-01-04 17:55:54Z lou $ */
2// Copyright (C) 2002, International Business Machines
3// Corporation and others. 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 ClpSimplexDual_H
12#define ClpSimplexDual_H
13
14#include "ClpSimplex.hpp"
15
16/** This solves LPs using the dual simplex method
17
18 It inherits from ClpSimplex. It has no data of its own and
19 is never created - only cast from a ClpSimplex object at algorithm time.
20
21*/
22
23class ClpSimplexDual : public ClpSimplex {
24
25public:
26
27 /**@name Description of algorithm */
28 //@{
29 /** Dual algorithm
30
31 Method
32
33 It tries to be a single phase approach with a weight of 1.0 being
34 given to getting optimal and a weight of updatedDualBound_ being
35 given to getting dual feasible. In this version I have used the
36 idea that this weight can be thought of as a fake bound. If the
37 distance between the lower and upper bounds on a variable is less
38 than the feasibility weight then we are always better off flipping
39 to other bound to make dual feasible. If the distance is greater
40 then we make up a fake bound updatedDualBound_ away from one bound.
41 If we end up optimal or primal infeasible, we check to see if
42 bounds okay. If so we have finished, if not we increase updatedDualBound_
43 and continue (after checking if unbounded). I am undecided about
44 free variables - there is coding but I am not sure about it. At
45 present I put them in basis anyway.
46
47 The code is designed to take advantage of sparsity so arrays are
48 seldom zeroed out from scratch or gone over in their entirety.
49 The only exception is a full scan to find outgoing variable for
50 Dantzig row choice. For steepest edge we keep an updated list
51 of infeasibilities (actually squares).
52 On easy problems we don't need full scan - just
53 pick first reasonable.
54
55 One problem is how to tackle degeneracy and accuracy. At present
56 I am using the modification of costs which I put in OSL and some
57 of what I think is the dual analog of Gill et al.
58 I am still not sure of the exact details.
59
60 The flow of dual is three while loops as follows:
61
62 while (not finished) {
63
64 while (not clean solution) {
65
66 Factorize and/or clean up solution by flipping variables so
67 dual feasible. If looks finished check fake dual bounds.
68 Repeat until status is iterating (-1) or finished (0,1,2)
69
70 }
71
72 while (status==-1) {
73
74 Iterate until no pivot in or out or time to re-factorize.
75
76 Flow is:
77
78 choose pivot row (outgoing variable). if none then
79 we are primal feasible so looks as if done but we need to
80 break and check bounds etc.
81
82 Get pivot row in tableau
83
84 Choose incoming column. If we don't find one then we look
85 primal infeasible so break and check bounds etc. (Also the
86 pivot tolerance is larger after any iterations so that may be
87 reason)
88
89 If we do find incoming column, we may have to adjust costs to
90 keep going forwards (anti-degeneracy). Check pivot will be stable
91 and if unstable throw away iteration and break to re-factorize.
92 If minor error re-factorize after iteration.
93
94 Update everything (this may involve flipping variables to stay
95 dual feasible.
96
97 }
98
99 }
100
101 TODO's (or maybe not)
102
103 At present we never check we are going forwards. I overdid that in
104 OSL so will try and make a last resort.
105
106 Needs partial scan pivot out option.
107
108 May need other anti-degeneracy measures, especially if we try and use
109 loose tolerances as a way to solve in fewer iterations.
110
111 I like idea of dynamic scaling. This gives opportunity to decouple
112 different implications of scaling for accuracy, iteration count and
113 feasibility tolerance.
114
115 for use of exotic parameter startFinishoptions see Clpsimplex.hpp
116 */
117
118 int dual(int ifValuesPass, int startFinishOptions = 0);
119 /** For strong branching. On input lower and upper are new bounds
120 while on output they are change in objective function values
121 (>1.0e50 infeasible).
122 Return code is 0 if nothing interesting, -1 if infeasible both
123 ways and +1 if infeasible one way (check values to see which one(s))
124 Solutions are filled in as well - even down, odd up - also
125 status and number of iterations
126 */
127 int strongBranching(int numberVariables, const int * variables,
128 double * newLower, double * newUpper,
129 double ** outputSolution,
130 int * outputStatus, int * outputIterations,
131 bool stopOnFirstInfeasible = true,
132 bool alwaysFinish = false,
133 int startFinishOptions = 0);
134 /// This does first part of StrongBranching
135 ClpFactorization * setupForStrongBranching(char * arrays, int numberRows,
136 int numberColumns, bool solveLp = false);
137 /// This cleans up after strong branching
138 void cleanupAfterStrongBranching(ClpFactorization * factorization);
139 //@}
140
141 /**@name Functions used in dual */
142 //@{
143 /** This has the flow between re-factorizations
144 Broken out for clarity and will be used by strong branching
145
146 Reasons to come out:
147 -1 iterations etc
148 -2 inaccuracy
149 -3 slight inaccuracy (and done iterations)
150 +0 looks optimal (might be unbounded - but we will investigate)
151 +1 looks infeasible
152 +3 max iterations
153
154 If givenPi not NULL then in values pass
155 */
156 int whileIterating(double * & givenPi, int ifValuesPass);
157 /** The duals are updated by the given arrays.
158 Returns number of infeasibilities.
159 After rowArray and columnArray will just have those which
160 have been flipped.
161 Variables may be flipped between bounds to stay dual feasible.
162 The output vector has movement of primal
163 solution (row length array) */
164 int updateDualsInDual(CoinIndexedVector * rowArray,
165 CoinIndexedVector * columnArray,
166 CoinIndexedVector * outputArray,
167 double theta,
168 double & objectiveChange,
169 bool fullRecompute);
170 /** The duals are updated by the given arrays.
171 This is in values pass - so no changes to primal is made
172 */
173 void updateDualsInValuesPass(CoinIndexedVector * rowArray,
174 CoinIndexedVector * columnArray,
175 double theta);
176 /** While updateDualsInDual sees what effect is of flip
177 this does actual flipping.
178 */
179 void flipBounds(CoinIndexedVector * rowArray,
180 CoinIndexedVector * columnArray);
181 /**
182 Row array has row part of pivot row
183 Column array has column part.
184 This chooses pivot column.
185 Spare arrays are used to save pivots which will go infeasible
186 We will check for basic so spare array will never overflow.
187 If necessary will modify costs
188 For speed, we may need to go to a bucket approach when many
189 variables are being flipped.
190 Returns best possible pivot value
191 */
192 double dualColumn(CoinIndexedVector * rowArray,
193 CoinIndexedVector * columnArray,
194 CoinIndexedVector * spareArray,
195 CoinIndexedVector * spareArray2,
196 double accpetablePivot,
197 CoinBigIndex * dubiousWeights);
198 /// Does first bit of dualColumn
199 int dualColumn0(const CoinIndexedVector * rowArray,
200 const CoinIndexedVector * columnArray,
201 CoinIndexedVector * spareArray,
202 double acceptablePivot,
203 double & upperReturn, double &bestReturn, double & badFree);
204 /**
205 Row array has row part of pivot row
206 Column array has column part.
207 This sees what is best thing to do in dual values pass
208 if sequenceIn==sequenceOut can change dual on chosen row and leave variable in basis
209 */
210 void checkPossibleValuesMove(CoinIndexedVector * rowArray,
211 CoinIndexedVector * columnArray,
212 double acceptablePivot);
213 /**
214 Row array has row part of pivot row
215 Column array has column part.
216 This sees what is best thing to do in branch and bound cleanup
217 If sequenceIn_ < 0 then can't do anything
218 */
219 void checkPossibleCleanup(CoinIndexedVector * rowArray,
220 CoinIndexedVector * columnArray,
221 double acceptablePivot);
222 /**
223 This sees if we can move duals in dual values pass.
224 This is done before any pivoting
225 */
226 void doEasyOnesInValuesPass(double * givenReducedCosts);
227 /**
228 Chooses dual pivot row
229 Would be faster with separate region to scan
230 and will have this (with square of infeasibility) when steepest
231 For easy problems we can just choose one of the first rows we look at
232
233 If alreadyChosen >=0 then in values pass and that row has been
234 selected
235 */
236 void dualRow(int alreadyChosen);
237 /** Checks if any fake bounds active - if so returns number and modifies
238 updatedDualBound_ and everything.
239 Free variables will be left as free
240 Returns number of bounds changed if >=0
241 Returns -1 if not initialize and no effect
242 Fills in changeVector which can be used to see if unbounded
243 and cost of change vector
244 If 2 sets to original (just changed)
245 */
246 int changeBounds(int initialize, CoinIndexedVector * outputArray,
247 double & changeCost);
248 /** As changeBounds but just changes new bounds for a single variable.
249 Returns true if change */
250 bool changeBound( int iSequence);
251 /// Restores bound to original bound
252 void originalBound(int iSequence);
253 /** Checks if tentative optimal actually means unbounded in dual
254 Returns -3 if not, 2 if is unbounded */
255 int checkUnbounded(CoinIndexedVector * ray, CoinIndexedVector * spare,
256 double changeCost);
257 /** Refactorizes if necessary
258 Checks if finished. Updates status.
259 lastCleaned refers to iteration at which some objective/feasibility
260 cleaning too place.
261
262 type - 0 initial so set up save arrays etc
263 - 1 normal -if good update save
264 - 2 restoring from saved
265 */
266 void statusOfProblemInDual(int & lastCleaned, int type,
267 double * givenDjs, ClpDataSave & saveData,
268 int ifValuesPass);
269 /** Perturbs problem (method depends on perturbation())
270 returns nonzero if should go to dual */
271 int perturb();
272 /** Fast iterations. Misses out a lot of initialization.
273 Normally stops on maximum iterations, first re-factorization
274 or tentative optimum. If looks interesting then continues as
275 normal. Returns 0 if finished properly, 1 otherwise.
276 */
277 int fastDual(bool alwaysFinish = false);
278 /** Checks number of variables at fake bounds. This is used by fastDual
279 so can exit gracefully before end */
280 int numberAtFakeBound();
281
282 /** Pivot in a variable and choose an outgoing one. Assumes dual
283 feasible - will not go through a reduced cost. Returns step length in theta
284 Returns ray in ray_ (or NULL if no pivot)
285 Return codes as before but -1 means no acceptable pivot
286 */
287 int pivotResult();
288 /** Get next free , -1 if none */
289 int nextSuperBasic();
290 /** Startup part of dual (may be extended to other algorithms)
291 returns 0 if good, 1 if bad */
292 int startupSolve(int ifValuesPass, double * saveDuals, int startFinishOptions);
293 void finishSolve(int startFinishOptions);
294 void gutsOfDual(int ifValuesPass, double * & saveDuals, int initialStatus,
295 ClpDataSave & saveData);
296 //int dual2(int ifValuesPass,int startFinishOptions=0);
297 void resetFakeBounds(int type);
298
299 //@}
300};
301#endif
302