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

Last change on this file since 16674 was 16674, checked in by gkronber, 4 months ago

#2994: worked on AutoDiff implementation based on BatchInterpreter

File size: 26.9 KB
Line 
1using System;
2using System.Collections.Generic;
3using System.Linq;
4using HeuristicLab.Common;
5using HeuristicLab.Encodings.SymbolicExpressionTreeEncoding;
6using HEAL.Attic;
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: 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
71          case OpCodes.Exp: {
72              instr.value.AssignExp(code[c].value);
73              break;
74            }
75
76          case OpCodes.Log: {
77              instr.value.AssignLog(code[c].value);
78              break;
79            }
80        }
81      }
82      return code[0].value;
83    }
84
85    protected Instruction[] Compile(ISymbolicExpressionTree tree) {
86      var root = tree.Root.GetSubtree(0).GetSubtree(0);
87      var code = new Instruction[root.GetLength()];
88      if (root.SubtreeCount > ushort.MaxValue) throw new ArgumentException("Number of subtrees is too big (>65.535)");
89      int c = 1, i = 0;
90      foreach (var node in root.IterateNodesBreadth()) {
91        if (node.SubtreeCount > ushort.MaxValue) throw new ArgumentException("Number of subtrees is too big (>65.535)");
92        code[i] = new Instruction {
93          opcode = OpCodes.MapSymbolToOpCode(node),
94          narg = (ushort)node.SubtreeCount,
95          childIndex = c
96        };
97        if (node is VariableTreeNode variable) {
98          InitializeTerminalInstruction(ref code[i], variable);
99        } else if (node is ConstantTreeNode constant) {
100          InitializeTerminalInstruction(ref code[i], constant);
101        } else {
102          InitializeInternalInstruction(ref code[i], node);
103        }
104        c += node.SubtreeCount;
105        ++i;
106      }
107      return code;
108    }
109
110    protected abstract void InitializeTerminalInstruction(ref Instruction instruction, ConstantTreeNode constant);
111    protected abstract void InitializeTerminalInstruction(ref Instruction instruction, VariableTreeNode variable);
112    protected abstract void InitializeInternalInstruction(ref Instruction instruction, ISymbolicExpressionTreeNode node);
113
114    protected abstract void LoadVariable(Instruction a);
115
116  }
117
118
119  public class VectorEvaluator : Interpreter<DoubleVector> {
120    private const int BATCHSIZE = 128;
121    [ThreadStatic]
122    private Dictionary<string, double[]> cachedData;
123
124    [ThreadStatic]
125    private IDataset dataset;
126
127    [ThreadStatic]
128    private int rowIndex;
129
130    [ThreadStatic]
131    private int[] rows;
132
133    private void InitCache(IDataset dataset) {
134      this.dataset = dataset;
135      cachedData = new Dictionary<string, double[]>();
136      foreach (var v in dataset.DoubleVariables) {
137        cachedData[v] = dataset.GetReadOnlyDoubleValues(v).ToArray();
138      }
139    }
140
141    public double[] Evaluate(ISymbolicExpressionTree tree, IDataset dataset, int[] rows) {
142      if (cachedData == null || this.dataset != dataset) {
143        InitCache(dataset);
144      }
145
146      this.rows = rows;
147      var code = Compile(tree);
148      var remainingRows = rows.Length % BATCHSIZE;
149      var roundedTotal = rows.Length - remainingRows;
150
151      var result = new double[rows.Length];
152
153      for (rowIndex = 0; rowIndex < roundedTotal; rowIndex += BATCHSIZE) {
154        Evaluate(code);
155        code[0].value.CopyTo(result, rowIndex, BATCHSIZE);
156      }
157
158      if (remainingRows > 0) {
159        Evaluate(code);
160        code[0].value.CopyTo(result, roundedTotal, remainingRows);
161      }
162
163      return result;
164    }
165
166    protected override void InitializeTerminalInstruction(ref Instruction instruction, ConstantTreeNode constant) {
167      instruction.dblVal = constant.Value;
168      instruction.value = new DoubleVector(BATCHSIZE);
169      instruction.value.AssignConstant(instruction.dblVal);
170    }
171
172    protected override void InitializeTerminalInstruction(ref Instruction instruction, VariableTreeNode variable) {
173      instruction.dblVal = variable.Weight;
174      instruction.value = new DoubleVector(BATCHSIZE);
175      if (cachedData.ContainsKey(variable.VariableName)) {
176        instruction.data = cachedData[variable.VariableName];
177      } else {
178        instruction.data = dataset.GetDoubleValues(variable.VariableName).ToArray();
179        cachedData[variable.VariableName] = (double[])instruction.data;
180      }
181    }
182
183    protected override void InitializeInternalInstruction(ref Instruction instruction, ISymbolicExpressionTreeNode node) {
184      instruction.value = new DoubleVector(BATCHSIZE);
185    }
186
187    protected override void LoadVariable(Instruction a) {
188      var data = (double[])a.data;
189      for (int i = rowIndex; i < rows.Length && i - rowIndex < BATCHSIZE; i++) a.value[i - rowIndex] = data[rows[i]];
190      a.value.Scale(a.dblVal);
191    }
192  }
193
194  public class VectorAutoDiffEvaluator : Interpreter<VectorDual<DoubleVector>> {
195    private const int BATCHSIZE = 128;
196    [ThreadStatic]
197    private Dictionary<string, double[]> cachedData;
198
199    [ThreadStatic]
200    private IDataset dataset;
201
202    [ThreadStatic]
203    private int rowIndex;
204
205    [ThreadStatic]
206    private int[] rows;
207
208    [ThreadStatic]
209    private Dictionary<ISymbolicExpressionTreeNode, int> node2paramIdx;
210
211    private void InitCache(IDataset dataset) {
212      this.dataset = dataset;
213      cachedData = new Dictionary<string, double[]>();
214      foreach (var v in dataset.DoubleVariables) {
215        cachedData[v] = dataset.GetDoubleValues(v).ToArray();
216      }
217    }
218
219    public void Evaluate(ISymbolicExpressionTree tree, IDataset dataset, int[] rows, ISymbolicExpressionTreeNode[] parameterNodes, out double[] fi, out double[,] jac) {
220      if (cachedData == null || this.dataset != dataset) {
221        InitCache(dataset);
222      }
223
224      int nParams = parameterNodes.Length;
225      node2paramIdx = new Dictionary<ISymbolicExpressionTreeNode, int>();
226      for (int i = 0; i < parameterNodes.Length; i++) node2paramIdx.Add(parameterNodes[i], i);
227
228      var code = Compile(tree);
229
230      var remainingRows = rows.Length % BATCHSIZE;
231      var roundedTotal = rows.Length - remainingRows;
232
233      fi = new double[rows.Length];
234      jac = new double[rows.Length, nParams];
235
236      this.rows = rows;
237
238      for (rowIndex = 0; rowIndex < roundedTotal; rowIndex += BATCHSIZE) {
239        Evaluate(code);
240        code[0].value.Value.CopyTo(fi, rowIndex, BATCHSIZE);
241
242        // TRANSPOSE into JAC
243        // here we assume that we usually calculate the gradient over _all_ parameters (otherwise it would be better to propagate only relevant vectors in evaluation)
244        for (int j = 0; j < nParams; ++j) 
245          code[0].value.Gradient[j].CopyColumnTo(jac, j, rowIndex, BATCHSIZE); XXX CONTINUE HERE
246      }
247
248      if (remainingRows > 0) {
249        Evaluate(code);
250        code[0].value.Value.CopyTo(fi, roundedTotal, remainingRows);
251        for (int j = 0; j < nParams; ++j)
252          code[0].value.Gradient[j].CopyColumnTo(jac, j, roundedTotal, remainingRows);
253      }
254    }
255
256    protected override void InitializeInternalInstruction(ref Instruction instruction, ISymbolicExpressionTreeNode node) {
257      var zero = new DoubleVector(BATCHSIZE);
258      //       var g = new Vector<DoubleVector>(node2paramIdx.Count);
259      //       for (int i = 0; i < BATCHSIZE; ++i) g[i] = new DoubleVector(BATCHSIZE);
260      instruction.value = new VectorDual<DoubleVector>(zero, null);
261    }
262
263    protected override void InitializeTerminalInstruction(ref Instruction instruction, ConstantTreeNode constant) {
264      instruction.dblVal = constant.Value;
265
266      var g_arr = new double[BATCHSIZE];
267      if (node2paramIdx.ContainsKey(constant)) {
268        for (int i = 0; i < BATCHSIZE; i++) g_arr[i] = 1.0;
269      }
270      var g = new DoubleVector(g_arr);
271
272      instruction.value = new VectorDual<DoubleVector>(new DoubleVector(BATCHSIZE), new Vector<DoubleVector>(new[] { g })); // only a single column for the gradient
273      instruction.value.Value.AssignConstant(instruction.dblVal);
274    }
275
276    protected override void InitializeTerminalInstruction(ref Instruction instruction, VariableTreeNode variable) {
277      double[] data;
278      if (cachedData.ContainsKey(variable.VariableName)) {
279        data = cachedData[variable.VariableName];
280      } else {
281        data = dataset.GetReadOnlyDoubleValues(variable.VariableName).ToArray();
282        cachedData[variable.VariableName] = (double[])instruction.data;
283      }
284
285      var paramIdx = -1;
286      if (node2paramIdx.ContainsKey(variable)) {
287        paramIdx = node2paramIdx[variable];
288      }
289
290      instruction.dblVal = variable.Weight;
291      instruction.data = new object[] { data, paramIdx };
292
293      var f = new DoubleVector(BATCHSIZE);
294      var g = new DoubleVector(BATCHSIZE);
295      instruction.value = new VectorDual<DoubleVector>(f, new Vector<DoubleVector>(new[] { g }));
296    }
297
298    protected override void LoadVariable(Instruction a) {
299      var paramIdx = (int)((object[])a.data)[1];
300      var data = (double[])((object[])a.data)[0];
301
302      for (int i = rowIndex; i < rows.Length && i - rowIndex < BATCHSIZE; i++) a.value.Value[i - rowIndex] = data[rows[i]];
303      a.value.Scale(a.dblVal);
304
305      if (paramIdx >= 0) {
306        // update gradient with variable values
307        var g = a.value.Gradient[0];
308        for (int i = rowIndex; i < rows.Length && i - rowIndex < BATCHSIZE; i++) {
309          g[i] = data[rows[i]];
310        }
311      }
312    }
313  }
314
315
316  public class IntervalEvaluator : Interpreter<AlgebraicInterval> {
317    [ThreadStatic]
318    private Dictionary<string, AlgebraicInterval> intervals;
319
320    public AlgebraicInterval Evaluate(ISymbolicExpressionTree tree, Dictionary<string, AlgebraicInterval> intervals) {
321      this.intervals = intervals;
322      var code = Compile(tree);
323      Evaluate(code);
324      return code[0].value;
325    }
326
327
328    protected override void InitializeInternalInstruction(ref Instruction instruction, ISymbolicExpressionTreeNode node) {
329      instruction.value = new AlgebraicInterval(0, 0);
330    }
331
332
333    protected override void InitializeTerminalInstruction(ref Instruction instruction, ConstantTreeNode constant) {
334      instruction.dblVal = constant.Value;
335      instruction.value = new AlgebraicInterval(constant.Value, constant.Value);
336    }
337
338    protected override void InitializeTerminalInstruction(ref Instruction instruction, VariableTreeNode variable) {
339      instruction.dblVal = variable.Weight;
340      instruction.value = intervals[variable.VariableName];
341    }
342
343    protected override void LoadVariable(Instruction a) {
344      // nothing to do
345    }
346  }
347
348  public interface IAlgebraicType<T> {
349    T Assign(T a); // assign this to same value as a (copy!)
350    T AssignNeg(T a); // set this to negative(a)
351    T AssignInv(T a); // set this to inv(a);
352    T Scale(double s); // scale this with s
353    T Add(T a); // add a to this
354    T Sub(T a); // subtract a from this
355    T Mul(T a); // multiply this with a
356    T Div(T a); // divide this by a
357    T AssignLog(T a); // set this to log a
358    T AssignExp(T a); // set this to exp(a)
359    T AssignSin(T a); // set this to sin(a)
360    T AssignCos(T a); // set this to cos(a)
361    T AssignIntPower(T a, int p);
362    T AssignIntRoot(T a, int r);
363    T Clone();
364  }
365
366
367  // a simple vector as an algebraic type
368  public class DoubleVector : IAlgebraicType<DoubleVector> {
369    private double[] arr;
370
371    public double this[int idx] { get { return arr[idx]; } set { arr[idx] = value; } }
372
373    public DoubleVector(int length) {
374      arr = new double[length];
375    }
376
377    /// <summary>
378    ///
379    /// </summary>
380    /// <param name="arr">array is not copied</param>
381    public DoubleVector(double[] arr) {
382      this.arr = arr;
383    }
384
385    public DoubleVector(int length, double constantValue) : this(length) {
386      for (int i = 0; i < length; ++i) arr[i] = constantValue;
387    }
388
389    public void CopyTo(double[] dest, int idx, int length) {
390      Array.Copy(arr, 0, dest, idx, length);
391    }
392
393    public void CopyRowTo(double[,] dest, int row) {
394      for (int j = 0; j < arr.Length; ++j) dest[row, j] = arr[j];
395    }
396    internal void CopyColumnTo(double[,] dest, int column, int row, int len) {
397      for (int j = 0; j < len; ++j) dest[row + j, column] = arr[j];
398    }
399
400
401    public DoubleVector Add(DoubleVector a) {
402      for (int i = 0; i < arr.Length; ++i) {
403        arr[i] += a.arr[i];
404      }
405      return this;
406    }
407
408    public void AssignConstant(double constantValue) {
409      for (int i = 0; i < arr.Length; ++i) {
410        arr[i] = constantValue;
411      }
412    }
413
414    public DoubleVector Assign(DoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = a.arr[i]; } return this; }
415    public DoubleVector AssignCos(DoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = Math.Cos(a.arr[i]); } return this; }
416    public DoubleVector Div(DoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] /= a.arr[i]; } return this; }
417    public DoubleVector AssignExp(DoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = Math.Exp(a.arr[i]); } return this; }
418    public DoubleVector AssignIntPower(DoubleVector a, int p) { throw new NotImplementedException(); }
419    public DoubleVector AssignIntRoot(DoubleVector a, int r) { throw new NotImplementedException(); }
420    public DoubleVector AssignInv(DoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = 1.0 / a.arr[i]; } return this; }
421    public DoubleVector AssignLog(DoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = Math.Log(a.arr[i]); } return this; }
422    public DoubleVector Mul(DoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] *= a.arr[i]; } return this; }
423    public DoubleVector AssignNeg(DoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = -a.arr[i]; } return this; }
424    public DoubleVector Scale(double s) { for (int i = 0; i < arr.Length; ++i) { arr[i] *= s; } return this; }
425    public DoubleVector AssignSin(DoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = Math.Sin(a.arr[i]); } return this; }
426    public DoubleVector Sub(DoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] -= a.arr[i]; } return this; }
427
428    public DoubleVector Clone() {
429      var v = new DoubleVector(this.arr.Length);
430      Array.Copy(arr, v.arr, v.arr.Length);
431      return v;
432    }
433
434    public void CopyFrom(double[] data, int rowIndex) {
435      Array.Copy(data, rowIndex, arr, 0, Math.Min(arr.Length, data.Length - rowIndex));
436    }
437
438  }
439
440
441  // vectors of algebraic types
442  public class Vector<T> : IAlgebraicType<Vector<T>> where T : IAlgebraicType<T> {
443    private T[] elems;
444
445    public T this[int idx] { get { return elems[idx]; } set { elems[idx] = value; } }
446
447    public int Length => elems.Length;
448
449    private Vector() { }
450
451    public Vector(int len) {
452      elems = new T[len];
453    }
454
455    /// <summary>
456    ///
457    /// </summary>
458    /// <param name="elems">The array is copied</param>
459    public Vector(T[] elems) {
460      this.elems = new T[elems.Length];
461      for (int i = 0; i < elems.Length; ++i) { this.elems[i] = elems[i].Clone(); }
462    }
463
464    /// <summary>
465    ///
466    /// </summary>
467    /// <param name="elems">Array is not copied!</param>
468    /// <returns></returns>
469    public Vector<T> FromArray(T[] elems) {
470      var v = new Vector<T>();
471      v.elems = elems;
472      return v;
473    }
474
475    public void CopyTo(T[] dest) {
476      if (dest.Length != elems.Length) throw new InvalidOperationException("arr lengths do not match in Vector<T>.Copy");
477      Array.Copy(elems, dest, dest.Length);
478    }
479
480    public Vector<T> Clone() {
481      return new Vector<T>(elems);
482    }
483
484    public Vector<T> Concat(Vector<T> other) {
485      var oldLen = Length;
486      Array.Resize(ref this.elems, oldLen + other.Length);
487      for (int i = oldLen; i < Length; i++) {
488        elems[i] = other.elems[i - oldLen].Clone();
489      }
490      return this;
491    }
492
493    public Vector<T> Add(Vector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].Add(a.elems[i]); } return this; }
494    public Vector<T> Assign(Vector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].Assign(a.elems[i]); } return this; }
495    public Vector<T> AssignCos(Vector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].AssignCos(a.elems[i]); } return this; }
496    public Vector<T> AssignExp(Vector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].AssignExp(a.elems[i]); } return this; }
497    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; }
498    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; }
499    public Vector<T> AssignInv(Vector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].AssignInv(a.elems[i]); } return this; }
500    public Vector<T> AssignLog(Vector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].AssignLog(a.elems[i]); } return this; }
501    public Vector<T> AssignNeg(Vector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].AssignNeg(a.elems[i]); } return this; }
502    public Vector<T> AssignSin(Vector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].AssignSin(a.elems[i]); } return this; }
503    public Vector<T> Div(Vector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].Div(a.elems[i]); } return this; }
504    public Vector<T> Mul(Vector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].Mul(a.elems[i]); } return this; }
505    public Vector<T> Scale(double s) { for (int i = 0; i < elems.Length; ++i) { elems[i].Scale(s); } return this; }
506    public Vector<T> Scale(T s) { for (int i = 0; i < elems.Length; ++i) { elems[i].Mul(s); } return this; }
507    public Vector<T> Sub(Vector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].Sub(a.elems[i]); } return this; }
508  }
509
510
511  public class AlgebraicInterval : IAlgebraicType<AlgebraicInterval> {
512    private double low;
513    private double high;
514
515    public double LowerBound => low;
516    public double UpperBound => high;
517
518    public AlgebraicInterval() : this(double.NegativeInfinity, double.PositiveInfinity) { }
519
520    public AlgebraicInterval(double low, double high) {
521      this.low = low;
522      this.high = high;
523    }
524
525    public AlgebraicInterval Add(AlgebraicInterval a) {
526      low += a.low;
527      high += a.high;
528      return this;
529    }
530
531    public AlgebraicInterval Assign(AlgebraicInterval a) {
532      low = a.low;
533      high = a.high;
534      return this;
535    }
536
537    public AlgebraicInterval AssignCos(AlgebraicInterval a) {
538      throw new NotImplementedException();
539    }
540
541    public AlgebraicInterval Div(AlgebraicInterval a) {
542      if (Contains(0.0)) {
543        if (a.low.IsAlmost(0.0))
544          Mul(new AlgebraicInterval(1.0 / a.high, double.PositiveInfinity));
545        else if (a.high.IsAlmost(0.0))
546          Mul(new AlgebraicInterval(double.NegativeInfinity, 1.0 / a.low));
547        else {
548          low = double.NegativeInfinity;
549          high = double.PositiveInfinity;
550        }
551      } else {
552        Mul(new AlgebraicInterval(1.0 / a.high, 1.0 / a.low));
553      }
554      return this;
555    }
556
557    public AlgebraicInterval AssignExp(AlgebraicInterval a) {
558      low = Math.Exp(a.low);
559      high = Math.Exp(a.high);
560      return this;
561    }
562
563    public AlgebraicInterval AssignIntPower(AlgebraicInterval a, int p) {
564      throw new NotImplementedException();
565    }
566
567    public AlgebraicInterval AssignIntRoot(AlgebraicInterval a, int r) {
568      throw new NotImplementedException();
569    }
570
571    public AlgebraicInterval AssignInv(AlgebraicInterval a) {
572      low = 1.0;
573      high = 1.0;
574      return Div(a);
575    }
576
577    public AlgebraicInterval AssignLog(AlgebraicInterval a) {
578      low = Math.Log(a.low);
579      high = Math.Log(a.high);
580      return this;
581    }
582
583    public AlgebraicInterval Mul(AlgebraicInterval a) {
584      double v1 = low * a.low;
585      double v2 = low * a.high;
586      double v3 = high * a.low;
587      double v4 = high * a.high;
588
589      low = Math.Min(Math.Min(v1, v2), Math.Min(v3, v4));
590      high = Math.Max(Math.Max(v1, v2), Math.Max(v3, v4));
591      return this;
592    }
593
594    public AlgebraicInterval AssignNeg(AlgebraicInterval a) {
595      low = 0;
596      high = 0;
597      return Sub(a);
598    }
599
600    public AlgebraicInterval Scale(double s) {
601      if (s < 0) {
602        low = s * high;
603        high = s * low;
604      } else {
605        low = s * low;
606        high = s * high;
607      }
608      return this;
609    }
610
611    public AlgebraicInterval AssignSin(AlgebraicInterval a) {
612      throw new NotImplementedException();
613    }
614
615    public AlgebraicInterval Sub(AlgebraicInterval a) {
616      low -= a.high;
617      high -= a.low;
618      return this;
619    }
620
621    public AlgebraicInterval Clone() {
622      return new AlgebraicInterval(low, high);
623    }
624
625    public bool Contains(double val) {
626      return low <= val && val <= high;
627    }
628  }
629
630  public class Dual<V> : IAlgebraicType<Dual<V>>
631    where V : IAlgebraicType<V> {
632    private V v;
633    private V dv;
634
635    public Dual(V v, V dv) {
636      this.v = v;
637      this.dv = dv;
638    }
639
640    public Dual<V> Clone() {
641      return new Dual<V>(v.Clone(), dv.Clone());
642    }
643
644    public Dual<V> Add(Dual<V> a) {
645      v.Add(a.v);
646      dv.Add(a.dv);
647      return this;
648    }
649
650    public Dual<V> Assign(Dual<V> a) {
651      v.Assign(a.v);
652      dv.Assign(a.dv);
653      return this;
654    }
655
656    public Dual<V> AssignCos(Dual<V> a) {
657      v.AssignCos(a.v);
658      dv.AssignNeg(dv.AssignSin(a.v));
659      return this;
660    }
661
662    public Dual<V> Div(Dual<V> a) {
663      throw new NotImplementedException();
664    }
665
666    public Dual<V> AssignExp(Dual<V> a) {
667      v.AssignExp(a.v);
668      dv.Assign(a.dv).Mul(v); // exp(f(x)) = exp(f(x))*f(x)'
669      return this;
670    }
671
672    public Dual<V> AssignIntPower(Dual<V> a, int p) {
673      throw new NotImplementedException();
674    }
675
676    public Dual<V> AssignIntRoot(Dual<V> a, int r) {
677      throw new NotImplementedException();
678    }
679
680    public Dual<V> AssignInv(Dual<V> a) {
681      throw new NotImplementedException();
682    }
683
684    public Dual<V> AssignLog(Dual<V> a) {
685      v.AssignLog(a.v);
686      dv.Assign(a.dv).Div(a.v); // log(x)' = 1/f(x) * f(x)'
687      return this;
688    }
689
690    public Dual<V> Mul(Dual<V> a) {
691      // (a(x) * b(x))' = b(x)*a(x)' + b(x)'*a(x);
692
693      V t1 = default(V);
694      t1.Assign(a.dv).Mul(v);
695
696      var t2 = default(V);
697      t2.Assign(dv).Mul(a.v);
698
699      dv.Assign(t1).Add(t2);
700
701      v.Mul(a.v);
702      return this;
703    }
704
705    public Dual<V> AssignNeg(Dual<V> a) {
706      throw new NotImplementedException();
707    }
708
709    public Dual<V> Scale(double s) {
710      v.Scale(s);
711      dv.Scale(s);
712      return this;
713    }
714
715    public Dual<V> AssignSin(Dual<V> a) {
716      throw new NotImplementedException();
717    }
718
719    public Dual<V> Sub(Dual<V> a) {
720      throw new NotImplementedException();
721    }
722  }
723
724  public class VectorDual<V> : IAlgebraicType<VectorDual<V>> where V : IAlgebraicType<V> {
725    private V v;
726    public V Value => v;
727
728    private Vector<V> dv;
729    public Vector<V> Gradient => dv; // column-oriented (dv[0] is the vector for the first parameter)
730
731
732    public VectorDual(V v, Vector<V> dv) {
733      this.v = v;
734      this.dv = dv;
735    }
736
737    public VectorDual<V> Clone() {
738      return new VectorDual<V>(v.Clone(), dv.Clone());
739    }
740
741    public VectorDual<V> Add(VectorDual<V> a) {
742      v.Add(a.v);
743      dv.Concat(a.dv);
744      return this;
745    }
746
747    public VectorDual<V> Assign(VectorDual<V> a) {
748      v.Assign(a.v);
749      dv = a.dv.Clone();
750      return this;
751    }
752
753    public VectorDual<V> AssignCos(VectorDual<V> a) {
754      v.AssignCos(a.v);
755      V t = default(V);
756      t.AssignNeg(t.AssignSin(a.v));
757      dv = a.dv.Clone().Scale(t);
758      return this;
759    }
760
761    public VectorDual<V> AssignExp(VectorDual<V> a) {
762      v.AssignExp(a.v);
763      dv = a.dv.Clone().Scale(v);
764      return this;
765    }
766
767    public VectorDual<V> AssignIntPower(VectorDual<V> a, int p) {
768      throw new NotImplementedException();
769    }
770
771    public VectorDual<V> AssignIntRoot(VectorDual<V> a, int r) {
772      throw new NotImplementedException();
773    }
774
775    public VectorDual<V> AssignInv(VectorDual<V> a) {
776      throw new NotImplementedException();
777    }
778
779    public VectorDual<V> AssignLog(VectorDual<V> a) {
780      throw new NotImplementedException();
781    }
782
783    public VectorDual<V> AssignNeg(VectorDual<V> a) {
784      throw new NotImplementedException();
785    }
786
787    public VectorDual<V> AssignSin(VectorDual<V> a) {
788      throw new NotImplementedException();
789    }
790
791    public VectorDual<V> Div(VectorDual<V> a) {
792      throw new NotImplementedException();
793    }
794
795    public VectorDual<V> Mul(VectorDual<V> a) {
796      // (a(x) * b(x))' = b(x)*a(x)' + b(x)'*a(x);
797
798      var t1 = a.dv.Clone().Scale(v);
799      var t2 = dv.Clone().Scale(a.v);
800      dv.Assign(t1).Concat(t2); // TODO: CHECK ORDER
801
802      v.Mul(a.v);
803      return this;
804    }
805
806    public VectorDual<V> Scale(double s) {
807      v.Scale(s);
808      dv.Scale(s);
809      return this;
810    }
811
812    public VectorDual<V> Sub(VectorDual<V> a) {
813      v.Sub(a.v);
814      dv.Concat(a.dv.AssignNeg(a.dv));
815      return this;
816    }
817  }
818}
Note: See TracBrowser for help on using the repository browser.