#![forbid(unsafe_code)]
pub mod bspline;
pub mod cyclic;
pub mod error;
pub mod fourier;
pub mod ispline;
mod local;
pub mod monotone;
pub mod mspline;
pub mod natural;
pub mod open_uniform;
pub mod order;
pub mod penalty;
pub mod periodic;
pub mod row_basis;
pub mod tensor;
pub mod truncated_power;
pub use bspline::{BSplineBasis, pspline_design};
pub use cyclic::{CyclicSplineDesign, CyclicSplineSpec};
pub use error::{FourierError, SplineError};
pub use fourier::FourierDesign;
pub use ispline::{ISplineBasis, ISplineDesign};
pub use monotone::{MonotoneDirection, MonotoneISplineDesign};
pub use mspline::{MSplineBasis, MSplineDesign};
pub use natural::{NaturalCubicSplineBasis, NaturalCubicSplineDesign};
pub use open_uniform::{OpenUniformSplineBasis, OpenUniformSplineDesign};
pub use order::SplineOrder;
pub use penalty::{
CyclicDifferencePenalty, DifferencePenalty, EdgeMonotonicPenalty,
PreparedCyclicDifferencePenalty, PreparedDifferencePenalty, SlopeLimitPenalty,
};
pub use periodic::{PeriodicSplineDesign, PeriodicSplineSpec};
pub use row_basis::SplineRowBasis;
pub use tensor::TensorSplineDesign;
pub use truncated_power::{TruncatedPowerBasis, TruncatedPowerDesign};
pub mod prelude {
pub use crate::{
BSplineBasis, CyclicDifferencePenalty, CyclicSplineDesign, CyclicSplineSpec,
DifferencePenalty, EdgeMonotonicPenalty, FourierDesign, FourierError, ISplineBasis,
ISplineDesign, MSplineBasis, MSplineDesign, MonotoneDirection, MonotoneISplineDesign,
NaturalCubicSplineBasis, NaturalCubicSplineDesign, OpenUniformSplineBasis,
OpenUniformSplineDesign, PeriodicSplineDesign, PeriodicSplineSpec,
PreparedCyclicDifferencePenalty, PreparedDifferencePenalty, SlopeLimitPenalty, SplineError,
SplineOrder, SplineRowBasis, TensorSplineDesign, TruncatedPowerBasis, TruncatedPowerDesign,
pspline_design,
};
}
#[cfg(test)]
mod tests {
use approx::assert_relative_eq;
use gamlss_core::{Penalty, PredictorBlock};
use super::{
BSplineBasis, CyclicDifferencePenalty, CyclicSplineDesign, CyclicSplineSpec,
DifferencePenalty, EdgeMonotonicPenalty, FourierDesign, FourierError, ISplineBasis,
MSplineBasis, MonotoneDirection, MonotoneISplineDesign, NaturalCubicSplineBasis,
OpenUniformSplineBasis, OpenUniformSplineDesign, PeriodicSplineDesign, PeriodicSplineSpec,
PreparedCyclicDifferencePenalty, PreparedDifferencePenalty, SlopeLimitPenalty, SplineError,
SplineOrder, TensorSplineDesign, TruncatedPowerBasis,
};
#[test]
fn open_uniform_bspline_partitions_unity_inside_range() {
let x = [0.0, 0.25, 0.5, 0.75, 1.0];
let basis = BSplineBasis::open_uniform_from_data(&x, 6, 3).unwrap();
for value in x {
let sum = basis.evaluate(value).iter().sum::<f64>();
assert_relative_eq!(sum, 1.0, epsilon = 1.0e-12);
}
}
#[test]
fn difference_penalty_gradient_matches_finite_difference() {
let penalty = DifferencePenalty::new(0.7, 2);
let beta = vec![0.2, -0.4, 0.9, 1.1, -0.3];
let eps = 1.0e-6;
let mut grad = vec![0.0; beta.len()];
penalty.add_gradient(&beta, &mut grad);
for index in 0..beta.len() {
let mut plus = beta.clone();
plus[index] += eps;
let mut minus = beta.clone();
minus[index] -= eps;
let finite_difference = (penalty.value(&plus) - penalty.value(&minus)) / (2.0 * eps);
assert_relative_eq!(grad[index], finite_difference, epsilon = 1.0e-6);
}
}
#[test]
fn prepared_difference_penalty_matches_unprepared() {
let beta = [0.2, -0.4, 0.9, 1.1, -0.3];
for order in 0..=2 {
let unprepared = DifferencePenalty::new(0.7, order);
let prepared = PreparedDifferencePenalty::new(0.7, order);
let mut unprepared_grad = vec![0.0; beta.len()];
let mut prepared_grad = vec![0.0; beta.len()];
unprepared.add_gradient(&beta, &mut unprepared_grad);
prepared.add_gradient(&beta, &mut prepared_grad);
assert_relative_eq!(
prepared.value(&beta),
unprepared.value(&beta),
epsilon = 1.0e-12
);
for (prepared, unprepared) in prepared_grad.iter().zip(&unprepared_grad) {
assert_relative_eq!(prepared, unprepared, epsilon = 1.0e-12);
}
}
}
#[test]
fn open_uniform_basis_visitor_matches_design_row_basis() {
let x = [-0.25, 0.0, 0.2, 0.75, 1.25];
let basis = OpenUniformSplineBasis::new(0.0, 1.0, 6, SplineOrder::Cubic).unwrap();
let design = basis.design(&x).unwrap();
for (row, value) in x.iter().copied().enumerate() {
let mut from_design = Vec::new();
let mut from_basis = Vec::new();
super::SplineRowBasis::for_each_row_basis(&design, row, |index, weight| {
from_design.push((index, weight));
});
basis
.for_each_value_basis(value, |index, weight| {
from_basis.push((index, weight));
})
.unwrap();
assert_eq!(from_basis.len(), from_design.len());
for ((basis_index, basis_weight), (design_index, design_weight)) in
from_basis.iter().zip(&from_design)
{
assert_eq!(basis_index, design_index);
assert_relative_eq!(basis_weight, design_weight, epsilon = 1.0e-12);
}
}
}
#[test]
fn open_uniform_basis_visitor_rejects_non_finite_input() {
let basis = OpenUniformSplineBasis::new(0.0, 1.0, 6, SplineOrder::Cubic).unwrap();
let err = basis.for_each_value_basis(f64::NAN, |_, _| {}).unwrap_err();
assert_eq!(err, SplineError::NonFiniteValue);
}
#[test]
fn allocating_and_buffer_basis_evaluation_match() {
let x = [0.0, 0.25, 0.5, 0.75, 1.0];
let b = BSplineBasis::open_uniform_from_data(&x, 6, 3).unwrap();
let m = MSplineBasis::open_uniform_from_data(&x, 6, 3).unwrap();
let i = ISplineBasis::open_uniform_from_data(&x, 6, 3).unwrap();
let natural = NaturalCubicSplineBasis::new(vec![0.0, 0.5, 1.0]).unwrap();
let mut buffer = vec![0.0; b.n_basis()];
b.evaluate_into(0.4, &mut buffer);
assert_eq!(buffer, b.evaluate(0.4));
buffer.resize(m.n_basis(), 0.0);
m.evaluate_into(0.4, &mut buffer);
assert_eq!(buffer, m.evaluate(0.4));
buffer.resize(i.n_basis(), 0.0);
i.evaluate_into(0.4, &mut buffer);
assert_eq!(buffer, i.evaluate(0.4));
buffer.resize(natural.n_basis(), 0.0);
natural.evaluate_into(0.4, &mut buffer);
assert_eq!(buffer, natural.evaluate(0.4));
natural.evaluate_derivative_into(0.4, &mut buffer);
assert_eq!(buffer, natural.evaluate_derivative(0.4));
}
#[test]
fn spline_row_basis_visitors_match_allocating_evaluate_values() {
let x = [0.0, 0.25, 0.5, 0.75, 1.0];
let m = MSplineBasis::open_uniform_from_data(&x, 6, 3).unwrap();
let i = ISplineBasis::open_uniform_from_data(&x, 6, 3).unwrap();
let natural = NaturalCubicSplineBasis::new(vec![0.0, 0.5, 1.0]).unwrap();
assert_row_basis_matches_evaluate(&m.design(&x).unwrap(), |row| m.evaluate(x[row]));
assert_row_basis_matches_evaluate(&i.design(&x).unwrap(), |row| i.evaluate(x[row]));
assert_row_basis_matches_evaluate(&natural.design(&x).unwrap(), |row| {
natural.evaluate(x[row])
});
}
#[test]
fn cyclic_spline_design_wraps_and_partitions_unity() {
let design =
CyclicSplineDesign::new(&[0.0, 0.25, 1.0, -0.25], 8, SplineOrder::Cubic).unwrap();
let beta = vec![1.0; design.nparams()];
for row in 0..design.nrows() {
assert_relative_eq!(design.eta_row(row, &beta), 1.0, epsilon = 1.0e-12);
}
let ramp = (0..8).map(|value| value as f64).collect::<Vec<_>>();
assert_relative_eq!(design.eta_row(0, &ramp), design.eta_row(2, &ramp));
}
#[test]
fn fourier_design_evaluates_harmonics_without_materialized_matrix() {
let design = FourierDesign::new(&[0.0, 0.25, 0.5, 1.25], 1.0, 2, true).unwrap();
let beta = [0.5, 1.0, 2.0, -0.25, 0.75];
assert_eq!(design.nrows(), 4);
assert_eq!(design.nparams(), 5);
assert_eq!(design.order(), 2);
assert_relative_eq!(design.period(), 1.0, epsilon = 1.0e-12);
assert!(design.include_intercept());
for row in 0..design.nrows() {
let x = design.x()[row];
let phase1 = std::f64::consts::TAU * x;
let phase2 = 2.0 * std::f64::consts::TAU * x;
let expected = beta[0]
+ beta[1] * phase1.sin()
+ beta[2] * phase1.cos()
+ beta[3] * phase2.sin()
+ beta[4] * phase2.cos();
assert_relative_eq!(design.eta_row(row, &beta), expected, epsilon = 1.0e-12);
}
assert_relative_eq!(
design.eta_row(1, &beta),
design.eta_row(3, &beta),
epsilon = 1.0e-12
);
}
#[test]
fn fourier_design_without_intercept_uses_two_coefficients_per_harmonic() {
let design = FourierDesign::new(&[1.0], 4.0, 1, false).unwrap();
let beta = [2.0, 3.0];
let phase = std::f64::consts::TAU * 1.0 / 4.0;
assert_eq!(design.nparams(), 2);
assert_relative_eq!(
design.eta_row(0, &beta),
beta[0] * phase.sin() + beta[1] * phase.cos(),
epsilon = 1.0e-12
);
}
#[test]
fn fourier_row_basis_matches_eta_unit_coefficients() {
let design = FourierDesign::new(&[0.0, 0.2, 0.7], 1.0, 3, true).unwrap();
for row in 0..design.nrows() {
let mut basis = vec![0.0; design.nparams()];
crate::SplineRowBasis::for_each_row_basis(&design, row, |index, weight| {
basis[index] = weight;
});
for index in 0..design.nparams() {
let mut beta = vec![0.0; design.nparams()];
beta[index] = 1.0;
assert_relative_eq!(basis[index], design.eta_row(row, &beta), epsilon = 1.0e-12);
}
}
}
#[test]
fn fourier_design_rejects_invalid_inputs() {
assert_eq!(
FourierDesign::new(&[f64::NAN], 1.0, 1, true).unwrap_err(),
FourierError::NonFiniteValue
);
assert_eq!(
FourierDesign::new(&[0.0], 0.0, 1, true).unwrap_err(),
FourierError::InvalidPeriod
);
assert_eq!(
FourierDesign::new(&[0.0], 1.0, 0, true).unwrap_err(),
FourierError::InvalidOrder
);
assert_eq!(
FourierDesign::new(&[0.0], 1.0, usize::MAX, true).unwrap_err(),
FourierError::CoefficientOverflow
);
}
#[test]
fn fourier_design_gradient_matches_finite_difference() {
let design = FourierDesign::new(&[0.0, 0.2, 0.7, 1.4], 1.0, 3, true).unwrap();
let beta = vec![0.3, -0.5, 1.2, 0.4, -0.8, 0.9, 0.1];
let scores = vec![0.7, -1.1, 0.2, 1.4];
let eps = 1.0e-6;
let mut grad = vec![0.0; design.nparams()];
design.add_gradient(&scores, &beta, &mut grad);
for index in 0..beta.len() {
let mut plus = beta.clone();
plus[index] += eps;
let mut minus = beta.clone();
minus[index] -= eps;
let objective = |candidate: &[f64]| {
(0..design.nrows())
.map(|row| scores[row] * design.eta_row(row, candidate))
.sum::<f64>()
};
let finite_difference = (objective(&plus) - objective(&minus)) / (2.0 * eps);
assert_relative_eq!(grad[index], finite_difference, epsilon = 1.0e-6);
}
}
#[test]
fn open_uniform_spline_design_handles_boundaries_and_extrapolation() {
let design = OpenUniformSplineDesign::with_range(
&[-0.5, 0.0, 0.5, 1.0, 1.5],
0.0,
1.0,
6,
SplineOrder::Cubic,
)
.unwrap();
let beta = vec![1.0; design.nparams()];
for row in 0..design.nrows() {
assert_relative_eq!(design.eta_row(row, &beta), 1.0, epsilon = 1.0e-12);
}
}
#[test]
fn open_uniform_spline_basis_reuses_training_range_for_new_data() {
let train = [0.0, 0.5, 1.0];
let basis = OpenUniformSplineBasis::from_data(&train, 6, SplineOrder::Cubic).unwrap();
let train_design = basis.design(&train).unwrap();
let new_design = basis.design(&[-0.25, 0.25, 1.25]).unwrap();
assert_eq!(train_design.basis(), basis);
assert_eq!(new_design.nparams(), train_design.nparams());
assert_relative_eq!(basis.min(), 0.0, epsilon = 1.0e-12);
assert_relative_eq!(basis.max(), 1.0, epsilon = 1.0e-12);
let beta = vec![1.0; new_design.nparams()];
for row in 0..new_design.nrows() {
assert_relative_eq!(new_design.eta_row(row, &beta), 1.0, epsilon = 1.0e-12);
}
}
#[test]
fn cyclic_spline_spec_reuses_basis_shape_for_new_data() {
let spec = CyclicSplineSpec::new(6, SplineOrder::Cubic).unwrap();
let train_design = spec.design(&[0.0, 0.5]).unwrap();
let new_design = spec.design(&[0.25, 1.25]).unwrap();
assert_eq!(train_design.spec(), spec);
assert_eq!(new_design.nparams(), train_design.nparams());
let beta = vec![1.0; new_design.nparams()];
for row in 0..new_design.nrows() {
assert_relative_eq!(new_design.eta_row(row, &beta), 1.0, epsilon = 1.0e-12);
}
}
#[test]
fn cyclic_difference_penalty_gradient_matches_finite_difference() {
let penalty = CyclicDifferencePenalty::new(0.7, 2);
let beta = vec![0.2, -0.4, 0.9, 1.1, -0.3];
assert_penalty_gradient_matches_finite_difference(&penalty, &beta);
}
#[test]
fn prepared_cyclic_difference_penalty_matches_unprepared() {
let beta = [0.2, -0.4, 0.9, 1.1, -0.3];
for order in 0..=2 {
let unprepared = CyclicDifferencePenalty::new(0.7, order);
let prepared = PreparedCyclicDifferencePenalty::new(0.7, order);
let mut unprepared_grad = vec![0.0; beta.len()];
let mut prepared_grad = vec![0.0; beta.len()];
unprepared.add_gradient(&beta, &mut unprepared_grad);
prepared.add_gradient(&beta, &mut prepared_grad);
assert_relative_eq!(
prepared.value(&beta),
unprepared.value(&beta),
epsilon = 1.0e-12
);
for (prepared, unprepared) in prepared_grad.iter().zip(&unprepared_grad) {
assert_relative_eq!(prepared, unprepared, epsilon = 1.0e-12);
}
}
}
#[test]
fn edge_monotonic_penalty_gradient_matches_finite_difference() {
let penalty = EdgeMonotonicPenalty::new(3.0);
let beta = vec![0.2, 0.8, 0.4, -0.1];
assert_penalty_gradient_matches_finite_difference(&penalty, &beta);
}
#[test]
fn slope_limit_penalty_gradient_matches_finite_difference() {
let penalty = SlopeLimitPenalty::new(5.0, 2.0, Some(0.4), Some(0.3));
let beta = vec![0.6, 0.1, -0.1, 0.4];
assert_penalty_gradient_matches_finite_difference(&penalty, &beta);
}
#[test]
fn slope_limit_penalty_ignores_invalid_inputs() {
let cases = [
SlopeLimitPenalty::new(f64::NAN, 2.0, Some(0.4), Some(0.3)),
SlopeLimitPenalty::new(5.0, f64::NAN, Some(0.4), Some(0.3)),
SlopeLimitPenalty::new(5.0, 2.0, Some(f64::NAN), Some(-0.3)),
];
let beta = vec![0.6, 0.1, -0.1, 0.4];
for penalty in cases {
let mut grad = vec![0.0; beta.len()];
assert_eq!(penalty.value(&beta), 0.0);
penalty.add_gradient(&beta, &mut grad);
assert_eq!(grad, vec![0.0; beta.len()]);
}
}
#[test]
fn natural_cubic_validates_knots_and_has_linear_edges() {
assert_eq!(
NaturalCubicSplineBasis::new(vec![0.0]).unwrap_err(),
SplineError::NotEnoughKnots { min: 2 }
);
assert_eq!(
NaturalCubicSplineBasis::new(vec![0.0, f64::NAN]).unwrap_err(),
SplineError::InvalidKnots
);
let basis = NaturalCubicSplineBasis::new(vec![0.0, 0.5, 1.0]).unwrap();
let left = basis.evaluate_derivative(-0.25);
let at_left = basis.evaluate_derivative(0.0);
let right = basis.evaluate_derivative(1.25);
let at_right = basis.evaluate_derivative(1.0);
for index in 0..basis.n_basis() {
assert_relative_eq!(left[index], at_left[index], epsilon = 1.0e-12);
assert_relative_eq!(right[index], at_right[index], epsilon = 1.0e-12);
}
}
#[test]
fn natural_cubic_design_gradient_matches_finite_difference() {
let design = NaturalCubicSplineBasis::new(vec![0.0, 0.4, 0.7, 1.0])
.unwrap()
.design(&[0.0, 0.2, 0.8, 1.1])
.unwrap();
assert_predictor_gradient_matches_finite_difference(
&design,
&[0.4, -0.2, 0.7, 1.1],
&[0.3, -0.8, 0.5, 1.0],
);
}
#[test]
fn truncated_power_basis_uses_documented_column_order() {
let basis = TruncatedPowerBasis::new(vec![0.25, 0.75], SplineOrder::Cubic, true).unwrap();
let values = basis.evaluate(1.25);
assert_eq!(basis.n_basis(), 6);
assert_eq!(basis.order(), SplineOrder::Cubic);
assert!(basis.include_intercept());
assert_relative_eq!(values[0], 1.0, epsilon = 1.0e-12);
assert_relative_eq!(values[1], 1.25, epsilon = 1.0e-12);
assert_relative_eq!(values[2], 1.25_f64.powi(2), epsilon = 1.0e-12);
assert_relative_eq!(values[3], 1.25_f64.powi(3), epsilon = 1.0e-12);
assert_relative_eq!(values[4], 1.0, epsilon = 1.0e-12);
assert_relative_eq!(values[5], 0.5_f64.powi(3), epsilon = 1.0e-12);
}
#[test]
fn truncated_power_basis_validates_inputs_and_builds_uniform_knots() {
assert_eq!(
TruncatedPowerBasis::new(vec![0.0, 0.0], SplineOrder::Linear, true).unwrap_err(),
SplineError::InvalidKnots
);
assert_eq!(
TruncatedPowerBasis::uniform_from_data(&[0.0, f64::NAN], 1, SplineOrder::Linear, true)
.unwrap_err(),
SplineError::NonFiniteValue
);
let basis =
TruncatedPowerBasis::uniform_from_data(&[0.0, 0.5, 1.0], 2, SplineOrder::Linear, false)
.unwrap();
assert_relative_eq!(basis.knots()[0], 1.0 / 3.0, epsilon = 1.0e-12);
assert_relative_eq!(basis.knots()[1], 2.0 / 3.0, epsilon = 1.0e-12);
assert_eq!(basis.n_basis(), 3);
}
#[test]
fn truncated_power_buffer_evaluation_clears_stale_values() {
let basis = TruncatedPowerBasis::new(vec![0.25, 0.75], SplineOrder::Cubic, true).unwrap();
let mut values = vec![42.0; basis.n_basis()];
basis.evaluate_into(0.0, &mut values);
assert_relative_eq!(values[0], 1.0, epsilon = 1.0e-12);
for value in &values[1..] {
assert_relative_eq!(*value, 0.0, epsilon = 1.0e-12);
}
}
#[test]
fn truncated_power_design_gradient_matches_finite_difference() {
let design = TruncatedPowerBasis::new(vec![0.25, 0.75], SplineOrder::Cubic, true)
.unwrap()
.design(&[-0.2, 0.1, 0.5, 1.2])
.unwrap();
assert_predictor_gradient_matches_finite_difference(
&design,
&[0.2, -0.4, 0.7, 0.1, -0.3, 0.9],
&[0.3, -0.8, 0.5, 1.0],
);
assert_row_basis_matches_evaluate(&design, |row| design.basis().evaluate(design.x()[row]));
let mut visited = Vec::new();
crate::SplineRowBasis::for_each_row_basis(&design, 1, |index, weight| {
visited.push((index, weight));
});
assert_eq!(visited.len(), 4);
for ((actual_index, actual_weight), (expected_index, expected_weight)) in visited
.iter()
.zip([(0, 1.0), (1, 0.1), (2, 0.01), (3, 0.001)])
{
assert_eq!(*actual_index, expected_index);
assert_relative_eq!(*actual_weight, expected_weight, epsilon = 1.0e-12);
}
}
#[test]
fn spline_derivatives_match_finite_difference() {
let x = [0.1, 0.25, 0.6, 0.9];
let beta = vec![0.1, -0.4, 0.7, 0.2, -0.1, 0.5];
assert_eta_derivative_matches_coordinate_difference(
|points| {
OpenUniformSplineDesign::with_range(points, 0.0, 1.0, 6, SplineOrder::Cubic)
.unwrap()
},
|design, row, beta| design.eta_derivative_row(row, beta),
&x,
&beta,
);
let cyclic = CyclicSplineDesign::new(&x, 6, SplineOrder::Cubic).unwrap();
assert_eta_derivative_matches_coordinate_difference(
|points| CyclicSplineDesign::new(points, 6, SplineOrder::Cubic).unwrap(),
|design, row, beta| design.eta_derivative_row(row, beta),
&x,
&beta,
);
let natural = NaturalCubicSplineBasis::new(vec![0.0, 0.4, 0.7, 1.0])
.unwrap()
.design(&x)
.unwrap();
let natural_beta = vec![0.1, -0.4, 0.7, 0.2];
assert_eta_derivative_matches_coordinate_difference(
|points| {
NaturalCubicSplineBasis::new(vec![0.0, 0.4, 0.7, 1.0])
.unwrap()
.design(points)
.unwrap()
},
|design, row, beta| design.eta_derivative_row(row, beta),
&x,
&natural_beta,
);
let truncated = TruncatedPowerBasis::new(vec![0.3, 0.7], SplineOrder::Cubic, true)
.unwrap()
.design(&x)
.unwrap();
let truncated_beta = vec![0.1, -0.4, 0.7, 0.2, -0.1, 0.5];
assert_eta_derivative_matches_coordinate_difference(
|points| {
TruncatedPowerBasis::new(vec![0.3, 0.7], SplineOrder::Cubic, true)
.unwrap()
.design(points)
.unwrap()
},
|design, row, beta| design.eta_derivative_row(row, beta),
&x,
&truncated_beta,
);
assert!(cyclic.eta_derivative_row(0, &beta).is_finite());
assert!(natural.eta_derivative_row(0, &natural_beta).is_finite());
assert!(truncated.eta_derivative_row(0, &truncated_beta).is_finite());
}
#[test]
fn periodic_spline_wraps_physical_coordinates_and_rejects_invalid_period() {
assert_eq!(
PeriodicSplineSpec::new(6, SplineOrder::Cubic, 0.0, 0.0).unwrap_err(),
SplineError::InvalidPeriod
);
let design =
PeriodicSplineDesign::new(&[0.0, 12.0, 3.0, 15.0], 6, SplineOrder::Cubic, 12.0, 0.0)
.unwrap();
let beta = vec![0.2, -0.1, 0.7, 0.4, -0.3, 0.9];
assert_relative_eq!(
design.eta_row(0, &beta),
design.eta_row(1, &beta),
epsilon = 1.0e-12
);
assert_relative_eq!(
design.eta_row(2, &beta),
design.eta_row(3, &beta),
epsilon = 1.0e-12
);
assert!(design.eta_derivative_row(2, &beta).is_finite());
}
#[test]
fn tensor_spline_matches_rowwise_kronecker_and_gradient() {
let x = [0.0, 0.3, 0.8];
let left =
OpenUniformSplineDesign::with_range(&x, 0.0, 1.0, 4, SplineOrder::Cubic).unwrap();
let right = CyclicSplineDesign::new(&x, 5, SplineOrder::Cubic).unwrap();
let tensor = TensorSplineDesign::new(left.clone(), right.clone()).unwrap();
let beta = (0..tensor.nparams())
.map(|index| index as f64 / 10.0)
.collect::<Vec<_>>();
let row = 1;
let mut expected = 0.0;
for left_index in 0..left.nparams() {
let mut left_unit = vec![0.0; left.nparams()];
left_unit[left_index] = 1.0;
let left_weight = left.eta_row(row, &left_unit);
for right_index in 0..right.nparams() {
let mut right_unit = vec![0.0; right.nparams()];
right_unit[right_index] = 1.0;
let right_weight = right.eta_row(row, &right_unit);
expected +=
beta[left_index * right.nparams() + right_index] * left_weight * right_weight;
}
}
assert_relative_eq!(tensor.eta_row(row, &beta), expected, epsilon = 1.0e-12);
assert_predictor_gradient_matches_finite_difference(&tensor, &beta, &[0.5, -0.2, 0.9]);
let short = CyclicSplineDesign::new(&[0.0], 5, SplineOrder::Cubic).unwrap();
assert_eq!(
TensorSplineDesign::new(left, short).unwrap_err(),
SplineError::RowMismatch {
expected: 3,
actual: 1
}
);
}
#[test]
fn m_and_i_splines_are_nonnegative_monotone_and_differentiable() {
let x = [0.0, 0.2, 0.5, 0.8, 1.0];
let m = MSplineBasis::open_uniform_from_data(&x, 6, 3).unwrap();
let i = ISplineBasis::open_uniform_from_data(&x, 6, 3).unwrap();
for value in [0.1, 0.35, 0.7] {
for weight in m.evaluate(value) {
assert!(weight >= -1.0e-12);
}
}
for basis_index in 0..i.n_basis() {
let a = i.evaluate(0.25)[basis_index];
let b = i.evaluate(0.75)[basis_index];
assert!(a <= b + 1.0e-12);
assert_relative_eq!(i.evaluate(-0.1)[basis_index], 0.0, epsilon = 1.0e-12);
assert_relative_eq!(i.evaluate(1.1)[basis_index], 1.0, epsilon = 1.0e-9);
}
let point = 0.43;
let eps = 1.0e-6;
for basis_index in 0..i.n_basis() {
let finite = (i.evaluate(point + eps)[basis_index]
- i.evaluate(point - eps)[basis_index])
/ (2.0 * eps);
assert_relative_eq!(
i.evaluate_derivative(point)[basis_index],
finite,
epsilon = 1.0e-5
);
}
}
#[test]
fn monotone_i_spline_is_hard_monotone_and_gradient_matches_finite_difference() {
let x = [0.0, 0.2, 0.5, 0.8, 1.0];
let basis = ISplineBasis::open_uniform_from_data(&x, 5, 2).unwrap();
let design =
MonotoneISplineDesign::new(&x, basis.clone(), MonotoneDirection::Increasing).unwrap();
let beta = vec![0.1, -1.0, 0.2, -0.4, 0.7, 1.0];
for row in 1..design.nrows() {
assert!(design.eta_row(row - 1, &beta) <= design.eta_row(row, &beta) + 1.0e-12);
assert!(design.eta_derivative_row(row, &beta) >= -1.0e-10);
}
assert_predictor_gradient_matches_finite_difference(
&design,
&beta,
&[0.2, -0.3, 0.5, 0.7, -0.1],
);
let decreasing =
MonotoneISplineDesign::new(&x, basis, MonotoneDirection::Decreasing).unwrap();
for row in 1..decreasing.nrows() {
assert!(decreasing.eta_row(row - 1, &beta) >= decreasing.eta_row(row, &beta) - 1.0e-12);
}
}
fn assert_penalty_gradient_matches_finite_difference<P>(penalty: &P, beta: &[f64])
where
P: Penalty,
{
let eps = 1.0e-6;
let mut grad = vec![0.0; beta.len()];
penalty.add_gradient(beta, &mut grad);
for index in 0..beta.len() {
let mut plus = beta.to_vec();
plus[index] += eps;
let mut minus = beta.to_vec();
minus[index] -= eps;
let finite_difference = (penalty.value(&plus) - penalty.value(&minus)) / (2.0 * eps);
assert_relative_eq!(grad[index], finite_difference, epsilon = 1.0e-6);
}
}
fn assert_predictor_gradient_matches_finite_difference<P>(
predictor: &P,
beta: &[f64],
scores: &[f64],
) where
P: PredictorBlock,
{
let eps = 1.0e-6;
let mut grad = vec![0.0; beta.len()];
predictor.add_gradient(scores, beta, &mut grad);
for index in 0..beta.len() {
let mut plus = beta.to_vec();
plus[index] += eps;
let mut minus = beta.to_vec();
minus[index] -= eps;
let objective = |candidate: &[f64]| {
(0..predictor.nrows())
.map(|row| scores[row] * predictor.eta_row(row, candidate))
.sum::<f64>()
};
let finite_difference = (objective(&plus) - objective(&minus)) / (2.0 * eps);
assert_relative_eq!(grad[index], finite_difference, epsilon = 1.0e-6);
}
}
fn assert_eta_derivative_matches_coordinate_difference<P>(
build: impl Fn(&[f64]) -> P,
derivative: impl Fn(&P, usize, &[f64]) -> f64,
x: &[f64],
beta: &[f64],
) where
P: PredictorBlock,
{
let eps = 1.0e-6;
let design = build(x);
for row in 0..x.len() {
let mut plus_x = x.to_vec();
plus_x[row] += eps;
let plus = build(&plus_x);
let mut minus_x = x.to_vec();
minus_x[row] -= eps;
let minus = build(&minus_x);
let finite_difference =
(plus.eta_row(row, beta) - minus.eta_row(row, beta)) / (2.0 * eps);
assert_relative_eq!(
derivative(&design, row, beta),
finite_difference,
epsilon = 1.0e-5
);
}
}
fn assert_row_basis_matches_evaluate<B>(basis: &B, evaluate: impl Fn(usize) -> Vec<f64>)
where
B: super::SplineRowBasis,
{
for row in 0..basis.nrows() {
let expected = evaluate(row);
let mut actual = vec![0.0; basis.nparams()];
basis.for_each_row_basis(row, |index, weight| {
actual[index] = weight;
});
for (actual, expected) in actual.iter().zip(expected) {
assert_relative_eq!(*actual, expected, epsilon = 1.0e-12);
}
}
}
}