1/* $Id: ClpNonLinearCost.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#ifndef ClpNonLinearCost_H
7#define ClpNonLinearCost_H
8
9
10#include "CoinPragma.hpp"
11
12class ClpSimplex;
13class CoinIndexedVector;
14
15/** Trivial class to deal with non linear costs
16
17 I don't make any explicit assumptions about convexity but I am
18 sure I do make implicit ones.
19
20 One interesting idea for normal LP's will be to allow non-basic
21 variables to come into basis as infeasible i.e. if variable at
22 lower bound has very large positive reduced cost (when problem
23 is infeasible) could it reduce overall problem infeasibility more
24 by bringing it into basis below its lower bound.
25
26 Another feature would be to automatically discover when problems
27 are convex piecewise linear and re-formulate to use non-linear.
28 I did some work on this many years ago on "grade" problems, but
29 while it improved primal interior point algorithms were much better
30 for that particular problem.
31*/
32/* status has original status and current status
33 0 - below lower so stored is upper
34 1 - in range
35 2 - above upper so stored is lower
36 4 - (for current) - same as original
37*/
38#define CLP_BELOW_LOWER 0
39#define CLP_FEASIBLE 1
40#define CLP_ABOVE_UPPER 2
41#define CLP_SAME 4
42inline int originalStatus(unsigned char status)
43{
44 return (status & 15);
45}
46inline int currentStatus(unsigned char status)
47{
48 return (status >> 4);
49}
50inline void setOriginalStatus(unsigned char & status, int value)
51{
52 status = static_cast<unsigned char>(status & ~15);
53 status = static_cast<unsigned char>(status | value);
54}
55inline void setCurrentStatus(unsigned char &status, int value)
56{
57 status = static_cast<unsigned char>(status & ~(15 << 4));
58 status = static_cast<unsigned char>(status | (value << 4));
59}
60inline void setInitialStatus(unsigned char &status)
61{
62 status = static_cast<unsigned char>(CLP_FEASIBLE | (CLP_SAME << 4));
63}
64inline void setSameStatus(unsigned char &status)
65{
66 status = static_cast<unsigned char>(status & ~(15 << 4));
67 status = static_cast<unsigned char>(status | (CLP_SAME << 4));
68}
69// Use second version to get more speed
70//#define FAST_CLPNON
71#ifndef FAST_CLPNON
72#define CLP_METHOD1 ((method_&1)!=0)
73#define CLP_METHOD2 ((method_&2)!=0)
74#else
75#define CLP_METHOD1 (false)
76#define CLP_METHOD2 (true)
77#endif
78class ClpNonLinearCost {
79
80public:
81
82public:
83
84 /**@name Constructors, destructor */
85 //@{
86 /// Default constructor.
87 ClpNonLinearCost();
88 /** Constructor from simplex.
89 This will just set up wasteful arrays for linear, but
90 later may do dual analysis and even finding duplicate columns .
91 */
92 ClpNonLinearCost(ClpSimplex * model, int method = 1);
93 /** Constructor from simplex and list of non-linearities (columns only)
94 First lower of each column has to match real lower
95 Last lower has to be <= upper (if == then cost ignored)
96 This could obviously be changed to make more user friendly
97 */
98 ClpNonLinearCost(ClpSimplex * model, const int * starts,
99 const double * lower, const double * cost);
100 /// Destructor
101 ~ClpNonLinearCost();
102 // Copy
103 ClpNonLinearCost(const ClpNonLinearCost&);
104 // Assignment
105 ClpNonLinearCost& operator=(const ClpNonLinearCost&);
106 //@}
107
108
109 /**@name Actual work in primal */
110 //@{
111 /** Changes infeasible costs and computes number and cost of infeas
112 Puts all non-basic (non free) variables to bounds
113 and all free variables to zero if oldTolerance is non-zero
114 - but does not move those <= oldTolerance away*/
115 void checkInfeasibilities(double oldTolerance = 0.0);
116 /** Changes infeasible costs for each variable
117 The indices are row indices and need converting to sequences
118 */
119 void checkInfeasibilities(int numberInArray, const int * index);
120 /** Puts back correct infeasible costs for each variable
121 The input indices are row indices and need converting to sequences
122 for costs.
123 On input array is empty (but indices exist). On exit just
124 changed costs will be stored as normal CoinIndexedVector
125 */
126 void checkChanged(int numberInArray, CoinIndexedVector * update);
127 /** Goes through one bound for each variable.
128 If multiplier*work[iRow]>0 goes down, otherwise up.
129 The indices are row indices and need converting to sequences
130 Temporary offsets may be set
131 Rhs entries are increased
132 */
133 void goThru(int numberInArray, double multiplier,
134 const int * index, const double * work,
135 double * rhs);
136 /** Takes off last iteration (i.e. offsets closer to 0)
137 */
138 void goBack(int numberInArray, const int * index,
139 double * rhs);
140 /** Puts back correct infeasible costs for each variable
141 The input indices are row indices and need converting to sequences
142 for costs.
143 At the end of this all temporary offsets are zero
144 */
145 void goBackAll(const CoinIndexedVector * update);
146 /// Temporary zeroing of feasible costs
147 void zapCosts();
148 /// Refreshes costs always makes row costs zero
149 void refreshCosts(const double * columnCosts);
150 /// Puts feasible bounds into lower and upper
151 void feasibleBounds();
152 /** Sets bounds and cost for one variable
153 Returns change in cost
154 May need to be inline for speed */
155 double setOne(int sequence, double solutionValue);
156 /** Sets bounds and infeasible cost and true cost for one variable
157 This is for gub and column generation etc */
158 void setOne(int sequence, double solutionValue, double lowerValue, double upperValue,
159 double costValue = 0.0);
160 /** Sets bounds and cost for outgoing variable
161 may change value
162 Returns direction */
163 int setOneOutgoing(int sequence, double &solutionValue);
164 /// Returns nearest bound
165 double nearest(int sequence, double solutionValue);
166 /** Returns change in cost - one down if alpha >0.0, up if <0.0
167 Value is current - new
168 */
169 inline double changeInCost(int sequence, double alpha) const {
170 double returnValue = 0.0;
171 if (CLP_METHOD1) {
172 int iRange = whichRange_[sequence] + offset_[sequence];
173 if (alpha > 0.0)
174 returnValue = cost_[iRange] - cost_[iRange-1];
175 else
176 returnValue = cost_[iRange] - cost_[iRange+1];
177 }
178 if (CLP_METHOD2) {
179 returnValue = (alpha > 0.0) ? infeasibilityWeight_ : -infeasibilityWeight_;
180 }
181 return returnValue;
182 }
183 inline double changeUpInCost(int sequence) const {
184 double returnValue = 0.0;
185 if (CLP_METHOD1) {
186 int iRange = whichRange_[sequence] + offset_[sequence];
187 if (iRange + 1 != start_[sequence+1] && !infeasible(iRange + 1))
188 returnValue = cost_[iRange] - cost_[iRange+1];
189 else
190 returnValue = -1.0e100;
191 }
192 if (CLP_METHOD2) {
193 returnValue = -infeasibilityWeight_;
194 }
195 return returnValue;
196 }
197 inline double changeDownInCost(int sequence) const {
198 double returnValue = 0.0;
199 if (CLP_METHOD1) {
200 int iRange = whichRange_[sequence] + offset_[sequence];
201 if (iRange != start_[sequence] && !infeasible(iRange - 1))
202 returnValue = cost_[iRange] - cost_[iRange-1];
203 else
204 returnValue = 1.0e100;
205 }
206 if (CLP_METHOD2) {
207 returnValue = infeasibilityWeight_;
208 }
209 return returnValue;
210 }
211 /// This also updates next bound
212 inline double changeInCost(int sequence, double alpha, double &rhs) {
213 double returnValue = 0.0;
214#ifdef NONLIN_DEBUG
215 double saveRhs = rhs;
216#endif
217 if (CLP_METHOD1) {
218 int iRange = whichRange_[sequence] + offset_[sequence];
219 if (alpha > 0.0) {
220 assert(iRange - 1 >= start_[sequence]);
221 offset_[sequence]--;
222 rhs += lower_[iRange] - lower_[iRange-1];
223 returnValue = alpha * (cost_[iRange] - cost_[iRange-1]);
224 } else {
225 assert(iRange + 1 < start_[sequence+1] - 1);
226 offset_[sequence]++;
227 rhs += lower_[iRange+2] - lower_[iRange+1];
228 returnValue = alpha * (cost_[iRange] - cost_[iRange+1]);
229 }
230 }
231 if (CLP_METHOD2) {
232#ifdef NONLIN_DEBUG
233 double saveRhs1 = rhs;
234 rhs = saveRhs;
235#endif
236 unsigned char iStatus = status_[sequence];
237 int iWhere = currentStatus(iStatus);
238 if (iWhere == CLP_SAME)
239 iWhere = originalStatus(iStatus);
240 // rhs always increases
241 if (iWhere == CLP_FEASIBLE) {
242 if (alpha > 0.0) {
243 // going below
244 iWhere = CLP_BELOW_LOWER;
245 rhs = COIN_DBL_MAX;
246 } else {
247 // going above
248 iWhere = CLP_ABOVE_UPPER;
249 rhs = COIN_DBL_MAX;
250 }
251 } else if (iWhere == CLP_BELOW_LOWER) {
252 assert (alpha < 0);
253 // going feasible
254 iWhere = CLP_FEASIBLE;
255 rhs += bound_[sequence] - model_->upperRegion()[sequence];
256 } else {
257 assert (iWhere == CLP_ABOVE_UPPER);
258 // going feasible
259 iWhere = CLP_FEASIBLE;
260 rhs += model_->lowerRegion()[sequence] - bound_[sequence];
261 }
262 setCurrentStatus(status_[sequence], iWhere);
263#ifdef NONLIN_DEBUG
264 assert(saveRhs1 == rhs);
265#endif
266 returnValue = fabs(alpha) * infeasibilityWeight_;
267 }
268 return returnValue;
269 }
270 /// Returns current lower bound
271 inline double lower(int sequence) const {
272 return lower_[whichRange_[sequence] + offset_[sequence]];
273 }
274 /// Returns current upper bound
275 inline double upper(int sequence) const {
276 return lower_[whichRange_[sequence] + offset_[sequence] + 1];
277 }
278 /// Returns current cost
279 inline double cost(int sequence) const {
280 return cost_[whichRange_[sequence] + offset_[sequence]];
281 }
282 //@}
283
284
285 /**@name Gets and sets */
286 //@{
287 /// Number of infeasibilities
288 inline int numberInfeasibilities() const {
289 return numberInfeasibilities_;
290 }
291 /// Change in cost
292 inline double changeInCost() const {
293 return changeCost_;
294 }
295 /// Feasible cost
296 inline double feasibleCost() const {
297 return feasibleCost_;
298 }
299 /// Feasible cost with offset and direction (i.e. for reporting)
300 double feasibleReportCost() const;
301 /// Sum of infeasibilities
302 inline double sumInfeasibilities() const {
303 return sumInfeasibilities_;
304 }
305 /// Largest infeasibility
306 inline double largestInfeasibility() const {
307 return largestInfeasibility_;
308 }
309 /// Average theta
310 inline double averageTheta() const {
311 return averageTheta_;
312 }
313 inline void setAverageTheta(double value) {
314 averageTheta_ = value;
315 }
316 inline void setChangeInCost(double value) {
317 changeCost_ = value;
318 }
319 inline void setMethod(int value) {
320 method_ = value;
321 }
322 /// See if may want to look both ways
323 inline bool lookBothWays() const {
324 return bothWays_;
325 }
326 //@}
327 ///@name Private functions to deal with infeasible regions
328 inline bool infeasible(int i) const {
329 return ((infeasible_[i>>5] >> (i & 31)) & 1) != 0;
330 }
331 inline void setInfeasible(int i, bool trueFalse) {
332 unsigned int & value = infeasible_[i>>5];
333 int bit = i & 31;
334 if (trueFalse)
335 value |= (1 << bit);
336 else
337 value &= ~(1 << bit);
338 }
339 inline unsigned char * statusArray() const {
340 return status_;
341 }
342 /// For debug
343 void validate();
344 //@}
345
346private:
347 /**@name Data members */
348 //@{
349 /// Change in cost because of infeasibilities
350 double changeCost_;
351 /// Feasible cost
352 double feasibleCost_;
353 /// Current infeasibility weight
354 double infeasibilityWeight_;
355 /// Largest infeasibility
356 double largestInfeasibility_;
357 /// Sum of infeasibilities
358 double sumInfeasibilities_;
359 /// Average theta - kept here as only for primal
360 double averageTheta_;
361 /// Number of rows (mainly for checking and copy)
362 int numberRows_;
363 /// Number of columns (mainly for checking and copy)
364 int numberColumns_;
365 /// Starts for each entry (columns then rows)
366 int * start_;
367 /// Range for each entry (columns then rows)
368 int * whichRange_;
369 /// Temporary range offset for each entry (columns then rows)
370 int * offset_;
371 /** Lower bound for each range (upper bound is next lower).
372 For various reasons there is always an infeasible range
373 at bottom - even if lower bound is - infinity */
374 double * lower_;
375 /// Cost for each range
376 double * cost_;
377 /// Model
378 ClpSimplex * model_;
379 // Array to say which regions are infeasible
380 unsigned int * infeasible_;
381 /// Number of infeasibilities found
382 int numberInfeasibilities_;
383 // new stuff
384 /// Contains status at beginning and current
385 unsigned char * status_;
386 /// Bound which has been replaced in lower_ or upper_
387 double * bound_;
388 /// Feasible cost array
389 double * cost2_;
390 /// Method 1 old, 2 new, 3 both!
391 int method_;
392 /// If all non-linear costs convex
393 bool convex_;
394 /// If we should look both ways for djs
395 bool bothWays_;
396 //@}
397};
398
399#endif
400