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
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143

use num_traits::Float;
use crate::woltring::basis::create_basis;
use crate::woltring::prep::create_weighted_matrix;
use crate::woltring::splc::fit_spline_coefficients_with_stats;
use crate::woltring::support::{check_increasing, check_order, check_vector_length, FittingError};

pub(crate) fn fit_gcv_spline<T: Float>(knots: &Vec<T>, data: &Vec<T>, weight_factors: &Vec<T>,
                      half_order: usize, error_variance: T)
        -> Result<Vec<T>, FittingError> {
    let num_knots = knots.len();
    check_order(half_order, num_knots)?;
    check_increasing(knots)?;
    check_vector_length(data, num_knots)?;
    check_vector_length(weight_factors, num_knots)?;

    let smoothing_ratio = T::from(2.).expect("Cannot convert to type from f64");
    let tau = T::from(1.618033983).expect("Cannot convert to type from f64");
    let epsilon = T::from(1E-15).expect("Cannot convert to type from f64");
    let tolerance = T::from(1E-6).expect("Cannot convert to type from f64");
    let mut coefficients = vec![T::from(0.).expect("Cannot convert to type from f64"); num_knots];
    let mut stats = vec![T::from(0.).expect("Cannot convert to type from f64"); 6];

    // Compute design matrices and norms
    let (spline_tableau, basis_l1_norm): (Vec<T>, T) = create_basis(half_order, &knots)?;
    let (weighted_matrix, mut weighted_matrix_norm): (Vec<T>, T)
        = create_weighted_matrix(half_order, &knots, &weight_factors)?;
    weighted_matrix_norm = weighted_matrix_norm / basis_l1_norm;

    // Store temporary GCV function values
    let (mut gcv_f1, mut gcv_f2, mut gcv_f3, mut gcv_f4): (T, T, T, T);
    let (mut smoothing_1, mut smoothing_2, mut smoothing_3, mut smoothing_4): (T, T, T, T);
    let mut traced_matrix = vec![T::from(0.).expect("Cannot convert to type from f64");
                                 weighted_matrix.len()];
    // Zero variance case
    if error_variance == T::from(0.).expect("Cannot convert to type from f64") {
        smoothing_1 = T::from(0.).expect("Cannot convert to type from f64");
        let _gcv_f1 = fit_spline_coefficients_with_stats(
            half_order, &data, &weight_factors, error_variance, smoothing_1, epsilon, &spline_tableau,
            &weighted_matrix, weighted_matrix_norm, &mut coefficients, &mut stats,
            &mut traced_matrix)?;
    } else {
        let mut solved = false;
        smoothing_1 = T::from(1.).expect("Cannot convert to type from f64") / weighted_matrix_norm;
        smoothing_2 = smoothing_1 * smoothing_ratio;
        gcv_f2 = fit_spline_coefficients_with_stats(
            half_order, &data, &weight_factors, error_variance, smoothing_2, epsilon, &spline_tableau,
            &weighted_matrix, weighted_matrix_norm, &mut coefficients, &mut stats,
            &mut traced_matrix)?;
        gcv_f1 = fit_spline_coefficients_with_stats(
            half_order, &data, &weight_factors, error_variance, smoothing_1, epsilon, &spline_tableau,
            &weighted_matrix, weighted_matrix_norm, &mut coefficients, &mut stats,
            &mut traced_matrix)?;
        while gcv_f1 <= gcv_f2 && !solved {
            if stats[3] <= T::from(0.).expect("Cannot convert to type from f64") {
                solved = true;
            } else {
                smoothing_2 = smoothing_1;
                gcv_f2 = gcv_f1;
                smoothing_1 = smoothing_1 / smoothing_ratio;
                gcv_f1 = fit_spline_coefficients_with_stats(
                    half_order, &data, &weight_factors, error_variance, smoothing_1, epsilon, &spline_tableau,
                    &weighted_matrix, weighted_matrix_norm, &mut coefficients, &mut stats,
                    &mut traced_matrix)?;
            }
        }
        if !solved {
            smoothing_3 = smoothing_2 * smoothing_ratio;
            gcv_f3 = fit_spline_coefficients_with_stats(
                half_order, &data, &weight_factors, error_variance, smoothing_3, epsilon, &spline_tableau,
                &weighted_matrix, weighted_matrix_norm, &mut coefficients, &mut stats,
                &mut traced_matrix)?;
            while gcv_f3 <= gcv_f2 && !solved {
                if stats[3] >= T::from(1.).expect("Cannot convert to type from f64") {
                    solved = true;
                } else {
                    gcv_f2 = gcv_f3;
                    smoothing_3 = smoothing_3 * smoothing_ratio;
                    gcv_f3 = fit_spline_coefficients_with_stats(
                        half_order, &data, &weight_factors, error_variance, smoothing_3, epsilon, &spline_tableau,
                        &weighted_matrix, weighted_matrix_norm, &mut coefficients, &mut stats,
                        &mut traced_matrix)?;
                }
            }
            if !solved {
                smoothing_2 = smoothing_3;
                let mut alpha = (smoothing_2 - smoothing_1) / tau;
                smoothing_4 = smoothing_1 + alpha;
                smoothing_3 = smoothing_2 - alpha;
                gcv_f3 = fit_spline_coefficients_with_stats(
                    half_order, &data, &weight_factors, error_variance, smoothing_3, epsilon, &spline_tableau,
                    &weighted_matrix, weighted_matrix_norm, &mut coefficients, &mut stats,
                    &mut traced_matrix)?;
                gcv_f4 = fit_spline_coefficients_with_stats(
                    half_order, &data, &weight_factors, error_variance, smoothing_4, epsilon, &spline_tableau,
                    &weighted_matrix, weighted_matrix_norm, &mut coefficients, &mut stats,
                    &mut traced_matrix)?;
                while !solved {
                    if gcv_f3 <= gcv_f4 {
                        smoothing_2 = smoothing_4;
                        let error = (smoothing_2 - smoothing_1) / (smoothing_1 + smoothing_2);
                        if error * error + T::from(1.).expect("Cannot convert to type from f64")
                            == T::from(1.).expect("Cannot convert to type from f64") || error <= tolerance {
                            solved = true;
                        } else {
                            smoothing_4 = smoothing_3;
                            gcv_f4 = gcv_f3;
                            alpha = alpha / tau;
                            smoothing_3 = smoothing_2 - alpha;
                            gcv_f3 = fit_spline_coefficients_with_stats(
                                half_order, &data, &weight_factors, error_variance, smoothing_3, epsilon, &spline_tableau,
                                &weighted_matrix, weighted_matrix_norm, &mut coefficients, &mut stats,
                                &mut traced_matrix)?;
                        }
                    } else {
                        smoothing_1 = smoothing_3;
                        let error = (smoothing_2 - smoothing_1) / (smoothing_1 + smoothing_2);
                        if error * error + T::from(1.).expect("Cannot convert to type from f64")
                            <= T::from(1.).expect("Cannot convert to type from f64") || error <= tolerance {
                            solved = true;
                        } else {
                            smoothing_3 = smoothing_4;
                            gcv_f3 = gcv_f4;
                            alpha = alpha / tau;
                            smoothing_4 = smoothing_1 + alpha;
                            gcv_f4 = fit_spline_coefficients_with_stats(half_order, &data, &weight_factors,
                                error_variance, smoothing_4, epsilon, &spline_tableau,
                                &weighted_matrix, weighted_matrix_norm, &mut coefficients, &mut stats,
                                &mut traced_matrix)?;
                        }
                    }
                }
                smoothing_1 = T::from(0.5).expect("Cannot convert to type from f64") *
                    (smoothing_1 + smoothing_2);
                let _gcv_f1 = fit_spline_coefficients_with_stats(
                    half_order, &data, &weight_factors, error_variance, smoothing_1, epsilon, &spline_tableau,
                    &weighted_matrix, weighted_matrix_norm, &mut coefficients, &mut stats,
                    &mut traced_matrix)?;
            }
        }
    }

    Ok(coefficients)
}