1 | using System;
|
---|
2 | using System.Diagnostics;
|
---|
3 | using System.Linq;
|
---|
4 |
|
---|
5 | namespace HeuristicLab.Problems.DataAnalysis.Symbolic {
|
---|
6 | // this is our own implementation of interval arithmetic
|
---|
7 | // for a well worked out definition of interval operations for IEEE reals see:
|
---|
8 | // Stahl: Interval Methods for Bounding the Range of Polynomials and Solving Systems of Nonlinear Equations, Dissertation, JKU, 1995
|
---|
9 | [DebuggerDisplay("[{low.Value}..{high.Value}]")]
|
---|
10 | public class AlgebraicInterval : IAlgebraicType<AlgebraicInterval> {
|
---|
11 | [DebuggerBrowsable(DebuggerBrowsableState.Never)]
|
---|
12 | private MultivariateDual<AlgebraicDouble> low;
|
---|
13 | public MultivariateDual<AlgebraicDouble> LowerBound => low.Clone();
|
---|
14 |
|
---|
15 | [DebuggerBrowsable(DebuggerBrowsableState.Never)]
|
---|
16 | private MultivariateDual<AlgebraicDouble> high;
|
---|
17 | public MultivariateDual<AlgebraicDouble> UpperBound => high.Clone();
|
---|
18 |
|
---|
19 |
|
---|
20 | public AlgebraicInterval() : this(double.NegativeInfinity, double.PositiveInfinity) { }
|
---|
21 |
|
---|
22 | public AlgebraicInterval(MultivariateDual<AlgebraicDouble> low, MultivariateDual<AlgebraicDouble> high) {
|
---|
23 | this.low = low.Clone();
|
---|
24 | this.high = high.Clone();
|
---|
25 | }
|
---|
26 |
|
---|
27 | public AlgebraicInterval(double low, double high) {
|
---|
28 | this.low = new MultivariateDual<AlgebraicDouble>(new AlgebraicDouble(low));
|
---|
29 | this.high = new MultivariateDual<AlgebraicDouble>(new AlgebraicDouble(high));
|
---|
30 | }
|
---|
31 |
|
---|
32 | [DebuggerBrowsable(DebuggerBrowsableState.Never)]
|
---|
33 | public AlgebraicInterval Zero => new AlgebraicInterval(0.0, 0.0);
|
---|
34 | [DebuggerBrowsable(DebuggerBrowsableState.Never)]
|
---|
35 | public AlgebraicInterval One => new AlgebraicInterval(1.0, 1.0);
|
---|
36 |
|
---|
37 | public AlgebraicInterval Add(AlgebraicInterval a) {
|
---|
38 | low.Add(a.low);
|
---|
39 | high.Add(a.high);
|
---|
40 | return this;
|
---|
41 | }
|
---|
42 |
|
---|
43 | public AlgebraicInterval Mul(AlgebraicInterval a) {
|
---|
44 | var v1 = low.Clone().Mul(a.low);
|
---|
45 | var v2 = low.Clone().Mul(a.high);
|
---|
46 | var v3 = high.Clone().Mul(a.low);
|
---|
47 | var v4 = high.Clone().Mul(a.high);
|
---|
48 |
|
---|
49 | AssignLowAndHigh(v1, v2, v3, v4);
|
---|
50 |
|
---|
51 | return this;
|
---|
52 | }
|
---|
53 |
|
---|
54 | public AlgebraicInterval Assign(AlgebraicInterval a) {
|
---|
55 | low = a.low;
|
---|
56 | high = a.high;
|
---|
57 | return this;
|
---|
58 | }
|
---|
59 |
|
---|
60 | public AlgebraicInterval AssignCos(AlgebraicInterval a) {
|
---|
61 | return AssignSin(a.Clone().Add(new AlgebraicInterval(Math.PI / 2, Math.PI / 2)));
|
---|
62 | }
|
---|
63 |
|
---|
64 | public AlgebraicInterval Div(AlgebraicInterval a) {
|
---|
65 | if (a.Contains(0.0)) {
|
---|
66 | if (a.low.Value.Value == 0 && a.high.Value.Value == 0) {
|
---|
67 | if (this.low.Value >= 0) {
|
---|
68 | // pos / 0
|
---|
69 | } else if (this.high.Value <= 0) {
|
---|
70 | // neg / 0
|
---|
71 | } else {
|
---|
72 | low = new MultivariateDual<AlgebraicDouble>(double.NegativeInfinity);
|
---|
73 | high = new MultivariateDual<AlgebraicDouble>(double.PositiveInfinity);
|
---|
74 | }
|
---|
75 | } else if (a.low.Value.Value >= 0) {
|
---|
76 | // a is positive
|
---|
77 | Mul(new AlgebraicInterval(a.Clone().high.Inv(), new MultivariateDual<AlgebraicDouble>(double.PositiveInfinity)));
|
---|
78 | } else if (a.high.Value <= 0) {
|
---|
79 | // a is negative
|
---|
80 | Mul(new AlgebraicInterval(new MultivariateDual<AlgebraicDouble>(double.NegativeInfinity), a.low.Clone().Inv()));
|
---|
81 | } else {
|
---|
82 | // a is interval over zero
|
---|
83 | low = new MultivariateDual<AlgebraicDouble>(double.NegativeInfinity);
|
---|
84 | high = new MultivariateDual<AlgebraicDouble>(double.PositiveInfinity);
|
---|
85 | }
|
---|
86 | } else {
|
---|
87 | Mul(new AlgebraicInterval(a.high.Clone().Inv(), a.low.Clone().Inv())); // inverting leads to inverse roles of high and low
|
---|
88 | }
|
---|
89 | return this;
|
---|
90 | }
|
---|
91 |
|
---|
92 | public AlgebraicInterval AssignExp(AlgebraicInterval a) {
|
---|
93 | low.AssignExp(a.low);
|
---|
94 | high.AssignExp(a.high);
|
---|
95 | return this;
|
---|
96 | }
|
---|
97 |
|
---|
98 | // tanh is a bijective function
|
---|
99 | public AlgebraicInterval AssignTanh(AlgebraicInterval a) {
|
---|
100 | low.AssignTanh(a.low);
|
---|
101 | high.AssignTanh(a.high);
|
---|
102 | return this;
|
---|
103 | }
|
---|
104 |
|
---|
105 | public AlgebraicInterval AssignIntPower(AlgebraicInterval a, int p) {
|
---|
106 | if (p < 0) { // x^-3 == 1/(x^3)
|
---|
107 | AssignIntPower(a, -p);
|
---|
108 | return AssignInv(this);
|
---|
109 | } else if (p == 0) {
|
---|
110 | if (a.Contains(0.0)) {
|
---|
111 | // => 0^0 = 0 ; might not be relevant
|
---|
112 | low = new MultivariateDual<AlgebraicDouble>(0.0);
|
---|
113 | high = new MultivariateDual<AlgebraicDouble>(1.0);
|
---|
114 | return this;
|
---|
115 | } else {
|
---|
116 | // => 1
|
---|
117 | low = new MultivariateDual<AlgebraicDouble>(1.0);
|
---|
118 | high = new MultivariateDual<AlgebraicDouble>(1.0);
|
---|
119 | return this;
|
---|
120 | }
|
---|
121 | } else if (p == 1) return this;
|
---|
122 | else if (p % 2 == 0) {
|
---|
123 | // p is even => interval must be positive
|
---|
124 | if (a.Contains(0.0)) {
|
---|
125 | low = new MultivariateDual<AlgebraicDouble>(0.0);
|
---|
126 | AssignMax(high, a.low.IntPower(p), a.high.IntPower(p));
|
---|
127 | } else {
|
---|
128 | AssignLowAndHigh(a.low.IntPower(p), a.high.IntPower(p));
|
---|
129 | }
|
---|
130 | } else {
|
---|
131 | // p is uneven
|
---|
132 | if (a.Contains(0.0)) {
|
---|
133 | low.AssignIntPower(a.low, p);
|
---|
134 | high.AssignIntPower(a.high, p);
|
---|
135 | } else {
|
---|
136 | var lowPower = a.low.IntPower(p);
|
---|
137 | var highPower = a.high.IntPower(p);
|
---|
138 | AssignMin(low, lowPower, highPower);
|
---|
139 | AssignMax(high, lowPower, highPower);
|
---|
140 | }
|
---|
141 | }
|
---|
142 | return this;
|
---|
143 | }
|
---|
144 |
|
---|
145 | public AlgebraicInterval AssignIntRoot(AlgebraicInterval a, int r) {
|
---|
146 | if (r == 0) { low = new MultivariateDual<AlgebraicDouble>(double.NaN); high = new MultivariateDual<AlgebraicDouble>(double.NaN); return this; }
|
---|
147 | if (r == 1) return this;
|
---|
148 | if (r < 0) {
|
---|
149 | // x^ (-1/2) = 1 / (x^(1/2))
|
---|
150 | AssignIntRoot(a, -r);
|
---|
151 | return AssignInv(this);
|
---|
152 | } else {
|
---|
153 | // root only defined for positive arguments for even roots
|
---|
154 | if (r % 2 == 0 && a.LowerBound.Value.Value < 0) {
|
---|
155 | low = new MultivariateDual<AlgebraicDouble>(double.NaN);
|
---|
156 | high = new MultivariateDual<AlgebraicDouble>(double.NaN);
|
---|
157 | return this;
|
---|
158 | } else {
|
---|
159 | low.AssignIntRoot(a.low, r);
|
---|
160 | high.AssignIntRoot(a.high, r);
|
---|
161 | return this;
|
---|
162 | }
|
---|
163 | }
|
---|
164 | }
|
---|
165 |
|
---|
166 | public AlgebraicInterval AssignInv(AlgebraicInterval a) {
|
---|
167 | low = new MultivariateDual<AlgebraicDouble>(1.0);
|
---|
168 | high = new MultivariateDual<AlgebraicDouble>(1.0);
|
---|
169 | return Div(a);
|
---|
170 | }
|
---|
171 |
|
---|
172 | public AlgebraicInterval AssignLog(AlgebraicInterval a) {
|
---|
173 | low.AssignLog(a.low);
|
---|
174 | high.AssignLog(a.high);
|
---|
175 | return this;
|
---|
176 | }
|
---|
177 |
|
---|
178 | public AlgebraicInterval AssignNeg(AlgebraicInterval a) {
|
---|
179 | low.AssignNeg(a.high);
|
---|
180 | high.AssignNeg(a.low);
|
---|
181 | return this;
|
---|
182 | }
|
---|
183 |
|
---|
184 | public AlgebraicInterval Scale(double s) {
|
---|
185 | low.Scale(s);
|
---|
186 | high.Scale(s);
|
---|
187 | if (s < 0) {
|
---|
188 | var t = low;
|
---|
189 | low = high;
|
---|
190 | high = t;
|
---|
191 | }
|
---|
192 | return this;
|
---|
193 | }
|
---|
194 |
|
---|
195 | public AlgebraicInterval AssignSin(AlgebraicInterval a) {
|
---|
196 | var lower = a.LowerBound.Value.Value;
|
---|
197 | var size = a.UpperBound.Value.Value - lower;
|
---|
198 | if (size < 0) throw new InvalidProgramException(); // ASSERT interval >= 0;
|
---|
199 |
|
---|
200 | if (size >= Math.PI * 2) {
|
---|
201 | low = new MultivariateDual<AlgebraicDouble>(-1.0); // zero gradient
|
---|
202 | high = new MultivariateDual<AlgebraicDouble>(1.0);
|
---|
203 | return this;
|
---|
204 | }
|
---|
205 |
|
---|
206 | // assume low and high are in the same quadrant
|
---|
207 | AssignLowAndHigh(a.LowerBound.Clone().Sin(), a.UpperBound.Clone().Sin()); // AssignLowAndHigh determines the lower and higher value
|
---|
208 |
|
---|
209 | // override min and max if necessary
|
---|
210 |
|
---|
211 | // shift interval 'a' into the range [-2pi .. 2pi] without changing the size of the interval to simplify the checks
|
---|
212 | lower = lower % (2 * Math.PI); // lower in [-2pi .. 2pi]
|
---|
213 |
|
---|
214 | // handle min = -1 and max = 1 cases explicitly
|
---|
215 | var pi_2 = Math.PI / 2.0;
|
---|
216 | var maxima = new double[] { -3 * pi_2, pi_2 };
|
---|
217 | var minima = new double[] { -pi_2, 3 * pi_2 };
|
---|
218 |
|
---|
219 | // override min and max if necessary
|
---|
220 | if (maxima.Any(m => lower < m && lower + size > m)) {
|
---|
221 | // max = 1
|
---|
222 | high = new MultivariateDual<AlgebraicDouble>(1.0); // zero gradient
|
---|
223 | }
|
---|
224 |
|
---|
225 | if (minima.Any(m => lower < m && lower + size > m)) {
|
---|
226 | // min = -1;
|
---|
227 | low = new MultivariateDual<AlgebraicDouble>(-1.0); // zero gradient
|
---|
228 | }
|
---|
229 | return this;
|
---|
230 | }
|
---|
231 |
|
---|
232 | public AlgebraicInterval Sub(AlgebraicInterval a) {
|
---|
233 | // [x1,x2] − [y1,y2] = [x1 − y2,x2 − y1]
|
---|
234 | low.Sub(a.high);
|
---|
235 | high.Sub(a.low);
|
---|
236 | return this;
|
---|
237 | }
|
---|
238 |
|
---|
239 | public AlgebraicInterval Clone() {
|
---|
240 | return new AlgebraicInterval(low, high);
|
---|
241 | }
|
---|
242 |
|
---|
243 | public bool Contains(double val) {
|
---|
244 | return LowerBound.Value.Value <= val && val <= UpperBound.Value.Value;
|
---|
245 | }
|
---|
246 |
|
---|
247 | public AlgebraicInterval AssignAbs(AlgebraicInterval a) {
|
---|
248 | if (a.Contains(0.0)) {
|
---|
249 | low.Assign(new MultivariateDual<AlgebraicDouble>(0.0)); // lost gradient for lower bound
|
---|
250 | AssignMax(high, a.low.Clone().Abs(), a.high.Clone().Abs());
|
---|
251 | } else {
|
---|
252 | var abslow = a.low.Clone().Abs();
|
---|
253 | var abshigh = a.high.Clone().Abs();
|
---|
254 | AssignLowAndHigh(abslow, abshigh);
|
---|
255 | }
|
---|
256 | return this;
|
---|
257 | }
|
---|
258 |
|
---|
259 | public AlgebraicInterval AssignSgn(AlgebraicInterval a) {
|
---|
260 | low.AssignSgn(a.low);
|
---|
261 | high.AssignSgn(a.high);
|
---|
262 | return this;
|
---|
263 | }
|
---|
264 |
|
---|
265 |
|
---|
266 | #region helper
|
---|
267 | private void AssignMin(MultivariateDual<AlgebraicDouble> dest, MultivariateDual<AlgebraicDouble> a, MultivariateDual<AlgebraicDouble> b) {
|
---|
268 | if (a.Value.CompareTo(b.Value) <= 0) {
|
---|
269 | dest.Assign(a);
|
---|
270 | } else {
|
---|
271 | dest.Assign(b);
|
---|
272 | }
|
---|
273 | }
|
---|
274 | private void AssignMax(MultivariateDual<AlgebraicDouble> dest, MultivariateDual<AlgebraicDouble> a, MultivariateDual<AlgebraicDouble> b) {
|
---|
275 | if (a.Value.CompareTo(b.Value) <= 0) {
|
---|
276 | dest.Assign(b);
|
---|
277 | } else {
|
---|
278 | dest.Assign(a);
|
---|
279 | }
|
---|
280 | }
|
---|
281 |
|
---|
282 | // determines the smaller and larger value and sets low and high accordingly
|
---|
283 | private void AssignLowAndHigh(MultivariateDual<AlgebraicDouble> a, MultivariateDual<AlgebraicDouble> b) {
|
---|
284 | // we must make sure that low and high are different objects when a == b
|
---|
285 | if (a.Value.CompareTo(b.Value) == 0) {
|
---|
286 | low.Assign(a);
|
---|
287 | high.Assign(b);
|
---|
288 | } else {
|
---|
289 | AssignMin(low, a, b);
|
---|
290 | AssignMax(high, a, b);
|
---|
291 | }
|
---|
292 | }
|
---|
293 | private void AssignLowAndHigh(params MultivariateDual<AlgebraicDouble>[] xs) {
|
---|
294 | if (xs.Length <= 2) throw new ArgumentException("need at least 3 arguments");
|
---|
295 | AssignLowAndHigh(xs[0], xs[1]);
|
---|
296 | for(int i=2;i<xs.Length;i++) {
|
---|
297 | // we must make sure that low and high are different objects when a == b
|
---|
298 | if (low.Value.CompareTo(xs[i].Value) > 0) low.Assign(xs[i]);
|
---|
299 | if (high.Value.CompareTo(xs[i].Value) < 0) high.Assign(xs[i]);
|
---|
300 | }
|
---|
301 | }
|
---|
302 | #endregion
|
---|
303 |
|
---|
304 | }
|
---|
305 | } |
---|