Changeset 16727 for branches/2994AutoDiffForIntervals/HeuristicLab.Problems.DataAnalysis.Symbolic/3.4/Interpreter/Interpreter.cs
 Timestamp:
 03/29/19 15:01:47 (5 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

branches/2994AutoDiffForIntervals/HeuristicLab.Problems.DataAnalysis.Symbolic/3.4/Interpreter/Interpreter.cs
r16695 r16727 72 72 break; 73 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 } 74 86 case OpCodes.Exp: { 75 87 instr.value.AssignExp(code[c].value); 76 88 break; 77 89 } 78 79 90 case OpCodes.Log: { 80 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; 110 t.Add(code[c + j].value.Clone().IntPower(2)); 111 instr.value.Div(t.IntRoot(2)); 112 } 81 113 break; 82 114 } … … 221 253 } 222 254 223 public void Evaluate(ISymbolicExpressionTree tree, IDataset dataset, int[] rows, ISymbolicExpressionTreeNode[] parameterNodes, out double[] fi, out double[,] jac) { 255 /// <summary> 256 /// 257 /// </summary> 258 /// <param name="tree"></param> 259 /// <param name="dataset"></param> 260 /// <param name="rows"></param> 261 /// <param name="parameterNodes"></param> 262 /// <param name="fi">Function output. Must be allocated by the caller.</param> 263 /// <param name="jac">Jacobian matrix. Must be allocated by the caller.</param> 264 public void Evaluate(ISymbolicExpressionTree tree, IDataset dataset, int[] rows, ISymbolicExpressionTreeNode[] parameterNodes, double[] fi, double[,] jac) { 224 265 if (cachedData == null  this.dataset != dataset) { 225 266 InitCache(dataset); … … 235 276 var roundedTotal = rows.Length  remainingRows; 236 277 237 fi = new double[rows.Length];238 jac = new double[rows.Length, nParams];239 240 278 this.rows = rows; 241 279 … … 247 285 var g = code[0].value.Gradient; 248 286 for (int j = 0; j < nParams; ++j) { 249 g.Elements[j].CopyColumnTo(jac, j, rowIndex, BATCHSIZE); 287 if(g.Elements.TryGetValue(j, out AlgebraicDoubleVector v)) { 288 v.CopyColumnTo(jac, j, rowIndex, BATCHSIZE); 289 } else { 290 for (int r = 0; r < BATCHSIZE; r++) jac[rowIndex + r, j] = 0.0; 291 } 250 292 } 251 293 } … … 257 299 var g = code[0].value.Gradient; 258 300 for (int j = 0; j < nParams; ++j) 259 g.Elements[j].CopyColumnTo(jac, j, roundedTotal, remainingRows); 301 if (g.Elements.TryGetValue(j, out AlgebraicDoubleVector v)) { 302 v.CopyColumnTo(jac, j, roundedTotal, remainingRows); 303 } else { 304 for (int r = 0; r < remainingRows; r++) jac[roundedTotal + r, j] = 0.0; 305 } 260 306 } 261 307 } … … 315 361 var g = a.value.Gradient.Elements[paramIdx]; 316 362 for (int i = rowIndex; i < rows.Length && i  rowIndex < BATCHSIZE; i++) { 317 g[i ] = data[rows[i]];363 g[i  rowIndex] = data[rows[i]]; 318 364 } 319 365 } … … 377 423 public interface IAlgebraicType<T> { 378 424 T Zero { get; } 425 T One { get; } 379 426 380 427 T AssignAbs(T a); // set this to assign abs(a) … … 429 476 [DebuggerBrowsable(DebuggerBrowsableState.Never)] 430 477 public AlgebraicDouble Zero => new AlgebraicDouble(0.0); 478 [DebuggerBrowsable(DebuggerBrowsableState.Never)] 479 public AlgebraicDouble One => new AlgebraicDouble(1.0); 431 480 public AlgebraicDouble() { } 432 481 public AlgebraicDouble(double value) { this.Value = value; } … … 441 490 public AlgebraicDouble AssignSin(AlgebraicDouble a) { Value = Math.Sin(a.Value); return this; } 442 491 public AlgebraicDouble AssignCos(AlgebraicDouble a) { Value = Math.Cos(a.Value); return this; } 443 public AlgebraicDouble AssignLog(AlgebraicDouble a) { Value = Math.Log(a.Value); return this; }492 public AlgebraicDouble AssignLog(AlgebraicDouble a) { Value = a.Value <= 0?double.NegativeInfinity : Math.Log(a.Value); return this; } // alternative definiton of log to prevent NaN 444 493 public AlgebraicDouble AssignExp(AlgebraicDouble a) { Value = Math.Exp(a.Value); return this; } 445 494 public AlgebraicDouble AssignIntPower(AlgebraicDouble a, int p) { Value = Math.Pow(a.Value, p); return this; } 446 495 public AlgebraicDouble AssignIntRoot(AlgebraicDouble a, int r) { Value = Math.Pow(a.Value, 1.0 / r); return this; } 447 496 public AlgebraicDouble AssignAbs(AlgebraicDouble a) { Value = Math.Abs(a.Value); return this; } 448 public AlgebraicDouble AssignSgn(AlgebraicDouble a) { Value = Math.Sign(a.Value); return this; }497 public AlgebraicDouble AssignSgn(AlgebraicDouble a) { Value = double.IsNaN(a.Value)? double.NaN : Math.Sign(a.Value); return this; } 449 498 public AlgebraicDouble Clone() { return new AlgebraicDouble(Value); } 450 499 … … 473 522 [DebuggerBrowsable(DebuggerBrowsableState.Never)] 474 523 public AlgebraicDoubleVector Zero => new AlgebraicDoubleVector(new double[this.Length]); // must return vector of same length as this (therefore Zero is not static) 524 [DebuggerBrowsable(DebuggerBrowsableState.Never)] 525 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) 475 526 public AlgebraicDoubleVector Assign(AlgebraicDoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = a.arr[i]; } return this; } 476 527 public AlgebraicDoubleVector Add(AlgebraicDoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] += a.arr[i]; } return this; } … … 496 547 } 497 548 498 public voidAssignConstant(double constantValue) {549 public AlgebraicDoubleVector AssignConstant(double constantValue) { 499 550 for (int i = 0; i < arr.Length; ++i) { 500 551 arr[i] = constantValue; 501 552 } 553 return this; 502 554 } 503 555 … … 522 574 } 523 575 576 577 /* 524 578 // vectors of algebraic types 525 public sealed class AlgebraicVector<T> : IAlgebraicType<AlgebraicVector<T>> where T : IAlgebraicType<T> {579 public sealed class AlgebraicVector<T> : IAlgebraicType<AlgebraicVector<T>> where T : IAlgebraicType<T>, new() { 526 580 private T[] elems; 527 581 … … 561 615 public AlgebraicVector<T> Clone() { return new AlgebraicVector<T>(elems); } 562 616 563 public AlgebraicVector<T> Concat(AlgebraicVector<T> other) {564 var oldLen = Length;565 Array.Resize(ref this.elems, oldLen + other.Length);566 for (int i = oldLen; i < Length; i++) {567 elems[i] = other.elems[i  oldLen].Clone();568 }569 return this;570 }571 617 572 618 [DebuggerBrowsable(DebuggerBrowsableState.Never)] 573 619 public AlgebraicVector<T> Zero => new AlgebraicVector<T>(Length); 620 [DebuggerBrowsable(DebuggerBrowsableState.Never)] 621 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; } } 574 622 public AlgebraicVector<T> Assign(AlgebraicVector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].Assign(a.elems[i]); } return this; } 575 623 public AlgebraicVector<T> Add(AlgebraicVector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].Add(a.elems[i]); } return this; } … … 591 639 } 592 640 641 */ 642 593 643 594 644 /// <summary> … … 625 675 } 626 676 627 628 629 // combined elements from both vectors630 private void UnionAssign(AlgebraicSparseVector<K, T> a, Func<T, T, T> mapAssign) {631 // elements from a632 foreach (var kvp in a.elems) {633 // this = f(a, this)634 if (elems.TryGetValue(kvp.Key, out T value))635 mapAssign(kvp.Value, value);636 else {637 // this = f(a, 0)638 var newValue = kvp.Value.Zero;639 elems.Add(kvp.Key, newValue);640 mapAssign(kvp.Value, newValue);641 }642 }643 // elements from this (without a)644 foreach (var kvp in elems) {645 if (a.elems.ContainsKey(kvp.Key)) continue; // already processed above646 // this = f(0, this)647 mapAssign(kvp.Value.Zero, kvp.Value);648 }649 }650 651 // keep only elements in both vectors652 private void IntersectAssign(AlgebraicSparseVector<K, T> a, Func<T, T, T> mapAssign) {653 List<K> keysToRemove = new List<K>();654 foreach (var kvp in elems) {655 if (a.elems.TryGetValue(kvp.Key, out T value))656 mapAssign(value, kvp.Value);657 else658 keysToRemove.Add(kvp.Key);659 }660 foreach (var o in keysToRemove) elems.Remove(o); // > zero661 }662 663 677 // keep only elements from a 664 678 private void AssignFromSource(AlgebraicSparseVector<K, T> a, Func<T, T, T> mapAssign) { … … 680 694 [DebuggerBrowsable(DebuggerBrowsableState.Never)] 681 695 public AlgebraicSparseVector<K, T> Zero => new AlgebraicSparseVector<K, T>(); 696 [DebuggerBrowsable(DebuggerBrowsableState.Never)] 697 public AlgebraicSparseVector<K, T> One => throw new NotSupportedException(); 682 698 683 699 public AlgebraicSparseVector<K, T> Scale(T s) { foreach (var kvp in elems) { kvp.Value.Mul(s); } return this; } … … 685 701 686 702 public AlgebraicSparseVector<K, T> Assign(AlgebraicSparseVector<K, T> a) { elems.Clear(); AssignFromSource(a, (src, dest) => dest.Assign(src)); return this; } 687 public AlgebraicSparseVector<K, T> Add(AlgebraicSparseVector<K, T> a) { UnionAssign(a, (src, dest) => dest.Add(src)); return this; }688 public AlgebraicSparseVector<K, T> Sub(AlgebraicSparseVector<K, T> a) { UnionAssign(a, (src, dest) => dest.Sub(src)); return this; }689 public AlgebraicSparseVector<K, T> Mul(AlgebraicSparseVector<K, T> a) { IntersectAssign(a, (src, dest) => dest.Mul(src)); return this; }690 public AlgebraicSparseVector<K, T> Div(AlgebraicSparseVector<K, T> a) { UnionAssign(a, (src, dest) => dest.Div(src)); return this; }691 703 public AlgebraicSparseVector<K, T> AssignInv(AlgebraicSparseVector<K, T> a) { AssignFromSource(a, (src, dest) => dest.AssignInv(src)); return this; } 692 704 public AlgebraicSparseVector<K, T> AssignNeg(AlgebraicSparseVector<K, T> a) { AssignFromSource(a, (src, dest) => dest.AssignNeg(src)); return this; } … … 699 711 public AlgebraicSparseVector<K, T> AssignAbs(AlgebraicSparseVector<K, T> a) { AssignFromSource(a, (src, dest) => dest.AssignAbs(src)); return this; } 700 712 public AlgebraicSparseVector<K, T> AssignSgn(AlgebraicSparseVector<K, T> a) { AssignFromSource(a, (src, dest) => dest.AssignSgn(src)); return this; } 713 public AlgebraicSparseVector<K, T> Add(AlgebraicSparseVector<K, T> a) { 714 foreach (var kvp in a.elems) { 715 if (elems.TryGetValue(kvp.Key, out T value)) 716 value.Add(kvp.Value); 717 else 718 elems.Add(kvp.Key, kvp.Value.Clone()); // 0 + a 719 } 720 return this; 721 } 722 723 public AlgebraicSparseVector<K, T> Sub(AlgebraicSparseVector<K, T> a) { 724 foreach (var kvp in a.elems) { 725 if (elems.TryGetValue(kvp.Key, out T value)) 726 value.Sub(kvp.Value); 727 else 728 elems.Add(kvp.Key, kvp.Value.Zero.Sub(kvp.Value)); // 0  a 729 } 730 return this; 731 } 732 733 public AlgebraicSparseVector<K, T> Mul(AlgebraicSparseVector<K, T> a) { 734 var keys = elems.Keys.ToArray(); 735 foreach (var k in keys) if (!a.elems.ContainsKey(k)) elems.Remove(k); // 0 * a => 0 736 foreach (var kvp in a.elems) { 737 if (elems.TryGetValue(kvp.Key, out T value)) 738 value.Mul(kvp.Value); // this * a 739 } 740 return this; 741 } 742 743 public AlgebraicSparseVector<K, T> Div(AlgebraicSparseVector<K, T> a) { 744 return Mul(a.Clone().Inv()); 745 } 701 746 702 747 public AlgebraicSparseVector<K, T> Clone() { … … 734 779 [DebuggerBrowsable(DebuggerBrowsableState.Never)] 735 780 public AlgebraicInterval Zero => new AlgebraicInterval(0.0, 0.0); 781 public AlgebraicInterval One => new AlgebraicInterval(1.0, 1.0); 736 782 public AlgebraicInterval Add(AlgebraicInterval a) { 737 783 low.Add(a.low); … … 799 845 if (a.Contains(0.0)) { 800 846 low = new MultivariateDual<AlgebraicDouble>(0.0); 801 high = Algebraic.Max( low.Clone().IntPower(p), high.Clone().IntPower(p));847 high = Algebraic.Max(a.low.IntPower(p), a.high.IntPower(p)); 802 848 } else { 803 var lowPower = low.Clone().IntPower(p);804 var highPower = high.Clone().IntPower(p);849 var lowPower = a.low.IntPower(p); 850 var highPower = a.high.IntPower(p); 805 851 low = Algebraic.Min(lowPower, highPower); 806 852 high = Algebraic.Max(lowPower, highPower); … … 808 854 } else { 809 855 // p is uneven 810 var lowPower = low.Clone().IntPower(p);811 var highPower = high.Clone().IntPower(p);856 var lowPower = a.low.IntPower(p); 857 var highPower = a.high.IntPower(p); 812 858 low = Algebraic.Min(lowPower, highPower); 813 859 high = Algebraic.Max(lowPower, highPower); … … 951 997 [DebuggerBrowsable(DebuggerBrowsableState.Never)] 952 998 public Dual<V> Zero => new Dual<V>(Value.Zero, Derivative.Zero); 999 [DebuggerBrowsable(DebuggerBrowsableState.Never)] 1000 public Dual<V> One => new Dual<V>(Value.One, Derivative.Zero); 953 1001 954 1002 public Dual<V> Assign(Dual<V> a) { v.Assign(a.v); dv.Assign(a.dv); return this; } … … 1034 1082 1035 1083 public MultivariateDual<V> Zero => new MultivariateDual<V>(Value.Zero, Gradient.Zero); 1084 public MultivariateDual<V> One => new MultivariateDual<V>(Value.One, Gradient.Zero); 1036 1085 1037 1086 public MultivariateDual<V> Scale(double s) { v.Scale(s); dv.Scale(s); return this; }
Note: See TracChangeset
for help on using the changeset viewer.