gcv_spline 0.2.0

A library for fitting and evaluating GCV splines.
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126

use num_traits::Float;
use crate::woltring::search::find_knot_interval;

pub(crate) fn evaluate_spline<T: Float>(derivative_order: usize, half_order: usize, point: T, knots: &Vec<T>,
                       coefficients: &Vec<T>, knot_guess: usize) -> T {
    let num_knots = knots.len();

    // Derivatives of order >= 2 * half_order are always zero
    let order = half_order as i32 * 2 - derivative_order as i32;
    if order < 1 {
        return T::from(0.).expect("Cannot convert to type from f64");
    }

    // Search for interval value
    let knot_interval = find_knot_interval(&knots, point, knot_guess);

    // Initialize parameters and first row of B-spline coefficients tableau
    let mut tableau = vec![T::from(0.).expect("Cannot convert to type from f64"); 2 * half_order];
    let mut lower_index = knot_interval as i32 + 1;
    let upper_index = knot_interval as i32 + half_order as i32 * 2;
    let mut inner_index = num_knots - 2 * half_order;

    for index in lower_index ..= upper_index {
        if index >= half_order as i32 + 1 && index <= num_knots as i32 + half_order as i32 {
            tableau[(index - knot_interval as i32 - 1) as usize] =
                coefficients[(index - half_order as i32 - 1) as usize];
        } else {
            tableau[(index - knot_interval as i32 - 1) as usize] =
                T::from(0.).expect("Cannot convert to type from f64");
        }
    }

    // The following loop computes differences of the B-spline coefficients. If the value of the
    // spline is required, differencing is not necessary.
    if derivative_order > 0 {
        lower_index -= half_order as i32 * 2;
        let index_bound = half_order as i32 * 2 - knot_interval as i32;
        for der_index in 1 ..= derivative_order {
            lower_index += 1;
            inner_index += 1;
            let idx_1 = 1.max(lower_index) as usize;
            let idx_2 = knot_interval.min(inner_index);
            let mut idx = idx_2 + 1;
            if idx_1 <= idx_2 {
                for _ in idx_1 ..= idx_2 {
                    idx -= 1;
                    let work_idx = (index_bound + idx as i32) as usize;
                    tableau[work_idx - 1] = (tableau[work_idx - 1] - tableau[work_idx - 2])
                        / (knots[idx + half_order * 2 - der_index - 1] - knots[idx - 1]);
                }
            }
            if lower_index < 1 {
                idx = (index_bound + 1) as usize;
                if der_index as i32 + 1 <= index_bound {
                    for _ in der_index as i32 + 1 ..= index_bound {
                        idx -= 1;
                        tableau[idx - 1] = -T::from(1.).expect("Cannot convert to type from f64")
                            * tableau[idx - 2];
                    }
                }
            }
        }
        for idx in 1 ..= order {
            tableau[idx as usize - 1] = tableau[idx as usize + derivative_order - 1];
        }
    }

    let mut solution;

    // Compute lower half of the evaluation tableau
    if order - 1 >= 1 { // Tableau is ready if derivative order == 2 * half_order - 1
        for idx in 1 ..= order as usize - 1 {
            let order_idx = (num_knots as i32 - order) as usize + idx;
            let mut working_idx = order as usize;
            let mut knot_idx = knot_interval;

            // Right hand splines
            if knot_interval >= order_idx + 1 {
                for _ in order_idx + 1 ..= knot_interval {
                    tableau[working_idx - 1] = tableau[working_idx - 2] + (point - knots[knot_idx - 1])
                        * tableau[working_idx - 1];
                    knot_idx -= 1;
                    working_idx -= 1;
                }
            }

            // Middle B-splines
            let idx_1 = 1.max(knot_interval as i32 - order + 1 + idx as i32);
            let idx_2 = knot_interval.min(order_idx) as i32;
            if idx_1 <= idx_2 {
                for _ in idx_1 ..= idx_2 {
                    let working_knot = knots[(knot_idx as i32 + order - idx as i32 - 1) as usize];
                    solution = tableau[working_idx - 1];
                    tableau[working_idx - 1] = solution + (working_knot - point)
                        * (tableau[working_idx - 2] - solution) / (working_knot - knots[knot_idx - 1]);
                    working_idx -= 1;
                    knot_idx -= 1;
                }
            }

            // Left hand B-splines
            if knot_interval as i32 - order + 1 + idx as i32 <= 0 {
                knot_idx = (order - idx as i32) as usize;
                for _ in 1 ..= 1 - (knot_interval as i32 - order + 1 + idx as i32) {
                    tableau[working_idx - 1] = tableau[working_idx - 1]
                        + (knots[knot_idx - 1] - point) * tableau[working_idx - 2];
                    knot_idx -= 1;
                    working_idx -= 1;
                }
            }
        }
    }

    // Compute the return value
    solution = tableau[order as usize - 1];

    // Multiply with factorial of derivative order > 0
    if derivative_order > 0 {
        for idx in order ..= 2 * half_order as i32 - 1 {
            let index = idx as usize;
            solution = solution * T::from(index).expect("Cannot convert to type from usize");
        }
    }

    solution
}