Free cookie consent management tool by TermsFeed Policy Generator

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

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

#2994: worked on implementations of remaining methods for IAlgebraicType

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