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

Last change on this file since 17294 was 17294, checked in by gkronber, 3 years ago

#2994: min() and max() are problematic in AutoDiff. We have used max(a,b) = (a+b + abs(b-a)) / 2 but this is problematic because of catastrophic cancellation. Min() and Max() are necessary for interval calculation. However, there we only need to determine min and max for double values (or duals). We do not need min() and max() in general.

File size: 56.3 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; } // 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
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
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) {
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) {
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
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))
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
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
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) {
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; }
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);
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>
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
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);
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  }
1259}
Note: See TracBrowser for help on using the repository browser.