@@ -492,3 +492,132 @@ mod tests_cubic_spline_interpolation {
492
492
) ;
493
493
}
494
494
}
495
+
496
+ #[ cfg( test) ]
497
+ mod tests_cubic_spline_helper_functions {
498
+ use super :: * ;
499
+
500
+ #[ test]
501
+ fn test_compute_lower_tri_and_diagonal_matrices ( ) {
502
+ let xs: Vec < f64 > = vec ! [ 1. , 2. , 3. , 4. , 5. ] ;
503
+ let ys: Vec < f64 > = vec ! [ 1. , 2. , 3. , 4. , 5. ] ;
504
+
505
+ let time_steps: Vec < f64 > = xs. windows ( 2 )
506
+ . map ( |x| ( x[ 1 ] - x[ 0 ] ) )
507
+ . collect ( ) ;
508
+
509
+ let interpolator: CubicSplineInterpolator < f64 , f64 > = CubicSplineInterpolator :: new ( xs, ys) . unwrap ( ) ;
510
+ let ( lower_triangular, diagonal) = interpolator. compute_lower_tri_and_diagonal_matrices ( & time_steps) ;
511
+
512
+ assert ! ( lower_triangular == vec![ 4.0 , 3.75 , 3.7333333333333334 ] ) ;
513
+ assert ! ( diagonal == vec![ 0.25 , 0.26666666666666666 ] ) ;
514
+ }
515
+
516
+ #[ test]
517
+ fn test_invert_lower_tri_matrix ( ) {
518
+ let xs: Vec < f64 > = vec ! [ 1. , 2. , 3. , 4. , 5. ] ;
519
+ let ys: Vec < f64 > = vec ! [ 1. , 2. , 3. , 4. , 5. ] ;
520
+
521
+ let interpolator: CubicSplineInterpolator < f64 , f64 > = CubicSplineInterpolator :: new ( xs, ys) . unwrap ( ) ;
522
+ let inv_lower_tri_matrix: Vec < Vec < f64 > > = interpolator. invert_lower_tri_matrix ( & [ 3.0 , 2.0 , 1.0 ] ) ;
523
+
524
+ assert ! (
525
+ inv_lower_tri_matrix == vec![
526
+ vec![ 1.0 ] ,
527
+ vec![ -3.0 , 1.0 ] ,
528
+ vec![ 6.0 , -2.0 , 1.0 ] ,
529
+ vec![ -6.0 , 2.0 , -1.0 , 1.0 ]
530
+ ]
531
+ )
532
+ }
533
+
534
+ #[ test]
535
+ fn test_transpose ( ) {
536
+ let xs: Vec < f64 > = vec ! [ 1. , 2. , 3. , 4. , 5. ] ;
537
+ let ys: Vec < f64 > = vec ! [ 1. , 2. , 3. , 4. , 5. ] ;
538
+
539
+ let interpolator: CubicSplineInterpolator < f64 , f64 > = CubicSplineInterpolator :: new ( xs, ys) . unwrap ( ) ;
540
+ let transposed_matrix = interpolator. transpose ( vec ! [
541
+ vec![ 1.0 ] ,
542
+ vec![ 2.0 , 3.0 ] ,
543
+ vec![ 4.0 , 5.0 , 6.0 ] ,
544
+ ] ) ;
545
+
546
+ assert ! (
547
+ transposed_matrix == vec![
548
+ vec![ 1.0 , 2.0 , 4.0 ] ,
549
+ vec![ 3.0 , 5.0 ] ,
550
+ vec![ 6.0 ] ,
551
+ ]
552
+ ) ;
553
+ }
554
+
555
+ #[ test]
556
+ fn lower_tri_transpose_times_diag_inv ( ) {
557
+ let xs: Vec < f64 > = vec ! [ 1. , 2. , 3. , 4. , 5. ] ;
558
+ let ys: Vec < f64 > = vec ! [ 1. , 2. , 3. , 4. , 5. ] ;
559
+
560
+ let interpolator: CubicSplineInterpolator < f64 , f64 > = CubicSplineInterpolator :: new ( xs, ys) . unwrap ( ) ;
561
+ let lower_tri_inverse: Vec < Vec < f64 > > = vec ! [ vec![ 1.0 ] , vec![ 2.0 , 3.0 ] , vec![ 4.0 , 5.0 , 6.0 ] ] ;
562
+ let diagonal: Vec < f64 > = vec ! [ 1.0 , 2.0 , 3.0 ] ;
563
+
564
+ assert ! (
565
+ interpolator. lower_tri_transpose_times_diag_inv( & lower_tri_inverse, & diagonal) == vec![
566
+ vec![ 1.0 , 2.0 , 4.0 ] ,
567
+ vec![ 0.5 , 2.5 ] ,
568
+ vec![ 0.3333333333333333 ]
569
+ ]
570
+ ) ;
571
+ }
572
+
573
+ #[ test]
574
+ fn test_upper_tri_times_lower_tri ( ) {
575
+ let xs: Vec < f64 > = vec ! [ 1. , 2. , 3. , 4. , 5. ] ;
576
+ let ys: Vec < f64 > = vec ! [ 1. , 2. , 3. , 4. , 5. ] ;
577
+
578
+ let interpolator: CubicSplineInterpolator < f64 , f64 > = CubicSplineInterpolator :: new ( xs, ys) . unwrap ( ) ;
579
+ let upper_tri_matrix = vec ! [
580
+ vec![ 1.0 , 2.0 , 3.0 ] ,
581
+ vec![ 4.0 , 5.0 ] ,
582
+ vec![ 6.0 ] ,
583
+ ] ;
584
+ let lower_tri_matrix_transpose = vec ! [
585
+ vec![ -1.0 , -2.0 , -3.0 ] ,
586
+ vec![ -4.0 , -5.0 ] ,
587
+ vec![ -6.0 ] ,
588
+ ] ;
589
+
590
+ let product: Vec < Vec < f64 > > = interpolator. upper_tri_times_lower_tri (
591
+ & upper_tri_matrix,
592
+ & lower_tri_matrix_transpose
593
+ ) ;
594
+
595
+ assert ! (
596
+ product == vec![
597
+ [ -14.0 , -23.0 , -18.0 ] ,
598
+ [ -23.0 , -41.0 , -30.0 ] ,
599
+ [ -18.0 , -30.0 , -36.0 ]
600
+ ]
601
+ ) ;
602
+ }
603
+
604
+ #[ test]
605
+ fn test_spline ( ) {
606
+ let xs: Vec < f64 > = vec ! [ 0. , 1. , 2. , 3. , 4. ] ;
607
+ let ys: Vec < f64 > = vec ! [ 0. , 1. , 16. , 81. , 256. ] ;
608
+
609
+ let interpolator: CubicSplineInterpolator < f64 , f64 > = CubicSplineInterpolator :: new ( xs, ys) . unwrap ( ) ;
610
+ let spline_result: f64 = interpolator. spline (
611
+ 2.5 ,
612
+ 1.0 ,
613
+ 2.0 ,
614
+ 3.0 ,
615
+ 16.0 ,
616
+ 81.0 ,
617
+ 32.75733333333333 ,
618
+ 156.85714285714286
619
+ ) ;
620
+
621
+ assert ! ( spline_result == 36.64909523809524 ) ;
622
+ }
623
+ }
0 commit comments