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

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

#2994: min() and max() are problematic in AutoDiff. We have used max(a,b) = (a+b + abs(b-a)) / 2 but this is problematic because of catastrophic cancellation. Min() and Max() are necessary for interval calculation. However, there we only need to determine min and max for double values (or duals). We do not need min() and max() in general.

File size: 56.3 KB
Line 
1using System;
2using System.Collections.Generic;
3using System.Diagnostics;
4using System.Linq;
5using HeuristicLab.Common;
6using HeuristicLab.Encodings.SymbolicExpressionTreeEncoding;
7
8namespace HeuristicLab.Problems.DataAnalysis.Symbolic {
9  public abstract class Interpreter<T> where T : IAlgebraicType<T> {
10    public struct Instruction {
11      public byte opcode;
12      public ushort narg;
13      public int childIndex;
14      public double dblVal;
15      public object data; // any kind of data you want to store in instructions
16      public T value;
17    }
18
19    public T Evaluate(Instruction[] code) {
20      for (int i = code.Length - 1; i >= 0; --i) {
21        var instr = code[i];
22        var c = instr.childIndex;
23        var n = instr.narg;
24
25        switch (instr.opcode) {
26          case OpCodes.Variable: {
27              LoadVariable(instr);
28              break;
29            }
30          case OpCodes.Constant: { break; }  // we initialize constants in Compile. The value never changes afterwards
31          case OpCodes.Add: {
32              instr.value.Assign(code[c].value);
33              for (int j = 1; j < n; ++j) {
34                instr.value.Add(code[c + j].value);
35              }
36              break;
37            }
38
39          case OpCodes.Sub: {
40              if (n == 1) {
41                instr.value.AssignNeg(code[c].value);
42              } else {
43                instr.value.Assign(code[c].value);
44                for (int j = 1; j < n; ++j) {
45                  instr.value.Sub(code[c + j].value);
46                }
47              }
48              break;
49            }
50
51          case OpCodes.Mul: {
52              instr.value.Assign(code[c].value);
53              for (int j = 1; j < n; ++j) {
54                instr.value.Mul(code[c + j].value);
55              }
56              break;
57            }
58
59          case OpCodes.Div: {
60              if (n == 1) {
61                instr.value.AssignInv(code[c].value);
62              } else {
63                instr.value.Assign(code[c].value);
64                for (int j = 1; j < n; ++j) {
65                  instr.value.Div(code[c + j].value);
66                }
67              }
68              break;
69            }
70          case OpCodes.Square: {
71              instr.value.AssignIntPower(code[c].value, 2);
72              break;
73            }
74          case OpCodes.SquareRoot: {
75              instr.value.AssignIntRoot(code[c].value, 2);
76              break;
77            }
78          case OpCodes.Cube: {
79              instr.value.AssignIntPower(code[c].value, 3);
80              break;
81            }
82          case OpCodes.CubeRoot: {
83              instr.value.AssignIntRoot(code[c].value, 3);
84              break;
85            }
86          case OpCodes.Exp: {
87              instr.value.AssignExp(code[c].value);
88              break;
89            }
90          case OpCodes.Log: {
91              instr.value.AssignLog(code[c].value);
92              break;
93            }
94          case OpCodes.Sin: {
95              instr.value.AssignSin(code[c].value);
96              break;
97            }
98          case OpCodes.Cos: {
99              instr.value.AssignCos(code[c].value);
100              break;
101            }
102          case OpCodes.Absolute: {
103              instr.value.AssignAbs(code[c].value);
104              break;
105            }
106          case OpCodes.AnalyticQuotient: {
107              instr.value.Assign(code[c].value);
108              for (int j = 1; j < n; ++j) {
109                var t = instr.value.One;
110                t.Add(code[c + j].value.Clone().IntPower(2));
111                instr.value.Div(t.IntRoot(2));
112              }
113              break;
114            }
115          case OpCodes.Tanh: {
116              instr.value.AssignTanh(code[c].value);
117              break;
118            }
119
120          default: throw new ArgumentException($"Unknown opcode {instr.opcode}");
121        }
122      }
123      return code[0].value;
124    }
125
126    protected Instruction[] Compile(ISymbolicExpressionTree tree) {
127      var root = tree.Root.GetSubtree(0).GetSubtree(0);
128      var code = new Instruction[root.GetLength()];
129      if (root.SubtreeCount > ushort.MaxValue) throw new ArgumentException("Number of subtrees is too big (>65.535)");
130      int c = 1, i = 0;
131      foreach (var node in root.IterateNodesBreadth()) {
132        if (node.SubtreeCount > ushort.MaxValue) throw new ArgumentException("Number of subtrees is too big (>65.535)");
133        code[i] = new Instruction {
134          opcode = OpCodes.MapSymbolToOpCode(node),
135          narg = (ushort)node.SubtreeCount,
136          childIndex = c
137        };
138        if (node is VariableTreeNode variable) {
139          InitializeTerminalInstruction(ref code[i], variable);
140        } else if (node is ConstantTreeNode constant) {
141          InitializeTerminalInstruction(ref code[i], constant);
142        } else {
143          InitializeInternalInstruction(ref code[i], node);
144        }
145        c += node.SubtreeCount;
146        ++i;
147      }
148      return code;
149    }
150
151    protected abstract void InitializeTerminalInstruction(ref Instruction instruction, ConstantTreeNode constant);
152    protected abstract void InitializeTerminalInstruction(ref Instruction instruction, VariableTreeNode variable);
153    protected abstract void InitializeInternalInstruction(ref Instruction instruction, ISymbolicExpressionTreeNode node);
154
155    protected abstract void LoadVariable(Instruction a);
156
157  }
158
159
160  public sealed class VectorEvaluator : Interpreter<AlgebraicDoubleVector> {
161    private const int BATCHSIZE = 128;
162    [ThreadStatic]
163    private Dictionary<string, double[]> cachedData;
164
165    [ThreadStatic]
166    private IDataset dataset;
167
168    [ThreadStatic]
169    private int rowIndex;
170
171    [ThreadStatic]
172    private int[] rows;
173
174    private void InitCache(IDataset dataset) {
175      this.dataset = dataset;
176      cachedData = new Dictionary<string, double[]>();
177      foreach (var v in dataset.DoubleVariables) {
178        cachedData[v] = dataset.GetReadOnlyDoubleValues(v).ToArray();
179      }
180    }
181
182    public double[] Evaluate(ISymbolicExpressionTree tree, IDataset dataset, int[] rows) {
183      if (cachedData == null || this.dataset != dataset) {
184        InitCache(dataset);
185      }
186
187      this.rows = rows;
188      var code = Compile(tree);
189      var remainingRows = rows.Length % BATCHSIZE;
190      var roundedTotal = rows.Length - remainingRows;
191
192      var result = new double[rows.Length];
193
194      for (rowIndex = 0; rowIndex < roundedTotal; rowIndex += BATCHSIZE) {
195        Evaluate(code);
196        code[0].value.CopyTo(result, rowIndex, BATCHSIZE);
197      }
198
199      if (remainingRows > 0) {
200        Evaluate(code);
201        code[0].value.CopyTo(result, roundedTotal, remainingRows);
202      }
203
204      return result;
205    }
206
207    protected override void InitializeTerminalInstruction(ref Instruction instruction, ConstantTreeNode constant) {
208      instruction.dblVal = constant.Value;
209      instruction.value = new AlgebraicDoubleVector(BATCHSIZE);
210      instruction.value.AssignConstant(instruction.dblVal);
211    }
212
213    protected override void InitializeTerminalInstruction(ref Instruction instruction, VariableTreeNode variable) {
214      instruction.dblVal = variable.Weight;
215      instruction.value = new AlgebraicDoubleVector(BATCHSIZE);
216      if (cachedData.ContainsKey(variable.VariableName)) {
217        instruction.data = cachedData[variable.VariableName];
218      } else {
219        instruction.data = dataset.GetDoubleValues(variable.VariableName).ToArray();
220        cachedData[variable.VariableName] = (double[])instruction.data;
221      }
222    }
223
224    protected override void InitializeInternalInstruction(ref Instruction instruction, ISymbolicExpressionTreeNode node) {
225      instruction.value = new AlgebraicDoubleVector(BATCHSIZE);
226    }
227
228    protected override void LoadVariable(Instruction a) {
229      var data = (double[])a.data;
230      for (int i = rowIndex; i < rows.Length && i - rowIndex < BATCHSIZE; i++) a.value[i - rowIndex] = data[rows[i]];
231      a.value.Scale(a.dblVal);
232    }
233  }
234
235  public sealed class VectorAutoDiffEvaluator : Interpreter<MultivariateDual<AlgebraicDoubleVector>> {
236    private const int BATCHSIZE = 128;
237    [ThreadStatic]
238    private Dictionary<string, double[]> cachedData;
239
240    [ThreadStatic]
241    private IDataset dataset;
242
243    [ThreadStatic]
244    private int rowIndex;
245
246    [ThreadStatic]
247    private int[] rows;
248
249    [ThreadStatic]
250    private Dictionary<ISymbolicExpressionTreeNode, int> node2paramIdx;
251
252    private void InitCache(IDataset dataset) {
253      this.dataset = dataset;
254      cachedData = new Dictionary<string, double[]>();
255      foreach (var v in dataset.DoubleVariables) {
256        cachedData[v] = dataset.GetDoubleValues(v).ToArray();
257      }
258    }
259
260    /// <summary>
261    ///
262    /// </summary>
263    /// <param name="tree"></param>
264    /// <param name="dataset"></param>
265    /// <param name="rows"></param>
266    /// <param name="parameterNodes"></param>
267    /// <param name="fi">Function output. Must be allocated by the caller.</param>
268    /// <param name="jac">Jacobian matrix. Must be allocated by the caller.</param>
269    public void Evaluate(ISymbolicExpressionTree tree, IDataset dataset, int[] rows, ISymbolicExpressionTreeNode[] parameterNodes, double[] fi, double[,] jac) {
270      if (cachedData == null || this.dataset != dataset) {
271        InitCache(dataset);
272      }
273
274      int nParams = parameterNodes.Length;
275      node2paramIdx = new Dictionary<ISymbolicExpressionTreeNode, int>();
276      for (int i = 0; i < parameterNodes.Length; i++) node2paramIdx.Add(parameterNodes[i], i);
277
278      var code = Compile(tree);
279
280      var remainingRows = rows.Length % BATCHSIZE;
281      var roundedTotal = rows.Length - remainingRows;
282
283      this.rows = rows;
284
285      for (rowIndex = 0; rowIndex < roundedTotal; rowIndex += BATCHSIZE) {
286        Evaluate(code);
287        code[0].value.Value.CopyTo(fi, rowIndex, BATCHSIZE);
288
289        // TRANSPOSE into JAC
290        var g = code[0].value.Gradient;
291        for (int j = 0; j < nParams; ++j) {
292          if (g.Elements.TryGetValue(j, out AlgebraicDoubleVector v)) {
293            v.CopyColumnTo(jac, j, rowIndex, BATCHSIZE);
294          } else {
295            for (int r = 0; r < BATCHSIZE; r++) jac[rowIndex + r, j] = 0.0;
296          }
297        }
298      }
299
300      if (remainingRows > 0) {
301        Evaluate(code);
302        code[0].value.Value.CopyTo(fi, roundedTotal, remainingRows);
303
304        var g = code[0].value.Gradient;
305        for (int j = 0; j < nParams; ++j)
306          if (g.Elements.TryGetValue(j, out AlgebraicDoubleVector v)) {
307            v.CopyColumnTo(jac, j, roundedTotal, remainingRows);
308          } else {
309            for (int r = 0; r < remainingRows; r++) jac[roundedTotal + r, j] = 0.0;
310          }
311      }
312    }
313
314    protected override void InitializeInternalInstruction(ref Instruction instruction, ISymbolicExpressionTreeNode node) {
315      var zero = new AlgebraicDoubleVector(BATCHSIZE);
316      instruction.value = new MultivariateDual<AlgebraicDoubleVector>(zero);
317    }
318
319    protected override void InitializeTerminalInstruction(ref Instruction instruction, ConstantTreeNode constant) {
320      var g_arr = new double[BATCHSIZE];
321      if (node2paramIdx.TryGetValue(constant, out var paramIdx)) {
322        for (int i = 0; i < BATCHSIZE; i++) g_arr[i] = 1.0;
323        var g = new AlgebraicDoubleVector(g_arr);
324        instruction.value = new MultivariateDual<AlgebraicDoubleVector>(new AlgebraicDoubleVector(BATCHSIZE), paramIdx, g); // only a single column for the gradient
325      } else {
326        instruction.value = new MultivariateDual<AlgebraicDoubleVector>(new AlgebraicDoubleVector(BATCHSIZE));
327      }
328
329      instruction.dblVal = constant.Value;
330      instruction.value.Value.AssignConstant(instruction.dblVal);
331    }
332
333    protected override void InitializeTerminalInstruction(ref Instruction instruction, VariableTreeNode variable) {
334      double[] data;
335      if (cachedData.ContainsKey(variable.VariableName)) {
336        data = cachedData[variable.VariableName];
337      } else {
338        data = dataset.GetReadOnlyDoubleValues(variable.VariableName).ToArray();
339        cachedData[variable.VariableName] = (double[])instruction.data;
340      }
341
342      var paramIdx = -1;
343      if (node2paramIdx.ContainsKey(variable)) {
344        paramIdx = node2paramIdx[variable];
345        var f = new AlgebraicDoubleVector(BATCHSIZE);
346        var g = new AlgebraicDoubleVector(BATCHSIZE);
347        instruction.value = new MultivariateDual<AlgebraicDoubleVector>(f, paramIdx, g);
348      } else {
349        var f = new AlgebraicDoubleVector(BATCHSIZE);
350        instruction.value = new MultivariateDual<AlgebraicDoubleVector>(f);
351      }
352
353      instruction.dblVal = variable.Weight;
354      instruction.data = new object[] { data, paramIdx };
355    }
356
357    protected override void LoadVariable(Instruction a) {
358      var paramIdx = (int)((object[])a.data)[1];
359      var data = (double[])((object[])a.data)[0];
360
361      for (int i = rowIndex; i < rows.Length && i - rowIndex < BATCHSIZE; i++) a.value.Value[i - rowIndex] = data[rows[i]];
362      a.value.Scale(a.dblVal);
363
364      if (paramIdx >= 0) {
365        // update gradient with variable values
366        var g = a.value.Gradient.Elements[paramIdx];
367        for (int i = rowIndex; i < rows.Length && i - rowIndex < BATCHSIZE; i++) {
368          g[i - rowIndex] = data[rows[i]];
369        }
370      }
371    }
372  }
373
374
375  public sealed class IntervalEvaluator : Interpreter<AlgebraicInterval> {
376    [ThreadStatic]
377    private IDictionary<string, Interval> intervals;
378
379    public Interval Evaluate(ISymbolicExpressionTree tree, IDictionary<string, Interval> intervals) {
380      this.intervals = intervals;
381      var code = Compile(tree);
382      Evaluate(code);
383      if (code[0].value.LowerBound.Value.Value > code[0].value.UpperBound.Value.Value) throw new InvalidProgramException($"lower: {code[0].value.LowerBound.Value.Value} > upper: {code[0].value.UpperBound.Value.Value}");
384      return new Interval(code[0].value.LowerBound.Value.Value, code[0].value.UpperBound.Value.Value);
385    }
386
387    public Interval Evaluate(ISymbolicExpressionTree tree, IDictionary<string, Interval> intervals, ISymbolicExpressionTreeNode[] paramNodes, out double[] lowerGradient, out double[] upperGradient) {
388      this.intervals = intervals;
389      var code = Compile(tree);
390      Evaluate(code);
391      lowerGradient = new double[paramNodes.Length];
392      upperGradient = new double[paramNodes.Length];
393      var l = code[0].value.LowerBound;
394      var u = code[0].value.UpperBound;
395      for (int i = 0; i < paramNodes.Length; ++i) {
396        if (paramNodes[i] == null) continue;
397        if (l.Gradient.Elements.TryGetValue(paramNodes[i], out AlgebraicDouble value)) lowerGradient[i] = value;
398        if (u.Gradient.Elements.TryGetValue(paramNodes[i], out value)) upperGradient[i] = value;
399      }
400      return new Interval(code[0].value.LowerBound.Value.Value, code[0].value.UpperBound.Value.Value);
401    }
402
403    protected override void InitializeInternalInstruction(ref Instruction instruction, ISymbolicExpressionTreeNode node) {
404      instruction.value = new AlgebraicInterval(0, 0);
405    }
406
407
408    protected override void InitializeTerminalInstruction(ref Instruction instruction, ConstantTreeNode constant) {
409      instruction.dblVal = constant.Value;
410      instruction.value = new AlgebraicInterval(
411        new MultivariateDual<AlgebraicDouble>(instruction.dblVal, constant, 1.0),
412        new MultivariateDual<AlgebraicDouble>(instruction.dblVal, constant, 1.0) // use node as key
413        );
414    }
415
416    protected override void InitializeTerminalInstruction(ref Instruction instruction, VariableTreeNode variable) {
417      instruction.dblVal = variable.Weight;
418      var v1 = instruction.dblVal * intervals[variable.VariableName].LowerBound;
419      var v2 = instruction.dblVal * intervals[variable.VariableName].UpperBound;
420      var lower = Math.Min(v1, v2);
421      var upper = Math.Max(v1, v2);
422      // we assume that the for variable nodes ( v(x,w) = w * x ) the gradient is returned for parameter w
423      instruction.value = new AlgebraicInterval(
424        low: new MultivariateDual<AlgebraicDouble>(v: lower, key: variable, dv: intervals[variable.VariableName].LowerBound),
425        high: new MultivariateDual<AlgebraicDouble>(v: upper, key: variable, dv: intervals[variable.VariableName].UpperBound)
426        );
427    }
428
429    protected override void LoadVariable(Instruction a) {
430      // nothing to do
431    }
432  }
433
434  public interface IAlgebraicType<T> {
435    T Zero { get; } // Zero and One must create new objects
436    T One { get; }
437
438    T AssignAbs(T a); // set this to assign abs(a)
439    T Assign(T a); // assign this to same value as a (copy!)
440    T AssignNeg(T a); // set this to negative(a)
441    T AssignInv(T a); // set this to inv(a);
442    T Scale(double s); // scale this with s
443    T Add(T a); // add a to this
444    T Sub(T a); // subtract a from this
445    T Mul(T a); // multiply this with a
446    T Div(T a); // divide this by a
447    T AssignLog(T a); // set this to log a
448    T AssignExp(T a); // set this to exp(a)
449    T AssignSin(T a); // set this to sin(a)
450    T AssignCos(T a); // set this to cos(a)
451    T AssignTanh(T a); // set this to tanh(a)
452    T AssignIntPower(T a, int p);
453    T AssignIntRoot(T a, int r);
454    T AssignSgn(T a); // set this to sign(a)
455    T AssignMin(T other); // set this min(this, other)
456    T AssignMax(T other); // set this max(this, other)
457    T Clone();
458  }
459
460  public static class Algebraic {
461    public static T Abs<T>(this T a) where T : IAlgebraicType<T> { a.AssignAbs(a.Clone()); return a; }
462    public static T Neg<T>(this T a) where T : IAlgebraicType<T> { a.AssignNeg(a.Clone()); return a; }
463    public static T Inv<T>(this T a) where T : IAlgebraicType<T> { a.AssignInv(a.Clone()); return a; }
464    public static T Log<T>(this T a) where T : IAlgebraicType<T> { a.AssignLog(a.Clone()); return a; }
465    public static T Exp<T>(this T a) where T : IAlgebraicType<T> { a.AssignExp(a.Clone()); return a; }
466    public static T Sin<T>(this T a) where T : IAlgebraicType<T> { a.AssignSin(a.Clone()); return a; }
467    public static T Cos<T>(this T a) where T : IAlgebraicType<T> { a.AssignCos(a.Clone()); return a; }
468    public static T Sgn<T>(this T a) where T : IAlgebraicType<T> { a.AssignSgn(a.Clone()); return a; }
469    public static T IntPower<T>(this T a, int p) where T : IAlgebraicType<T> { a.AssignIntPower(a.Clone(), p); return a; }
470    public static T IntRoot<T>(this T a, int r) where T : IAlgebraicType<T> { a.AssignIntRoot(a.Clone(), r); return a; }
471
472    internal static T Min<T>(T a, T b) where T : IAlgebraicType<T> { return a.Clone().AssignMin(b); }
473    internal static T Max<T>(T a, T b) where T : IAlgebraicType<T> { return a.Clone().AssignMax(b); }
474
475    // public static T Max<T>(T a, T b) where T : IAlgebraicType<T> {
476    //   // ((a + b) + abs(b - a)) / 2
477    //   return a.Clone().Add(b).Add(b.Clone().Sub(a).Abs()).Scale(1.0 / 2.0);
478    // }
479    // public static T Min<T>(T a, T b) where T : IAlgebraicType<T> {
480    //   // ((a + b) - abs(a - b)) / 2
481    //   return a.Clone().Add(b).Sub(a.Clone().Sub(b).Abs()).Scale(1.0 / 2.0);
482    // }
483  }
484
485
486  // algebraic type wrapper for a double value
487  [DebuggerDisplay("{Value}")]
488  public sealed class AlgebraicDouble : IAlgebraicType<AlgebraicDouble> {
489    public static implicit operator AlgebraicDouble(double value) { return new AlgebraicDouble(value); }
490    public static implicit operator double(AlgebraicDouble value) { return value.Value; }
491    public double Value;
492
493    [DebuggerBrowsable(DebuggerBrowsableState.Never)]
494    public AlgebraicDouble Zero => new AlgebraicDouble(0.0);
495    [DebuggerBrowsable(DebuggerBrowsableState.Never)]
496    public AlgebraicDouble One => new AlgebraicDouble(1.0);
497
498    public bool IsInfinity => IsNegativeInfinity || IsPositiveInfinity;
499    public bool IsNegativeInfinity => double.IsNegativeInfinity(Value);
500    public bool IsPositiveInfinity => double.IsPositiveInfinity(Value);
501    public AlgebraicDouble() { }
502    public AlgebraicDouble(double value) { this.Value = value; }
503    public AlgebraicDouble Assign(AlgebraicDouble a) { Value = a.Value; return this; }
504    public AlgebraicDouble Add(AlgebraicDouble a) { Value += a.Value; return this; }
505    public AlgebraicDouble Sub(AlgebraicDouble a) { Value -= a.Value; return this; }
506    public AlgebraicDouble Mul(AlgebraicDouble a) { Value *= a.Value; return this; }
507    public AlgebraicDouble Div(AlgebraicDouble a) { Value /= a.Value; return this; }
508    public AlgebraicDouble Scale(double s) { Value *= s; return this; }
509    public AlgebraicDouble AssignInv(AlgebraicDouble a) { Value = 1.0 / a.Value; return this; }
510    public AlgebraicDouble AssignNeg(AlgebraicDouble a) { Value = -a.Value; return this; }
511    public AlgebraicDouble AssignSin(AlgebraicDouble a) { Value = Math.Sin(a.Value); return this; }
512    public AlgebraicDouble AssignCos(AlgebraicDouble a) { Value = Math.Cos(a.Value); return this; }
513    public AlgebraicDouble AssignTanh(AlgebraicDouble a) { Value = Math.Tanh(a.Value); return this; }
514    public AlgebraicDouble AssignLog(AlgebraicDouble a) { Value = Math.Log(a.Value); return this; }
515    public AlgebraicDouble AssignExp(AlgebraicDouble a) { Value = Math.Exp(a.Value); return this; }
516    public AlgebraicDouble AssignIntPower(AlgebraicDouble a, int p) { Value = Math.Pow(a.Value, p); return this; }
517    public AlgebraicDouble AssignIntRoot(AlgebraicDouble a, int r) { Value = IntRoot(a.Value, r); return this; }
518    public AlgebraicDouble AssignMin(AlgebraicDouble other) { Value = Math.Min(Value, other.Value); return this; }
519    public AlgebraicDouble AssignMax(AlgebraicDouble other) { Value = Math.Max(Value, other.Value); return this; }
520
521    // helper
522    private static double IntRoot(double value, int r) {
523      if (r % 2 == 0) return Math.Pow(value, 1.0 / r);
524      else return value < 0 ? -Math.Pow(-value, 1.0 / r) : Math.Pow(value, 1.0 / r);
525    }
526
527    public AlgebraicDouble AssignAbs(AlgebraicDouble a) { Value = Math.Abs(a.Value); return this; }
528    public AlgebraicDouble AssignSgn(AlgebraicDouble a) { Value = double.IsNaN(a.Value) ? double.NaN : Math.Sign(a.Value); return this; }
529    public AlgebraicDouble Clone() { return new AlgebraicDouble(Value); }
530
531    public override string ToString() {
532      return Value.ToString();
533    }
534  }
535
536  // a simple vector as an algebraic type
537  [DebuggerDisplay("DoubleVector(len={Length}): {string.}")]
538  public class AlgebraicDoubleVector : IAlgebraicType<AlgebraicDoubleVector> {
539    private double[] arr;
540    public double this[int idx] { get { return arr[idx]; } set { arr[idx] = value; } }
541    public int Length => arr.Length;
542
543    public AlgebraicDoubleVector(int length) { arr = new double[length]; }
544
545    public AlgebraicDoubleVector() { }
546
547    /// <summary>
548    ///
549    /// </summary>
550    /// <param name="arr">array is not copied</param>
551    public AlgebraicDoubleVector(double[] arr) { this.arr = arr; }
552
553    [DebuggerBrowsable(DebuggerBrowsableState.Never)]
554    public AlgebraicDoubleVector Zero => new AlgebraicDoubleVector(new double[this.Length]); // must return vector of same length as this (therefore Zero is not static)
555    [DebuggerBrowsable(DebuggerBrowsableState.Never)]
556    public AlgebraicDoubleVector One => new AlgebraicDoubleVector(new double[this.Length]).AssignConstant(1.0); // must return vector of same length as this (therefore Zero is not static)
557    public AlgebraicDoubleVector Assign(AlgebraicDoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = a.arr[i]; } return this; }
558    public AlgebraicDoubleVector Add(AlgebraicDoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] += a.arr[i]; } return this; }
559    public AlgebraicDoubleVector Sub(AlgebraicDoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] -= a.arr[i]; } return this; }
560    public AlgebraicDoubleVector Mul(AlgebraicDoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] *= a.arr[i]; } return this; }
561    public AlgebraicDoubleVector Div(AlgebraicDoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] /= a.arr[i]; } return this; }
562    public AlgebraicDoubleVector AssignNeg(AlgebraicDoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = -a.arr[i]; } return this; }
563    public AlgebraicDoubleVector AssignInv(AlgebraicDoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = 1.0 / a.arr[i]; } return this; }
564    public AlgebraicDoubleVector Scale(double s) { for (int i = 0; i < arr.Length; ++i) { arr[i] *= s; } return this; }
565    public AlgebraicDoubleVector AssignLog(AlgebraicDoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = Math.Log(a.arr[i]); } return this; }
566    public AlgebraicDoubleVector AssignSin(AlgebraicDoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = Math.Sin(a.arr[i]); } return this; }
567    public AlgebraicDoubleVector AssignExp(AlgebraicDoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = Math.Exp(a.arr[i]); } return this; }
568    public AlgebraicDoubleVector AssignCos(AlgebraicDoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = Math.Cos(a.arr[i]); } return this; }
569    public AlgebraicDoubleVector AssignTanh(AlgebraicDoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = Math.Tanh(a.arr[i]); } return this; }
570    public AlgebraicDoubleVector AssignIntPower(AlgebraicDoubleVector a, int p) { for (int i = 0; i < arr.Length; ++i) { arr[i] = Math.Pow(a.arr[i], p); } return this; }
571    public AlgebraicDoubleVector AssignIntRoot(AlgebraicDoubleVector a, int r) { for (int i = 0; i < arr.Length; ++i) { arr[i] = IntRoot(a.arr[i], r); } return this; }
572    public AlgebraicDoubleVector AssignMin(AlgebraicDoubleVector other) { for (int i = 0; i < arr.Length; ++i) { arr[i] = Math.Min(arr[i], other.arr[i]); } return this; }
573    public AlgebraicDoubleVector AssignMax(AlgebraicDoubleVector other) { for (int i = 0; i < arr.Length; ++i) { arr[i] = Math.Max(arr[i], other.arr[i]); } return this; }
574
575    // helper
576    private double IntRoot(double v, int r) {
577      if (r % 2 == 0) return Math.Pow(v, 1.0 / r);
578      else return v < 0 ? -Math.Pow(-v, 1.0 / r) : Math.Pow(v, 1.0 / r);
579    }
580
581    public AlgebraicDoubleVector AssignAbs(AlgebraicDoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = Math.Abs(a.arr[i]); } return this; }
582    public AlgebraicDoubleVector AssignSgn(AlgebraicDoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = Math.Sign(a.arr[i]); } return this; }
583
584    public AlgebraicDoubleVector Clone() {
585      var v = new AlgebraicDoubleVector(this.arr.Length);
586      Array.Copy(arr, v.arr, v.arr.Length);
587      return v;
588    }
589
590    public AlgebraicDoubleVector AssignConstant(double constantValue) {
591      for (int i = 0; i < arr.Length; ++i) {
592        arr[i] = constantValue;
593      }
594      return this;
595    }
596
597    public void CopyTo(double[] dest, int idx, int length) {
598      Array.Copy(arr, 0, dest, idx, length);
599    }
600
601    public void CopyFrom(double[] data, int rowIndex) {
602      Array.Copy(data, rowIndex, arr, 0, Math.Min(arr.Length, data.Length - rowIndex));
603    }
604    public void CopyRowTo(double[,] dest, int row) {
605      for (int j = 0; j < arr.Length; ++j) dest[row, j] = arr[j];
606    }
607
608    internal void CopyColumnTo(double[,] dest, int column, int row, int len) {
609      for (int j = 0; j < len; ++j) dest[row + j, column] = arr[j];
610    }
611
612    public override string ToString() {
613      return "{" + string.Join(", ", arr.Take(Math.Max(5, arr.Length))) + (arr.Length > 5 ? "..." : string.Empty) + "}";
614    }
615  }
616
617
618  /*
619  // vectors of algebraic types
620  public sealed class AlgebraicVector<T> : IAlgebraicType<AlgebraicVector<T>> where T : IAlgebraicType<T>, new() {
621    private T[] elems;
622
623    public T this[int idx] { get { return elems[idx]; } set { elems[idx] = value; } }
624
625    public int Length => elems.Length;
626
627    private AlgebraicVector() { }
628
629    public AlgebraicVector(int len) { elems = new T[len]; }
630
631    /// <summary>
632    ///
633    /// </summary>
634    /// <param name="elems">The array is copied (element-wise clone)</param>
635    public AlgebraicVector(T[] elems) {
636      this.elems = new T[elems.Length];
637      for (int i = 0; i < elems.Length; ++i) { this.elems[i] = elems[i].Clone(); }
638    }
639
640    /// <summary>
641    ///
642    /// </summary>
643    /// <param name="elems">Array is not copied!</param>
644    /// <returns></returns>
645    public AlgebraicVector<T> FromArray(T[] elems) {
646      var v = new AlgebraicVector<T>();
647      v.elems = elems;
648      return v;
649    }
650
651    public void CopyTo(T[] dest) {
652      if (dest.Length != elems.Length) throw new InvalidOperationException("arr lengths do not match in Vector<T>.Copy");
653      Array.Copy(elems, dest, dest.Length);
654    }
655
656    public AlgebraicVector<T> Clone() { return new AlgebraicVector<T>(elems); }
657
658
659    [DebuggerBrowsable(DebuggerBrowsableState.Never)]
660    public AlgebraicVector<T> Zero => new AlgebraicVector<T>(Length);
661    [DebuggerBrowsable(DebuggerBrowsableState.Never)]
662    public AlgebraicVector<T> One { get { var v = new AlgebraicVector<T>(Length); for (int i = 0; i < elems.Length; ++i) elems[i] = new T().One; return v; } }
663    public AlgebraicVector<T> Assign(AlgebraicVector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].Assign(a.elems[i]); } return this; }
664    public AlgebraicVector<T> Add(AlgebraicVector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].Add(a.elems[i]); } return this; }
665    public AlgebraicVector<T> Sub(AlgebraicVector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].Sub(a.elems[i]); } return this; }
666    public AlgebraicVector<T> Mul(AlgebraicVector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].Mul(a.elems[i]); } return this; }
667    public AlgebraicVector<T> Div(AlgebraicVector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].Div(a.elems[i]); } return this; }
668    public AlgebraicVector<T> AssignNeg(AlgebraicVector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].AssignNeg(a.elems[i]); } return this; }
669    public AlgebraicVector<T> Scale(double s) { for (int i = 0; i < elems.Length; ++i) { elems[i].Scale(s); } return this; }
670    public AlgebraicVector<T> Scale(T s) { for (int i = 0; i < elems.Length; ++i) { elems[i].Mul(s); } return this; }
671    public AlgebraicVector<T> AssignInv(AlgebraicVector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].AssignInv(a.elems[i]); } return this; }
672    public AlgebraicVector<T> AssignLog(AlgebraicVector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].AssignLog(a.elems[i]); } return this; }
673    public AlgebraicVector<T> AssignExp(AlgebraicVector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].AssignExp(a.elems[i]); } return this; }
674    public AlgebraicVector<T> AssignSin(AlgebraicVector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].AssignSin(a.elems[i]); } return this; }
675    public AlgebraicVector<T> AssignCos(AlgebraicVector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].AssignCos(a.elems[i]); } return this; }
676    public AlgebraicVector<T> AssignIntPower(AlgebraicVector<T> a, int p) { for (int i = 0; i < elems.Length; ++i) { elems[i].AssignIntPower(a.elems[i], p); } return this; }
677    public AlgebraicVector<T> AssignIntRoot(AlgebraicVector<T> a, int r) { for (int i = 0; i < elems.Length; ++i) { elems[i].AssignIntRoot(a.elems[i], r); } return this; }
678    public AlgebraicVector<T> AssignAbs(AlgebraicVector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].AssignAbs(a.elems[i]); } return this; }
679    public AlgebraicVector<T> AssignSgn(AlgebraicVector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].AssignSgn(a.elems[i]); } return this; }
680  }
681
682  */
683
684
685  /// <summary>
686  /// A sparse vector of algebraic types. Elements are accessed via a key of type K
687  /// </summary>
688  /// <typeparam name="K">Key type</typeparam>
689  /// <typeparam name="T">Element type</typeparam>
690  [DebuggerDisplay("SparseVector: {ToString()}")]
691  public sealed class AlgebraicSparseVector<K, T> : IAlgebraicType<AlgebraicSparseVector<K, T>> where T : IAlgebraicType<T> {
692    [DebuggerBrowsable(DebuggerBrowsableState.Never)]
693    private Dictionary<K, T> elems;
694    public IReadOnlyDictionary<K, T> Elements => elems;
695
696
697    public AlgebraicSparseVector(AlgebraicSparseVector<K, T> original) {
698      elems = original.elems.ToDictionary(kvp => kvp.Key, kvp => kvp.Value.Clone());
699    }
700
701    /// <summary>
702    ///
703    /// </summary>
704    /// <param name="keys"></param>
705    /// <param name="values">values are cloned</param>
706    public AlgebraicSparseVector(K[] keys, T[] values) {
707      if (keys.Length != values.Length) throw new ArgumentException("lengths of keys and values doesn't match in SparseVector");
708      elems = new Dictionary<K, T>(keys.Length);
709      for (int i = 0; i < keys.Length; ++i) {
710        elems.Add(keys[i], values[i].Clone());
711      }
712    }
713
714    public AlgebraicSparseVector() {
715      this.elems = new Dictionary<K, T>();
716    }
717
718    // keep only elements from a
719    private void AssignFromSource(AlgebraicSparseVector<K, T> a, Func<T, T, T> mapAssign) {
720      // remove elems from this which don't occur in a
721      List<K> keysToRemove = new List<K>();
722      foreach (var kvp in elems) {
723        if (!a.elems.ContainsKey(kvp.Key)) keysToRemove.Add(kvp.Key);
724      }
725      foreach (var o in keysToRemove) elems.Remove(o); // -> zero
726
727      foreach (var kvp in a.elems) {
728        if (elems.TryGetValue(kvp.Key, out T value))
729          mapAssign(kvp.Value, value);
730        else
731          elems.Add(kvp.Key, mapAssign(kvp.Value, kvp.Value.Zero));
732      }
733    }
734
735    [DebuggerBrowsable(DebuggerBrowsableState.Never)]
736    public AlgebraicSparseVector<K, T> Zero => new AlgebraicSparseVector<K, T>();
737    [DebuggerBrowsable(DebuggerBrowsableState.Never)]
738    public AlgebraicSparseVector<K, T> One => throw new NotSupportedException();
739
740    public AlgebraicSparseVector<K, T> Scale(T s) { foreach (var kvp in elems) { kvp.Value.Mul(s); } return this; }
741    public AlgebraicSparseVector<K, T> Scale(double s) { foreach (var kvp in elems) { kvp.Value.Scale(s); } return this; }
742
743    public AlgebraicSparseVector<K, T> Assign(AlgebraicSparseVector<K, T> a) { elems.Clear(); AssignFromSource(a, (src, dest) => dest.Assign(src)); return this; }
744    public AlgebraicSparseVector<K, T> AssignInv(AlgebraicSparseVector<K, T> a) { AssignFromSource(a, (src, dest) => dest.AssignInv(src)); return this; }
745    public AlgebraicSparseVector<K, T> AssignNeg(AlgebraicSparseVector<K, T> a) { AssignFromSource(a, (src, dest) => dest.AssignNeg(src)); return this; }
746    public AlgebraicSparseVector<K, T> AssignLog(AlgebraicSparseVector<K, T> a) { AssignFromSource(a, (src, dest) => dest.AssignLog(src)); return this; }
747    public AlgebraicSparseVector<K, T> AssignExp(AlgebraicSparseVector<K, T> a) { AssignFromSource(a, (src, dest) => dest.AssignExp(src)); return this; }
748    public AlgebraicSparseVector<K, T> AssignIntPower(AlgebraicSparseVector<K, T> a, int p) { AssignFromSource(a, (src, dest) => dest.AssignIntPower(src, p)); return this; }
749    public AlgebraicSparseVector<K, T> AssignIntRoot(AlgebraicSparseVector<K, T> a, int r) { AssignFromSource(a, (src, dest) => dest.AssignIntRoot(src, r)); return this; }
750    public AlgebraicSparseVector<K, T> AssignSin(AlgebraicSparseVector<K, T> a) { AssignFromSource(a, (src, dest) => dest.AssignSin(src)); return this; }
751    public AlgebraicSparseVector<K, T> AssignCos(AlgebraicSparseVector<K, T> a) { AssignFromSource(a, (src, dest) => dest.AssignCos(src)); return this; }
752    public AlgebraicSparseVector<K, T> AssignTanh(AlgebraicSparseVector<K, T> a) { AssignFromSource(a, (src, dest) => dest.AssignTanh(src)); return this; }
753    public AlgebraicSparseVector<K, T> AssignAbs(AlgebraicSparseVector<K, T> a) { AssignFromSource(a, (src, dest) => dest.AssignAbs(src)); return this; }
754    public AlgebraicSparseVector<K, T> AssignSgn(AlgebraicSparseVector<K, T> a) { AssignFromSource(a, (src, dest) => dest.AssignSgn(src)); return this; }
755    public AlgebraicSparseVector<K, T> Add(AlgebraicSparseVector<K, T> a) {
756      foreach (var kvp in a.elems) {
757        if (elems.TryGetValue(kvp.Key, out T value))
758          value.Add(kvp.Value);
759        else
760          elems.Add(kvp.Key, kvp.Value.Clone()); // 0 + a
761      }
762      return this;
763    }
764
765    public AlgebraicSparseVector<K, T> Sub(AlgebraicSparseVector<K, T> a) {
766      foreach (var kvp in a.elems) {
767        if (elems.TryGetValue(kvp.Key, out T value))
768          value.Sub(kvp.Value);
769        else
770          elems.Add(kvp.Key, kvp.Value.Zero.Sub(kvp.Value)); // 0 - a
771      }
772      return this;
773    }
774
775    public AlgebraicSparseVector<K, T> Mul(AlgebraicSparseVector<K, T> a) {
776      var keys = elems.Keys.ToArray();
777      foreach (var k in keys) if (!a.elems.ContainsKey(k)) elems.Remove(k); // 0 * a => 0
778      foreach (var kvp in a.elems) {
779        if (elems.TryGetValue(kvp.Key, out T value))
780          value.Mul(kvp.Value); // this * a
781      }
782      return this;
783    }
784
785    public AlgebraicSparseVector<K, T> Div(AlgebraicSparseVector<K, T> a) {
786      return Mul(a.Clone().Inv());
787    }
788
789    public AlgebraicSparseVector<K, T> AssignMin(AlgebraicSparseVector<K, T> other) {
790      // assumes that keys without a matching key in other are zero and vice versa
791      foreach (var kvp in elems) if (!other.elems.ContainsKey(kvp.Key)) kvp.Value.AssignMin(kvp.Value.Zero); // min(v, 0)
792      foreach (var kvp in other.elems) {
793        if (elems.TryGetValue(kvp.Key, out T value))
794          value.AssignMin(kvp.Value);
795        else
796          elems.Add(kvp.Key, kvp.Value.Zero.AssignMin(kvp.Value));
797      }
798      return this;
799    }
800
801    public AlgebraicSparseVector<K, T> AssignMax(AlgebraicSparseVector<K, T> other) {
802      // assumes that keys without a matching key in other are zero and vice versa
803      foreach (var kvp in elems) if (!other.elems.ContainsKey(kvp.Key)) kvp.Value.AssignMax(kvp.Value.Zero); // max(v, 0)
804      foreach (var kvp in other.elems) {
805        if (elems.TryGetValue(kvp.Key, out T value))
806          value.AssignMax(kvp.Value);
807        else
808          elems.Add(kvp.Key, kvp.Value.Zero.AssignMax(kvp.Value));
809      }
810      return this;
811    }
812
813
814    public AlgebraicSparseVector<K, T> Clone() {
815      return new AlgebraicSparseVector<K, T>(this);
816    }
817
818    public override string ToString() {
819      return "[" + string.Join(" ", elems.Select(kvp => kvp.Key + ": " + kvp.Value)) + "]";
820    }
821  }
822
823  // this is our own implementation of interval arithmetic
824  // for a well worked out definition of interval operations for IEEE reals see:
825  // Stahl: Interval Methods for Bounding the Range of Polynomials and Solving Systems of Nonlinear Equations, Dissertation, JKU, 1995
826  [DebuggerDisplay("[{low.Value}..{high.Value}]")]
827  public class AlgebraicInterval : IAlgebraicType<AlgebraicInterval> {
828    [DebuggerBrowsable(DebuggerBrowsableState.Never)]
829    private MultivariateDual<AlgebraicDouble> low;
830    public MultivariateDual<AlgebraicDouble> LowerBound => low.Clone();
831
832    [DebuggerBrowsable(DebuggerBrowsableState.Never)]
833    private MultivariateDual<AlgebraicDouble> high;
834    public MultivariateDual<AlgebraicDouble> UpperBound => high.Clone();
835
836
837    public AlgebraicInterval() : this(double.NegativeInfinity, double.PositiveInfinity) { }
838
839    public AlgebraicInterval(MultivariateDual<AlgebraicDouble> low, MultivariateDual<AlgebraicDouble> high) {
840      this.low = low.Clone();
841      this.high = high.Clone();
842    }
843
844    public AlgebraicInterval(double low, double high) {
845      this.low = new MultivariateDual<AlgebraicDouble>(new AlgebraicDouble(low));
846      this.high = new MultivariateDual<AlgebraicDouble>(new AlgebraicDouble(high));
847    }
848
849    [DebuggerBrowsable(DebuggerBrowsableState.Never)]
850    public AlgebraicInterval Zero => new AlgebraicInterval(0.0, 0.0);
851    [DebuggerBrowsable(DebuggerBrowsableState.Never)]
852    public AlgebraicInterval One => new AlgebraicInterval(1.0, 1.0);
853
854    public AlgebraicInterval Add(AlgebraicInterval a) {
855      low.Add(a.low);
856      high.Add(a.high);
857      return this;
858    }
859
860    public AlgebraicInterval Mul(AlgebraicInterval a) {
861      var v1 = low.Clone().Mul(a.low);
862      var v2 = low.Clone().Mul(a.high);
863      var v3 = high.Clone().Mul(a.low);
864      var v4 = high.Clone().Mul(a.high);
865
866      low = Min(Min(v1, v2), Min(v3, v4));
867      high = Max(Max(v1, v2), Max(v3, v4));
868
869      return this;
870    }
871
872
873    private static MultivariateDual<AlgebraicDouble> Min(MultivariateDual<AlgebraicDouble> a, MultivariateDual<AlgebraicDouble> b) {
874      return a.Value < b.Value ? a : b;
875    }
876    private static MultivariateDual<AlgebraicDouble> Max(MultivariateDual<AlgebraicDouble> a, MultivariateDual<AlgebraicDouble> b) {
877      return a.Value > b.Value ? a : b;
878    }
879
880    public AlgebraicInterval Assign(AlgebraicInterval a) {
881      low = a.low;
882      high = a.high;
883      return this;
884    }
885
886    public AlgebraicInterval AssignCos(AlgebraicInterval a) {
887      return AssignSin(a.Clone().Add(new AlgebraicInterval(Math.PI / 2, Math.PI / 2)));
888    }
889
890    public AlgebraicInterval Div(AlgebraicInterval a) {
891      if (a.Contains(0.0)) {
892        if (a.low.Value.Value == 0 && a.high.Value.Value == 0) {
893          if (this.low.Value >= 0) {
894            // pos / 0
895          } else if (this.high.Value <= 0) {
896            // neg / 0
897          } else {
898            low = new MultivariateDual<AlgebraicDouble>(double.NegativeInfinity);
899            high = new MultivariateDual<AlgebraicDouble>(double.PositiveInfinity);
900          }
901        } else if (a.low.Value.Value >= 0) {
902          // a is positive
903          Mul(new AlgebraicInterval(a.Clone().high.Inv(), new MultivariateDual<AlgebraicDouble>(double.PositiveInfinity)));
904        } else if (a.high.Value <= 0) {
905          // a is negative
906          Mul(new AlgebraicInterval(new MultivariateDual<AlgebraicDouble>(double.NegativeInfinity), a.low.Clone().Inv()));
907        } else {
908          // a is interval over zero
909          low = new MultivariateDual<AlgebraicDouble>(double.NegativeInfinity);
910          high = new MultivariateDual<AlgebraicDouble>(double.PositiveInfinity);
911        }
912      } else {
913        Mul(new AlgebraicInterval(a.high.Clone().Inv(), a.low.Clone().Inv())); // inverting leads to inverse roles of high and low
914      }
915      return this;
916    }
917
918    public AlgebraicInterval AssignExp(AlgebraicInterval a) {
919      low.AssignExp(a.low);
920      high.AssignExp(a.high);
921      return this;
922    }
923
924    // tanh is a bijective function
925    public AlgebraicInterval AssignTanh(AlgebraicInterval a) {
926      low.AssignTanh(a.low);
927      high.AssignTanh(a.high);
928      return this;
929    }
930
931    public AlgebraicInterval AssignIntPower(AlgebraicInterval a, int p) {
932      if (p < 0) {  // x^-3 == 1/(x^3)
933        AssignIntPower(a, -p);
934        return AssignInv(this);
935      } else if (p == 0) {
936        if (a.Contains(0.0)) {
937          // => 0^0 = 0 ; might not be relevant
938          low = new MultivariateDual<AlgebraicDouble>(0.0);
939          high = new MultivariateDual<AlgebraicDouble>(1.0);
940          return this;
941        } else {
942          // => 1
943          low = new MultivariateDual<AlgebraicDouble>(1.0);
944          high = new MultivariateDual<AlgebraicDouble>(1.0);
945          return this;
946        }
947      } else if (p == 1) return this;
948      else if (p % 2 == 0) {
949        // p is even => interval must be positive
950        if (a.Contains(0.0)) {
951          low = new MultivariateDual<AlgebraicDouble>(0.0);
952          high = a.low.IntPower(p).AssignMax(a.high.IntPower(p));
953        } else {
954          var lowPower = a.low.IntPower(p);
955          var highPower = a.high.IntPower(p);
956          low = lowPower.AssignMin(highPower);
957          high = lowPower.AssignMax(highPower);
958        }
959      } else {
960        // p is uneven
961        if (a.Contains(0.0)) {
962          low.AssignIntPower(a.low, p);
963          high.AssignIntPower(a.high, p);
964        } else {
965          var lowPower = a.low.IntPower(p);
966          var highPower = a.high.IntPower(p);
967          low = lowPower.AssignMin(highPower);
968          high = lowPower.AssignMax(highPower);
969        }
970      }
971      return this;
972    }
973
974    public AlgebraicInterval AssignIntRoot(AlgebraicInterval a, int r) {
975      if (r == 0) { low = new MultivariateDual<AlgebraicDouble>(double.NaN); high = new MultivariateDual<AlgebraicDouble>(double.NaN); return this; }
976      if (r == 1) return this;
977      if (r < 0) {
978        // x^ (-1/2) = 1 / (x^(1/2))
979        AssignIntRoot(a, -r);
980        return AssignInv(this);
981      } else {
982        // root only defined for positive arguments for even roots
983        if (r % 2 == 0 && a.LowerBound.Value.Value < 0) {
984          low = new MultivariateDual<AlgebraicDouble>(double.NaN);
985          high = new MultivariateDual<AlgebraicDouble>(double.NaN);
986          return this;
987        } else {
988          low.AssignIntRoot(a.low, r);
989          high.AssignIntRoot(a.high, r);
990          return this;
991        }
992      }
993    }
994
995    public AlgebraicInterval AssignInv(AlgebraicInterval a) {
996      low = new MultivariateDual<AlgebraicDouble>(1.0);
997      high = new MultivariateDual<AlgebraicDouble>(1.0);
998      return Div(a);
999    }
1000
1001    public AlgebraicInterval AssignLog(AlgebraicInterval a) {
1002      low.AssignLog(a.low);
1003      high.AssignLog(a.high);
1004      return this;
1005    }
1006
1007    public AlgebraicInterval AssignNeg(AlgebraicInterval a) {
1008      low.AssignNeg(a.high);
1009      high.AssignNeg(a.low);
1010      return this;
1011    }
1012
1013    public AlgebraicInterval Scale(double s) {
1014      low.Scale(s);
1015      high.Scale(s);
1016      if (s < 0) {
1017        var t = low;
1018        low = high;
1019        high = t;
1020      }
1021      return this;
1022    }
1023
1024    public AlgebraicInterval AssignSin(AlgebraicInterval a) {
1025      var lower = a.LowerBound.Value.Value;
1026      var size = a.UpperBound.Value.Value - lower;
1027      if (size < 0) throw new InvalidProgramException(); // ASSERT interval >= 0;
1028
1029      if (size >= Math.PI * 2) {
1030        low = new MultivariateDual<AlgebraicDouble>(-1.0); // zero gradient
1031        high = new MultivariateDual<AlgebraicDouble>(1.0);
1032        return this;
1033      }
1034
1035      // assume low and high are in the same quadrant
1036      low = Algebraic.Min(a.LowerBound.Clone().Sin(), a.UpperBound.Clone().Sin());
1037      high = Algebraic.Max(a.LowerBound.Clone().Sin(), a.UpperBound.Clone().Sin());
1038
1039      // override min and max if necessary
1040
1041      // shift interval 'a' into the range [-2pi .. 2pi] without changing the size of the interval to simplify the checks
1042      lower = lower % (2 * Math.PI); // lower in [-2pi .. 2pi]     
1043
1044      // handle min = -1 and max = 1 cases explicitly
1045      var pi_2 = Math.PI / 2.0;
1046      var maxima = new double[] { -3 * pi_2, pi_2 };
1047      var minima = new double[] { -pi_2, 3 * pi_2 };
1048
1049      // override min and max if necessary
1050      if (maxima.Any(m => lower < m && lower + size > m)) {
1051        // max = 1
1052        high = new MultivariateDual<AlgebraicDouble>(1.0); // zero gradient
1053      }
1054
1055      if (minima.Any(m => lower < m && lower + size > m)) {
1056        // min = -1;
1057        low = new MultivariateDual<AlgebraicDouble>(-1.0); // zero gradient
1058      }
1059      return this;
1060    }
1061
1062    public AlgebraicInterval Sub(AlgebraicInterval a) {
1063      // [x1,x2] − [y1,y2] = [x1 − y2,x2 − y1]
1064      low.Sub(a.high);
1065      high.Sub(a.low);
1066      return this;
1067    }
1068
1069    public AlgebraicInterval Clone() {
1070      return new AlgebraicInterval(low, high);
1071    }
1072
1073    public bool Contains(double val) {
1074      return LowerBound.Value.Value <= val && val <= UpperBound.Value.Value;
1075    }
1076
1077    public AlgebraicInterval AssignAbs(AlgebraicInterval a) {
1078      if (a.Contains(0.0)) {
1079        var abslow = a.low.Clone().Abs();
1080        var abshigh = a.high.Clone().Abs();
1081        a.high.Assign(Algebraic.Max(abslow, abshigh));
1082        a.low.Assign(new MultivariateDual<AlgebraicDouble>(0.0)); // lost gradient for lower bound
1083      } else {
1084        var abslow = a.low.Clone().Abs();
1085        var abshigh = a.high.Clone().Abs();
1086        a.low.Assign(Algebraic.Min(abslow, abshigh));
1087        a.high.Assign(Algebraic.Max(abslow, abshigh));
1088      }
1089      return this;
1090    }
1091
1092    public AlgebraicInterval AssignSgn(AlgebraicInterval a) {
1093      low.AssignSgn(a.low);
1094      high.AssignSgn(a.high);
1095      return this;
1096    }
1097
1098    public AlgebraicInterval AssignMin(AlgebraicInterval other) {
1099      low.AssignMin(other.low);
1100      high.AssignMin(other.high);
1101      return this;
1102    }
1103
1104    public AlgebraicInterval AssignMax(AlgebraicInterval other) {
1105      low.AssignMax(other.low);
1106      high.AssignMax(other.high);
1107      return this;
1108    }
1109  }
1110
1111  public class Dual<V> : IAlgebraicType<Dual<V>>
1112    where V : IAlgebraicType<V> {
1113    [DebuggerBrowsable(DebuggerBrowsableState.Never)]
1114    private V v;
1115    public V Value => v;
1116
1117    [DebuggerBrowsable(DebuggerBrowsableState.Never)]
1118    private V dv;
1119    public V Derivative => dv;
1120
1121    public Dual(V v, V dv) { this.v = v; this.dv = dv; }
1122
1123    [DebuggerBrowsable(DebuggerBrowsableState.Never)]
1124    public Dual<V> Zero => new Dual<V>(Value.Zero, Derivative.Zero);
1125    [DebuggerBrowsable(DebuggerBrowsableState.Never)]
1126    public Dual<V> One => new Dual<V>(Value.One, Derivative.Zero);
1127
1128    public Dual<V> Assign(Dual<V> a) { v.Assign(a.v); dv.Assign(a.dv); return this; }
1129    public Dual<V> Scale(double s) { v.Scale(s); dv.Scale(s); return this; }
1130    public Dual<V> Add(Dual<V> a) { v.Add(a.v); dv.Add(a.dv); return this; }
1131    public Dual<V> Sub(Dual<V> a) { v.Sub(a.v); dv.Sub(a.dv); return this; }
1132    public Dual<V> AssignNeg(Dual<V> a) { v.AssignNeg(a.v); dv.AssignNeg(a.dv); return this; }
1133    public Dual<V> AssignInv(Dual<V> a) { v.AssignInv(a.v); dv.AssignNeg(a.dv).Mul(v).Mul(v); return this; } // (1/f(x))' = - f(x)' / f(x)^2
1134
1135    // (a(x) * b(x))' = b(x)*a(x)' + b(x)'*a(x);
1136    public Dual<V> Mul(Dual<V> a) {
1137      var t1 = a.dv.Clone().Mul(v);
1138      var t2 = dv.Clone().Mul(a.v);
1139      dv.Assign(t1).Add(t2);
1140
1141      v.Mul(a.v);
1142      return this;
1143    }
1144    public Dual<V> Div(Dual<V> a) { Mul(a.Inv()); return this; }
1145
1146    public Dual<V> AssignExp(Dual<V> a) { v.AssignExp(a.v); dv.Assign(a.dv).Mul(v); return this; } // exp(f(x)) = exp(f(x))*f(x)'
1147    public Dual<V> AssignLog(Dual<V> a) { v.AssignLog(a.v); dv.Assign(a.dv).Div(a.v); return this; }     // log(x)' = 1/f(x) * f(x)'
1148
1149    public Dual<V> AssignIntPower(Dual<V> a, int p) { v.AssignIntPower(a.v, p); dv.Assign(a.dv).Scale(p).Mul(a.v.Clone().IntPower(p - 1)); return this; }
1150    public Dual<V> AssignIntRoot(Dual<V> a, int r) { v.AssignIntRoot(a.v, r); dv.Assign(a.dv).Scale(1.0 / r).Mul(a.v.IntRoot(r - 1)); return this; }
1151
1152    public Dual<V> AssignSin(Dual<V> a) { v.AssignSin(a.v); dv.Assign(a.dv).Mul(a.v.Clone().Cos()); return this; }
1153    public Dual<V> AssignCos(Dual<V> a) { v.AssignCos(a.v); dv.AssignNeg(a.dv).Mul(a.v.Clone().Sin()); return this; }
1154    public Dual<V> AssignTanh(Dual<V> a) { v.AssignTanh(a.v); dv.Assign(a.dv.Mul(v.Clone().IntPower(2).Neg().Add(Value.One))); return this; }
1155
1156    public Dual<V> AssignAbs(Dual<V> a) { v.AssignAbs(a.v); dv.Assign(a.dv).Mul(a.v.Clone().Sgn()); return this; }       // abs(f(x))' = f(x)*f'(x) / |f(x)|     
1157    public Dual<V> AssignSgn(Dual<V> a) { v.AssignSgn(a.v); dv.Assign(a.dv.Zero); return this; }
1158
1159    public Dual<V> Clone() { return new Dual<V>(v.Clone(), dv.Clone()); }
1160
1161    public Dual<V> AssignMin(Dual<V> other) {
1162      throw new NotImplementedException();
1163    }
1164
1165    public Dual<V> AssignMax(Dual<V> other) {
1166      throw new NotImplementedException();
1167    }
1168  }
1169
1170  /// <summary>
1171  /// An algebraic type which has a value as well as the partial derivatives of the value over multiple variables.
1172  /// </summary>
1173  /// <typeparam name="V"></typeparam>
1174  [DebuggerDisplay("v={Value}; dv={dv}")]
1175  public class MultivariateDual<V> : IAlgebraicType<MultivariateDual<V>> where V : IAlgebraicType<V>, new() {
1176    [DebuggerBrowsable(DebuggerBrowsableState.Never)]
1177    private V v;
1178    public V Value => v;
1179
1180    [DebuggerBrowsable(DebuggerBrowsableState.Never)]
1181    private AlgebraicSparseVector<object, V> dv;
1182    public AlgebraicSparseVector<object, V> Gradient => dv; // <key,value> partial derivative identified via the key
1183
1184    private MultivariateDual(MultivariateDual<V> orig) { this.v = orig.v.Clone(); this.dv = orig.dv.Clone(); }
1185
1186    /// <summary>
1187    /// Constructor without partial derivative
1188    /// </summary>
1189    /// <param name="v"></param>
1190    public MultivariateDual(V v) { this.v = v.Clone(); this.dv = new AlgebraicSparseVector<object, V>(); }
1191
1192    /// <summary>
1193    /// Constructor for multiple partial derivatives
1194    /// </summary>
1195    /// <param name="v"></param>
1196    /// <param name="keys"></param>
1197    /// <param name="dv"></param>
1198    public MultivariateDual(V v, object[] keys, V[] dv) { this.v = v.Clone(); this.dv = new AlgebraicSparseVector<object, V>(keys, dv); }
1199
1200    /// <summary>
1201    /// Constructor for a single partial derivative
1202    /// </summary>
1203    /// <param name="v"></param>
1204    /// <param name="key"></param>
1205    /// <param name="dv"></param>
1206    public MultivariateDual(V v, object key, V dv) { this.v = v.Clone(); this.dv = new AlgebraicSparseVector<object, V>(new[] { key }, new[] { dv }); }
1207
1208    /// <summary>
1209    /// Constructor with a given value and gradient. For internal use.
1210    /// </summary>
1211    /// <param name="v">The value (not cloned).</param>
1212    /// <param name="gradient">The gradient (not cloned).</param>
1213    internal MultivariateDual(V v, AlgebraicSparseVector<object, V> gradient) { this.v = v; this.dv = gradient; }
1214
1215    public MultivariateDual<V> Clone() { return new MultivariateDual<V>(this); }
1216
1217    public MultivariateDual<V> Zero => new MultivariateDual<V>(Value.Zero, Gradient.Zero);
1218    public MultivariateDual<V> One => new MultivariateDual<V>(Value.One, Gradient.Zero);
1219
1220    public MultivariateDual<V> Scale(double s) { v.Scale(s); dv.Scale(s); return this; }
1221
1222    public MultivariateDual<V> Add(MultivariateDual<V> a) { v.Add(a.v); dv.Add(a.dv); return this; }
1223    public MultivariateDual<V> Sub(MultivariateDual<V> a) { v.Sub(a.v); dv.Sub(a.dv); return this; }
1224    public MultivariateDual<V> Assign(MultivariateDual<V> a) { v.Assign(a.v); dv.Assign(a.dv); return this; }
1225    public MultivariateDual<V> Mul(MultivariateDual<V> a) {
1226      // (a(x) * b(x))' = b(x)*a(x)' + b(x)'*a(x);
1227      var t1 = a.dv.Clone().Scale(v);
1228      var t2 = dv.Clone().Scale(a.v);
1229      dv.Assign(t1).Add(t2);
1230
1231      v.Mul(a.v);
1232      return this;
1233    }
1234    public MultivariateDual<V> Div(MultivariateDual<V> a) { v.Div(a.v); dv.Mul(a.dv.Inv()); return this; }
1235    public MultivariateDual<V> AssignNeg(MultivariateDual<V> a) { v.AssignNeg(a.v); dv.AssignNeg(a.dv); return this; }
1236    public MultivariateDual<V> AssignInv(MultivariateDual<V> a) { v.AssignInv(a.v); dv.AssignNeg(a.dv).Scale(v).Scale(v); return this; }   // (1/f(x))' = - f(x)' / f(x)^2
1237
1238    public MultivariateDual<V> AssignSin(MultivariateDual<V> a) { v.AssignSin(a.v); dv.Assign(a.dv).Scale(a.v.Clone().Cos()); return this; }
1239    public MultivariateDual<V> AssignCos(MultivariateDual<V> a) { v.AssignCos(a.v); dv.AssignNeg(a.dv).Scale(a.v.Clone().Sin()); return this; }
1240    public MultivariateDual<V> AssignTanh(MultivariateDual<V> a) { v.AssignTanh(a.v); dv.Assign(a.dv.Scale(v.Clone().IntPower(2).Neg().Add(Value.One))); return this; }     // tanh(f(x))' = f(x)'sech²(f(x)) = f(x)'(1 - tanh²(f(x)))
1241
1242    public MultivariateDual<V> AssignIntPower(MultivariateDual<V> a, int p) { v.AssignIntPower(a.v, p); dv.Assign(a.dv).Scale(p).Scale(a.v.Clone().IntPower(p - 1)); return this; }
1243    public MultivariateDual<V> AssignIntRoot(MultivariateDual<V> a, int r) { v.AssignIntRoot(a.v, r); dv.Assign(a.dv).Scale(1.0 / r).Scale(a.v.IntRoot(r - 1)); return this; }
1244
1245    public MultivariateDual<V> AssignExp(MultivariateDual<V> a) { v.AssignExp(a.v); dv.Assign(a.dv).Scale(v); return this; } // exp(f(x)) = exp(f(x))*f(x)'     
1246    public MultivariateDual<V> AssignLog(MultivariateDual<V> a) { v.AssignLog(a.v); dv.Assign(a.dv).Scale(a.v.Clone().Inv()); return this; }     // log(x)' = 1/f(x) * f(x)'
1247
1248    public MultivariateDual<V> AssignAbs(MultivariateDual<V> a) { v.AssignAbs(a.v); dv.Assign(a.dv).Scale(a.v.Clone().Sgn()); return this; }      // abs(f(x))' = f(x)*f'(x) / |f(x)|  doesn't work for intervals     
1249    public MultivariateDual<V> AssignSgn(MultivariateDual<V> a) { v.AssignSgn(a.v); dv = a.dv.Zero; return this; } // sign(f(x))' = 0;     
1250
1251    public MultivariateDual<V> AssignMin(MultivariateDual<V> other) {
1252      XXX
1253    }
1254
1255    public MultivariateDual<V> AssignMax(MultivariateDual<V> other) {
1256      XXX
1257    }
1258  }
1259}
Note: See TracBrowser for help on using the repository browser.