Free cookie consent management tool by TermsFeed Policy Generator

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

Last change on this file since 16744 was 16744, checked in by gkronber, 5 years ago

#2994: slightly changed calculation of integer powers for intervals and added unit tests.

File size: 50.4 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          default: throw new ArgumentException($"Unknown opcode {code[c].opcode}");
116        }
117      }
118      return code[0].value;
119    }
120
121    protected Instruction[] Compile(ISymbolicExpressionTree tree) {
122      var root = tree.Root.GetSubtree(0).GetSubtree(0);
123      var code = new Instruction[root.GetLength()];
124      if (root.SubtreeCount > ushort.MaxValue) throw new ArgumentException("Number of subtrees is too big (>65.535)");
125      int c = 1, i = 0;
126      foreach (var node in root.IterateNodesBreadth()) {
127        if (node.SubtreeCount > ushort.MaxValue) throw new ArgumentException("Number of subtrees is too big (>65.535)");
128        code[i] = new Instruction {
129          opcode = OpCodes.MapSymbolToOpCode(node),
130          narg = (ushort)node.SubtreeCount,
131          childIndex = c
132        };
133        if (node is VariableTreeNode variable) {
134          InitializeTerminalInstruction(ref code[i], variable);
135        } else if (node is ConstantTreeNode constant) {
136          InitializeTerminalInstruction(ref code[i], constant);
137        } else {
138          InitializeInternalInstruction(ref code[i], node);
139        }
140        c += node.SubtreeCount;
141        ++i;
142      }
143      return code;
144    }
145
146    protected abstract void InitializeTerminalInstruction(ref Instruction instruction, ConstantTreeNode constant);
147    protected abstract void InitializeTerminalInstruction(ref Instruction instruction, VariableTreeNode variable);
148    protected abstract void InitializeInternalInstruction(ref Instruction instruction, ISymbolicExpressionTreeNode node);
149
150    protected abstract void LoadVariable(Instruction a);
151
152  }
153
154
155  public sealed class VectorEvaluator : Interpreter<AlgebraicDoubleVector> {
156    private const int BATCHSIZE = 128;
157    [ThreadStatic]
158    private Dictionary<string, double[]> cachedData;
159
160    [ThreadStatic]
161    private IDataset dataset;
162
163    [ThreadStatic]
164    private int rowIndex;
165
166    [ThreadStatic]
167    private int[] rows;
168
169    private void InitCache(IDataset dataset) {
170      this.dataset = dataset;
171      cachedData = new Dictionary<string, double[]>();
172      foreach (var v in dataset.DoubleVariables) {
173        cachedData[v] = dataset.GetReadOnlyDoubleValues(v).ToArray();
174      }
175    }
176
177    public double[] Evaluate(ISymbolicExpressionTree tree, IDataset dataset, int[] rows) {
178      if (cachedData == null || this.dataset != dataset) {
179        InitCache(dataset);
180      }
181
182      this.rows = rows;
183      var code = Compile(tree);
184      var remainingRows = rows.Length % BATCHSIZE;
185      var roundedTotal = rows.Length - remainingRows;
186
187      var result = new double[rows.Length];
188
189      for (rowIndex = 0; rowIndex < roundedTotal; rowIndex += BATCHSIZE) {
190        Evaluate(code);
191        code[0].value.CopyTo(result, rowIndex, BATCHSIZE);
192      }
193
194      if (remainingRows > 0) {
195        Evaluate(code);
196        code[0].value.CopyTo(result, roundedTotal, remainingRows);
197      }
198
199      return result;
200    }
201
202    protected override void InitializeTerminalInstruction(ref Instruction instruction, ConstantTreeNode constant) {
203      instruction.dblVal = constant.Value;
204      instruction.value = new AlgebraicDoubleVector(BATCHSIZE);
205      instruction.value.AssignConstant(instruction.dblVal);
206    }
207
208    protected override void InitializeTerminalInstruction(ref Instruction instruction, VariableTreeNode variable) {
209      instruction.dblVal = variable.Weight;
210      instruction.value = new AlgebraicDoubleVector(BATCHSIZE);
211      if (cachedData.ContainsKey(variable.VariableName)) {
212        instruction.data = cachedData[variable.VariableName];
213      } else {
214        instruction.data = dataset.GetDoubleValues(variable.VariableName).ToArray();
215        cachedData[variable.VariableName] = (double[])instruction.data;
216      }
217    }
218
219    protected override void InitializeInternalInstruction(ref Instruction instruction, ISymbolicExpressionTreeNode node) {
220      instruction.value = new AlgebraicDoubleVector(BATCHSIZE);
221    }
222
223    protected override void LoadVariable(Instruction a) {
224      var data = (double[])a.data;
225      for (int i = rowIndex; i < rows.Length && i - rowIndex < BATCHSIZE; i++) a.value[i - rowIndex] = data[rows[i]];
226      a.value.Scale(a.dblVal);
227    }
228  }
229
230  public sealed class VectorAutoDiffEvaluator : Interpreter<MultivariateDual<AlgebraicDoubleVector>> {
231    private const int BATCHSIZE = 128;
232    [ThreadStatic]
233    private Dictionary<string, double[]> cachedData;
234
235    [ThreadStatic]
236    private IDataset dataset;
237
238    [ThreadStatic]
239    private int rowIndex;
240
241    [ThreadStatic]
242    private int[] rows;
243
244    [ThreadStatic]
245    private Dictionary<ISymbolicExpressionTreeNode, int> node2paramIdx;
246
247    private void InitCache(IDataset dataset) {
248      this.dataset = dataset;
249      cachedData = new Dictionary<string, double[]>();
250      foreach (var v in dataset.DoubleVariables) {
251        cachedData[v] = dataset.GetDoubleValues(v).ToArray();
252      }
253    }
254
255    /// <summary>
256    ///
257    /// </summary>
258    /// <param name="tree"></param>
259    /// <param name="dataset"></param>
260    /// <param name="rows"></param>
261    /// <param name="parameterNodes"></param>
262    /// <param name="fi">Function output. Must be allocated by the caller.</param>
263    /// <param name="jac">Jacobian matrix. Must be allocated by the caller.</param>
264    public void Evaluate(ISymbolicExpressionTree tree, IDataset dataset, int[] rows, ISymbolicExpressionTreeNode[] parameterNodes, double[] fi, double[,] jac) {
265      if (cachedData == null || this.dataset != dataset) {
266        InitCache(dataset);
267      }
268
269      int nParams = parameterNodes.Length;
270      node2paramIdx = new Dictionary<ISymbolicExpressionTreeNode, int>();
271      for (int i = 0; i < parameterNodes.Length; i++) node2paramIdx.Add(parameterNodes[i], i);
272
273      var code = Compile(tree);
274
275      var remainingRows = rows.Length % BATCHSIZE;
276      var roundedTotal = rows.Length - remainingRows;
277
278      this.rows = rows;
279
280      for (rowIndex = 0; rowIndex < roundedTotal; rowIndex += BATCHSIZE) {
281        Evaluate(code);
282        code[0].value.Value.CopyTo(fi, rowIndex, BATCHSIZE);
283
284        // TRANSPOSE into JAC
285        var g = code[0].value.Gradient;
286        for (int j = 0; j < nParams; ++j) {
287          if (g.Elements.TryGetValue(j, out AlgebraicDoubleVector v)) {
288            v.CopyColumnTo(jac, j, rowIndex, BATCHSIZE);
289          } else {
290            for (int r = 0; r < BATCHSIZE; r++) jac[rowIndex + r, j] = 0.0;
291          }
292        }
293      }
294
295      if (remainingRows > 0) {
296        Evaluate(code);
297        code[0].value.Value.CopyTo(fi, roundedTotal, remainingRows);
298
299        var g = code[0].value.Gradient;
300        for (int j = 0; j < nParams; ++j)
301          if (g.Elements.TryGetValue(j, out AlgebraicDoubleVector v)) {
302            v.CopyColumnTo(jac, j, roundedTotal, remainingRows);
303          } else {
304            for (int r = 0; r < remainingRows; r++) jac[roundedTotal + r, j] = 0.0;
305          }
306      }
307    }
308
309    protected override void InitializeInternalInstruction(ref Instruction instruction, ISymbolicExpressionTreeNode node) {
310      var zero = new AlgebraicDoubleVector(BATCHSIZE);
311      instruction.value = new MultivariateDual<AlgebraicDoubleVector>(zero);
312    }
313
314    protected override void InitializeTerminalInstruction(ref Instruction instruction, ConstantTreeNode constant) {
315      var g_arr = new double[BATCHSIZE];
316      if (node2paramIdx.TryGetValue(constant, out var paramIdx)) {
317        for (int i = 0; i < BATCHSIZE; i++) g_arr[i] = 1.0;
318        var g = new AlgebraicDoubleVector(g_arr);
319        instruction.value = new MultivariateDual<AlgebraicDoubleVector>(new AlgebraicDoubleVector(BATCHSIZE), paramIdx, g); // only a single column for the gradient
320      } else {
321        instruction.value = new MultivariateDual<AlgebraicDoubleVector>(new AlgebraicDoubleVector(BATCHSIZE));
322      }
323
324      instruction.dblVal = constant.Value;
325      instruction.value.Value.AssignConstant(instruction.dblVal);
326    }
327
328    protected override void InitializeTerminalInstruction(ref Instruction instruction, VariableTreeNode variable) {
329      double[] data;
330      if (cachedData.ContainsKey(variable.VariableName)) {
331        data = cachedData[variable.VariableName];
332      } else {
333        data = dataset.GetReadOnlyDoubleValues(variable.VariableName).ToArray();
334        cachedData[variable.VariableName] = (double[])instruction.data;
335      }
336
337      var paramIdx = -1;
338      if (node2paramIdx.ContainsKey(variable)) {
339        paramIdx = node2paramIdx[variable];
340        var f = new AlgebraicDoubleVector(BATCHSIZE);
341        var g = new AlgebraicDoubleVector(BATCHSIZE);
342        instruction.value = new MultivariateDual<AlgebraicDoubleVector>(f, paramIdx, g);
343      } else {
344        var f = new AlgebraicDoubleVector(BATCHSIZE);
345        instruction.value = new MultivariateDual<AlgebraicDoubleVector>(f);
346      }
347
348      instruction.dblVal = variable.Weight;
349      instruction.data = new object[] { data, paramIdx };
350    }
351
352    protected override void LoadVariable(Instruction a) {
353      var paramIdx = (int)((object[])a.data)[1];
354      var data = (double[])((object[])a.data)[0];
355
356      for (int i = rowIndex; i < rows.Length && i - rowIndex < BATCHSIZE; i++) a.value.Value[i - rowIndex] = data[rows[i]];
357      a.value.Scale(a.dblVal);
358
359      if (paramIdx >= 0) {
360        // update gradient with variable values
361        var g = a.value.Gradient.Elements[paramIdx];
362        for (int i = rowIndex; i < rows.Length && i - rowIndex < BATCHSIZE; i++) {
363          g[i - rowIndex] = data[rows[i]];
364        }
365      }
366    }
367  }
368
369
370  public sealed class IntervalEvaluator : Interpreter<AlgebraicInterval> {
371    [ThreadStatic]
372    private Dictionary<string, Interval> intervals;
373
374    public Interval Evaluate(ISymbolicExpressionTree tree, Dictionary<string, Interval> intervals) {
375      this.intervals = intervals;
376      var code = Compile(tree);
377      Evaluate(code);
378      return new Interval(code[0].value.LowerBound.Value.Value, code[0].value.UpperBound.Value.Value);
379    }
380
381    public Interval Evaluate(ISymbolicExpressionTree tree, Dictionary<string, Interval> intervals, ISymbolicExpressionTreeNode[] paramNodes, out double[] lowerGradient, out double[] upperGradient) {
382      this.intervals = intervals;
383      var code = Compile(tree);
384      Evaluate(code);
385      lowerGradient = new double[paramNodes.Length];
386      upperGradient = new double[paramNodes.Length];
387      var l = code[0].value.LowerBound;
388      var u = code[0].value.UpperBound;
389      for (int i = 0; i < paramNodes.Length; ++i) {
390        if (paramNodes[i] == null) continue;
391        if (l.Gradient.Elements.TryGetValue(paramNodes[i], out AlgebraicDouble value)) lowerGradient[i] = value;
392        if (u.Gradient.Elements.TryGetValue(paramNodes[i], out value)) upperGradient[i] = value;
393      }
394      return new Interval(code[0].value.LowerBound.Value.Value, code[0].value.UpperBound.Value.Value);
395    }
396
397    protected override void InitializeInternalInstruction(ref Instruction instruction, ISymbolicExpressionTreeNode node) {
398      instruction.value = new AlgebraicInterval(0, 0);
399    }
400
401
402    protected override void InitializeTerminalInstruction(ref Instruction instruction, ConstantTreeNode constant) {
403      instruction.dblVal = constant.Value;
404      instruction.value = new AlgebraicInterval(
405        new MultivariateDual<AlgebraicDouble>(constant.Value, constant, 1.0),
406        new MultivariateDual<AlgebraicDouble>(constant.Value, constant, 1.0) // use node as key
407        );
408    }
409
410    protected override void InitializeTerminalInstruction(ref Instruction instruction, VariableTreeNode variable) {
411      instruction.dblVal = variable.Weight;
412      instruction.value = new AlgebraicInterval(
413        low: new MultivariateDual<AlgebraicDouble>(intervals[variable.VariableName].LowerBound, variable, intervals[variable.VariableName].LowerBound),  // bounds change by variable value d/dc (c I(var)) = I(var)
414        high: new MultivariateDual<AlgebraicDouble>(intervals[variable.VariableName].UpperBound, variable, intervals[variable.VariableName].UpperBound)
415        );
416    }
417
418    protected override void LoadVariable(Instruction a) {
419      // nothing to do
420    }
421  }
422
423  public interface IAlgebraicType<T> {
424    T Zero { get; }
425    T One { get; }
426
427    T AssignAbs(T a); // set this to assign abs(a)
428    T Assign(T a); // assign this to same value as a (copy!)
429    T AssignNeg(T a); // set this to negative(a)
430    T AssignInv(T a); // set this to inv(a);
431    T Scale(double s); // scale this with s
432    T Add(T a); // add a to this
433    T Sub(T a); // subtract a from this
434    T Mul(T a); // multiply this with a
435    T Div(T a); // divide this by a
436    T AssignLog(T a); // set this to log a
437    T AssignExp(T a); // set this to exp(a)
438    T AssignSin(T a); // set this to sin(a)
439    T AssignCos(T a); // set this to cos(a)
440    T AssignIntPower(T a, int p);
441    T AssignIntRoot(T a, int r);
442    T AssignSgn(T a); // set this to sign(a)
443    T Clone();
444  }
445
446  public static class Algebraic {
447    public static T Abs<T>(this T a) where T : IAlgebraicType<T> { a.AssignAbs(a.Clone()); return a; }
448    public static T Neg<T>(this T a) where T : IAlgebraicType<T> { a.AssignNeg(a.Clone()); return a; }
449    public static T Inv<T>(this T a) where T : IAlgebraicType<T> { a.AssignInv(a.Clone()); return a; }
450    public static T Log<T>(this T a) where T : IAlgebraicType<T> { a.AssignLog(a.Clone()); return a; }
451    public static T Exp<T>(this T a) where T : IAlgebraicType<T> { a.AssignExp(a.Clone()); return a; }
452    public static T Sin<T>(this T a) where T : IAlgebraicType<T> { a.AssignSin(a.Clone()); return a; }
453    public static T Cos<T>(this T a) where T : IAlgebraicType<T> { a.AssignCos(a.Clone()); return a; }
454    public static T Sgn<T>(this T a) where T : IAlgebraicType<T> { a.AssignSgn(a.Clone()); return a; }
455    public static T IntPower<T>(this T a, int p) where T : IAlgebraicType<T> { a.AssignIntPower(a.Clone(), p); return a; }
456    public static T IntRoot<T>(this T a, int r) where T : IAlgebraicType<T> { a.AssignIntRoot(a.Clone(), r); return a; }
457
458    public static T Max<T>(T a, T b) where T : IAlgebraicType<T> {
459      // ((a + b) + abs(b - a)) / 2
460      return a.Clone().Add(b).Add(b.Clone().Sub(a).Abs()).Scale(1.0 / 2.0);
461    }
462    public static T Min<T>(T a, T b) where T : IAlgebraicType<T> {
463      // ((a + b) - abs(a - b)) / 2
464      return a.Clone().Add(b).Sub(a.Clone().Sub(b).Abs()).Scale(1.0 / 2.0);
465    }
466  }
467
468
469  // algebraic type wrapper for a double value
470  [DebuggerDisplay("{Value}")]
471  public sealed class AlgebraicDouble : IAlgebraicType<AlgebraicDouble> {
472    public static implicit operator AlgebraicDouble(double value) { return new AlgebraicDouble(value); }
473    public static implicit operator double(AlgebraicDouble value) { return value.Value; }
474    public double Value;
475
476    [DebuggerBrowsable(DebuggerBrowsableState.Never)]
477    public AlgebraicDouble Zero => new AlgebraicDouble(0.0);
478    [DebuggerBrowsable(DebuggerBrowsableState.Never)]
479    public AlgebraicDouble One => new AlgebraicDouble(1.0);
480    public AlgebraicDouble() { }
481    public AlgebraicDouble(double value) { this.Value = value; }
482    public AlgebraicDouble Assign(AlgebraicDouble a) { Value = a.Value; return this; }
483    public AlgebraicDouble Add(AlgebraicDouble a) { Value += a.Value; return this; }
484    public AlgebraicDouble Sub(AlgebraicDouble a) { Value -= a.Value; return this; }
485    public AlgebraicDouble Mul(AlgebraicDouble a) { Value *= a.Value; return this; }
486    public AlgebraicDouble Div(AlgebraicDouble a) { Value /= a.Value; return this; }
487    public AlgebraicDouble Scale(double s) { Value *= s; return this; }
488    public AlgebraicDouble AssignInv(AlgebraicDouble a) { Value = 1.0 / a.Value; return this; }
489    public AlgebraicDouble AssignNeg(AlgebraicDouble a) { Value = -a.Value; return this; }
490    public AlgebraicDouble AssignSin(AlgebraicDouble a) { Value = Math.Sin(a.Value); return this; }
491    public AlgebraicDouble AssignCos(AlgebraicDouble a) { Value = Math.Cos(a.Value); return this; }
492    public AlgebraicDouble AssignLog(AlgebraicDouble a) { Value = a.Value <= 0 ? double.NegativeInfinity : Math.Log(a.Value); return this; } // alternative definiton of log to prevent NaN
493    public AlgebraicDouble AssignExp(AlgebraicDouble a) { Value = Math.Exp(a.Value); return this; }
494    public AlgebraicDouble AssignIntPower(AlgebraicDouble a, int p) { Value = Math.Pow(a.Value, p); return this; }
495    public AlgebraicDouble AssignIntRoot(AlgebraicDouble a, int r) { Value = Math.Pow(a.Value, 1.0 / r); return this; }
496    public AlgebraicDouble AssignAbs(AlgebraicDouble a) { Value = Math.Abs(a.Value); return this; }
497    public AlgebraicDouble AssignSgn(AlgebraicDouble a) { Value = double.IsNaN(a.Value) ? double.NaN : Math.Sign(a.Value); return this; }
498    public AlgebraicDouble Clone() { return new AlgebraicDouble(Value); }
499
500    public override string ToString() {
501      return Value.ToString();
502    }
503  }
504
505  // a simple vector as an algebraic type
506  [DebuggerDisplay("DoubleVector(len={Length}): {string.}")]
507  public class AlgebraicDoubleVector : IAlgebraicType<AlgebraicDoubleVector> {
508    private double[] arr;
509    public double this[int idx] { get { return arr[idx]; } set { arr[idx] = value; } }
510    public int Length => arr.Length;
511
512    public AlgebraicDoubleVector(int length) { arr = new double[length]; }
513
514    public AlgebraicDoubleVector() { }
515
516    /// <summary>
517    ///
518    /// </summary>
519    /// <param name="arr">array is not copied</param>
520    public AlgebraicDoubleVector(double[] arr) { this.arr = arr; }
521
522    [DebuggerBrowsable(DebuggerBrowsableState.Never)]
523    public AlgebraicDoubleVector Zero => new AlgebraicDoubleVector(new double[this.Length]); // must return vector of same length as this (therefore Zero is not static)
524    [DebuggerBrowsable(DebuggerBrowsableState.Never)]
525    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)
526    public AlgebraicDoubleVector Assign(AlgebraicDoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = a.arr[i]; } return this; }
527    public AlgebraicDoubleVector Add(AlgebraicDoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] += a.arr[i]; } return this; }
528    public AlgebraicDoubleVector Sub(AlgebraicDoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] -= a.arr[i]; } return this; }
529    public AlgebraicDoubleVector Mul(AlgebraicDoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] *= a.arr[i]; } return this; }
530    public AlgebraicDoubleVector Div(AlgebraicDoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] /= a.arr[i]; } return this; }
531    public AlgebraicDoubleVector AssignNeg(AlgebraicDoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = -a.arr[i]; } return this; }
532    public AlgebraicDoubleVector AssignInv(AlgebraicDoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = 1.0 / a.arr[i]; } return this; }
533    public AlgebraicDoubleVector Scale(double s) { for (int i = 0; i < arr.Length; ++i) { arr[i] *= s; } return this; }
534    public AlgebraicDoubleVector AssignLog(AlgebraicDoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = Math.Log(a.arr[i]); } return this; }
535    public AlgebraicDoubleVector AssignSin(AlgebraicDoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = Math.Sin(a.arr[i]); } return this; }
536    public AlgebraicDoubleVector AssignExp(AlgebraicDoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = Math.Exp(a.arr[i]); } return this; }
537    public AlgebraicDoubleVector AssignCos(AlgebraicDoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = Math.Cos(a.arr[i]); } return this; }
538    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; }
539    public AlgebraicDoubleVector AssignIntRoot(AlgebraicDoubleVector a, int r) { for (int i = 0; i < arr.Length; ++i) { arr[i] = Math.Pow(a.arr[i], 1.0 / r); } return this; }
540    public AlgebraicDoubleVector AssignAbs(AlgebraicDoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = Math.Abs(a.arr[i]); } return this; }
541    public AlgebraicDoubleVector AssignSgn(AlgebraicDoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = Math.Sign(a.arr[i]); } return this; }
542
543    public AlgebraicDoubleVector Clone() {
544      var v = new AlgebraicDoubleVector(this.arr.Length);
545      Array.Copy(arr, v.arr, v.arr.Length);
546      return v;
547    }
548
549    public AlgebraicDoubleVector AssignConstant(double constantValue) {
550      for (int i = 0; i < arr.Length; ++i) {
551        arr[i] = constantValue;
552      }
553      return this;
554    }
555
556    public void CopyTo(double[] dest, int idx, int length) {
557      Array.Copy(arr, 0, dest, idx, length);
558    }
559
560    public void CopyFrom(double[] data, int rowIndex) {
561      Array.Copy(data, rowIndex, arr, 0, Math.Min(arr.Length, data.Length - rowIndex));
562    }
563    public void CopyRowTo(double[,] dest, int row) {
564      for (int j = 0; j < arr.Length; ++j) dest[row, j] = arr[j];
565    }
566
567    internal void CopyColumnTo(double[,] dest, int column, int row, int len) {
568      for (int j = 0; j < len; ++j) dest[row + j, column] = arr[j];
569    }
570
571    public override string ToString() {
572      return "{" + string.Join(", ", arr.Take(Math.Max(5, arr.Length))) + (arr.Length > 5 ? "..." : string.Empty) + "}";
573    }
574  }
575
576
577  /*
578  // vectors of algebraic types
579  public sealed class AlgebraicVector<T> : IAlgebraicType<AlgebraicVector<T>> where T : IAlgebraicType<T>, new() {
580    private T[] elems;
581
582    public T this[int idx] { get { return elems[idx]; } set { elems[idx] = value; } }
583
584    public int Length => elems.Length;
585
586    private AlgebraicVector() { }
587
588    public AlgebraicVector(int len) { elems = new T[len]; }
589
590    /// <summary>
591    ///
592    /// </summary>
593    /// <param name="elems">The array is copied (element-wise clone)</param>
594    public AlgebraicVector(T[] elems) {
595      this.elems = new T[elems.Length];
596      for (int i = 0; i < elems.Length; ++i) { this.elems[i] = elems[i].Clone(); }
597    }
598
599    /// <summary>
600    ///
601    /// </summary>
602    /// <param name="elems">Array is not copied!</param>
603    /// <returns></returns>
604    public AlgebraicVector<T> FromArray(T[] elems) {
605      var v = new AlgebraicVector<T>();
606      v.elems = elems;
607      return v;
608    }
609
610    public void CopyTo(T[] dest) {
611      if (dest.Length != elems.Length) throw new InvalidOperationException("arr lengths do not match in Vector<T>.Copy");
612      Array.Copy(elems, dest, dest.Length);
613    }
614
615    public AlgebraicVector<T> Clone() { return new AlgebraicVector<T>(elems); }
616
617
618    [DebuggerBrowsable(DebuggerBrowsableState.Never)]
619    public AlgebraicVector<T> Zero => new AlgebraicVector<T>(Length);
620    [DebuggerBrowsable(DebuggerBrowsableState.Never)]
621    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; } }
622    public AlgebraicVector<T> Assign(AlgebraicVector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].Assign(a.elems[i]); } return this; }
623    public AlgebraicVector<T> Add(AlgebraicVector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].Add(a.elems[i]); } return this; }
624    public AlgebraicVector<T> Sub(AlgebraicVector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].Sub(a.elems[i]); } return this; }
625    public AlgebraicVector<T> Mul(AlgebraicVector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].Mul(a.elems[i]); } return this; }
626    public AlgebraicVector<T> Div(AlgebraicVector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].Div(a.elems[i]); } return this; }
627    public AlgebraicVector<T> AssignNeg(AlgebraicVector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].AssignNeg(a.elems[i]); } return this; }
628    public AlgebraicVector<T> Scale(double s) { for (int i = 0; i < elems.Length; ++i) { elems[i].Scale(s); } return this; }
629    public AlgebraicVector<T> Scale(T s) { for (int i = 0; i < elems.Length; ++i) { elems[i].Mul(s); } return this; }
630    public AlgebraicVector<T> AssignInv(AlgebraicVector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].AssignInv(a.elems[i]); } return this; }
631    public AlgebraicVector<T> AssignLog(AlgebraicVector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].AssignLog(a.elems[i]); } return this; }
632    public AlgebraicVector<T> AssignExp(AlgebraicVector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].AssignExp(a.elems[i]); } return this; }
633    public AlgebraicVector<T> AssignSin(AlgebraicVector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].AssignSin(a.elems[i]); } return this; }
634    public AlgebraicVector<T> AssignCos(AlgebraicVector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].AssignCos(a.elems[i]); } return this; }
635    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; }
636    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; }
637    public AlgebraicVector<T> AssignAbs(AlgebraicVector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].AssignAbs(a.elems[i]); } return this; }
638    public AlgebraicVector<T> AssignSgn(AlgebraicVector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].AssignSgn(a.elems[i]); } return this; }
639  }
640
641  */
642
643
644  /// <summary>
645  /// A sparse vector of algebraic types. Elements are accessed via a key of type K
646  /// </summary>
647  /// <typeparam name="K">Key type</typeparam>
648  /// <typeparam name="T">Element type</typeparam>
649  [DebuggerDisplay("SparseVector: {ToString()}")]
650  public sealed class AlgebraicSparseVector<K, T> : IAlgebraicType<AlgebraicSparseVector<K, T>> where T : IAlgebraicType<T> {
651    [DebuggerBrowsable(DebuggerBrowsableState.Never)]
652    private Dictionary<K, T> elems;
653    public IReadOnlyDictionary<K, T> Elements => elems;
654
655
656    public AlgebraicSparseVector(AlgebraicSparseVector<K, T> original) {
657      elems = original.elems.ToDictionary(kvp => kvp.Key, kvp => kvp.Value.Clone());
658    }
659
660    /// <summary>
661    ///
662    /// </summary>
663    /// <param name="keys"></param>
664    /// <param name="values">values are cloned</param>
665    public AlgebraicSparseVector(K[] keys, T[] values) {
666      if (keys.Length != values.Length) throw new ArgumentException("lengths of keys and values doesn't match in SparseVector");
667      elems = new Dictionary<K, T>(keys.Length);
668      for (int i = 0; i < keys.Length; ++i) {
669        elems.Add(keys[i], values[i].Clone());
670      }
671    }
672
673    public AlgebraicSparseVector() {
674      this.elems = new Dictionary<K, T>();
675    }
676
677    // keep only elements from a
678    private void AssignFromSource(AlgebraicSparseVector<K, T> a, Func<T, T, T> mapAssign) {
679      // remove elems from this which don't occur in a
680      List<K> keysToRemove = new List<K>();
681      foreach (var kvp in elems) {
682        if (!a.elems.ContainsKey(kvp.Key)) keysToRemove.Add(kvp.Key);
683      }
684      foreach (var o in keysToRemove) elems.Remove(o); // -> zero
685
686      foreach (var kvp in a.elems) {
687        if (elems.TryGetValue(kvp.Key, out T value))
688          mapAssign(kvp.Value, value);
689        else
690          elems.Add(kvp.Key, mapAssign(kvp.Value, kvp.Value.Zero));
691      }
692    }
693
694    [DebuggerBrowsable(DebuggerBrowsableState.Never)]
695    public AlgebraicSparseVector<K, T> Zero => new AlgebraicSparseVector<K, T>();
696    [DebuggerBrowsable(DebuggerBrowsableState.Never)]
697    public AlgebraicSparseVector<K, T> One => throw new NotSupportedException();
698
699    public AlgebraicSparseVector<K, T> Scale(T s) { foreach (var kvp in elems) { kvp.Value.Mul(s); } return this; }
700    public AlgebraicSparseVector<K, T> Scale(double s) { foreach (var kvp in elems) { kvp.Value.Scale(s); } return this; }
701
702    public AlgebraicSparseVector<K, T> Assign(AlgebraicSparseVector<K, T> a) { elems.Clear(); AssignFromSource(a, (src, dest) => dest.Assign(src)); return this; }
703    public AlgebraicSparseVector<K, T> AssignInv(AlgebraicSparseVector<K, T> a) { AssignFromSource(a, (src, dest) => dest.AssignInv(src)); return this; }
704    public AlgebraicSparseVector<K, T> AssignNeg(AlgebraicSparseVector<K, T> a) { AssignFromSource(a, (src, dest) => dest.AssignNeg(src)); return this; }
705    public AlgebraicSparseVector<K, T> AssignLog(AlgebraicSparseVector<K, T> a) { AssignFromSource(a, (src, dest) => dest.AssignLog(src)); return this; }
706    public AlgebraicSparseVector<K, T> AssignExp(AlgebraicSparseVector<K, T> a) { AssignFromSource(a, (src, dest) => dest.AssignExp(src)); return this; }
707    public AlgebraicSparseVector<K, T> AssignIntPower(AlgebraicSparseVector<K, T> a, int p) { AssignFromSource(a, (src, dest) => dest.AssignIntPower(src, p)); return this; }
708    public AlgebraicSparseVector<K, T> AssignIntRoot(AlgebraicSparseVector<K, T> a, int r) { AssignFromSource(a, (src, dest) => dest.AssignIntRoot(src, r)); return this; }
709    public AlgebraicSparseVector<K, T> AssignSin(AlgebraicSparseVector<K, T> a) { AssignFromSource(a, (src, dest) => dest.AssignSin(src)); return this; }
710    public AlgebraicSparseVector<K, T> AssignCos(AlgebraicSparseVector<K, T> a) { AssignFromSource(a, (src, dest) => dest.AssignCos(src)); return this; }
711    public AlgebraicSparseVector<K, T> AssignAbs(AlgebraicSparseVector<K, T> a) { AssignFromSource(a, (src, dest) => dest.AssignAbs(src)); return this; }
712    public AlgebraicSparseVector<K, T> AssignSgn(AlgebraicSparseVector<K, T> a) { AssignFromSource(a, (src, dest) => dest.AssignSgn(src)); return this; }
713    public AlgebraicSparseVector<K, T> Add(AlgebraicSparseVector<K, T> a) {
714      foreach (var kvp in a.elems) {
715        if (elems.TryGetValue(kvp.Key, out T value))
716          value.Add(kvp.Value);
717        else
718          elems.Add(kvp.Key, kvp.Value.Clone()); // 0 + a
719      }
720      return this;
721    }
722
723    public AlgebraicSparseVector<K, T> Sub(AlgebraicSparseVector<K, T> a) {
724      foreach (var kvp in a.elems) {
725        if (elems.TryGetValue(kvp.Key, out T value))
726          value.Sub(kvp.Value);
727        else
728          elems.Add(kvp.Key, kvp.Value.Zero.Sub(kvp.Value)); // 0 - a
729      }
730      return this;
731    }
732
733    public AlgebraicSparseVector<K, T> Mul(AlgebraicSparseVector<K, T> a) {
734      var keys = elems.Keys.ToArray();
735      foreach (var k in keys) if (!a.elems.ContainsKey(k)) elems.Remove(k); // 0 * a => 0
736      foreach (var kvp in a.elems) {
737        if (elems.TryGetValue(kvp.Key, out T value))
738          value.Mul(kvp.Value); // this * a
739      }
740      return this;
741    }
742
743    public AlgebraicSparseVector<K, T> Div(AlgebraicSparseVector<K, T> a) {
744      return Mul(a.Clone().Inv());
745    }
746
747    public AlgebraicSparseVector<K, T> Clone() {
748      return new AlgebraicSparseVector<K, T>(this);
749    }
750
751    public override string ToString() {
752      return "[" + string.Join(" ", elems.Select(kvp => kvp.Key + ": " + kvp.Value)) + "]";
753    }
754  }
755
756  [DebuggerDisplay("[{low.Value}..{high.Value}]")]
757  public class AlgebraicInterval : IAlgebraicType<AlgebraicInterval> {
758    [DebuggerBrowsable(DebuggerBrowsableState.Never)]
759    private MultivariateDual<AlgebraicDouble> low;
760    public MultivariateDual<AlgebraicDouble> LowerBound => low.Clone();
761
762    [DebuggerBrowsable(DebuggerBrowsableState.Never)]
763    private MultivariateDual<AlgebraicDouble> high;
764    public MultivariateDual<AlgebraicDouble> UpperBound => high.Clone();
765
766
767    public AlgebraicInterval() : this(double.NegativeInfinity, double.PositiveInfinity) { }
768
769    public AlgebraicInterval(MultivariateDual<AlgebraicDouble> low, MultivariateDual<AlgebraicDouble> high) {
770      this.low = low.Clone();
771      this.high = high.Clone();
772    }
773
774    public AlgebraicInterval(double low, double high) {
775      this.low = new MultivariateDual<AlgebraicDouble>(new AlgebraicDouble(low));
776      this.high = new MultivariateDual<AlgebraicDouble>(new AlgebraicDouble(high));
777    }
778
779    [DebuggerBrowsable(DebuggerBrowsableState.Never)]
780    public AlgebraicInterval Zero => new AlgebraicInterval(0.0, 0.0);
781    public AlgebraicInterval One => new AlgebraicInterval(1.0, 1.0);
782    public AlgebraicInterval Add(AlgebraicInterval a) {
783      low.Add(a.low);
784      high.Add(a.high);
785      return this;
786    }
787
788    public AlgebraicInterval Mul(AlgebraicInterval a) {
789      var v1 = low.Clone().Mul(a.low);
790      var v2 = low.Clone().Mul(a.high);
791      var v3 = high.Clone().Mul(a.low);
792      var v4 = high.Clone().Mul(a.high);
793
794      low = Algebraic.Min(Algebraic.Min(v1, v2), Algebraic.Min(v3, v4));
795      high = Algebraic.Max(Algebraic.Max(v1, v2), Algebraic.Max(v3, v4));
796      return this;
797    }
798
799    public AlgebraicInterval Assign(AlgebraicInterval a) {
800      low = a.low;
801      high = a.high;
802      return this;
803    }
804
805    public AlgebraicInterval AssignCos(AlgebraicInterval a) {
806      return AssignSin(a.Clone().Sub(new AlgebraicInterval(Math.PI / 2, Math.PI / 2)));
807    }
808
809    public AlgebraicInterval Div(AlgebraicInterval a) {
810      if (a.Contains(0.0)) {
811        if (a.low.Value.Value.IsAlmost(0.0) && a.high.Value.Value.IsAlmost(0.0)) {
812          low = new MultivariateDual<AlgebraicDouble>(double.NegativeInfinity);
813          high = new MultivariateDual<AlgebraicDouble>(double.PositiveInfinity);
814        } else if (a.low.Value.Value.IsAlmost(0.0))
815          Mul(new AlgebraicInterval(a.Clone().high.Inv(), new MultivariateDual<AlgebraicDouble>(double.PositiveInfinity)));
816        else
817          Mul(new AlgebraicInterval(new MultivariateDual<AlgebraicDouble>(double.NegativeInfinity), a.low.Clone().Inv()));
818      } else {
819        Mul(new AlgebraicInterval(a.high.Clone().Inv(), a.low.Clone().Inv())); // inverting leads to inverse roles of high and low
820      }
821      return this;
822    }
823
824    public AlgebraicInterval AssignExp(AlgebraicInterval a) {
825      low.AssignExp(a.low);
826      high.AssignExp(a.high);
827      return this;
828    }
829
830    public AlgebraicInterval AssignIntPower(AlgebraicInterval a, int p) {
831      if (p < 0) {  // x^-3 == 1/(x^3)
832        AssignIntPower(a, -p);
833        return AssignInv(this);
834      } else if (p == 0) {
835        if (a.Contains(0.0)) {
836          // => 0^0 = 0 ; might not be relevant
837          low = new MultivariateDual<AlgebraicDouble>(0.0);
838          high = new MultivariateDual<AlgebraicDouble>(1.0);
839          return this;
840        } else {
841          // => 1
842          low = new MultivariateDual<AlgebraicDouble>(1.0);
843          high = new MultivariateDual<AlgebraicDouble>(1.0);
844          return this;
845        }
846      } else if (p == 1) return this;
847      else if (p % 2 == 0) {
848        // p is even => interval must be positive
849        if (a.Contains(0.0)) {
850          low = new MultivariateDual<AlgebraicDouble>(0.0);
851          high = Algebraic.Max(a.low.IntPower(p), a.high.IntPower(p));
852        } else {
853          var lowPower = a.low.IntPower(p);
854          var highPower = a.high.IntPower(p);
855          low = Algebraic.Min(lowPower, highPower);
856          high = Algebraic.Max(lowPower, highPower);
857        }
858      } else {
859        // p is uneven
860        if (a.Contains(0.0)) {
861          low.AssignIntPower(a.low, p);
862          high.AssignIntPower(a.high, p);
863        } else {
864          var lowPower = a.low.IntPower(p);
865          var highPower = a.high.IntPower(p);
866          low = Algebraic.Min(lowPower, highPower);
867          high = Algebraic.Max(lowPower, highPower);
868        }
869      }
870      return this;
871    }
872
873    public AlgebraicInterval AssignIntRoot(AlgebraicInterval a, int r) {
874      if (r == 0) { low = new MultivariateDual<AlgebraicDouble>(double.NaN); high = new MultivariateDual<AlgebraicDouble>(double.NaN); return this; }
875      if (r == 1) return this;
876      if (r < 0) {
877        // x^ (-1/2) = 1 / (x^(1/2))
878        AssignIntRoot(a, -r);
879        return AssignInv(this);
880      } else {
881        // root only defined for positive arguments
882        if (a.LowerBound.Value.Value < 0) {
883          low = new MultivariateDual<AlgebraicDouble>(double.NaN);
884          high = new MultivariateDual<AlgebraicDouble>(double.NaN);
885          return this;
886        } else {
887          low.AssignIntRoot(a.low, r);
888          high.AssignIntRoot(a.high, r);
889          return this;
890        }
891      }
892    }
893
894    public AlgebraicInterval AssignInv(AlgebraicInterval a) {
895      low = new MultivariateDual<AlgebraicDouble>(1.0);
896      high = new MultivariateDual<AlgebraicDouble>(1.0);
897      return Div(a);
898    }
899
900    public AlgebraicInterval AssignLog(AlgebraicInterval a) {
901      low.AssignLog(a.low);
902      high.AssignLog(a.high);
903      return this;
904    }
905
906    public AlgebraicInterval AssignNeg(AlgebraicInterval a) {
907      low.AssignNeg(a.high);
908      high.AssignNeg(a.low);
909      return this;
910    }
911
912    public AlgebraicInterval Scale(double s) {
913      low.Scale(s);
914      high.Scale(s);
915      if (s < 0) {
916        var t = low;
917        low = high;
918        high = t;
919      }
920      return this;
921    }
922
923    public AlgebraicInterval AssignSin(AlgebraicInterval a) {
924      if (Math.Abs(a.UpperBound.Value.Value - a.LowerBound.Value.Value) >= Math.PI * 2) {
925        low = new MultivariateDual<AlgebraicDouble>(-1.0);
926        high = new MultivariateDual<AlgebraicDouble>(1.0);
927      }
928
929      //divide the interval by PI/2 so that the optima lie at x element of N (0,1,2,3,4,...)
930      double Pihalf = Math.PI / 2;
931      var scaled = this.Clone().Scale(1.0 / Pihalf);
932      //move to positive scale
933      if (scaled.LowerBound.Value.Value < 0) {
934        int periodsToMove = Math.Abs((int)scaled.LowerBound.Value.Value / 4) + 1;
935        scaled.Add(new AlgebraicInterval(periodsToMove * 4, periodsToMove * 4));
936      }
937
938      double scaledLowerBound = scaled.LowerBound.Value.Value % 4.0;
939      double scaledUpperBound = scaled.UpperBound.Value.Value % 4.0;
940      if (scaledUpperBound < scaledLowerBound) scaledUpperBound += 4.0;
941      List<double> sinValues = new List<double>();
942      sinValues.Add(Math.Sin(scaledLowerBound * Pihalf));
943      sinValues.Add(Math.Sin(scaledUpperBound * Pihalf));
944
945      int startValue = (int)Math.Ceiling(scaledLowerBound);
946      while (startValue < scaledUpperBound) {
947        sinValues.Add(Math.Sin(startValue * Pihalf));
948        startValue += 1;
949      }
950
951      low = new MultivariateDual<AlgebraicDouble>(sinValues.Min());
952      high = new MultivariateDual<AlgebraicDouble>(sinValues.Max());
953      return this;
954    }
955
956    public AlgebraicInterval Sub(AlgebraicInterval a) {
957      // [x1,x2] − [y1,y2] = [x1 − y2,x2 − y1]
958      low.Sub(a.high);
959      high.Sub(a.low);
960      return this;
961    }
962
963    public AlgebraicInterval Clone() {
964      return new AlgebraicInterval(low, high);
965    }
966
967    public bool Contains(double val) {
968      return LowerBound.Value.Value <= val && val <= UpperBound.Value.Value;
969    }
970
971    public AlgebraicInterval AssignAbs(AlgebraicInterval a) {
972      if (a.Contains(0.0)) {
973        var abslow = a.low.Clone().Abs();
974        var abshigh = a.high.Clone().Abs();
975        a.high.Assign(Algebraic.Max(abslow, abshigh));
976        a.low.Assign(new MultivariateDual<AlgebraicDouble>(0.0)); // lost gradient for lower bound
977      } else {
978        var abslow = a.low.Clone().Abs();
979        var abshigh = a.high.Clone().Abs();
980        a.low.Assign(Algebraic.Min(abslow, abshigh));
981        a.high.Assign(Algebraic.Max(abslow, abshigh));
982      }
983      return this;
984    }
985
986    public AlgebraicInterval AssignSgn(AlgebraicInterval a) {
987      low.AssignSgn(a.low);
988      high.AssignSgn(a.high);
989      return this;
990    }
991  }
992
993  public class Dual<V> : IAlgebraicType<Dual<V>>
994    where V : IAlgebraicType<V> {
995    [DebuggerBrowsable(DebuggerBrowsableState.Never)]
996    private V v;
997    public V Value => v;
998
999    [DebuggerBrowsable(DebuggerBrowsableState.Never)]
1000    private V dv;
1001    public V Derivative => dv;
1002
1003    public Dual(V v, V dv) { this.v = v; this.dv = dv; }
1004
1005    [DebuggerBrowsable(DebuggerBrowsableState.Never)]
1006    public Dual<V> Zero => new Dual<V>(Value.Zero, Derivative.Zero);
1007    [DebuggerBrowsable(DebuggerBrowsableState.Never)]
1008    public Dual<V> One => new Dual<V>(Value.One, Derivative.Zero);
1009
1010    public Dual<V> Assign(Dual<V> a) { v.Assign(a.v); dv.Assign(a.dv); return this; }
1011    public Dual<V> Scale(double s) { v.Scale(s); dv.Scale(s); return this; }
1012    public Dual<V> Add(Dual<V> a) { v.Add(a.v); dv.Add(a.dv); return this; }
1013    public Dual<V> Sub(Dual<V> a) { v.Sub(a.v); dv.Sub(a.dv); return this; }
1014    public Dual<V> AssignNeg(Dual<V> a) { v.AssignNeg(a.v); dv.AssignNeg(a.dv); return this; }
1015    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
1016
1017    // (a(x) * b(x))' = b(x)*a(x)' + b(x)'*a(x);
1018    public Dual<V> Mul(Dual<V> a) {
1019      var t1 = a.dv.Clone().Mul(v);
1020      var t2 = dv.Clone().Mul(a.v);
1021      dv.Assign(t1).Add(t2);
1022
1023      v.Mul(a.v);
1024      return this;
1025    }
1026    public Dual<V> Div(Dual<V> a) { Mul(a.Inv()); return this; }
1027
1028    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)'
1029    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)'
1030
1031    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; }
1032    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; }
1033
1034    public Dual<V> AssignSin(Dual<V> a) { v.AssignSin(a.v); dv.Assign(a.dv).Mul(a.v.Clone().Cos()); return this; }
1035    public Dual<V> AssignCos(Dual<V> a) { v.AssignCos(a.v); dv.AssignNeg(a.dv).Mul(a.v.Clone().Sin()); return this; }
1036
1037    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)|     
1038    public Dual<V> AssignSgn(Dual<V> a) { v.AssignSgn(a.v); dv.Assign(a.dv.Zero); return this; }
1039
1040    public Dual<V> Clone() { return new Dual<V>(v.Clone(), dv.Clone()); }
1041
1042  }
1043
1044  /// <summary>
1045  /// An algebraic type which has a value as well as the partial derivatives of the value over multiple variables.
1046  /// </summary>
1047  /// <typeparam name="V"></typeparam>
1048  [DebuggerDisplay("v={Value}; dv={dv}")]
1049  public class MultivariateDual<V> : IAlgebraicType<MultivariateDual<V>> where V : IAlgebraicType<V>, new() {
1050    [DebuggerBrowsable(DebuggerBrowsableState.Never)]
1051    private V v;
1052    public V Value => v;
1053
1054    [DebuggerBrowsable(DebuggerBrowsableState.Never)]
1055    private AlgebraicSparseVector<object, V> dv;
1056    public AlgebraicSparseVector<object, V> Gradient => dv; // <key,value> partial derivative identified via the key
1057
1058    private MultivariateDual(MultivariateDual<V> orig) { this.v = orig.v.Clone(); this.dv = orig.dv.Clone(); }
1059
1060    /// <summary>
1061    /// Constructor without partial derivative
1062    /// </summary>
1063    /// <param name="v"></param>
1064    public MultivariateDual(V v) { this.v = v.Clone(); this.dv = new AlgebraicSparseVector<object, V>(); }
1065
1066    /// <summary>
1067    /// Constructor for multiple partial derivatives
1068    /// </summary>
1069    /// <param name="v"></param>
1070    /// <param name="keys"></param>
1071    /// <param name="dv"></param>
1072    public MultivariateDual(V v, object[] keys, V[] dv) { this.v = v.Clone(); this.dv = new AlgebraicSparseVector<object, V>(keys, dv); }
1073
1074    /// <summary>
1075    /// Constructor for a single partial derivative
1076    /// </summary>
1077    /// <param name="v"></param>
1078    /// <param name="key"></param>
1079    /// <param name="dv"></param>
1080    public MultivariateDual(V v, object key, V dv) { this.v = v.Clone(); this.dv = new AlgebraicSparseVector<object, V>(new[] { key }, new[] { dv }); }
1081
1082    /// <summary>
1083    /// Constructor with a given value and gradient. For internal use.
1084    /// </summary>
1085    /// <param name="v">The value (not cloned).</param>
1086    /// <param name="gradient">The gradient (not cloned).</param>
1087    internal MultivariateDual(V v, AlgebraicSparseVector<object, V> gradient) { this.v = v; this.dv = gradient; }
1088
1089    public MultivariateDual<V> Clone() { return new MultivariateDual<V>(this); }
1090
1091    public MultivariateDual<V> Zero => new MultivariateDual<V>(Value.Zero, Gradient.Zero);
1092    public MultivariateDual<V> One => new MultivariateDual<V>(Value.One, Gradient.Zero);
1093
1094    public MultivariateDual<V> Scale(double s) { v.Scale(s); dv.Scale(s); return this; }
1095
1096    public MultivariateDual<V> Add(MultivariateDual<V> a) { v.Add(a.v); dv.Add(a.dv); return this; }
1097    public MultivariateDual<V> Sub(MultivariateDual<V> a) { v.Sub(a.v); dv.Sub(a.dv); return this; }
1098    public MultivariateDual<V> Assign(MultivariateDual<V> a) { v.Assign(a.v); dv.Assign(a.dv); return this; }
1099    public MultivariateDual<V> Mul(MultivariateDual<V> a) {
1100      // (a(x) * b(x))' = b(x)*a(x)' + b(x)'*a(x);
1101      var t1 = a.dv.Clone().Scale(v);
1102      var t2 = dv.Clone().Scale(a.v);
1103      dv.Assign(t1).Add(t2);
1104
1105      v.Mul(a.v);
1106      return this;
1107    }
1108    public MultivariateDual<V> Div(MultivariateDual<V> a) { v.Div(a.v); dv.Mul(a.dv.Inv()); return this; }
1109    public MultivariateDual<V> AssignNeg(MultivariateDual<V> a) { v.AssignNeg(a.v); dv.AssignNeg(a.dv); return this; }
1110    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
1111
1112    public MultivariateDual<V> AssignSin(MultivariateDual<V> a) { v.AssignSin(a.v); dv.Assign(a.dv).Scale(a.v.Clone().Cos()); return this; }
1113    public MultivariateDual<V> AssignCos(MultivariateDual<V> a) { v.AssignCos(a.v); dv.AssignNeg(a.dv).Scale(a.v.Clone().Sin()); return this; }
1114
1115    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; }
1116    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; }
1117
1118    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)'     
1119    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)'
1120
1121    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     
1122    public MultivariateDual<V> AssignSgn(MultivariateDual<V> a) { v.AssignSgn(a.v); dv = a.dv.Zero; return this; } // sign(f(x))' = 0;     
1123  }
1124}
Note: See TracBrowser for help on using the repository browser.