Changeset 16693 for branches/2994AutoDiffForIntervals/HeuristicLab.Problems.DataAnalysis.Symbolic/3.4/Interpreter/Interpreter.cs
 Timestamp:
 03/19/19 10:19:40 (5 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

branches/2994AutoDiffForIntervals/HeuristicLab.Problems.DataAnalysis.Symbolic/3.4/Interpreter/Interpreter.cs
r16682 r16693 4 4 using HeuristicLab.Common; 5 5 using HeuristicLab.Encodings.SymbolicExpressionTreeEncoding; 6 using HEAL.Attic;7 6 8 7 namespace HeuristicLab.Problems.DataAnalysis.Symbolic { … … 117 116 118 117 119 public class VectorEvaluator : Interpreter<DoubleVector> {118 public sealed class VectorEvaluator : Interpreter<DoubleVector> { 120 119 private const int BATCHSIZE = 128; 121 120 [ThreadStatic] … … 192 191 } 193 192 194 public class VectorAutoDiffEvaluator : Interpreter<MultivariateDual<DoubleVector>> {193 public sealed class VectorAutoDiffEvaluator : Interpreter<MultivariateDual<DoubleVector>> { 195 194 private const int BATCHSIZE = 128; 196 195 [ThreadStatic] … … 318 317 319 318 320 public class IntervalEvaluator : Interpreter<AlgebraicInterval> {319 public sealed class IntervalEvaluator : Interpreter<AlgebraicInterval> { 321 320 [ThreadStatic] 322 321 private Dictionary<string, Interval> intervals; … … 352 351 instruction.dblVal = constant.Value; 353 352 instruction.value = new AlgebraicInterval( 354 new MultivariateDual< Double>(constant.Value, constant, 1.0),355 new MultivariateDual< Double>(constant.Value, constant, 1.0) // use node as key353 new MultivariateDual<AlgebraicDouble>(constant.Value, constant, 1.0), 354 new MultivariateDual<AlgebraicDouble>(constant.Value, constant, 1.0) // use node as key 356 355 ); 357 356 } … … 360 359 instruction.dblVal = variable.Weight; 361 360 instruction.value = new AlgebraicInterval( 362 low: new MultivariateDual< Double>(intervals[variable.VariableName].LowerBound, variable, intervals[variable.VariableName].LowerBound), // bounds change by variable value d/dc (c I(var)) = I(var)363 high: new MultivariateDual< Double>(intervals[variable.VariableName].UpperBound, variable, intervals[variable.VariableName].UpperBound)361 low: new MultivariateDual<AlgebraicDouble>(intervals[variable.VariableName].LowerBound, variable, intervals[variable.VariableName].LowerBound), // bounds change by variable value d/dc (c I(var)) = I(var) 362 high: new MultivariateDual<AlgebraicDouble>(intervals[variable.VariableName].UpperBound, variable, intervals[variable.VariableName].UpperBound) 364 363 ); 365 364 } … … 371 370 372 371 public interface IAlgebraicType<T> { 372 T Zero { get; } 373 373 374 T AssignAbs(T a); // set this to assign abs(a) 374 375 T Assign(T a); // assign this to same value as a (copy!) … … 414 415 415 416 // algebraic type wrapper for a double value 416 public class Double : IAlgebraicType<Double> {417 public static implicit operator Double(double value) { return newDouble(value); }418 public static implicit operator double( Double value) { return value.Value; }417 public sealed class AlgebraicDouble : IAlgebraicType<AlgebraicDouble> { 418 public static implicit operator AlgebraicDouble(double value) { return new AlgebraicDouble(value); } 419 public static implicit operator double(AlgebraicDouble value) { return value.Value; } 419 420 public double Value; 420 public Double() { } 421 public Double(double value) { this.Value = value; }422 public Double Add(Double a) { Value += a.Value; return this;}423 public Double Assign(Double a) { Value = a.Value; return this; }424 425 public Double AssignAbs(Double a) { Value = Math.Abs(a.Value); return this; }426 public Double AssignCos(Double a) { Value = Math.Cos(a.Value); return this; }427 public Double AssignExp(Double a) { Value = Math.Exp(a.Value); return this; }428 public Double AssignIntPower(Double a, int p) { Value = Math.Pow(a.Value, p); return this; }429 public Double AssignIntRoot(Double a, int r) { Value = Math.Pow(a.Value, 1.0 / r); return this; }430 public Double AssignInv(Double a) { Value = 1.0 / a.Value; return this; }431 public Double AssignLog(Double a) { Value = Math.Log(a.Value); return this; }432 public Double AssignNeg(Double a) { Value = a.Value; return this; }433 public Double AssignSin(Double a) { Value = Math.Sin(a.Value); return this; }434 public Double AssignSgn(Double a) { Value = Math.Sign(a.Value); return this; }435 public Double Clone() { return new Double(Value); }436 public Double Div(Double a) { Value /= a.Value; return this; }437 public Double Mul(Double a) { Value *= a.Value; return this; }438 public Double Scale(double s) { Value *= s; return this; }439 public Double Sub(Double a) { Value = a.Value; return this; }440 421 422 public AlgebraicDouble Zero => new AlgebraicDouble(0.0); 423 public AlgebraicDouble() { } 424 public AlgebraicDouble(double value) { this.Value = value; } 425 public AlgebraicDouble Assign(AlgebraicDouble a) { Value = a.Value; return this; } 426 public AlgebraicDouble Add(AlgebraicDouble a) { Value += a.Value; return this; } 427 public AlgebraicDouble Sub(AlgebraicDouble a) { Value = a.Value; return this; } 428 public AlgebraicDouble Mul(AlgebraicDouble a) { Value *= a.Value; return this; } 429 public AlgebraicDouble Div(AlgebraicDouble a) { Value /= a.Value; return this; } 430 public AlgebraicDouble Scale(double s) { Value *= s; return this; } 431 public AlgebraicDouble AssignInv(AlgebraicDouble a) { Value = 1.0 / a.Value; return this; } 432 public AlgebraicDouble AssignNeg(AlgebraicDouble a) { Value = a.Value; return this; } 433 public AlgebraicDouble AssignSin(AlgebraicDouble a) { Value = Math.Sin(a.Value); return this; } 434 public AlgebraicDouble AssignCos(AlgebraicDouble a) { Value = Math.Cos(a.Value); return this; } 435 public AlgebraicDouble AssignLog(AlgebraicDouble a) { Value = Math.Log(a.Value); return this; } 436 public AlgebraicDouble AssignExp(AlgebraicDouble a) { Value = Math.Exp(a.Value); return this; } 437 public AlgebraicDouble AssignIntPower(AlgebraicDouble a, int p) { Value = Math.Pow(a.Value, p); return this; } 438 public AlgebraicDouble AssignIntRoot(AlgebraicDouble a, int r) { Value = Math.Pow(a.Value, 1.0 / r); return this; } 439 public AlgebraicDouble AssignAbs(AlgebraicDouble a) { Value = Math.Abs(a.Value); return this; } 440 public AlgebraicDouble AssignSgn(AlgebraicDouble a) { Value = Math.Sign(a.Value); return this; } 441 public AlgebraicDouble Clone() { return new AlgebraicDouble(Value); } 441 442 } 442 443 … … 445 446 private double[] arr; 446 447 448 447 449 public double this[int idx] { get { return arr[idx]; } set { arr[idx] = value; } } 450 451 public int Length => arr.Length; 448 452 449 453 public DoubleVector(int length) { … … 461 465 } 462 466 463 public DoubleVector(int length, double constantValue) : this(length) { 464 for (int i = 0; i < length; ++i) arr[i] = constantValue; 465 } 466 467 public void CopyTo(double[] dest, int idx, int length) { 468 Array.Copy(arr, 0, dest, idx, length); 469 } 470 471 public void CopyRowTo(double[,] dest, int row) { 472 for (int j = 0; j < arr.Length; ++j) dest[row, j] = arr[j]; 473 } 474 internal void CopyColumnTo(double[,] dest, int column, int row, int len) { 475 for (int j = 0; j < len; ++j) dest[row + j, column] = arr[j]; 476 } 477 478 public void AssignConstant(double constantValue) { 479 for (int i = 0; i < arr.Length; ++i) { 480 arr[i] = constantValue; 481 } 482 } 483 484 public DoubleVector Assign(DoubleVector a) { if (arr == null) arr = new double[a.arr.Length]; for (int i = 0; i < arr.Length; ++i) { arr[i] = a.arr[i]; } return this; } 485 public DoubleVector AssignCos(DoubleVector a) { if (arr == null) arr = new double[a.arr.Length]; for (int i = 0; i < arr.Length; ++i) { arr[i] = Math.Cos(a.arr[i]); } return this; } 486 public DoubleVector Div(DoubleVector a) { if (arr == null) arr = new double[a.arr.Length]; for (int i = 0; i < arr.Length; ++i) { arr[i] /= a.arr[i]; } return this; } 487 public DoubleVector AssignExp(DoubleVector a) { if (arr == null) arr = new double[a.arr.Length]; for (int i = 0; i < arr.Length; ++i) { arr[i] = Math.Exp(a.arr[i]); } return this; } 488 public DoubleVector AssignIntPower(DoubleVector a, int p) { throw new NotImplementedException(); } 489 public DoubleVector AssignIntRoot(DoubleVector a, int r) { throw new NotImplementedException(); } 490 public DoubleVector AssignInv(DoubleVector a) { if (arr == null) arr = new double[a.arr.Length]; for (int i = 0; i < arr.Length; ++i) { arr[i] = 1.0 / a.arr[i]; } return this; } 491 public DoubleVector AssignLog(DoubleVector a) { if (arr == null) arr = new double[a.arr.Length]; for (int i = 0; i < arr.Length; ++i) { arr[i] = Math.Log(a.arr[i]); } return this; } 467 public DoubleVector Zero => new DoubleVector(new double[this.Length]); // must return vector of same length as this (therefore Zero is not static) 468 public DoubleVector Assign(DoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = a.arr[i]; } return this; } 492 469 public DoubleVector Add(DoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] += a.arr[i]; } return this; } 470 public DoubleVector Sub(DoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = a.arr[i]; } return this; } 493 471 public DoubleVector Mul(DoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] *= a.arr[i]; } return this; } 494 public DoubleVector AssignNeg(DoubleVector a) { if (arr == null) arr = new double[a.arr.Length]; for (int i = 0; i < arr.Length; ++i) { arr[i] = a.arr[i]; } return this; } 472 public DoubleVector Div(DoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] /= a.arr[i]; } return this; } 473 public DoubleVector AssignNeg(DoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = a.arr[i]; } return this; } 474 public DoubleVector AssignInv(DoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = 1.0 / a.arr[i]; } return this; } 495 475 public DoubleVector Scale(double s) { for (int i = 0; i < arr.Length; ++i) { arr[i] *= s; } return this; } 496 public DoubleVector AssignSin(DoubleVector a) { if (arr == null) arr = new double[a.arr.Length]; for (int i = 0; i < arr.Length; ++i) { arr[i] = Math.Sin(a.arr[i]); } return this; } 497 public DoubleVector Sub(DoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = a.arr[i]; } return this; } 498 public DoubleVector AssignAbs(DoubleVector a) { if (arr == null) arr = new double[a.arr.Length]; for (int i = 0; i < arr.Length; ++i) { arr[i] = Math.Abs(a.arr[i]); } return this; } 499 public DoubleVector AssignSgn(DoubleVector a) { if (arr == null) arr = new double[a.arr.Length]; for (int i = 0; i < arr.Length; ++i) { arr[i] = Math.Sign(a.arr[i]); } return this; } 476 public DoubleVector AssignLog(DoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = Math.Log(a.arr[i]); } return this; } 477 public DoubleVector AssignSin(DoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = Math.Sin(a.arr[i]); } return this; } 478 public DoubleVector AssignExp(DoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = Math.Exp(a.arr[i]); } return this; } 479 public DoubleVector AssignCos(DoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = Math.Cos(a.arr[i]); } return this; } 480 public DoubleVector AssignIntPower(DoubleVector a, int p) { for (int i = 0; i < arr.Length; ++i) { arr[i] = Math.Pow(a.arr[i], p); } return this; } 481 public DoubleVector AssignIntRoot(DoubleVector a, int r) { for (int i = 0; i < arr.Length; ++i) { arr[i] = Math.Pow(a.arr[i], 1.0 / r); } return this; } 482 public DoubleVector AssignAbs(DoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = Math.Abs(a.arr[i]); } return this; } 483 public DoubleVector AssignSgn(DoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = Math.Sign(a.arr[i]); } return this; } 484 500 485 public DoubleVector Clone() { 501 486 var v = new DoubleVector(this.arr.Length); … … 504 489 } 505 490 491 public void AssignConstant(double constantValue) { 492 for (int i = 0; i < arr.Length; ++i) { 493 arr[i] = constantValue; 494 } 495 } 496 497 public void CopyTo(double[] dest, int idx, int length) { 498 Array.Copy(arr, 0, dest, idx, length); 499 } 500 506 501 public void CopyFrom(double[] data, int rowIndex) { 507 502 Array.Copy(data, rowIndex, arr, 0, Math.Min(arr.Length, data.Length  rowIndex)); 508 503 } 509 504 505 public void CopyRowTo(double[,] dest, int row) { 506 for (int j = 0; j < arr.Length; ++j) dest[row, j] = arr[j]; 507 } 508 509 internal void CopyColumnTo(double[,] dest, int column, int row, int len) { 510 for (int j = 0; j < len; ++j) dest[row + j, column] = arr[j]; 511 } 510 512 } 511 513 512 514 // vectors of algebraic types 513 public class Vector<T> : IAlgebraicType<Vector<T>> where T : IAlgebraicType<T> {515 public sealed class Vector<T> : IAlgebraicType<Vector<T>> where T : IAlgebraicType<T> { 514 516 private T[] elems; 515 517 … … 517 519 518 520 public int Length => elems.Length; 521 519 522 520 523 private Vector() { } … … 527 530 /// 528 531 /// </summary> 529 /// <param name="elems">The array is copied </param>532 /// <param name="elems">The array is copied (elementwise clone)</param> 530 533 public Vector(T[] elems) { 531 534 this.elems = new T[elems.Length]; … … 562 565 } 563 566 567 public Vector<T> Zero => new Vector<T>(Length); 568 public Vector<T> Assign(Vector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].Assign(a.elems[i]); } return this; } 564 569 public Vector<T> Add(Vector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].Add(a.elems[i]); } return this; } 565 public Vector<T> Assign(Vector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].Assign(a.elems[i]); } return this; } 570 public Vector<T> Sub(Vector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].Sub(a.elems[i]); } return this; } 571 public Vector<T> Mul(Vector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].Mul(a.elems[i]); } return this; } 572 public Vector<T> Div(Vector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].Div(a.elems[i]); } return this; } 573 public Vector<T> AssignNeg(Vector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].AssignNeg(a.elems[i]); } return this; } 574 public Vector<T> Scale(double s) { for (int i = 0; i < elems.Length; ++i) { elems[i].Scale(s); } return this; } 575 public Vector<T> Scale(T s) { for (int i = 0; i < elems.Length; ++i) { elems[i].Mul(s); } return this; } 576 public Vector<T> AssignInv(Vector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].AssignInv(a.elems[i]); } return this; } 577 public Vector<T> AssignLog(Vector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].AssignLog(a.elems[i]); } return this; } 578 public Vector<T> AssignExp(Vector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].AssignExp(a.elems[i]); } return this; } 579 public Vector<T> AssignSin(Vector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].AssignSin(a.elems[i]); } return this; } 566 580 public Vector<T> AssignCos(Vector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].AssignCos(a.elems[i]); } return this; } 567 public Vector<T> AssignExp(Vector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].AssignExp(a.elems[i]); } return this; }568 581 public Vector<T> AssignIntPower(Vector<T> a, int p) { for (int i = 0; i < elems.Length; ++i) { elems[i].AssignIntPower(a.elems[i], p); } return this; } 569 582 public Vector<T> AssignIntRoot(Vector<T> a, int r) { for (int i = 0; i < elems.Length; ++i) { elems[i].AssignIntRoot(a.elems[i], r); } return this; } 570 public Vector<T> AssignInv(Vector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].AssignInv(a.elems[i]); } return this; }571 public Vector<T> AssignLog(Vector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].AssignLog(a.elems[i]); } return this; }572 public Vector<T> AssignNeg(Vector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].AssignNeg(a.elems[i]); } return this; }573 public Vector<T> AssignSin(Vector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].AssignSin(a.elems[i]); } return this; }574 public Vector<T> Div(Vector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].Div(a.elems[i]); } return this; }575 public Vector<T> Mul(Vector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].Mul(a.elems[i]); } return this; }576 public Vector<T> Scale(double s) { for (int i = 0; i < elems.Length; ++i) { elems[i].Scale(s); } return this; }577 public Vector<T> Scale(T s) { for (int i = 0; i < elems.Length; ++i) { elems[i].Mul(s); } return this; }578 public Vector<T> Sub(Vector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].Sub(a.elems[i]); } return this; }579 583 public Vector<T> AssignAbs(Vector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].AssignAbs(a.elems[i]); } return this; } 580 584 public Vector<T> AssignSgn(Vector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].AssignSgn(a.elems[i]); } return this; } … … 587 591 /// <typeparam name="K">Key type</typeparam> 588 592 /// <typeparam name="T">Element type</typeparam> 589 public class SparseVector<K, T> : IAlgebraicType<SparseVector<K, T>> where T : IAlgebraicType<T> {593 public sealed class SparseVector<K, T> : IAlgebraicType<SparseVector<K, T>> where T : IAlgebraicType<T> { 590 594 591 595 private Dictionary<K, T> elems; 592 593 596 public IReadOnlyDictionary<K, T> Elements => elems; 597 594 598 595 599 public SparseVector(SparseVector<K, T> original) { … … 597 601 } 598 602 603 /// <summary> 604 /// 605 /// </summary> 606 /// <param name="keys"></param> 607 /// <param name="values">values are cloned</param> 599 608 public SparseVector(K[] keys, T[] values) { 600 609 if (keys.Length != values.Length) throw new ArgumentException("lengths of keys and values doesn't match in SparseVector"); 601 610 elems = new Dictionary<K, T>(keys.Length); 602 611 for (int i = 0; i < keys.Length; ++i) { 603 elems.Add(keys[i], values[i] );612 elems.Add(keys[i], values[i].Clone()); 604 613 } 605 614 } … … 609 618 } 610 619 611 public SparseVector<K, T> Add(SparseVector<K, T> a) { 620 621 622 private void AssignTransformed(SparseVector<K, T> a, Func<T, T, T> mapAssign) { 612 623 foreach (var kvp in a.elems) { 613 if (elems.TryGetValue(kvp.Key, out var value)) { 614 value.Add(kvp.Value); 615 } else { 616 elems.Add(kvp.Key, kvp.Value.Clone()); 624 if (elems.TryGetValue(kvp.Key, out T value)) 625 mapAssign(kvp.Value, value); 626 else { 627 var newValue = kvp.Value.Zero; 628 elems.Add(kvp.Key, newValue); 629 mapAssign(kvp.Value, newValue); 617 630 } 618 631 } 619 return this; 620 } 621 622 public SparseVector<K, T> Scale(T s) { 623 foreach (var kvp in elems) { 624 kvp.Value.Mul(s); 625 } 626 return this; 627 } 628 629 public SparseVector<K, T> Scale(double s) { 630 foreach (var kvp in elems) { 631 kvp.Value.Scale(s); 632 } 633 return this; 634 } 635 636 637 public SparseVector<K, T> Assign(SparseVector<K, T> a) { 638 elems.Clear(); 639 elems = a.elems.ToDictionary(kvp => kvp.Key, kvp => kvp.Value.Clone()); 640 return this; 641 } 642 643 public SparseVector<K, T> AssignCos(SparseVector<K, T> a) { 644 throw new NotImplementedException(); 645 } 646 647 public SparseVector<K, T> AssignExp(SparseVector<K, T> a) { 648 throw new NotImplementedException(); 649 } 650 651 public SparseVector<K, T> AssignIntPower(SparseVector<K, T> a, int p) { 652 throw new NotImplementedException(); 653 } 654 655 public SparseVector<K, T> AssignIntRoot(SparseVector<K, T> a, int r) { 656 throw new NotImplementedException(); 657 } 658 659 public SparseVector<K, T> AssignInv(SparseVector<K, T> a) { 660 throw new NotImplementedException(); 661 } 662 663 public SparseVector<K, T> AssignLog(SparseVector<K, T> a) { 664 throw new NotImplementedException(); 665 } 666 667 public SparseVector<K, T> AssignNeg(SparseVector<K, T> a) { 668 throw new NotImplementedException(); 669 } 670 671 public SparseVector<K, T> AssignSin(SparseVector<K, T> a) { 672 throw new NotImplementedException(); 673 } 632 } 633 634 public SparseVector<K, T> Zero => new SparseVector<K, T>(); 635 636 public SparseVector<K, T> Scale(T s) { foreach (var kvp in elems) { kvp.Value.Mul(s); } return this; } 637 public SparseVector<K, T> Scale(double s) { foreach (var kvp in elems) { kvp.Value.Scale(s); } return this; } 638 639 public SparseVector<K, T> Assign(SparseVector<K, T> a) { AssignTransformed(a, (src, dest) => dest.Assign(src)); return this; } 640 public SparseVector<K, T> Add(SparseVector<K, T> a) { AssignTransformed(a, (src, dest) => dest.Add(src)); return this; } 641 public SparseVector<K, T> Mul(SparseVector<K, T> a) { AssignTransformed(a, (src, dest) => dest.Mul(src)); return this; } 642 public SparseVector<K, T> Sub(SparseVector<K, T> a) { AssignTransformed(a, (src, dest) => dest.Sub(src)); return this; } 643 public SparseVector<K, T> Div(SparseVector<K, T> a) { AssignTransformed(a, (src, dest) => dest.Div(src)); return this; } 644 public SparseVector<K, T> AssignInv(SparseVector<K, T> a) { AssignTransformed(a, (src, dest) => dest.AssignInv(src)); return this; } 645 public SparseVector<K, T> AssignNeg(SparseVector<K, T> a) { AssignTransformed(a, (src, dest) => dest.AssignNeg(src)); return this; } 646 public SparseVector<K, T> AssignLog(SparseVector<K, T> a) { AssignTransformed(a, (src, dest) => dest.AssignLog(src)); return this; } 647 public SparseVector<K, T> AssignExp(SparseVector<K, T> a) { AssignTransformed(a, (src, dest) => dest.AssignExp(src)); return this; } 648 public SparseVector<K, T> AssignIntPower(SparseVector<K, T> a, int p) { AssignTransformed(a, (src, dest) => dest.AssignIntPower(src, p)); return this; } 649 public SparseVector<K, T> AssignIntRoot(SparseVector<K, T> a, int r) { AssignTransformed(a, (src, dest) => dest.AssignIntRoot(src, r)); return this; } 650 public SparseVector<K, T> AssignSin(SparseVector<K, T> a) { AssignTransformed(a, (src, dest) => dest.AssignSin(src)); return this; } 651 public SparseVector<K, T> AssignCos(SparseVector<K, T> a) { AssignTransformed(a, (src, dest) => dest.AssignCos(src)); return this; } 652 public SparseVector<K, T> AssignAbs(SparseVector<K, T> a) { AssignTransformed(a, (src, dest) => dest.AssignAbs(src)); return this; } 653 public SparseVector<K, T> AssignSgn(SparseVector<K, T> a) { AssignTransformed(a, (src, dest) => dest.AssignSgn(src)); return this; } 674 654 675 655 public SparseVector<K, T> Clone() { 676 656 return new SparseVector<K, T>(this); 677 657 } 678 679 public SparseVector<K, T> Div(SparseVector<K, T> a) {680 throw new NotImplementedException();681 }682 683 public SparseVector<K, T> Mul(SparseVector<K, T> a) {684 throw new NotImplementedException();685 }686 687 public SparseVector<K, T> Sub(SparseVector<K, T> a) {688 foreach (var kvp in a.elems) {689 if (elems.TryGetValue(kvp.Key, out var value)) {690 value.Sub(kvp.Value);691 } else {692 elems.Add(kvp.Key, kvp.Value.Clone().Neg());693 }694 }695 return this;696 }697 698 public SparseVector<K, T> AssignAbs(SparseVector<K, T> a) {699 elems.Clear();700 foreach (var kvp in a.elems) {701 elems.Add(kvp.Key, kvp.Value.Clone().Abs());702 }703 return this;704 }705 706 public SparseVector<K, T> AssignSgn(SparseVector<K, T> a) {707 elems.Clear();708 foreach (var kvp in a.elems) {709 elems.Add(kvp.Key, kvp.Value.Clone().Sgn());710 }711 return this;712 }713 658 } 714 659 715 660 716 661 public class AlgebraicInterval : IAlgebraicType<AlgebraicInterval> { 717 private MultivariateDual<Double> low; 718 private MultivariateDual<Double> high; 719 720 public MultivariateDual<Double> LowerBound => low.Clone(); 721 public MultivariateDual<Double> UpperBound => high.Clone(); 662 private MultivariateDual<AlgebraicDouble> low; 663 private MultivariateDual<AlgebraicDouble> high; 664 665 public MultivariateDual<AlgebraicDouble> LowerBound => low.Clone(); 666 public MultivariateDual<AlgebraicDouble> UpperBound => high.Clone(); 667 722 668 723 669 public AlgebraicInterval() : this(double.NegativeInfinity, double.PositiveInfinity) { } 724 670 725 public AlgebraicInterval(MultivariateDual< Double> low, MultivariateDual<Double> high) {671 public AlgebraicInterval(MultivariateDual<AlgebraicDouble> low, MultivariateDual<AlgebraicDouble> high) { 726 672 this.low = low.Clone(); 727 673 this.high = high.Clone(); … … 729 675 730 676 public AlgebraicInterval(double low, double high) { 731 this.low = new MultivariateDual<Double>(new Double(low)); 732 this.high = new MultivariateDual<Double>(new Double(high)); 733 } 734 677 this.low = new MultivariateDual<AlgebraicDouble>(new AlgebraicDouble(low)); 678 this.high = new MultivariateDual<AlgebraicDouble>(new AlgebraicDouble(high)); 679 } 680 681 public AlgebraicInterval Zero => new AlgebraicInterval(0.0, 0.0); 735 682 public AlgebraicInterval Add(AlgebraicInterval a) { 736 683 low.Add(a.low); 737 684 high.Add(a.high); 738 return this;739 }740 741 public AlgebraicInterval Assign(AlgebraicInterval a) {742 low = a.low;743 high = a.high;744 return this;745 }746 747 public AlgebraicInterval AssignCos(AlgebraicInterval a) {748 throw new NotImplementedException();749 }750 751 public AlgebraicInterval Div(AlgebraicInterval a) {752 if (a.Contains(0.0)) {753 if (a.low.Value.Value.IsAlmost(0.0) && a.high.Value.Value.IsAlmost(0.0)) {754 low = new MultivariateDual<Double>(double.NegativeInfinity);755 high = new MultivariateDual<Double>(double.PositiveInfinity);756 } else if (a.low.Value.Value.IsAlmost(0.0))757 Mul(new AlgebraicInterval(a.Clone().high.Inv(), new MultivariateDual<Double>(double.PositiveInfinity)));758 else759 Mul(new AlgebraicInterval(new MultivariateDual<Double>(double.NegativeInfinity), a.low.Clone().Inv()));760 } else {761 Mul(new AlgebraicInterval(a.high.Clone().Inv(), a.low.Clone().Inv())); // inverting leads to inverse roles of high and low762 }763 return this;764 }765 766 public AlgebraicInterval AssignExp(AlgebraicInterval a) {767 low.AssignExp(a.low);768 high.AssignExp(a.high);769 return this;770 }771 772 public AlgebraicInterval AssignIntPower(AlgebraicInterval a, int p) {773 throw new NotImplementedException();774 }775 776 public AlgebraicInterval AssignIntRoot(AlgebraicInterval a, int r) {777 throw new NotImplementedException();778 }779 780 public AlgebraicInterval AssignInv(AlgebraicInterval a) {781 low = new MultivariateDual<Double>(1.0);782 high = new MultivariateDual<Double>(1.0);783 return Div(a);784 }785 786 public AlgebraicInterval AssignLog(AlgebraicInterval a) {787 low.AssignLog(a.low);788 high.AssignLog(a.high);789 685 return this; 790 686 } … … 801 697 } 802 698 699 public AlgebraicInterval Assign(AlgebraicInterval a) { 700 low = a.low; 701 high = a.high; 702 return this; 703 } 704 705 public AlgebraicInterval AssignCos(AlgebraicInterval a) { 706 return AssignSin(a.Clone().Sub(new AlgebraicInterval(Math.PI / 2, Math.PI / 2))); 707 } 708 709 public AlgebraicInterval Div(AlgebraicInterval a) { 710 if (a.Contains(0.0)) { 711 if (a.low.Value.Value.IsAlmost(0.0) && a.high.Value.Value.IsAlmost(0.0)) { 712 low = new MultivariateDual<AlgebraicDouble>(double.NegativeInfinity); 713 high = new MultivariateDual<AlgebraicDouble>(double.PositiveInfinity); 714 } else if (a.low.Value.Value.IsAlmost(0.0)) 715 Mul(new AlgebraicInterval(a.Clone().high.Inv(), new MultivariateDual<AlgebraicDouble>(double.PositiveInfinity))); 716 else 717 Mul(new AlgebraicInterval(new MultivariateDual<AlgebraicDouble>(double.NegativeInfinity), a.low.Clone().Inv())); 718 } else { 719 Mul(new AlgebraicInterval(a.high.Clone().Inv(), a.low.Clone().Inv())); // inverting leads to inverse roles of high and low 720 } 721 return this; 722 } 723 724 public AlgebraicInterval AssignExp(AlgebraicInterval a) { 725 low.AssignExp(a.low); 726 high.AssignExp(a.high); 727 return this; 728 } 729 730 public AlgebraicInterval AssignIntPower(AlgebraicInterval a, int p) { 731 if (p == 0) { 732 // => 1 733 low = new MultivariateDual<AlgebraicDouble>(1.0); 734 high = new MultivariateDual<AlgebraicDouble>(1.0); 735 return this; 736 } 737 if (p == 1) return this; 738 739 if (p < 0) { // x^3 == 1/(x^3) 740 AssignIntPower(a, p); 741 return AssignInv(this); 742 } else { 743 // p is even => interval must be positive 744 if (p % 2 == 0) { 745 if (a.Contains(0.0)) { 746 low = new MultivariateDual<AlgebraicDouble>(0.0); 747 high = Algebraic.Max(low.Clone().IntPower(p), high.Clone().IntPower(p)); 748 } else { 749 var lowPower = low.Clone().IntPower(p); 750 var highPower = high.Clone().IntPower(p); 751 low = Algebraic.Min(lowPower, highPower); 752 high = Algebraic.Max(lowPower, highPower); 753 } 754 } else { 755 // p is uneven 756 var lowPower = low.Clone().IntPower(p); 757 var highPower = high.Clone().IntPower(p); 758 low = Algebraic.Min(lowPower, highPower); 759 high = Algebraic.Max(lowPower, highPower); 760 } 761 return this; 762 } 763 } 764 765 public AlgebraicInterval AssignIntRoot(AlgebraicInterval a, int r) { 766 if (r == 0) { low = new MultivariateDual<AlgebraicDouble>(double.NaN); high = new MultivariateDual<AlgebraicDouble>(double.NaN); return this; } 767 if (r == 1) return this; 768 if (r < 0) { 769 // x^ (1/2) = 1 / (x^(1/2)) 770 AssignIntRoot(a, r); 771 return AssignInv(this); 772 } else { 773 // root only defined for positive arguments 774 if (a.LowerBound.Value.Value < 0) { 775 low = new MultivariateDual<AlgebraicDouble>(double.NaN); 776 high = new MultivariateDual<AlgebraicDouble>(double.NaN); 777 return this; 778 } else { 779 low.AssignIntRoot(a.low, r); 780 high.AssignIntRoot(a.high, r); 781 return this; 782 } 783 } 784 } 785 786 public AlgebraicInterval AssignInv(AlgebraicInterval a) { 787 low = new MultivariateDual<AlgebraicDouble>(1.0); 788 high = new MultivariateDual<AlgebraicDouble>(1.0); 789 return Div(a); 790 } 791 792 public AlgebraicInterval AssignLog(AlgebraicInterval a) { 793 low.AssignLog(a.low); 794 high.AssignLog(a.high); 795 return this; 796 } 797 803 798 public AlgebraicInterval AssignNeg(AlgebraicInterval a) { 804 throw new NotImplementedException(); 799 low.AssignNeg(a.high); 800 high.AssignNeg(a.low); 801 return this; 805 802 } 806 803 … … 817 814 818 815 public AlgebraicInterval AssignSin(AlgebraicInterval a) { 819 throw new NotImplementedException(); 816 if (Math.Abs(a.UpperBound.Value.Value  a.LowerBound.Value.Value) >= Math.PI * 2) { 817 low = new MultivariateDual<AlgebraicDouble>(1.0); 818 high = new MultivariateDual<AlgebraicDouble>(1.0); 819 } 820 821 //divide the interval by PI/2 so that the optima lie at x element of N (0,1,2,3,4,...) 822 double Pihalf = Math.PI / 2; 823 var scaled = this.Clone().Scale(1.0 / Pihalf); 824 //move to positive scale 825 if (scaled.LowerBound.Value.Value < 0) { 826 int periodsToMove = Math.Abs((int)scaled.LowerBound.Value.Value / 4) + 1; 827 scaled.Add(new AlgebraicInterval(periodsToMove * 4, periodsToMove * 4)); 828 } 829 830 double scaledLowerBound = scaled.LowerBound.Value.Value % 4.0; 831 double scaledUpperBound = scaled.UpperBound.Value.Value % 4.0; 832 if (scaledUpperBound < scaledLowerBound) scaledUpperBound += 4.0; 833 List<double> sinValues = new List<double>(); 834 sinValues.Add(Math.Sin(scaledLowerBound * Pihalf)); 835 sinValues.Add(Math.Sin(scaledUpperBound * Pihalf)); 836 837 int startValue = (int)Math.Ceiling(scaledLowerBound); 838 while (startValue < scaledUpperBound) { 839 sinValues.Add(Math.Sin(startValue * Pihalf)); 840 startValue += 1; 841 } 842 843 low = new MultivariateDual<AlgebraicDouble>(sinValues.Min()); 844 high = new MultivariateDual<AlgebraicDouble>(sinValues.Max()); 845 return this; 820 846 } 821 847 … … 840 866 var abshigh = a.high.Clone().Abs(); 841 867 a.high.Assign(Algebraic.Max(abslow, abshigh)); 842 a.low.Assign(new MultivariateDual< Double>(0.0)); // lost gradient for lower bound868 a.low.Assign(new MultivariateDual<AlgebraicDouble>(0.0)); // lost gradient for lower bound 843 869 } else { 844 870 var abslow = a.low.Clone().Abs(); … … 851 877 852 878 public AlgebraicInterval AssignSgn(AlgebraicInterval a) { 853 throw new NotImplementedException(); 879 low.AssignSgn(a.low); 880 high.AssignSgn(a.high); 881 return this; 854 882 } 855 883 } … … 864 892 public V Derivative => dv; 865 893 894 866 895 public Dual(V v, V dv) { 867 896 this.v = v; … … 873 902 } 874 903 904 public Dual<V> Zero => new Dual<V>(Value.Zero, Derivative.Zero); 905 875 906 public Dual<V> Add(Dual<V> a) { 876 907 v.Add(a.v); … … 892 923 893 924 public Dual<V> Div(Dual<V> a) { 894 throw new NotImplementedException(); 925 Mul(a.Inv()); 926 return this; 895 927 } 896 928 … … 902 934 903 935 public Dual<V> AssignIntPower(Dual<V> a, int p) { 904 throw new NotImplementedException(); 936 v.AssignIntPower(a.v, p); 937 dv.Assign(a.dv).Scale(p).Mul(a.v.IntPower(p  1)); 938 return this; 905 939 } 906 940 907 941 public Dual<V> AssignIntRoot(Dual<V> a, int r) { 908 throw new NotImplementedException(); 942 v.AssignIntRoot(a.v, r); 943 dv.Assign(a.dv).Scale(1.0 / r).Mul(a.v.IntRoot(r  1)); 944 return this; 909 945 } 910 946 911 947 public Dual<V> AssignInv(Dual<V> a) { 912 throw new NotImplementedException(); 948 v.AssignInv(a.v); 949 dv.AssignNeg(a.dv).Mul(v).Mul(v); // (1/f(x))' =  f(x)' / f(x)^2 950 return this; 913 951 } 914 952 … … 935 973 936 974 public Dual<V> AssignNeg(Dual<V> a) { 937 throw new NotImplementedException(); 975 v.AssignNeg(a.v); 976 dv.AssignNeg(a.dv); 977 return this; 938 978 } 939 979 … … 945 985 946 986 public Dual<V> AssignSin(Dual<V> a) { 947 throw new NotImplementedException(); 987 v.AssignSin(a.v); 988 dv.Assign(a.dv).Mul(a.v.Clone().Cos()); 989 return this; 948 990 } 949 991 950 992 public Dual<V> Sub(Dual<V> a) { 951 throw new NotImplementedException(); 993 v.Sub(a.v); 994 dv.Sub(a.dv); 995 return this; 952 996 } 953 997 … … 1011 1055 } 1012 1056 1057 /// <summary> 1058 /// 1059 /// </summary> 1060 /// <param name="v">not cloned</param> 1061 /// <param name="gradient">not cloned</param> 1062 internal MultivariateDual(V v, SparseVector<object, V> gradient) { 1063 this.v = v; 1064 this.dv = gradient; 1065 } 1066 1013 1067 public MultivariateDual<V> Clone() { 1014 1068 return new MultivariateDual<V>(this); 1015 1069 } 1070 1071 1072 public MultivariateDual<V> Zero => new MultivariateDual<V>(Value.Zero, Gradient.Zero); 1073 1016 1074 1017 1075 public MultivariateDual<V> Add(MultivariateDual<V> a) { … … 1097 1155 v.AssignSgn(a.v); 1098 1156 // sign(f(x)) = 0; 1099 dv = new SparseVector<object, V>();1157 dv = a.dv.Zero; 1100 1158 return this; 1101 1159 }
Note: See TracChangeset
for help on using the changeset viewer.