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

Last change on this file since 16682 was 16682, checked in by gkronber, 3 months ago

#2994: worked on auto diff for intervals and vectors

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