source: trunk/Clp/src/ClpNonLinearCost.hpp @ 2030

Last change on this file since 2030 was 1769, checked in by forrest, 8 years ago

changes for advanced use of Clp

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 14.1 KB
Line 
1/* $Id: ClpNonLinearCost.hpp 1769 2011-07-26 09:31:51Z forrest $ */
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     /// Refresh - assuming regions OK
153     void refresh();
154     /** Sets bounds and cost for one variable
155         Returns change in cost
156      May need to be inline for speed */
157     double setOne(int sequence, double solutionValue);
158     /** Sets bounds and infeasible cost and true cost for one variable
159         This is for gub and column generation etc */
160     void setOne(int sequence, double solutionValue, double lowerValue, double upperValue,
161                 double costValue = 0.0);
162     /** Sets bounds and cost for outgoing variable
163         may change value
164         Returns direction */
165     int setOneOutgoing(int sequence, double &solutionValue);
166     /// Returns nearest bound
167     double nearest(int sequence, double solutionValue);
168     /** Returns change in cost - one down if alpha >0.0, up if <0.0
169         Value is current - new
170      */
171     inline double changeInCost(int sequence, double alpha) const {
172          double returnValue = 0.0;
173          if (CLP_METHOD1) {
174               int iRange = whichRange_[sequence] + offset_[sequence];
175               if (alpha > 0.0)
176                    returnValue = cost_[iRange] - cost_[iRange-1];
177               else
178                    returnValue = cost_[iRange] - cost_[iRange+1];
179          }
180          if (CLP_METHOD2) {
181               returnValue = (alpha > 0.0) ? infeasibilityWeight_ : -infeasibilityWeight_;
182          }
183          return returnValue;
184     }
185     inline double changeUpInCost(int sequence) const {
186          double returnValue = 0.0;
187          if (CLP_METHOD1) {
188               int iRange = whichRange_[sequence] + offset_[sequence];
189               if (iRange + 1 != start_[sequence+1] && !infeasible(iRange + 1))
190                    returnValue = cost_[iRange] - cost_[iRange+1];
191               else
192                    returnValue = -1.0e100;
193          }
194          if (CLP_METHOD2) {
195               returnValue = -infeasibilityWeight_;
196          }
197          return returnValue;
198     }
199     inline double changeDownInCost(int sequence) const {
200          double returnValue = 0.0;
201          if (CLP_METHOD1) {
202               int iRange = whichRange_[sequence] + offset_[sequence];
203               if (iRange != start_[sequence] && !infeasible(iRange - 1))
204                    returnValue = cost_[iRange] - cost_[iRange-1];
205               else
206                    returnValue = 1.0e100;
207          }
208          if (CLP_METHOD2) {
209               returnValue = infeasibilityWeight_;
210          }
211          return returnValue;
212     }
213     /// This also updates next bound
214     inline double changeInCost(int sequence, double alpha, double &rhs) {
215          double returnValue = 0.0;
216#ifdef NONLIN_DEBUG
217          double saveRhs = rhs;
218#endif
219          if (CLP_METHOD1) {
220               int iRange = whichRange_[sequence] + offset_[sequence];
221               if (alpha > 0.0) {
222                    assert(iRange - 1 >= start_[sequence]);
223                    offset_[sequence]--;
224                    rhs += lower_[iRange] - lower_[iRange-1];
225                    returnValue = alpha * (cost_[iRange] - cost_[iRange-1]);
226               } else {
227                    assert(iRange + 1 < start_[sequence+1] - 1);
228                    offset_[sequence]++;
229                    rhs += lower_[iRange+2] - lower_[iRange+1];
230                    returnValue = alpha * (cost_[iRange] - cost_[iRange+1]);
231               }
232          }
233          if (CLP_METHOD2) {
234#ifdef NONLIN_DEBUG
235               double saveRhs1 = rhs;
236               rhs = saveRhs;
237#endif
238               unsigned char iStatus = status_[sequence];
239               int iWhere = currentStatus(iStatus);
240               if (iWhere == CLP_SAME)
241                    iWhere = originalStatus(iStatus);
242               // rhs always increases
243               if (iWhere == CLP_FEASIBLE) {
244                    if (alpha > 0.0) {
245                         // going below
246                         iWhere = CLP_BELOW_LOWER;
247                         rhs = COIN_DBL_MAX;
248                    } else {
249                         // going above
250                         iWhere = CLP_ABOVE_UPPER;
251                         rhs = COIN_DBL_MAX;
252                    }
253               } else if (iWhere == CLP_BELOW_LOWER) {
254                    assert (alpha < 0);
255                    // going feasible
256                    iWhere = CLP_FEASIBLE;
257                    rhs += bound_[sequence] - model_->upperRegion()[sequence];
258               } else {
259                    assert (iWhere == CLP_ABOVE_UPPER);
260                    // going feasible
261                    iWhere = CLP_FEASIBLE;
262                    rhs += model_->lowerRegion()[sequence] - bound_[sequence];
263               }
264               setCurrentStatus(status_[sequence], iWhere);
265#ifdef NONLIN_DEBUG
266               assert(saveRhs1 == rhs);
267#endif
268               returnValue = fabs(alpha) * infeasibilityWeight_;
269          }
270          return returnValue;
271     }
272     /// Returns current lower bound
273     inline double lower(int sequence) const {
274          return lower_[whichRange_[sequence] + offset_[sequence]];
275     }
276     /// Returns current upper bound
277     inline double upper(int sequence) const {
278          return lower_[whichRange_[sequence] + offset_[sequence] + 1];
279     }
280     /// Returns current cost
281     inline double cost(int sequence) const {
282          return cost_[whichRange_[sequence] + offset_[sequence]];
283     }
284     //@}
285
286
287     /**@name Gets and sets */
288     //@{
289     /// Number of infeasibilities
290     inline int numberInfeasibilities() const {
291          return numberInfeasibilities_;
292     }
293     /// Change in cost
294     inline double changeInCost() const {
295          return changeCost_;
296     }
297     /// Feasible cost
298     inline double feasibleCost() const {
299          return feasibleCost_;
300     }
301     /// Feasible cost with offset and direction (i.e. for reporting)
302     double feasibleReportCost() const;
303     /// Sum of infeasibilities
304     inline double sumInfeasibilities() const {
305          return sumInfeasibilities_;
306     }
307     /// Largest infeasibility
308     inline double largestInfeasibility() const {
309          return largestInfeasibility_;
310     }
311     /// Average theta
312     inline double averageTheta() const {
313          return averageTheta_;
314     }
315     inline void setAverageTheta(double value) {
316          averageTheta_ = value;
317     }
318     inline void setChangeInCost(double value) {
319          changeCost_ = value;
320     }
321     inline void setMethod(int value) {
322          method_ = value;
323     }
324     /// See if may want to look both ways
325     inline bool lookBothWays() const {
326          return bothWays_;
327     }
328     //@}
329     ///@name Private functions to deal with infeasible regions
330     inline bool infeasible(int i) const {
331          return ((infeasible_[i>>5] >> (i & 31)) & 1) != 0;
332     }
333     inline void setInfeasible(int i, bool trueFalse) {
334          unsigned int & value = infeasible_[i>>5];
335          int bit = i & 31;
336          if (trueFalse)
337               value |= (1 << bit);
338          else
339               value &= ~(1 << bit);
340     }
341     inline unsigned char * statusArray() const {
342          return status_;
343     }
344     /// For debug
345     void validate();
346     //@}
347
348private:
349     /**@name Data members */
350     //@{
351     /// Change in cost because of infeasibilities
352     double changeCost_;
353     /// Feasible cost
354     double feasibleCost_;
355     /// Current infeasibility weight
356     double infeasibilityWeight_;
357     /// Largest infeasibility
358     double largestInfeasibility_;
359     /// Sum of infeasibilities
360     double sumInfeasibilities_;
361     /// Average theta - kept here as only for primal
362     double averageTheta_;
363     /// Number of rows (mainly for checking and copy)
364     int numberRows_;
365     /// Number of columns (mainly for checking and copy)
366     int numberColumns_;
367     /// Starts for each entry (columns then rows)
368     int * start_;
369     /// Range for each entry (columns then rows)
370     int * whichRange_;
371     /// Temporary range offset for each entry (columns then rows)
372     int * offset_;
373     /** Lower bound for each range (upper bound is next lower).
374         For various reasons there is always an infeasible range
375         at bottom - even if lower bound is - infinity */
376     double * lower_;
377     /// Cost for each range
378     double * cost_;
379     /// Model
380     ClpSimplex * model_;
381     // Array to say which regions are infeasible
382     unsigned int * infeasible_;
383     /// Number of infeasibilities found
384     int numberInfeasibilities_;
385     // new stuff
386     /// Contains status at beginning and current
387     unsigned char * status_;
388     /// Bound which has been replaced in lower_ or upper_
389     double * bound_;
390     /// Feasible cost array
391     double * cost2_;
392     /// Method 1 old, 2 new, 3 both!
393     int method_;
394     /// If all non-linear costs convex
395     bool convex_;
396     /// If we should look both ways for djs
397     bool bothWays_;
398     //@}
399};
400
401#endif
Note: See TracBrowser for help on using the repository browser.