Free cookie consent management tool by TermsFeed Policy Generator

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

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

#2994: fixed a bug in the interval calculation for variables

File size: 54.2 KB
Line
1using System;
2using System.Collections.Generic;
3using System.Diagnostics;
4using System.Linq;
5using HeuristicLab.Common;
6using HeuristicLab.Encodings.SymbolicExpressionTreeEncoding;
7
8namespace HeuristicLab.Problems.DataAnalysis.Symbolic {
9  public abstract class Interpreter<T> where T : IAlgebraicType<T> {
10    public struct Instruction {
11      public byte opcode;
12      public ushort narg;
13      public int childIndex;
14      public double dblVal;
15      public object data; // any kind of data you want to store in instructions
16      public T value;
17    }
18
19    public T Evaluate(Instruction[] code) {
20      for (int i = code.Length - 1; i >= 0; --i) {
21        var instr = code[i];
22        var c = instr.childIndex;
23        var n = instr.narg;
24
25        switch (instr.opcode) {
26          case OpCodes.Variable: {
28              break;
29            }
30          case OpCodes.Constant: { break; }  // we initialize constants in Compile. The value never changes afterwards
32              instr.value.Assign(code[c].value);
33              for (int j = 1; j < n; ++j) {
35              }
36              break;
37            }
38
39          case OpCodes.Sub: {
40              if (n == 1) {
41                instr.value.AssignNeg(code[c].value);
42              } else {
43                instr.value.Assign(code[c].value);
44                for (int j = 1; j < n; ++j) {
45                  instr.value.Sub(code[c + j].value);
46                }
47              }
48              break;
49            }
50
51          case OpCodes.Mul: {
52              instr.value.Assign(code[c].value);
53              for (int j = 1; j < n; ++j) {
54                instr.value.Mul(code[c + j].value);
55              }
56              break;
57            }
58
59          case OpCodes.Div: {
60              if (n == 1) {
61                instr.value.AssignInv(code[c].value);
62              } else {
63                instr.value.Assign(code[c].value);
64                for (int j = 1; j < n; ++j) {
65                  instr.value.Div(code[c + j].value);
66                }
67              }
68              break;
69            }
70          case OpCodes.Square: {
71              instr.value.AssignIntPower(code[c].value, 2);
72              break;
73            }
74          case OpCodes.SquareRoot: {
75              instr.value.AssignIntRoot(code[c].value, 2);
76              break;
77            }
78          case OpCodes.Cube: {
79              instr.value.AssignIntPower(code[c].value, 3);
80              break;
81            }
82          case OpCodes.CubeRoot: {
83              instr.value.AssignIntRoot(code[c].value, 3);
84              break;
85            }
86          case OpCodes.Exp: {
87              instr.value.AssignExp(code[c].value);
88              break;
89            }
90          case OpCodes.Log: {
91              instr.value.AssignLog(code[c].value);
92              break;
93            }
94          case OpCodes.Sin: {
95              instr.value.AssignSin(code[c].value);
96              break;
97            }
98          case OpCodes.Cos: {
99              instr.value.AssignCos(code[c].value);
100              break;
101            }
102          case OpCodes.Absolute: {
103              instr.value.AssignAbs(code[c].value);
104              break;
105            }
106          case OpCodes.AnalyticQuotient: {
107              instr.value.Assign(code[c].value);
108              for (int j = 1; j < n; ++j) {
109                var t = instr.value.One;
111                instr.value.Div(t.IntRoot(2));
112              }
113              break;
114            }
115          case OpCodes.Tanh: {
116              instr.value.AssignTanh(code[c].value);
117              break;
118            }
119
120          default: throw new ArgumentException($"Unknown opcode {instr.opcode}"); 121 } 122 } 123 return code[0].value; 124 } 125 126 protected Instruction[] Compile(ISymbolicExpressionTree tree) { 127 var root = tree.Root.GetSubtree(0).GetSubtree(0); 128 var code = new Instruction[root.GetLength()]; 129 if (root.SubtreeCount > ushort.MaxValue) throw new ArgumentException("Number of subtrees is too big (>65.535)"); 130 int c = 1, i = 0; 131 foreach (var node in root.IterateNodesBreadth()) { 132 if (node.SubtreeCount > ushort.MaxValue) throw new ArgumentException("Number of subtrees is too big (>65.535)"); 133 code[i] = new Instruction { 134 opcode = OpCodes.MapSymbolToOpCode(node), 135 narg = (ushort)node.SubtreeCount, 136 childIndex = c 137 }; 138 if (node is VariableTreeNode variable) { 139 InitializeTerminalInstruction(ref code[i], variable); 140 } else if (node is ConstantTreeNode constant) { 141 InitializeTerminalInstruction(ref code[i], constant); 142 } else { 143 InitializeInternalInstruction(ref code[i], node); 144 } 145 c += node.SubtreeCount; 146 ++i; 147 } 148 return code; 149 } 150 151 protected abstract void InitializeTerminalInstruction(ref Instruction instruction, ConstantTreeNode constant); 152 protected abstract void InitializeTerminalInstruction(ref Instruction instruction, VariableTreeNode variable); 153 protected abstract void InitializeInternalInstruction(ref Instruction instruction, ISymbolicExpressionTreeNode node); 154 155 protected abstract void LoadVariable(Instruction a); 156 157 } 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);
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;
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; }
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
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 Clone();
456  }
457
458  public static class Algebraic {
459    public static T Abs<T>(this T a) where T : IAlgebraicType<T> { a.AssignAbs(a.Clone()); return a; }
460    public static T Neg<T>(this T a) where T : IAlgebraicType<T> { a.AssignNeg(a.Clone()); return a; }
461    public static T Inv<T>(this T a) where T : IAlgebraicType<T> { a.AssignInv(a.Clone()); return a; }
462    public static T Log<T>(this T a) where T : IAlgebraicType<T> { a.AssignLog(a.Clone()); return a; }
463    public static T Exp<T>(this T a) where T : IAlgebraicType<T> { a.AssignExp(a.Clone()); return a; }
464    public static T Sin<T>(this T a) where T : IAlgebraicType<T> { a.AssignSin(a.Clone()); return a; }
465    public static T Cos<T>(this T a) where T : IAlgebraicType<T> { a.AssignCos(a.Clone()); return a; }
466    public static T Sgn<T>(this T a) where T : IAlgebraicType<T> { a.AssignSgn(a.Clone()); return a; }
467    public static T IntPower<T>(this T a, int p) where T : IAlgebraicType<T> { a.AssignIntPower(a.Clone(), p); return a; }
468    public static T IntRoot<T>(this T a, int r) where T : IAlgebraicType<T> { a.AssignIntRoot(a.Clone(), r); return a; }
469
470    public static T Max<T>(T a, T b) where T : IAlgebraicType<T> {
471      // ((a + b) + abs(b - a)) / 2
473    }
474    public static T Min<T>(T a, T b) where T : IAlgebraicType<T> {
475      // ((a + b) - abs(a - b)) / 2
477    }
478  }
479
480
481  // algebraic type wrapper for a double value
482  [DebuggerDisplay("{Value}")]
483  public sealed class AlgebraicDouble : IAlgebraicType<AlgebraicDouble> {
484    public static implicit operator AlgebraicDouble(double value) { return new AlgebraicDouble(value); }
485    public static implicit operator double(AlgebraicDouble value) { return value.Value; }
486    public double Value;
487
488    [DebuggerBrowsable(DebuggerBrowsableState.Never)]
489    public AlgebraicDouble Zero => new AlgebraicDouble(0.0);
490    [DebuggerBrowsable(DebuggerBrowsableState.Never)]
491    public AlgebraicDouble One => new AlgebraicDouble(1.0);
492
493    public bool IsInfinity => IsNegativeInfinity || IsPositiveInfinity;
494    public bool IsNegativeInfinity => double.IsNegativeInfinity(Value);
495    public bool IsPositiveInfinity => double.IsPositiveInfinity(Value);
496    public AlgebraicDouble() { }
497    public AlgebraicDouble(double value) { this.Value = value; }
498    public AlgebraicDouble Assign(AlgebraicDouble a) { Value = a.Value; return this; }
499    public AlgebraicDouble Add(AlgebraicDouble a) { Value += a.Value; return this; }
500    public AlgebraicDouble Sub(AlgebraicDouble a) { Value -= a.Value; return this; }
501    public AlgebraicDouble Mul(AlgebraicDouble a) { Value *= a.Value; return this; }
502    public AlgebraicDouble Div(AlgebraicDouble a) { Value /= a.Value; return this; }
503    public AlgebraicDouble Scale(double s) { Value *= s; return this; }
504    public AlgebraicDouble AssignInv(AlgebraicDouble a) { Value = 1.0 / a.Value; return this; }
505    public AlgebraicDouble AssignNeg(AlgebraicDouble a) { Value = -a.Value; return this; }
506    public AlgebraicDouble AssignSin(AlgebraicDouble a) { Value = Math.Sin(a.Value); return this; }
507    public AlgebraicDouble AssignCos(AlgebraicDouble a) { Value = Math.Cos(a.Value); return this; }
508    public AlgebraicDouble AssignTanh(AlgebraicDouble a) { Value = Math.Tanh(a.Value); return this; }
509    public AlgebraicDouble AssignLog(AlgebraicDouble a) { Value = Math.Log(a.Value); return this; }
510    public AlgebraicDouble AssignExp(AlgebraicDouble a) { Value = Math.Exp(a.Value); return this; }
511    public AlgebraicDouble AssignIntPower(AlgebraicDouble a, int p) { Value = Math.Pow(a.Value, p); return this; }
512    public AlgebraicDouble AssignIntRoot(AlgebraicDouble a, int r) { Value = IntRoot(a.Value, r); return this; }
513
514    // helper
515    private static double IntRoot(double value, int r) {
516      if (r % 2 == 0) return Math.Pow(value, 1.0 / r);
517      else return value < 0 ? -Math.Pow(-value, 1.0 / r) : Math.Pow(value, 1.0 / r);
518    }
519
520    public AlgebraicDouble AssignAbs(AlgebraicDouble a) { Value = Math.Abs(a.Value); return this; }
521    public AlgebraicDouble AssignSgn(AlgebraicDouble a) { Value = double.IsNaN(a.Value) ? double.NaN : Math.Sign(a.Value); return this; }
522    public AlgebraicDouble Clone() { return new AlgebraicDouble(Value); }
523
524    public override string ToString() {
525      return Value.ToString();
526    }
527  }
528
529  // a simple vector as an algebraic type
530  [DebuggerDisplay("DoubleVector(len={Length}): {string.}")]
531  public class AlgebraicDoubleVector : IAlgebraicType<AlgebraicDoubleVector> {
532    private double[] arr;
533    public double this[int idx] { get { return arr[idx]; } set { arr[idx] = value; } }
534    public int Length => arr.Length;
535
536    public AlgebraicDoubleVector(int length) { arr = new double[length]; }
537
538    public AlgebraicDoubleVector() { }
539
540    /// <summary>
541    ///
542    /// </summary>
543    /// <param name="arr">array is not copied</param>
544    public AlgebraicDoubleVector(double[] arr) { this.arr = arr; }
545
546    [DebuggerBrowsable(DebuggerBrowsableState.Never)]
547    public AlgebraicDoubleVector Zero => new AlgebraicDoubleVector(new double[this.Length]); // must return vector of same length as this (therefore Zero is not static)
548    [DebuggerBrowsable(DebuggerBrowsableState.Never)]
549    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)
550    public AlgebraicDoubleVector Assign(AlgebraicDoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = a.arr[i]; } return this; }
551    public AlgebraicDoubleVector Add(AlgebraicDoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] += a.arr[i]; } return this; }
552    public AlgebraicDoubleVector Sub(AlgebraicDoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] -= a.arr[i]; } return this; }
553    public AlgebraicDoubleVector Mul(AlgebraicDoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] *= a.arr[i]; } return this; }
554    public AlgebraicDoubleVector Div(AlgebraicDoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] /= a.arr[i]; } return this; }
555    public AlgebraicDoubleVector AssignNeg(AlgebraicDoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = -a.arr[i]; } return this; }
556    public AlgebraicDoubleVector AssignInv(AlgebraicDoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = 1.0 / a.arr[i]; } return this; }
557    public AlgebraicDoubleVector Scale(double s) { for (int i = 0; i < arr.Length; ++i) { arr[i] *= s; } return this; }
558    public AlgebraicDoubleVector AssignLog(AlgebraicDoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = Math.Log(a.arr[i]); } return this; }
559    public AlgebraicDoubleVector AssignSin(AlgebraicDoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = Math.Sin(a.arr[i]); } return this; }
560    public AlgebraicDoubleVector AssignExp(AlgebraicDoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = Math.Exp(a.arr[i]); } return this; }
561    public AlgebraicDoubleVector AssignCos(AlgebraicDoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = Math.Cos(a.arr[i]); } return this; }
562    public AlgebraicDoubleVector AssignTanh(AlgebraicDoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = Math.Tanh(a.arr[i]); } return this; }
563    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; }
564    public AlgebraicDoubleVector AssignIntRoot(AlgebraicDoubleVector a, int r) { for (int i = 0; i < arr.Length; ++i) { arr[i] = IntRoot(a.arr[i], r); } return this; }
565
566    // helper
567    private double IntRoot(double v, int r) {
568      if (r % 2 == 0) return Math.Pow(v, 1.0 / r);
569      else return v < 0 ? -Math.Pow(-v, 1.0 / r) : Math.Pow(v, 1.0 / r);
570    }
571
572    public AlgebraicDoubleVector AssignAbs(AlgebraicDoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = Math.Abs(a.arr[i]); } return this; }
573    public AlgebraicDoubleVector AssignSgn(AlgebraicDoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = Math.Sign(a.arr[i]); } return this; }
574
575    public AlgebraicDoubleVector Clone() {
576      var v = new AlgebraicDoubleVector(this.arr.Length);
577      Array.Copy(arr, v.arr, v.arr.Length);
578      return v;
579    }
580
581    public AlgebraicDoubleVector AssignConstant(double constantValue) {
582      for (int i = 0; i < arr.Length; ++i) {
583        arr[i] = constantValue;
584      }
585      return this;
586    }
587
588    public void CopyTo(double[] dest, int idx, int length) {
589      Array.Copy(arr, 0, dest, idx, length);
590    }
591
592    public void CopyFrom(double[] data, int rowIndex) {
593      Array.Copy(data, rowIndex, arr, 0, Math.Min(arr.Length, data.Length - rowIndex));
594    }
595    public void CopyRowTo(double[,] dest, int row) {
596      for (int j = 0; j < arr.Length; ++j) dest[row, j] = arr[j];
597    }
598
599    internal void CopyColumnTo(double[,] dest, int column, int row, int len) {
600      for (int j = 0; j < len; ++j) dest[row + j, column] = arr[j];
601    }
602
603    public override string ToString() {
604      return "{" + string.Join(", ", arr.Take(Math.Max(5, arr.Length))) + (arr.Length > 5 ? "..." : string.Empty) + "}";
605    }
606  }
607
608
609  /*
610  // vectors of algebraic types
611  public sealed class AlgebraicVector<T> : IAlgebraicType<AlgebraicVector<T>> where T : IAlgebraicType<T>, new() {
612    private T[] elems;
613
614    public T this[int idx] { get { return elems[idx]; } set { elems[idx] = value; } }
615
616    public int Length => elems.Length;
617
618    private AlgebraicVector() { }
619
620    public AlgebraicVector(int len) { elems = new T[len]; }
621
622    /// <summary>
623    ///
624    /// </summary>
625    /// <param name="elems">The array is copied (element-wise clone)</param>
626    public AlgebraicVector(T[] elems) {
627      this.elems = new T[elems.Length];
628      for (int i = 0; i < elems.Length; ++i) { this.elems[i] = elems[i].Clone(); }
629    }
630
631    /// <summary>
632    ///
633    /// </summary>
634    /// <param name="elems">Array is not copied!</param>
635    /// <returns></returns>
636    public AlgebraicVector<T> FromArray(T[] elems) {
637      var v = new AlgebraicVector<T>();
638      v.elems = elems;
639      return v;
640    }
641
642    public void CopyTo(T[] dest) {
643      if (dest.Length != elems.Length) throw new InvalidOperationException("arr lengths do not match in Vector<T>.Copy");
644      Array.Copy(elems, dest, dest.Length);
645    }
646
647    public AlgebraicVector<T> Clone() { return new AlgebraicVector<T>(elems); }
648
649
650    [DebuggerBrowsable(DebuggerBrowsableState.Never)]
651    public AlgebraicVector<T> Zero => new AlgebraicVector<T>(Length);
652    [DebuggerBrowsable(DebuggerBrowsableState.Never)]
653    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; } }
654    public AlgebraicVector<T> Assign(AlgebraicVector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].Assign(a.elems[i]); } return this; }
655    public AlgebraicVector<T> Add(AlgebraicVector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].Add(a.elems[i]); } return this; }
656    public AlgebraicVector<T> Sub(AlgebraicVector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].Sub(a.elems[i]); } return this; }
657    public AlgebraicVector<T> Mul(AlgebraicVector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].Mul(a.elems[i]); } return this; }
658    public AlgebraicVector<T> Div(AlgebraicVector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].Div(a.elems[i]); } return this; }
659    public AlgebraicVector<T> AssignNeg(AlgebraicVector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].AssignNeg(a.elems[i]); } return this; }
660    public AlgebraicVector<T> Scale(double s) { for (int i = 0; i < elems.Length; ++i) { elems[i].Scale(s); } return this; }
661    public AlgebraicVector<T> Scale(T s) { for (int i = 0; i < elems.Length; ++i) { elems[i].Mul(s); } return this; }
662    public AlgebraicVector<T> AssignInv(AlgebraicVector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].AssignInv(a.elems[i]); } return this; }
663    public AlgebraicVector<T> AssignLog(AlgebraicVector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].AssignLog(a.elems[i]); } return this; }
664    public AlgebraicVector<T> AssignExp(AlgebraicVector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].AssignExp(a.elems[i]); } return this; }
665    public AlgebraicVector<T> AssignSin(AlgebraicVector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].AssignSin(a.elems[i]); } return this; }
666    public AlgebraicVector<T> AssignCos(AlgebraicVector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].AssignCos(a.elems[i]); } return this; }
667    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; }
668    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; }
669    public AlgebraicVector<T> AssignAbs(AlgebraicVector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].AssignAbs(a.elems[i]); } return this; }
670    public AlgebraicVector<T> AssignSgn(AlgebraicVector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].AssignSgn(a.elems[i]); } return this; }
671  }
672
673  */
674
675
676  /// <summary>
677  /// A sparse vector of algebraic types. Elements are accessed via a key of type K
678  /// </summary>
679  /// <typeparam name="K">Key type</typeparam>
680  /// <typeparam name="T">Element type</typeparam>
681  [DebuggerDisplay("SparseVector: {ToString()}")]
682  public sealed class AlgebraicSparseVector<K, T> : IAlgebraicType<AlgebraicSparseVector<K, T>> where T : IAlgebraicType<T> {
683    [DebuggerBrowsable(DebuggerBrowsableState.Never)]
684    private Dictionary<K, T> elems;
685    public IReadOnlyDictionary<K, T> Elements => elems;
686
687
688    public AlgebraicSparseVector(AlgebraicSparseVector<K, T> original) {
689      elems = original.elems.ToDictionary(kvp => kvp.Key, kvp => kvp.Value.Clone());
690    }
691
692    /// <summary>
693    ///
694    /// </summary>
695    /// <param name="keys"></param>
696    /// <param name="values">values are cloned</param>
697    public AlgebraicSparseVector(K[] keys, T[] values) {
698      if (keys.Length != values.Length) throw new ArgumentException("lengths of keys and values doesn't match in SparseVector");
699      elems = new Dictionary<K, T>(keys.Length);
700      for (int i = 0; i < keys.Length; ++i) {
702      }
703    }
704
705    public AlgebraicSparseVector() {
706      this.elems = new Dictionary<K, T>();
707    }
708
709    // keep only elements from a
710    private void AssignFromSource(AlgebraicSparseVector<K, T> a, Func<T, T, T> mapAssign) {
711      // remove elems from this which don't occur in a
712      List<K> keysToRemove = new List<K>();
713      foreach (var kvp in elems) {
715      }
716      foreach (var o in keysToRemove) elems.Remove(o); // -> zero
717
718      foreach (var kvp in a.elems) {
719        if (elems.TryGetValue(kvp.Key, out T value))
720          mapAssign(kvp.Value, value);
721        else
723      }
724    }
725
726    [DebuggerBrowsable(DebuggerBrowsableState.Never)]
727    public AlgebraicSparseVector<K, T> Zero => new AlgebraicSparseVector<K, T>();
728    [DebuggerBrowsable(DebuggerBrowsableState.Never)]
729    public AlgebraicSparseVector<K, T> One => throw new NotSupportedException();
730
731    public AlgebraicSparseVector<K, T> Scale(T s) { foreach (var kvp in elems) { kvp.Value.Mul(s); } return this; }
732    public AlgebraicSparseVector<K, T> Scale(double s) { foreach (var kvp in elems) { kvp.Value.Scale(s); } return this; }
733
734    public AlgebraicSparseVector<K, T> Assign(AlgebraicSparseVector<K, T> a) { elems.Clear(); AssignFromSource(a, (src, dest) => dest.Assign(src)); return this; }
735    public AlgebraicSparseVector<K, T> AssignInv(AlgebraicSparseVector<K, T> a) { AssignFromSource(a, (src, dest) => dest.AssignInv(src)); return this; }
736    public AlgebraicSparseVector<K, T> AssignNeg(AlgebraicSparseVector<K, T> a) { AssignFromSource(a, (src, dest) => dest.AssignNeg(src)); return this; }
737    public AlgebraicSparseVector<K, T> AssignLog(AlgebraicSparseVector<K, T> a) { AssignFromSource(a, (src, dest) => dest.AssignLog(src)); return this; }
738    public AlgebraicSparseVector<K, T> AssignExp(AlgebraicSparseVector<K, T> a) { AssignFromSource(a, (src, dest) => dest.AssignExp(src)); return this; }
739    public AlgebraicSparseVector<K, T> AssignIntPower(AlgebraicSparseVector<K, T> a, int p) { AssignFromSource(a, (src, dest) => dest.AssignIntPower(src, p)); return this; }
740    public AlgebraicSparseVector<K, T> AssignIntRoot(AlgebraicSparseVector<K, T> a, int r) { AssignFromSource(a, (src, dest) => dest.AssignIntRoot(src, r)); return this; }
741    public AlgebraicSparseVector<K, T> AssignSin(AlgebraicSparseVector<K, T> a) { AssignFromSource(a, (src, dest) => dest.AssignSin(src)); return this; }
742    public AlgebraicSparseVector<K, T> AssignCos(AlgebraicSparseVector<K, T> a) { AssignFromSource(a, (src, dest) => dest.AssignCos(src)); return this; }
743    public AlgebraicSparseVector<K, T> AssignTanh(AlgebraicSparseVector<K, T> a) { AssignFromSource(a, (src, dest) => dest.AssignTanh(src)); return this; }
744    public AlgebraicSparseVector<K, T> AssignAbs(AlgebraicSparseVector<K, T> a) { AssignFromSource(a, (src, dest) => dest.AssignAbs(src)); return this; }
745    public AlgebraicSparseVector<K, T> AssignSgn(AlgebraicSparseVector<K, T> a) { AssignFromSource(a, (src, dest) => dest.AssignSgn(src)); return this; }
746    public AlgebraicSparseVector<K, T> Add(AlgebraicSparseVector<K, T> a) {
747      foreach (var kvp in a.elems) {
748        if (elems.TryGetValue(kvp.Key, out T value))
750        else
751          elems.Add(kvp.Key, kvp.Value.Clone()); // 0 + a
752      }
753      return this;
754    }
755
756    public AlgebraicSparseVector<K, T> Sub(AlgebraicSparseVector<K, T> a) {
757      foreach (var kvp in a.elems) {
758        if (elems.TryGetValue(kvp.Key, out T value))
759          value.Sub(kvp.Value);
760        else
761          elems.Add(kvp.Key, kvp.Value.Zero.Sub(kvp.Value)); // 0 - a
762      }
763      return this;
764    }
765
766    public AlgebraicSparseVector<K, T> Mul(AlgebraicSparseVector<K, T> a) {
767      var keys = elems.Keys.ToArray();
768      foreach (var k in keys) if (!a.elems.ContainsKey(k)) elems.Remove(k); // 0 * a => 0
769      foreach (var kvp in a.elems) {
770        if (elems.TryGetValue(kvp.Key, out T value))
771          value.Mul(kvp.Value); // this * a
772      }
773      return this;
774    }
775
776    public AlgebraicSparseVector<K, T> Div(AlgebraicSparseVector<K, T> a) {
777      return Mul(a.Clone().Inv());
778    }
779
780    public AlgebraicSparseVector<K, T> Clone() {
781      return new AlgebraicSparseVector<K, T>(this);
782    }
783
784    public override string ToString() {
785      return "[" + string.Join(" ", elems.Select(kvp => kvp.Key + ": " + kvp.Value)) + "]";
786    }
787  }
788
789  // this is our own implementation of interval arithmetic
790  // for a well worked out definition of interval operations for IEEE reals see:
791  // Stahl: Interval Methods for Bounding the Range of Polynomials and Solving Systems of Nonlinear Equations, Dissertation, JKU, 1995
792  [DebuggerDisplay("[{low.Value}..{high.Value}]")]
793  public class AlgebraicInterval : IAlgebraicType<AlgebraicInterval> {
794    [DebuggerBrowsable(DebuggerBrowsableState.Never)]
795    private MultivariateDual<AlgebraicDouble> low;
796    public MultivariateDual<AlgebraicDouble> LowerBound => low.Clone();
797
798    [DebuggerBrowsable(DebuggerBrowsableState.Never)]
799    private MultivariateDual<AlgebraicDouble> high;
800    public MultivariateDual<AlgebraicDouble> UpperBound => high.Clone();
801
802
803    public AlgebraicInterval() : this(double.NegativeInfinity, double.PositiveInfinity) { }
804
805    public AlgebraicInterval(MultivariateDual<AlgebraicDouble> low, MultivariateDual<AlgebraicDouble> high) {
806      this.low = low.Clone();
807      this.high = high.Clone();
808    }
809
810    public AlgebraicInterval(double low, double high) {
811      this.low = new MultivariateDual<AlgebraicDouble>(new AlgebraicDouble(low));
812      this.high = new MultivariateDual<AlgebraicDouble>(new AlgebraicDouble(high));
813    }
814
815    [DebuggerBrowsable(DebuggerBrowsableState.Never)]
816    public AlgebraicInterval Zero => new AlgebraicInterval(0.0, 0.0);
817    [DebuggerBrowsable(DebuggerBrowsableState.Never)]
818    public AlgebraicInterval One => new AlgebraicInterval(1.0, 1.0);
819
820    public AlgebraicInterval Add(AlgebraicInterval a) {
823      return this;
824    }
825
826    public AlgebraicInterval Mul(AlgebraicInterval a) {
827      var v1 = low.Clone().Mul(a.low);
828      var v2 = low.Clone().Mul(a.high);
829      var v3 = high.Clone().Mul(a.low);
830      var v4 = high.Clone().Mul(a.high);
831
832      low = Min(Min(v1, v2), Min(v3, v4));
833      high = Max(Max(v1, v2), Max(v3, v4));
834
835      return this;
836    }
837
838    // algebraic min() / max() do not work for infinities
839    // detect and handle infinite values explicitly
840    private static MultivariateDual<AlgebraicDouble> Min(MultivariateDual<AlgebraicDouble> a, MultivariateDual<AlgebraicDouble> b) {
841      if (a.Value.IsInfinity || b.Value.IsInfinity) return a.Value < b.Value ? a : b; // NOTE: this is not differentiable but for infinite values we do not care
842      else return Algebraic.Min(a, b); // differentiable statement
843    }
844    private static MultivariateDual<AlgebraicDouble> Max(MultivariateDual<AlgebraicDouble> a, MultivariateDual<AlgebraicDouble> b) {
845      if (a.Value.IsInfinity || b.Value.IsInfinity) return a.Value >= b.Value ? a : b; // NOTE: this is not differentiable but for infinite values we do not care
846      else return Algebraic.Max(a, b); // differentiable statement
847    }
848
849    public AlgebraicInterval Assign(AlgebraicInterval a) {
850      low = a.low;
851      high = a.high;
852      return this;
853    }
854
855    public AlgebraicInterval AssignCos(AlgebraicInterval a) {
856      return AssignSin(a.Clone().Add(new AlgebraicInterval(Math.PI / 2, Math.PI / 2)));
857    }
858
859    public AlgebraicInterval Div(AlgebraicInterval a) {
860      if (a.Contains(0.0)) {
861        if (a.low.Value.Value == 0 && a.high.Value.Value == 0) {
862          if (this.low.Value >= 0) {
863            // pos / 0
864          } else if (this.high.Value <= 0) {
865            // neg / 0
866          } else {
867            low = new MultivariateDual<AlgebraicDouble>(double.NegativeInfinity);
868            high = new MultivariateDual<AlgebraicDouble>(double.PositiveInfinity);
869          }
870        } else if (a.low.Value.Value >= 0) {
871          // a is positive
872          Mul(new AlgebraicInterval(a.Clone().high.Inv(), new MultivariateDual<AlgebraicDouble>(double.PositiveInfinity)));
873        } else if (a.high.Value <= 0) {
874          // a is negative
875          Mul(new AlgebraicInterval(new MultivariateDual<AlgebraicDouble>(double.NegativeInfinity), a.low.Clone().Inv()));
876        } else {
877          // a is interval over zero
878          low = new MultivariateDual<AlgebraicDouble>(double.NegativeInfinity);
879          high = new MultivariateDual<AlgebraicDouble>(double.PositiveInfinity);
880        }
881      } else {
882        Mul(new AlgebraicInterval(a.high.Clone().Inv(), a.low.Clone().Inv())); // inverting leads to inverse roles of high and low
883      }
884      return this;
885    }
886
887    public AlgebraicInterval AssignExp(AlgebraicInterval a) {
888      low.AssignExp(a.low);
889      high.AssignExp(a.high);
890      return this;
891    }
892
893    // tanh is a bijective function
894    public AlgebraicInterval AssignTanh(AlgebraicInterval a) {
895      low.AssignTanh(a.low);
896      high.AssignTanh(a.high);
897      return this;
898    }
899
900    public AlgebraicInterval AssignIntPower(AlgebraicInterval a, int p) {
901      if (p < 0) {  // x^-3 == 1/(x^3)
902        AssignIntPower(a, -p);
903        return AssignInv(this);
904      } else if (p == 0) {
905        if (a.Contains(0.0)) {
906          // => 0^0 = 0 ; might not be relevant
907          low = new MultivariateDual<AlgebraicDouble>(0.0);
908          high = new MultivariateDual<AlgebraicDouble>(1.0);
909          return this;
910        } else {
911          // => 1
912          low = new MultivariateDual<AlgebraicDouble>(1.0);
913          high = new MultivariateDual<AlgebraicDouble>(1.0);
914          return this;
915        }
916      } else if (p == 1) return this;
917      else if (p % 2 == 0) {
918        // p is even => interval must be positive
919        if (a.Contains(0.0)) {
920          low = new MultivariateDual<AlgebraicDouble>(0.0);
921          high = Algebraic.Max(a.low.IntPower(p), a.high.IntPower(p));
922        } else {
923          var lowPower = a.low.IntPower(p);
924          var highPower = a.high.IntPower(p);
925          low = Algebraic.Min(lowPower, highPower);
926          high = Algebraic.Max(lowPower, highPower);
927        }
928      } else {
929        // p is uneven
930        if (a.Contains(0.0)) {
931          low.AssignIntPower(a.low, p);
932          high.AssignIntPower(a.high, p);
933        } else {
934          var lowPower = a.low.IntPower(p);
935          var highPower = a.high.IntPower(p);
936          low = Algebraic.Min(lowPower, highPower);
937          high = Algebraic.Max(lowPower, highPower);
938        }
939      }
940      return this;
941    }
942
943    public AlgebraicInterval AssignIntRoot(AlgebraicInterval a, int r) {
944      if (r == 0) { low = new MultivariateDual<AlgebraicDouble>(double.NaN); high = new MultivariateDual<AlgebraicDouble>(double.NaN); return this; }
945      if (r == 1) return this;
946      if (r < 0) {
947        // x^ (-1/2) = 1 / (x^(1/2))
948        AssignIntRoot(a, -r);
949        return AssignInv(this);
950      } else {
951        // root only defined for positive arguments for even roots
952        if (r % 2 == 0 && a.LowerBound.Value.Value < 0) {
953          low = new MultivariateDual<AlgebraicDouble>(double.NaN);
954          high = new MultivariateDual<AlgebraicDouble>(double.NaN);
955          return this;
956        } else {
957          low.AssignIntRoot(a.low, r);
958          high.AssignIntRoot(a.high, r);
959          return this;
960        }
961      }
962    }
963
964    public AlgebraicInterval AssignInv(AlgebraicInterval a) {
965      low = new MultivariateDual<AlgebraicDouble>(1.0);
966      high = new MultivariateDual<AlgebraicDouble>(1.0);
967      return Div(a);
968    }
969
970    public AlgebraicInterval AssignLog(AlgebraicInterval a) {
971      low.AssignLog(a.low);
972      high.AssignLog(a.high);
973      return this;
974    }
975
976    public AlgebraicInterval AssignNeg(AlgebraicInterval a) {
977      low.AssignNeg(a.high);
978      high.AssignNeg(a.low);
979      return this;
980    }
981
982    public AlgebraicInterval Scale(double s) {
983      low.Scale(s);
984      high.Scale(s);
985      if (s < 0) {
986        var t = low;
987        low = high;
988        high = t;
989      }
990      return this;
991    }
992
993    public AlgebraicInterval AssignSin(AlgebraicInterval a) {
994      var lower = a.LowerBound.Value.Value;
995      var size = a.UpperBound.Value.Value - lower;
996      if (size < 0) throw new InvalidProgramException(); // ASSERT interval >= 0;
997
998      if (size >= Math.PI * 2) {
999        low = new MultivariateDual<AlgebraicDouble>(-1.0); // zero gradient
1000        high = new MultivariateDual<AlgebraicDouble>(1.0);
1001        return this;
1002      }
1003
1004      // assume low and high are in the same quadrant
1005      low = Algebraic.Min(a.LowerBound.Clone().Sin(), a.UpperBound.Clone().Sin());
1006      high = Algebraic.Max(a.LowerBound.Clone().Sin(), a.UpperBound.Clone().Sin());
1007
1008      // override min and max if necessary
1009
1010      // shift interval 'a' into the range [-2pi .. 2pi] without changing the size of the interval to simplify the checks
1011      lower = lower % (2 * Math.PI); // lower in [-2pi .. 2pi]
1012
1013      // handle min = -1 and max = 1 cases explicitly
1014      var pi_2 = Math.PI / 2.0;
1015      var maxima = new double[] { -3 * pi_2, pi_2 };
1016      var minima = new double[] { -pi_2, 3 * pi_2 };
1017
1018      // override min and max if necessary
1019      if (maxima.Any(m => lower < m && lower + size > m)) {
1020        // max = 1
1021        high = new MultivariateDual<AlgebraicDouble>(1.0); // zero gradient
1022      }
1023
1024      if (minima.Any(m => lower < m && lower + size > m)) {
1025        // min = -1;
1026        low = new MultivariateDual<AlgebraicDouble>(-1.0); // zero gradient
1027      }
1028      return this;
1029    }
1030
1031    public AlgebraicInterval Sub(AlgebraicInterval a) {
1032      // [x1,x2] − [y1,y2] = [x1 − y2,x2 − y1]
1033      low.Sub(a.high);
1034      high.Sub(a.low);
1035      return this;
1036    }
1037
1038    public AlgebraicInterval Clone() {
1039      return new AlgebraicInterval(low, high);
1040    }
1041
1042    public bool Contains(double val) {
1043      return LowerBound.Value.Value <= val && val <= UpperBound.Value.Value;
1044    }
1045
1046    public AlgebraicInterval AssignAbs(AlgebraicInterval a) {
1047      if (a.Contains(0.0)) {
1048        var abslow = a.low.Clone().Abs();
1049        var abshigh = a.high.Clone().Abs();
1050        a.high.Assign(Algebraic.Max(abslow, abshigh));
1051        a.low.Assign(new MultivariateDual<AlgebraicDouble>(0.0)); // lost gradient for lower bound
1052      } else {
1053        var abslow = a.low.Clone().Abs();
1054        var abshigh = a.high.Clone().Abs();
1055        a.low.Assign(Algebraic.Min(abslow, abshigh));
1056        a.high.Assign(Algebraic.Max(abslow, abshigh));
1057      }
1058      return this;
1059    }
1060
1061    public AlgebraicInterval AssignSgn(AlgebraicInterval a) {
1062      low.AssignSgn(a.low);
1063      high.AssignSgn(a.high);
1064      return this;
1065    }
1066  }
1067
1068  public class Dual<V> : IAlgebraicType<Dual<V>>
1069    where V : IAlgebraicType<V> {
1070    [DebuggerBrowsable(DebuggerBrowsableState.Never)]
1071    private V v;
1072    public V Value => v;
1073
1074    [DebuggerBrowsable(DebuggerBrowsableState.Never)]
1075    private V dv;
1076    public V Derivative => dv;
1077
1078    public Dual(V v, V dv) { this.v = v; this.dv = dv; }
1079
1080    [DebuggerBrowsable(DebuggerBrowsableState.Never)]
1081    public Dual<V> Zero => new Dual<V>(Value.Zero, Derivative.Zero);
1082    [DebuggerBrowsable(DebuggerBrowsableState.Never)]
1083    public Dual<V> One => new Dual<V>(Value.One, Derivative.Zero);
1084
1085    public Dual<V> Assign(Dual<V> a) { v.Assign(a.v); dv.Assign(a.dv); return this; }
1086    public Dual<V> Scale(double s) { v.Scale(s); dv.Scale(s); return this; }
1088    public Dual<V> Sub(Dual<V> a) { v.Sub(a.v); dv.Sub(a.dv); return this; }
1089    public Dual<V> AssignNeg(Dual<V> a) { v.AssignNeg(a.v); dv.AssignNeg(a.dv); return this; }
1090    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
1091
1092    // (a(x) * b(x))' = b(x)*a(x)' + b(x)'*a(x);
1093    public Dual<V> Mul(Dual<V> a) {
1094      var t1 = a.dv.Clone().Mul(v);
1095      var t2 = dv.Clone().Mul(a.v);
1097
1098      v.Mul(a.v);
1099      return this;
1100    }
1101    public Dual<V> Div(Dual<V> a) { Mul(a.Inv()); return this; }
1102
1103    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)'
1104    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)'
1105
1106    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; }
1107    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; }
1108
1109    public Dual<V> AssignSin(Dual<V> a) { v.AssignSin(a.v); dv.Assign(a.dv).Mul(a.v.Clone().Cos()); return this; }
1110    public Dual<V> AssignCos(Dual<V> a) { v.AssignCos(a.v); dv.AssignNeg(a.dv).Mul(a.v.Clone().Sin()); return this; }
1111    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; }
1112
1113    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)|
1114    public Dual<V> AssignSgn(Dual<V> a) { v.AssignSgn(a.v); dv.Assign(a.dv.Zero); return this; }
1115
1116    public Dual<V> Clone() { return new Dual<V>(v.Clone(), dv.Clone()); }
1117
1118  }
1119
1120  /// <summary>
1121  /// An algebraic type which has a value as well as the partial derivatives of the value over multiple variables.
1122  /// </summary>
1123  /// <typeparam name="V"></typeparam>
1124  [DebuggerDisplay("v={Value}; dv={dv}")]
1125  public class MultivariateDual<V> : IAlgebraicType<MultivariateDual<V>> where V : IAlgebraicType<V>, new() {
1126    [DebuggerBrowsable(DebuggerBrowsableState.Never)]
1127    private V v;
1128    public V Value => v;
1129
1130    [DebuggerBrowsable(DebuggerBrowsableState.Never)]
1131    private AlgebraicSparseVector<object, V> dv;
1132    public AlgebraicSparseVector<object, V> Gradient => dv; // <key,value> partial derivative identified via the key
1133
1134    private MultivariateDual(MultivariateDual<V> orig) { this.v = orig.v.Clone(); this.dv = orig.dv.Clone(); }
1135
1136    /// <summary>
1137    /// Constructor without partial derivative
1138    /// </summary>
1139    /// <param name="v"></param>
1140    public MultivariateDual(V v) { this.v = v.Clone(); this.dv = new AlgebraicSparseVector<object, V>(); }
1141
1142    /// <summary>
1143    /// Constructor for multiple partial derivatives
1144    /// </summary>
1145    /// <param name="v"></param>
1146    /// <param name="keys"></param>
1147    /// <param name="dv"></param>
1148    public MultivariateDual(V v, object[] keys, V[] dv) { this.v = v.Clone(); this.dv = new AlgebraicSparseVector<object, V>(keys, dv); }
1149
1150    /// <summary>
1151    /// Constructor for a single partial derivative
1152    /// </summary>
1153    /// <param name="v"></param>
1154    /// <param name="key"></param>
1155    /// <param name="dv"></param>
1156    public MultivariateDual(V v, object key, V dv) { this.v = v.Clone(); this.dv = new AlgebraicSparseVector<object, V>(new[] { key }, new[] { dv }); }
1157
1158    /// <summary>
1159    /// Constructor with a given value and gradient. For internal use.
1160    /// </summary>
1161    /// <param name="v">The value (not cloned).</param>
1163    internal MultivariateDual(V v, AlgebraicSparseVector<object, V> gradient) { this.v = v; this.dv = gradient; }
1164
1165    public MultivariateDual<V> Clone() { return new MultivariateDual<V>(this); }
1166
1167    public MultivariateDual<V> Zero => new MultivariateDual<V>(Value.Zero, Gradient.Zero);
1168    public MultivariateDual<V> One => new MultivariateDual<V>(Value.One, Gradient.Zero);
1169
1170    public MultivariateDual<V> Scale(double s) { v.Scale(s); dv.Scale(s); return this; }
1171
1173    public MultivariateDual<V> Sub(MultivariateDual<V> a) { v.Sub(a.v); dv.Sub(a.dv); return this; }
1174    public MultivariateDual<V> Assign(MultivariateDual<V> a) { v.Assign(a.v); dv.Assign(a.dv); return this; }
1175    public MultivariateDual<V> Mul(MultivariateDual<V> a) {
1176      // (a(x) * b(x))' = b(x)*a(x)' + b(x)'*a(x);
1177      var t1 = a.dv.Clone().Scale(v);
1178      var t2 = dv.Clone().Scale(a.v);
1180
1181      v.Mul(a.v);
1182      return this;
1183    }
1184    public MultivariateDual<V> Div(MultivariateDual<V> a) { v.Div(a.v); dv.Mul(a.dv.Inv()); return this; }
1185    public MultivariateDual<V> AssignNeg(MultivariateDual<V> a) { v.AssignNeg(a.v); dv.AssignNeg(a.dv); return this; }
1186    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
1187
1188    public MultivariateDual<V> AssignSin(MultivariateDual<V> a) { v.AssignSin(a.v); dv.Assign(a.dv).Scale(a.v.Clone().Cos()); return this; }
1189    public MultivariateDual<V> AssignCos(MultivariateDual<V> a) { v.AssignCos(a.v); dv.AssignNeg(a.dv).Scale(a.v.Clone().Sin()); return this; }
1190    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)))
1191
1192    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; }
1193    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; }
1194
1195    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)'
1196    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)'
1197
1198    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
1199    public MultivariateDual<V> AssignSgn(MultivariateDual<V> a) { v.AssignSgn(a.v); dv = a.dv.Zero; return this; } // sign(f(x))' = 0;
1200
1201  }
1202}
Note: See TracBrowser for help on using the repository browser.