Changeset 16682 for branches/2994AutoDiffForIntervals/HeuristicLab.Problems.DataAnalysis.Symbolic/3.4/Interpreter/Interpreter.cs
 Timestamp:
 03/13/19 23:55:38 (5 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

branches/2994AutoDiffForIntervals/HeuristicLab.Problems.DataAnalysis.Symbolic/3.4/Interpreter/Interpreter.cs
r16674 r16682 192 192 } 193 193 194 public class VectorAutoDiffEvaluator : Interpreter< VectorDual<DoubleVector>> {194 public class VectorAutoDiffEvaluator : Interpreter<MultivariateDual<DoubleVector>> { 195 195 private const int BATCHSIZE = 128; 196 196 [ThreadStatic] … … 241 241 242 242 // TRANSPOSE into JAC 243 // here we assume that we usually calculate the gradient over _all_ parameters (otherwise it would be better to propagate only relevant vectors in evaluation) 244 for (int j = 0; j < nParams; ++j) 245 code[0].value.Gradient[j].CopyColumnTo(jac, j, rowIndex, BATCHSIZE); XXX CONTINUE HERE 243 var g = code[0].value.Gradient; 244 for (int j = 0; j < nParams; ++j) { 245 g.Elements[j].CopyColumnTo(jac, j, rowIndex, BATCHSIZE); 246 } 246 247 } 247 248 … … 249 250 Evaluate(code); 250 251 code[0].value.Value.CopyTo(fi, roundedTotal, remainingRows); 252 253 var g = code[0].value.Gradient; 251 254 for (int j = 0; j < nParams; ++j) 252 code[0].value.Gradient[j].CopyColumnTo(jac, j, roundedTotal, remainingRows);255 g.Elements[j].CopyColumnTo(jac, j, roundedTotal, remainingRows); 253 256 } 254 257 } … … 256 259 protected override void InitializeInternalInstruction(ref Instruction instruction, ISymbolicExpressionTreeNode node) { 257 260 var zero = new DoubleVector(BATCHSIZE); 258 // var g = new Vector<DoubleVector>(node2paramIdx.Count); 259 // for (int i = 0; i < BATCHSIZE; ++i) g[i] = new DoubleVector(BATCHSIZE); 260 instruction.value = new VectorDual<DoubleVector>(zero, null); 261 instruction.value = new MultivariateDual<DoubleVector>(zero); 261 262 } 262 263 263 264 protected override void InitializeTerminalInstruction(ref Instruction instruction, ConstantTreeNode constant) { 265 var g_arr = new double[BATCHSIZE]; 266 if (node2paramIdx.TryGetValue(constant, out var paramIdx)) { 267 for (int i = 0; i < BATCHSIZE; i++) g_arr[i] = 1.0; 268 var g = new DoubleVector(g_arr); 269 instruction.value = new MultivariateDual<DoubleVector>(new DoubleVector(BATCHSIZE), paramIdx, g); // only a single column for the gradient 270 } else { 271 instruction.value = new MultivariateDual<DoubleVector>(new DoubleVector(BATCHSIZE)); 272 } 273 264 274 instruction.dblVal = constant.Value; 265 266 var g_arr = new double[BATCHSIZE];267 if (node2paramIdx.ContainsKey(constant)) {268 for (int i = 0; i < BATCHSIZE; i++) g_arr[i] = 1.0;269 }270 var g = new DoubleVector(g_arr);271 272 instruction.value = new VectorDual<DoubleVector>(new DoubleVector(BATCHSIZE), new Vector<DoubleVector>(new[] { g })); // only a single column for the gradient273 275 instruction.value.Value.AssignConstant(instruction.dblVal); 274 276 } … … 286 288 if (node2paramIdx.ContainsKey(variable)) { 287 289 paramIdx = node2paramIdx[variable]; 290 var f = new DoubleVector(BATCHSIZE); 291 var g = new DoubleVector(BATCHSIZE); 292 instruction.value = new MultivariateDual<DoubleVector>(f, paramIdx, g); 293 } else { 294 var f = new DoubleVector(BATCHSIZE); 295 instruction.value = new MultivariateDual<DoubleVector>(f); 288 296 } 289 297 290 298 instruction.dblVal = variable.Weight; 291 299 instruction.data = new object[] { data, paramIdx }; 292 293 var f = new DoubleVector(BATCHSIZE);294 var g = new DoubleVector(BATCHSIZE);295 instruction.value = new VectorDual<DoubleVector>(f, new Vector<DoubleVector>(new[] { g }));296 300 } 297 301 … … 305 309 if (paramIdx >= 0) { 306 310 // update gradient with variable values 307 var g = a.value.Gradient [0];311 var g = a.value.Gradient.Elements[paramIdx]; 308 312 for (int i = rowIndex; i < rows.Length && i  rowIndex < BATCHSIZE; i++) { 309 313 g[i] = data[rows[i]]; … … 316 320 public class IntervalEvaluator : Interpreter<AlgebraicInterval> { 317 321 [ThreadStatic] 318 private Dictionary<string, AlgebraicInterval> intervals;319 320 public AlgebraicInterval Evaluate(ISymbolicExpressionTree tree, Dictionary<string, AlgebraicInterval> intervals) {322 private Dictionary<string, Interval> intervals; 323 324 public Interval Evaluate(ISymbolicExpressionTree tree, Dictionary<string, Interval> intervals) { 321 325 this.intervals = intervals; 322 326 var code = Compile(tree); 323 327 Evaluate(code); 324 return code[0].value; 325 } 326 328 return new Interval(code[0].value.LowerBound.Value.Value, code[0].value.UpperBound.Value.Value); 329 } 330 331 public Interval Evaluate(ISymbolicExpressionTree tree, Dictionary<string, Interval> intervals, ISymbolicExpressionTreeNode[] paramNodes, out double[] lowerGradient, out double[] upperGradient) { 332 this.intervals = intervals; 333 var code = Compile(tree); 334 Evaluate(code); 335 lowerGradient = new double[paramNodes.Length]; 336 upperGradient = new double[paramNodes.Length]; 337 var l = code[0].value.LowerBound; 338 var u = code[0].value.UpperBound; 339 for (int i = 0; i < paramNodes.Length; ++i) { 340 lowerGradient[i] = l.Gradient.Elements[paramNodes[i]]; 341 upperGradient[i] = u.Gradient.Elements[paramNodes[i]]; 342 } 343 return new Interval(code[0].value.LowerBound.Value.Value, code[0].value.UpperBound.Value.Value); 344 } 327 345 328 346 protected override void InitializeInternalInstruction(ref Instruction instruction, ISymbolicExpressionTreeNode node) { … … 333 351 protected override void InitializeTerminalInstruction(ref Instruction instruction, ConstantTreeNode constant) { 334 352 instruction.dblVal = constant.Value; 335 instruction.value = new AlgebraicInterval(constant.Value, constant.Value); 353 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 key 356 ); 336 357 } 337 358 338 359 protected override void InitializeTerminalInstruction(ref Instruction instruction, VariableTreeNode variable) { 339 360 instruction.dblVal = variable.Weight; 340 instruction.value = intervals[variable.VariableName]; 361 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) 364 ); 341 365 } 342 366 … … 347 371 348 372 public interface IAlgebraicType<T> { 373 T AssignAbs(T a); // set this to assign abs(a) 349 374 T Assign(T a); // assign this to same value as a (copy!) 350 375 T AssignNeg(T a); // set this to negative(a) … … 361 386 T AssignIntPower(T a, int p); 362 387 T AssignIntRoot(T a, int r); 388 T AssignSgn(T a); // set this to sign(a) 363 389 T Clone(); 364 390 } 365 391 392 public static class Algebraic { 393 public static T Abs<T>(this T a) where T : IAlgebraicType<T> { a.AssignAbs(a); return a; } 394 public static T Neg<T>(this T a) where T : IAlgebraicType<T> { a.AssignNeg(a); return a; } 395 public static T Inv<T>(this T a) where T : IAlgebraicType<T> { a.AssignInv(a); return a; } 396 public static T Log<T>(this T a) where T : IAlgebraicType<T> { a.AssignLog(a); return a; } 397 public static T Exp<T>(this T a) where T : IAlgebraicType<T> { a.AssignExp(a); return a; } 398 public static T Sin<T>(this T a) where T : IAlgebraicType<T> { a.AssignSin(a); return a; } 399 public static T Cos<T>(this T a) where T : IAlgebraicType<T> { a.AssignCos(a); return a; } 400 public static T Sgn<T>(this T a) where T : IAlgebraicType<T> { a.AssignSgn(a); return a; } 401 public static T IntPower<T>(this T a, int p) where T : IAlgebraicType<T> { a.AssignIntPower(a, p); return a; } 402 public static T IntRoot<T>(this T a, int r) where T : IAlgebraicType<T> { a.AssignIntRoot(a, r); return a; } 403 404 public static T Max<T>(T a, T b) where T : IAlgebraicType<T> { 405 // ((a + b) + abs(b  a)) / 2 406 return a.Clone().Add(b).Add(b.Clone().Sub(a).Abs()).Scale(1.0 / 2.0); 407 } 408 public static T Min<T>(T a, T b) where T : IAlgebraicType<T> { 409 // ((a + b)  abs(a  b)) / 2 410 return a.Clone().Add(b).Sub(a.Clone().Sub(b).Abs()).Scale(1.0 / 2.0); 411 } 412 } 413 414 415 // algebraic type wrapper for a double value 416 public class Double : IAlgebraicType<Double> { 417 public static implicit operator Double(double value) { return new Double(value); } 418 public static implicit operator double(Double value) { return value.Value; } 419 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 441 } 366 442 367 443 // a simple vector as an algebraic type … … 374 450 arr = new double[length]; 375 451 } 452 453 public DoubleVector() { } 376 454 377 455 /// <summary> … … 398 476 } 399 477 400 401 public DoubleVector Add(DoubleVector a) {402 for (int i = 0; i < arr.Length; ++i) {403 arr[i] += a.arr[i];404 }405 return this;406 }407 408 478 public void AssignConstant(double constantValue) { 409 479 for (int i = 0; i < arr.Length; ++i) { … … 412 482 } 413 483 414 public DoubleVector Assign(DoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = a.arr[i]; } return this; }415 public DoubleVector AssignCos(DoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = Math.Cos(a.arr[i]); } return this; }416 public DoubleVector Div(DoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] /= a.arr[i]; } return this; }417 public DoubleVector AssignExp(DoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = Math.Exp(a.arr[i]); } return this; }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; } 418 488 public DoubleVector AssignIntPower(DoubleVector a, int p) { throw new NotImplementedException(); } 419 489 public DoubleVector AssignIntRoot(DoubleVector a, int r) { throw new NotImplementedException(); } 420 public DoubleVector AssignInv(DoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = 1.0 / a.arr[i]; } return this; } 421 public DoubleVector AssignLog(DoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = Math.Log(a.arr[i]); } return this; } 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; } 492 public DoubleVector Add(DoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] += a.arr[i]; } return this; } 422 493 public DoubleVector Mul(DoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] *= a.arr[i]; } return this; } 423 public DoubleVector AssignNeg(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; } 424 495 public DoubleVector Scale(double s) { for (int i = 0; i < arr.Length; ++i) { arr[i] *= s; } return this; } 425 public DoubleVector AssignSin(DoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = Math.Sin(a.arr[i]); } 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; } 426 497 public DoubleVector Sub(DoubleVector a) { for (int i = 0; i < arr.Length; ++i) { arr[i] = a.arr[i]; } return this; } 427 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; } 428 500 public DoubleVector Clone() { 429 501 var v = new DoubleVector(this.arr.Length); … … 437 509 438 510 } 439 440 511 441 512 // vectors of algebraic types … … 506 577 public Vector<T> Scale(T s) { for (int i = 0; i < elems.Length; ++i) { elems[i].Mul(s); } return this; } 507 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 public Vector<T> AssignAbs(Vector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].AssignAbs(a.elems[i]); } return this; } 580 public Vector<T> AssignSgn(Vector<T> a) { for (int i = 0; i < elems.Length; ++i) { elems[i].AssignSgn(a.elems[i]); } return this; } 581 } 582 583 584 /// <summary> 585 /// A sparse vector of algebraic types. Elements are accessed via a key of type K 586 /// </summary> 587 /// <typeparam name="K">Key type</typeparam> 588 /// <typeparam name="T">Element type</typeparam> 589 public class SparseVector<K, T> : IAlgebraicType<SparseVector<K, T>> where T : IAlgebraicType<T> { 590 591 private Dictionary<K, T> elems; 592 593 public IReadOnlyDictionary<K, T> Elements => elems; 594 595 public SparseVector(SparseVector<K, T> original) { 596 elems = original.elems.ToDictionary(kvp => kvp.Key, kvp => kvp.Value.Clone()); 597 } 598 599 public SparseVector(K[] keys, T[] values) { 600 if (keys.Length != values.Length) throw new ArgumentException("lengths of keys and values doesn't match in SparseVector"); 601 elems = new Dictionary<K, T>(keys.Length); 602 for (int i = 0; i < keys.Length; ++i) { 603 elems.Add(keys[i], values[i]); 604 } 605 } 606 607 public SparseVector() { 608 this.elems = new Dictionary<K, T>(); 609 } 610 611 public SparseVector<K, T> Add(SparseVector<K, T> a) { 612 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()); 617 } 618 } 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 } 674 675 public SparseVector<K, T> Clone() { 676 return new SparseVector<K, T>(this); 677 } 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 } 508 713 } 509 714 510 715 511 716 public class AlgebraicInterval : IAlgebraicType<AlgebraicInterval> { 512 private doublelow;513 private doublehigh;514 515 public double LowerBound => low;516 public double UpperBound => high;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(); 517 722 518 723 public AlgebraicInterval() : this(double.NegativeInfinity, double.PositiveInfinity) { } 519 724 725 public AlgebraicInterval(MultivariateDual<Double> low, MultivariateDual<Double> high) { 726 this.low = low.Clone(); 727 this.high = high.Clone(); 728 } 729 520 730 public AlgebraicInterval(double low, double high) { 521 this.low = low;522 this.high = high;731 this.low = new MultivariateDual<Double>(new Double(low)); 732 this.high = new MultivariateDual<Double>(new Double(high)); 523 733 } 524 734 525 735 public AlgebraicInterval Add(AlgebraicInterval a) { 526 low += a.low;527 high += a.high;736 low.Add(a.low); 737 high.Add(a.high); 528 738 return this; 529 739 } … … 540 750 541 751 public AlgebraicInterval Div(AlgebraicInterval a) { 542 if (Contains(0.0)) { 543 if (a.low.IsAlmost(0.0)) 544 Mul(new AlgebraicInterval(1.0 / a.high, double.PositiveInfinity)); 545 else if (a.high.IsAlmost(0.0)) 546 Mul(new AlgebraicInterval(double.NegativeInfinity, 1.0 / a.low)); 547 else { 548 low = double.NegativeInfinity; 549 high = double.PositiveInfinity; 550 } 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 else 759 Mul(new AlgebraicInterval(new MultivariateDual<Double>(double.NegativeInfinity), a.low.Clone().Inv())); 551 760 } else { 552 Mul(new AlgebraicInterval( 1.0 / a.high, 1.0 / a.low));761 Mul(new AlgebraicInterval(a.high.Clone().Inv(), a.low.Clone().Inv())); // inverting leads to inverse roles of high and low 553 762 } 554 763 return this; … … 556 765 557 766 public AlgebraicInterval AssignExp(AlgebraicInterval a) { 558 low = Math.Exp(a.low);559 high = Math.Exp(a.high);767 low.AssignExp(a.low); 768 high.AssignExp(a.high); 560 769 return this; 561 770 } … … 570 779 571 780 public AlgebraicInterval AssignInv(AlgebraicInterval a) { 572 low = 1.0;573 high = 1.0;781 low = new MultivariateDual<Double>(1.0); 782 high = new MultivariateDual<Double>(1.0); 574 783 return Div(a); 575 784 } 576 785 577 786 public AlgebraicInterval AssignLog(AlgebraicInterval a) { 578 low = Math.Log(a.low);579 high = Math.Log(a.high);787 low.AssignLog(a.low); 788 high.AssignLog(a.high); 580 789 return this; 581 790 } 582 791 583 792 public AlgebraicInterval Mul(AlgebraicInterval a) { 584 double v1 = low * a.low;585 double v2 = low * a.high;586 double v3 = high * a.low;587 double v4 = high * a.high;588 589 low = Math.Min(Math.Min(v1, v2), Math.Min(v3, v4));590 high = Math.Max(Math.Max(v1, v2), Math.Max(v3, v4));793 var v1 = low.Clone().Mul(a.low); 794 var v2 = low.Clone().Mul(a.high); 795 var v3 = high.Clone().Mul(a.low); 796 var v4 = high.Clone().Mul(a.high); 797 798 low = Algebraic.Min(Algebraic.Min(v1, v2), Algebraic.Min(v3, v4)); 799 high = Algebraic.Max(Algebraic.Max(v1, v2), Algebraic.Max(v3, v4)); 591 800 return this; 592 801 } 593 802 594 803 public AlgebraicInterval AssignNeg(AlgebraicInterval a) { 595 low = 0; 596 high = 0; 597 return Sub(a); 804 throw new NotImplementedException(); 598 805 } 599 806 600 807 public AlgebraicInterval Scale(double s) { 808 low.Scale(s); 809 high.Scale(s); 601 810 if (s < 0) { 602 low = s * high; 603 high = s * low; 604 } else { 605 low = s * low; 606 high = s * high; 811 var t = low; 812 low = high; 813 high = t; 607 814 } 608 815 return this; … … 614 821 615 822 public AlgebraicInterval Sub(AlgebraicInterval a) { 616 low = a.high; 617 high = a.low; 823 // [x1,x2] − [y1,y2] = [x1 − y2,x2 − y1] 824 low.Sub(a.high); 825 high.Sub(a.low); 618 826 return this; 619 827 } … … 624 832 625 833 public bool Contains(double val) { 626 return low <= val && val <= high; 627 } 628 } 834 return LowerBound.Value.Value <= val && val <= UpperBound.Value.Value; 835 } 836 837 public AlgebraicInterval AssignAbs(AlgebraicInterval a) { 838 if (a.Contains(0.0)) { 839 var abslow = a.low.Clone().Abs(); 840 var abshigh = a.high.Clone().Abs(); 841 a.high.Assign(Algebraic.Max(abslow, abshigh)); 842 a.low.Assign(new MultivariateDual<Double>(0.0)); // lost gradient for lower bound 843 } else { 844 var abslow = a.low.Clone().Abs(); 845 var abshigh = a.high.Clone().Abs(); 846 a.low.Assign(Algebraic.Min(abslow, abshigh)); 847 a.high.Assign(Algebraic.Max(abslow, abshigh)); 848 } 849 return this; 850 } 851 852 public AlgebraicInterval AssignSgn(AlgebraicInterval a) { 853 throw new NotImplementedException(); 854 } 855 } 856 629 857 630 858 public class Dual<V> : IAlgebraicType<Dual<V>> … … 633 861 private V dv; 634 862 863 public V Value => v; 864 public V Derivative => dv; 865 635 866 public Dual(V v, V dv) { 636 867 this.v = v; … … 720 951 throw new NotImplementedException(); 721 952 } 722 } 723 724 public class VectorDual<V> : IAlgebraicType<VectorDual<V>> where V : IAlgebraicType<V> { 953 954 public Dual<V> AssignAbs(Dual<V> a) { 955 v.AssignAbs(a.v); 956 // abs(f(x))' = f(x)*f'(x) / f(x) 957 dv.Assign(a.dv).Mul(a.v.Clone().Sgn()); 958 return this; 959 } 960 961 public Dual<V> AssignSgn(Dual<V> a) { 962 throw new NotImplementedException(); 963 } 964 } 965 966 /// <summary> 967 /// An algebraic type which has a value as well as the partial derivatives of the value over multiple variables. 968 /// </summary> 969 /// <typeparam name="V"></typeparam> 970 public class MultivariateDual<V> : IAlgebraicType<MultivariateDual<V>> where V : IAlgebraicType<V>, new() { 725 971 private V v; 726 972 public V Value => v; 727 973 728 private Vector<V> dv; 729 public Vector<V> Gradient => dv; // columnoriented (dv[0] is the vector for the first parameter) 730 731 732 public VectorDual(V v, Vector<V> dv) { 733 this.v = v; 734 this.dv = dv; 735 } 736 737 public VectorDual<V> Clone() { 738 return new VectorDual<V>(v.Clone(), dv.Clone()); 739 } 740 741 public VectorDual<V> Add(VectorDual<V> a) { 974 private SparseVector<object, V> dv; 975 public SparseVector<object, V> Gradient => dv; // <key,value> partial derivative identified via the key 976 977 private MultivariateDual(MultivariateDual<V> orig) { 978 this.v = orig.v.Clone(); 979 this.dv = orig.dv.Clone(); 980 } 981 982 /// <summary> 983 /// Constructor without partial derivative 984 /// </summary> 985 /// <param name="v"></param> 986 public MultivariateDual(V v) { 987 this.v = v.Clone(); 988 this.dv = new SparseVector<object, V>(); 989 } 990 991 /// <summary> 992 /// Constructor for multiple partial derivatives 993 /// </summary> 994 /// <param name="v"></param> 995 /// <param name="keys"></param> 996 /// <param name="dv"></param> 997 public MultivariateDual(V v, object[] keys, V[] dv) { 998 this.v = v.Clone(); 999 this.dv = new SparseVector<object, V>(keys, dv); 1000 } 1001 1002 /// <summary> 1003 /// Constructor for a single partial derivative 1004 /// </summary> 1005 /// <param name="v"></param> 1006 /// <param name="key"></param> 1007 /// <param name="dv"></param> 1008 public MultivariateDual(V v, object key, V dv) { 1009 this.v = v.Clone(); 1010 this.dv = new SparseVector<object, V>(new[] { key }, new[] { dv }); 1011 } 1012 1013 public MultivariateDual<V> Clone() { 1014 return new MultivariateDual<V>(this); 1015 } 1016 1017 public MultivariateDual<V> Add(MultivariateDual<V> a) { 742 1018 v.Add(a.v); 743 dv. Concat(a.dv);744 return this; 745 } 746 747 public VectorDual<V> Assign(VectorDual<V> a) {1019 dv.Add(a.dv); 1020 return this; 1021 } 1022 1023 public MultivariateDual<V> Assign(MultivariateDual<V> a) { 748 1024 v.Assign(a.v); 749 dv = a.dv.Clone(); 750 return this; 751 } 752 753 public VectorDual<V> AssignCos(VectorDual<V> a) { 754 v.AssignCos(a.v); 755 V t = default(V); 756 t.AssignNeg(t.AssignSin(a.v)); 757 dv = a.dv.Clone().Scale(t); 758 return this; 759 } 760 761 public VectorDual<V> AssignExp(VectorDual<V> a) { 762 v.AssignExp(a.v); 763 dv = a.dv.Clone().Scale(v); 764 return this; 765 } 766 767 public VectorDual<V> AssignIntPower(VectorDual<V> a, int p) { 768 throw new NotImplementedException(); 769 } 770 771 public VectorDual<V> AssignIntRoot(VectorDual<V> a, int r) { 772 throw new NotImplementedException(); 773 } 774 775 public VectorDual<V> AssignInv(VectorDual<V> a) { 776 throw new NotImplementedException(); 777 } 778 779 public VectorDual<V> AssignLog(VectorDual<V> a) { 780 throw new NotImplementedException(); 781 } 782 783 public VectorDual<V> AssignNeg(VectorDual<V> a) { 784 throw new NotImplementedException(); 785 } 786 787 public VectorDual<V> AssignSin(VectorDual<V> a) { 788 throw new NotImplementedException(); 789 } 790 791 public VectorDual<V> Div(VectorDual<V> a) { 792 throw new NotImplementedException(); 793 } 794 795 public VectorDual<V> Mul(VectorDual<V> a) { 1025 dv.Assign(a.dv); 1026 return this; 1027 } 1028 1029 public MultivariateDual<V> AssignCos(MultivariateDual<V> a) { 1030 throw new NotImplementedException(); 1031 } 1032 1033 public MultivariateDual<V> AssignExp(MultivariateDual<V> a) { 1034 throw new NotImplementedException(); 1035 } 1036 1037 public MultivariateDual<V> AssignIntPower(MultivariateDual<V> a, int p) { 1038 throw new NotImplementedException(); 1039 } 1040 1041 public MultivariateDual<V> AssignIntRoot(MultivariateDual<V> a, int r) { 1042 throw new NotImplementedException(); 1043 } 1044 1045 public MultivariateDual<V> AssignInv(MultivariateDual<V> a) { 1046 throw new NotImplementedException(); 1047 } 1048 1049 public MultivariateDual<V> AssignLog(MultivariateDual<V> a) { 1050 throw new NotImplementedException(); 1051 } 1052 1053 public MultivariateDual<V> AssignNeg(MultivariateDual<V> a) { 1054 throw new NotImplementedException(); 1055 } 1056 1057 public MultivariateDual<V> AssignSin(MultivariateDual<V> a) { 1058 throw new NotImplementedException(); 1059 } 1060 1061 public MultivariateDual<V> Div(MultivariateDual<V> a) { 1062 throw new NotImplementedException(); 1063 } 1064 1065 public MultivariateDual<V> Mul(MultivariateDual<V> a) { 796 1066 // (a(x) * b(x))' = b(x)*a(x)' + b(x)'*a(x); 797 1067 798 1068 var t1 = a.dv.Clone().Scale(v); 799 1069 var t2 = dv.Clone().Scale(a.v); 800 dv.Assign(t1). Concat(t2); // TODO: CHECK ORDER1070 dv.Assign(t1).Add(t2); 801 1071 802 1072 v.Mul(a.v); … … 804 1074 } 805 1075 806 public VectorDual<V> Scale(double s) {1076 public MultivariateDual<V> Scale(double s) { 807 1077 v.Scale(s); 808 1078 dv.Scale(s); … … 810 1080 } 811 1081 812 public VectorDual<V> Sub(VectorDual<V> a) {1082 public MultivariateDual<V> Sub(MultivariateDual<V> a) { 813 1083 v.Sub(a.v); 814 dv.Concat(a.dv.AssignNeg(a.dv)); 1084 dv.Sub(a.dv); 1085 return this; 1086 } 1087 1088 public MultivariateDual<V> AssignAbs(MultivariateDual<V> a) { 1089 v.AssignAbs(a.v); 1090 // abs(f(x))' = f(x)*f'(x) / f(x) doesn't work for intervals 1091 1092 dv.Assign(a.dv).Scale(a.v.Clone().Sgn()); 1093 return this; 1094 } 1095 1096 public MultivariateDual<V> AssignSgn(MultivariateDual<V> a) { 1097 v.AssignSgn(a.v); 1098 // sign(f(x)) = 0; 1099 dv = new SparseVector<object, V>(); 815 1100 return this; 816 1101 }
Note: See TracChangeset
for help on using the changeset viewer.