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

Last change on this file since 2470 was 2385, checked in by unxusr, 11 months ago

formatting

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