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

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

#2994 continued refactoring and extended unit tests. Interval calculation still fails for some edge cases (mainly for undefined behaviour). VectorEvaluator and VectorAutoDiffEvaluator produce the same results as the LinearInterpreter. TODO: check gradient calculation

File size: 10.8 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      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}
Note: See TracBrowser for help on using the repository browser.