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

Last change on this file since 1665 was 1665, checked in by lou, 10 years ago

Add EPL license notice in src.

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 14.0 KB
Line 
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
Note: See TracBrowser for help on using the repository browser.