Free cookie consent management tool by TermsFeed Policy Generator

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

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

#2994: renamed classes Algebraic..., added DebuggerDisplay attributes to simplify debugging. Fixed bugs in static methods for Algebraic

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