source: branches/2994-AutoDiffForIntervals/HeuristicLab.Problems.DataAnalysis.Symbolic/3.4/Interpreter/AlgebraicInterval.cs @ 17295

Last change on this file since 17295 was 17295, checked in by gkronber, 3 years ago

#2994: refactoring: moved types into separate files

File size: 10.3 KB
Line 
1using System;
2using System.Diagnostics;
3using System.Linq;
4
5namespace 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      low = Min(Min(v1, v2), Min(v3, v4));
50      high = Max(Max(v1, v2), Max(v3, v4));
51
52      return this;
53    }
54
55
56    private static MultivariateDual<AlgebraicDouble> Min(MultivariateDual<AlgebraicDouble> a, MultivariateDual<AlgebraicDouble> b) {
57      return a.Value < b.Value ? a : b;
58    }
59    private static MultivariateDual<AlgebraicDouble> Max(MultivariateDual<AlgebraicDouble> a, MultivariateDual<AlgebraicDouble> b) {
60      return a.Value > b.Value ? a : b;
61    }
62
63    public AlgebraicInterval Assign(AlgebraicInterval a) {
64      low = a.low;
65      high = a.high;
66      return this;
67    }
68
69    public AlgebraicInterval AssignCos(AlgebraicInterval a) {
70      return AssignSin(a.Clone().Add(new AlgebraicInterval(Math.PI / 2, Math.PI / 2)));
71    }
72
73    public AlgebraicInterval Div(AlgebraicInterval a) {
74      if (a.Contains(0.0)) {
75        if (a.low.Value.Value == 0 && a.high.Value.Value == 0) {
76          if (this.low.Value >= 0) {
77            // pos / 0
78          } else if (this.high.Value <= 0) {
79            // neg / 0
80          } else {
81            low = new MultivariateDual<AlgebraicDouble>(double.NegativeInfinity);
82            high = new MultivariateDual<AlgebraicDouble>(double.PositiveInfinity);
83          }
84        } else if (a.low.Value.Value >= 0) {
85          // a is positive
86          Mul(new AlgebraicInterval(a.Clone().high.Inv(), new MultivariateDual<AlgebraicDouble>(double.PositiveInfinity)));
87        } else if (a.high.Value <= 0) {
88          // a is negative
89          Mul(new AlgebraicInterval(new MultivariateDual<AlgebraicDouble>(double.NegativeInfinity), a.low.Clone().Inv()));
90        } else {
91          // a is interval over zero
92          low = new MultivariateDual<AlgebraicDouble>(double.NegativeInfinity);
93          high = new MultivariateDual<AlgebraicDouble>(double.PositiveInfinity);
94        }
95      } else {
96        Mul(new AlgebraicInterval(a.high.Clone().Inv(), a.low.Clone().Inv())); // inverting leads to inverse roles of high and low
97      }
98      return this;
99    }
100
101    public AlgebraicInterval AssignExp(AlgebraicInterval a) {
102      low.AssignExp(a.low);
103      high.AssignExp(a.high);
104      return this;
105    }
106
107    // tanh is a bijective function
108    public AlgebraicInterval AssignTanh(AlgebraicInterval a) {
109      low.AssignTanh(a.low);
110      high.AssignTanh(a.high);
111      return this;
112    }
113
114    public AlgebraicInterval AssignIntPower(AlgebraicInterval a, int p) {
115      if (p < 0) {  // x^-3 == 1/(x^3)
116        AssignIntPower(a, -p);
117        return AssignInv(this);
118      } else if (p == 0) {
119        if (a.Contains(0.0)) {
120          // => 0^0 = 0 ; might not be relevant
121          low = new MultivariateDual<AlgebraicDouble>(0.0);
122          high = new MultivariateDual<AlgebraicDouble>(1.0);
123          return this;
124        } else {
125          // => 1
126          low = new MultivariateDual<AlgebraicDouble>(1.0);
127          high = new MultivariateDual<AlgebraicDouble>(1.0);
128          return this;
129        }
130      } else if (p == 1) return this;
131      else if (p % 2 == 0) {
132        // p is even => interval must be positive
133        if (a.Contains(0.0)) {
134          low = new MultivariateDual<AlgebraicDouble>(0.0);
135          high = a.low.IntPower(p).AssignMax(a.high.IntPower(p));
136        } else {
137          var lowPower = a.low.IntPower(p);
138          var highPower = a.high.IntPower(p);
139          low = lowPower.AssignMin(highPower);
140          high = lowPower.AssignMax(highPower);
141        }
142      } else {
143        // p is uneven
144        if (a.Contains(0.0)) {
145          low.AssignIntPower(a.low, p);
146          high.AssignIntPower(a.high, p);
147        } else {
148          var lowPower = a.low.IntPower(p);
149          var highPower = a.high.IntPower(p);
150          low = lowPower.AssignMin(highPower);
151          high = lowPower.AssignMax(highPower);
152        }
153      }
154      return this;
155    }
156
157    public AlgebraicInterval AssignIntRoot(AlgebraicInterval a, int r) {
158      if (r == 0) { low = new MultivariateDual<AlgebraicDouble>(double.NaN); high = new MultivariateDual<AlgebraicDouble>(double.NaN); return this; }
159      if (r == 1) return this;
160      if (r < 0) {
161        // x^ (-1/2) = 1 / (x^(1/2))
162        AssignIntRoot(a, -r);
163        return AssignInv(this);
164      } else {
165        // root only defined for positive arguments for even roots
166        if (r % 2 == 0 && a.LowerBound.Value.Value < 0) {
167          low = new MultivariateDual<AlgebraicDouble>(double.NaN);
168          high = new MultivariateDual<AlgebraicDouble>(double.NaN);
169          return this;
170        } else {
171          low.AssignIntRoot(a.low, r);
172          high.AssignIntRoot(a.high, r);
173          return this;
174        }
175      }
176    }
177
178    public AlgebraicInterval AssignInv(AlgebraicInterval a) {
179      low = new MultivariateDual<AlgebraicDouble>(1.0);
180      high = new MultivariateDual<AlgebraicDouble>(1.0);
181      return Div(a);
182    }
183
184    public AlgebraicInterval AssignLog(AlgebraicInterval a) {
185      low.AssignLog(a.low);
186      high.AssignLog(a.high);
187      return this;
188    }
189
190    public AlgebraicInterval AssignNeg(AlgebraicInterval a) {
191      low.AssignNeg(a.high);
192      high.AssignNeg(a.low);
193      return this;
194    }
195
196    public AlgebraicInterval Scale(double s) {
197      low.Scale(s);
198      high.Scale(s);
199      if (s < 0) {
200        var t = low;
201        low = high;
202        high = t;
203      }
204      return this;
205    }
206
207    public AlgebraicInterval AssignSin(AlgebraicInterval a) {
208      var lower = a.LowerBound.Value.Value;
209      var size = a.UpperBound.Value.Value - lower;
210      if (size < 0) throw new InvalidProgramException(); // ASSERT interval >= 0;
211
212      if (size >= Math.PI * 2) {
213        low = new MultivariateDual<AlgebraicDouble>(-1.0); // zero gradient
214        high = new MultivariateDual<AlgebraicDouble>(1.0);
215        return this;
216      }
217
218      // assume low and high are in the same quadrant
219      low = Algebraic.Min(a.LowerBound.Clone().Sin(), a.UpperBound.Clone().Sin());
220      high = Algebraic.Max(a.LowerBound.Clone().Sin(), a.UpperBound.Clone().Sin());
221
222      // override min and max if necessary
223
224      // shift interval 'a' into the range [-2pi .. 2pi] without changing the size of the interval to simplify the checks
225      lower = lower % (2 * Math.PI); // lower in [-2pi .. 2pi]     
226
227      // handle min = -1 and max = 1 cases explicitly
228      var pi_2 = Math.PI / 2.0;
229      var maxima = new double[] { -3 * pi_2, pi_2 };
230      var minima = new double[] { -pi_2, 3 * pi_2 };
231
232      // override min and max if necessary
233      if (maxima.Any(m => lower < m && lower + size > m)) {
234        // max = 1
235        high = new MultivariateDual<AlgebraicDouble>(1.0); // zero gradient
236      }
237
238      if (minima.Any(m => lower < m && lower + size > m)) {
239        // min = -1;
240        low = new MultivariateDual<AlgebraicDouble>(-1.0); // zero gradient
241      }
242      return this;
243    }
244
245    public AlgebraicInterval Sub(AlgebraicInterval a) {
246      // [x1,x2] − [y1,y2] = [x1 − y2,x2 − y1]
247      low.Sub(a.high);
248      high.Sub(a.low);
249      return this;
250    }
251
252    public AlgebraicInterval Clone() {
253      return new AlgebraicInterval(low, high);
254    }
255
256    public bool Contains(double val) {
257      return LowerBound.Value.Value <= val && val <= UpperBound.Value.Value;
258    }
259
260    public AlgebraicInterval AssignAbs(AlgebraicInterval a) {
261      if (a.Contains(0.0)) {
262        var abslow = a.low.Clone().Abs();
263        var abshigh = a.high.Clone().Abs();
264        a.high.Assign(Algebraic.Max(abslow, abshigh));
265        a.low.Assign(new MultivariateDual<AlgebraicDouble>(0.0)); // lost gradient for lower bound
266      } else {
267        var abslow = a.low.Clone().Abs();
268        var abshigh = a.high.Clone().Abs();
269        a.low.Assign(Algebraic.Min(abslow, abshigh));
270        a.high.Assign(Algebraic.Max(abslow, abshigh));
271      }
272      return this;
273    }
274
275    public AlgebraicInterval AssignSgn(AlgebraicInterval a) {
276      low.AssignSgn(a.low);
277      high.AssignSgn(a.high);
278      return this;
279    }
280
281    public AlgebraicInterval AssignMin(AlgebraicInterval other) {
282      low.AssignMin(other.low);
283      high.AssignMin(other.high);
284      return this;
285    }
286
287    public AlgebraicInterval AssignMax(AlgebraicInterval other) {
288      low.AssignMax(other.low);
289      high.AssignMax(other.high);
290      return this;
291    }
292  }
293}
Note: See TracBrowser for help on using the repository browser.