Free cookie consent management tool by TermsFeed Policy Generator

Ignore:
Timestamp:
10/02/19 17:09:31 (5 years ago)
Author:
gkronber
Message:

#2994: refactoring: moved types into separate files

Location:
branches/2994-AutoDiffForIntervals/HeuristicLab.Problems.DataAnalysis.Symbolic/3.4/Interpreter
Files:
11 added
1 edited

Legend:

Unmodified
Added
Removed

  • branches/2994-AutoDiffForIntervals/HeuristicLab.Problems.DataAnalysis.Symbolic/3.4/Interpreter/Interpreter.cs

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