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

Last change on this file since 754 was 754, checked in by andreasw, 14 years ago

first version

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