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

Last change on this file since 2030 was 2024, checked in by forrest, 6 years ago

changes to abc

  • Property svn:keywords set to Id
File size: 8.8 KB
Line 
1/* $Id: AbcNonLinearCost.hpp 2024 2014-03-07 17:18:15Z forrest $ */
2// Copyright (C) 2002, International Business Machines
3// Corporation and others, Copyright (C) 2012, FasterCoin.  All Rights Reserved.
4// This code is licensed under the terms of the Eclipse Public License (EPL).
5
6#ifndef AbcNonLinearCost_H
7#define AbcNonLinearCost_H
8
9
10#include "CoinPragma.hpp"
11#include "AbcCommon.hpp"
12
13class AbcSimplex;
14class CoinIndexedVector;
15
16/** Trivial class to deal with non linear costs
17   
18    I don't make any explicit assumptions about convexity but I am
19    sure I do make implicit ones.
20   
21    One interesting idea for normal LP's will be to allow non-basic
22    variables to come into basis as infeasible i.e. if variable at
23    lower bound has very large positive reduced cost (when problem
24    is infeasible) could it reduce overall problem infeasibility more
25    by bringing it into basis below its lower bound.
26   
27    Another feature would be to automatically discover when problems
28    are convex piecewise linear and re-formulate to use non-linear.
29    I did some work on this many years ago on "grade" problems, but
30    while it improved primal interior point algorithms were much better
31    for that particular problem.
32*/
33/* status has original status and current status
34   0 - below lower so stored is upper
35   1 - in range
36   2 - above upper so stored is lower
37   4 - (for current) - same as original
38*/
39#define CLP_BELOW_LOWER 0
40#define CLP_FEASIBLE 1
41#define CLP_ABOVE_UPPER 2
42#define CLP_SAME 4
43#ifndef ClpNonLinearCost_H
44inline int originalStatus(unsigned char status)
45{
46  return (status & 15);
47}
48inline int currentStatus(unsigned char status)
49{
50  return (status >> 4);
51}
52inline void setOriginalStatus(unsigned char & status, int value)
53{
54  status = static_cast<unsigned char>(status & ~15);
55  status = static_cast<unsigned char>(status | value);
56}
57inline void setCurrentStatus(unsigned char &status, int value)
58{
59  status = static_cast<unsigned char>(status & ~(15 << 4));
60  status = static_cast<unsigned char>(status | (value << 4));
61}
62inline void setInitialStatus(unsigned char &status)
63{
64  status = static_cast<unsigned char>(CLP_FEASIBLE | (CLP_SAME << 4));
65}
66inline void setSameStatus(unsigned char &status)
67{
68  status = static_cast<unsigned char>(status & ~(15 << 4));
69  status = static_cast<unsigned char>(status | (CLP_SAME << 4));
70}
71#endif
72class AbcNonLinearCost  {
73 
74public:
75 
76  /**@name Constructors, destructor */
77  //@{
78  /// Default constructor.
79  AbcNonLinearCost();
80  /** Constructor from simplex.
81      This will just set up wasteful arrays for linear, but
82      later may do dual analysis and even finding duplicate columns .
83  */
84  AbcNonLinearCost(AbcSimplex * model);
85  /// Destructor
86  ~AbcNonLinearCost();
87  // Copy
88  AbcNonLinearCost(const AbcNonLinearCost&);
89  // Assignment
90  AbcNonLinearCost& operator=(const AbcNonLinearCost&);
91  //@}
92 
93 
94  /**@name Actual work in primal */
95  //@{
96  /** Changes infeasible costs and computes number and cost of infeas
97      Puts all non-basic (non free) variables to bounds
98      and all free variables to zero if oldTolerance is non-zero
99      - but does not move those <= oldTolerance away*/
100  void checkInfeasibilities(double oldTolerance = 0.0);
101  /** Changes infeasible costs for each variable
102      The indices are row indices and need converting to sequences
103  */
104  void checkInfeasibilities(int numberInArray, const int * index);
105  /** Puts back correct infeasible costs for each variable
106      The input indices are row indices and need converting to sequences
107      for costs.
108      On input array is empty (but indices exist).  On exit just
109      changed costs will be stored as normal CoinIndexedVector
110  */
111  void checkChanged(int numberInArray, CoinIndexedVector * update);
112  /** Goes through one bound for each variable.
113      If multiplier*work[iRow]>0 goes down, otherwise up.
114      The indices are row indices and need converting to sequences
115      Temporary offsets may be set
116      Rhs entries are increased
117  */
118  void goThru(int numberInArray, double multiplier,
119              const int * index, const double * work,
120              double * rhs);
121  /** Takes off last iteration (i.e. offsets closer to 0)
122   */
123  void goBack(int numberInArray, const int * index,
124              double * rhs);
125  /** Puts back correct infeasible costs for each variable
126      The input indices are row indices and need converting to sequences
127      for costs.
128      At the end of this all temporary offsets are zero
129  */
130  void goBackAll(const CoinIndexedVector * update);
131  /// Temporary zeroing of feasible costs
132  void zapCosts();
133  /// Refreshes costs always makes row costs zero
134  void refreshCosts(const double * columnCosts);
135  /// Puts feasible bounds into lower and upper
136  void feasibleBounds();
137  /// Refresh - assuming regions OK
138  void refresh();
139  /// Refresh - from original
140  void refreshFromPerturbed(double tolerance);
141  /** Sets bounds and cost for one variable
142      Returns change in cost
143      May need to be inline for speed */
144  double setOne(int sequence, double solutionValue);
145  /** Sets bounds and cost for one variable
146      Returns change in cost
147      May need to be inline for speed */
148  double setOneBasic(int iRow, double solutionValue);
149  /** Sets bounds and cost for outgoing variable
150      may change value
151      Returns direction */
152  int setOneOutgoing(int sequence, double &solutionValue);
153  /// Returns nearest bound
154  double nearest(int iRow, double solutionValue);
155  /** Returns change in cost - one down if alpha >0.0, up if <0.0
156      Value is current - new
157  */
158  inline double changeInCost(int /*sequence*/, double alpha) const {
159    return (alpha > 0.0) ? infeasibilityWeight_ : -infeasibilityWeight_;
160  }
161  inline double changeUpInCost(int /*sequence*/) const {
162    return -infeasibilityWeight_;
163  }
164  inline double changeDownInCost(int /*sequence*/) const {
165    return infeasibilityWeight_;
166  }
167  /// This also updates next bound
168  inline double changeInCost(int iRow, double alpha, double &rhs) {
169    int sequence=model_->pivotVariable()[iRow];
170    double returnValue = 0.0;
171    unsigned char iStatus = status_[sequence];
172    int iWhere = currentStatus(iStatus);
173    if (iWhere == CLP_SAME)
174      iWhere = originalStatus(iStatus);
175    // rhs always increases
176    if (iWhere == CLP_FEASIBLE) {
177      if (alpha > 0.0) {
178        // going below
179        iWhere = CLP_BELOW_LOWER;
180        rhs = COIN_DBL_MAX;
181      } else {
182        // going above
183        iWhere = CLP_ABOVE_UPPER;
184        rhs = COIN_DBL_MAX;
185      }
186    } else if (iWhere == CLP_BELOW_LOWER) {
187      assert (alpha < 0);
188      // going feasible
189      iWhere = CLP_FEASIBLE;
190      rhs += bound_[sequence] - model_->upperRegion()[sequence];
191    } else {
192      assert (iWhere == CLP_ABOVE_UPPER);
193      // going feasible
194      iWhere = CLP_FEASIBLE;
195      rhs += model_->lowerRegion()[sequence] - bound_[sequence];
196    }
197    setCurrentStatus(status_[sequence], iWhere);
198    returnValue = fabs(alpha) * infeasibilityWeight_;
199    return returnValue;
200  }
201  //@}
202 
203 
204  /**@name Gets and sets */
205  //@{
206  /// Number of infeasibilities
207  inline int numberInfeasibilities() const {
208    return numberInfeasibilities_;
209  }
210  /// Change in cost
211  inline double changeInCost() const {
212    return changeCost_;
213  }
214  /// Feasible cost
215  inline double feasibleCost() const {
216    return feasibleCost_;
217  }
218  /// Feasible cost with offset and direction (i.e. for reporting)
219  double feasibleReportCost() const;
220  /// Sum of infeasibilities
221  inline double sumInfeasibilities() const {
222    return sumInfeasibilities_;
223  }
224  /// Largest infeasibility
225  inline double largestInfeasibility() const {
226    return largestInfeasibility_;
227  }
228  /// Average theta
229  inline double averageTheta() const {
230    return averageTheta_;
231  }
232  inline void setAverageTheta(double value) {
233    averageTheta_ = value;
234  }
235  inline void setChangeInCost(double value) {
236    changeCost_ = value;
237  }
238  //@}
239  ///@name Private functions to deal with infeasible regions
240  inline unsigned char * statusArray() const {
241    return status_;
242  }
243  inline int getCurrentStatus(int sequence)
244  {return (status_[sequence] >> 4);}
245  /// For debug
246  void validate();
247  //@}
248 
249private:
250  /**@name Data members */
251  //@{
252  /// Change in cost because of infeasibilities
253  double changeCost_;
254  /// Feasible cost
255  double feasibleCost_;
256  /// Current infeasibility weight
257  double infeasibilityWeight_;
258  /// Largest infeasibility
259  double largestInfeasibility_;
260  /// Sum of infeasibilities
261  double sumInfeasibilities_;
262  /// Average theta - kept here as only for primal
263  double averageTheta_;
264  /// Number of rows (mainly for checking and copy)
265  int numberRows_;
266  /// Number of columns (mainly for checking and copy)
267  int numberColumns_;
268  /// Model
269  AbcSimplex * model_;
270  /// Number of infeasibilities found
271  int numberInfeasibilities_;
272  // new stuff
273  /// Contains status at beginning and current
274  unsigned char * status_;
275  /// Bound which has been replaced in lower_ or upper_
276  double * bound_;
277  /// Feasible cost array
278  double * cost_;
279  //@}
280};
281
282#endif
Note: See TracBrowser for help on using the repository browser.