[17295] | 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 | 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 | } |
---|