cheby 0.4.0

Unit-safe Chebyshev approximation and spectral numerics for Rust.
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

use approx::assert_abs_diff_eq;
use cheby::approx::interpolation::BarycentricInterpolator;
use cheby::core::{basis, nodes, ChebySeries, Domain, NodeKind};
use cheby::piecewise::ChebySegmentTable;
use proptest::prelude::*;

proptest! {
    #[test]
    fn domain_normalization_properties(start in -1e3_f64..1e3, width in 1e-3_f64..1e3) {
        let domain = Domain::try_new(start, start + width).unwrap();
        // Maximum floating-point error in normalize is O(eps * |start| / half).
        // With these ranges (|start| ≤ 1e3, half ≥ 5e-4), that is < 1e-9.
        assert_abs_diff_eq!(domain.normalize(domain.start()), -1.0, epsilon = 1e-9);
        assert_abs_diff_eq!(domain.normalize(domain.midpoint()), 0.0, epsilon = 1e-9);
        assert_abs_diff_eq!(domain.normalize(domain.end()), 1.0, epsilon = 1e-9);
    }

    #[test]
    fn chebyshev_recurrence_holds(x in -1.0_f64..1.0, n in 1_usize..20) {
        let lhs = basis::t(n + 1, x);
        let rhs = 2.0 * x * basis::t(n, x) - basis::t(n - 1, x);
        assert!((lhs - rhs).abs() < 1e-12);
    }
}

#[test]
fn interpolation_reproduces_samples() {
    const N: usize = 9;
    let xs = nodes::<N>(NodeKind::Roots);
    let ys = std::array::from_fn(|k| xs[k].sin());
    let interp = BarycentricInterpolator::new(xs, ys).unwrap();
    for k in 0..N {
        assert_abs_diff_eq!(interp.evaluate(xs[k]).unwrap(), ys[k], epsilon = 1e-12);
    }
}

#[test]
fn interpolation_single_node_is_constant_even_at_zero() {
    // Regression: a single node at zero used to make `scale == 0`, producing
    // a NaN division during evaluation instead of the constant interpolant.
    let interp = BarycentricInterpolator::new([0.0_f64], [42.0_f64]).unwrap();
    for &x in &[-2.5_f64, -0.1, 0.0, 0.7, 17.0] {
        assert_abs_diff_eq!(interp.evaluate(x).unwrap(), 42.0, epsilon = 0.0);
    }
}

#[test]
fn interpolation_chebyshev_roots_constructor_matches_generic() {
    use cheby::core::{nodes_mapped, Domain, NodeKind};
    const N: usize = 12;
    let domain = Domain::new(-1.0_f64, 3.0);
    let xs = nodes_mapped::<f64, N>(domain, NodeKind::Roots);
    let ys: [f64; N] = std::array::from_fn(|k| (xs[k] + 0.3).cos());
    let generic = BarycentricInterpolator::new(xs, ys).unwrap();
    let fast = BarycentricInterpolator::on_chebyshev_roots(domain, ys).unwrap();
    for tau in [-0.95_f64, -0.4, 0.0, 0.2, 0.85] {
        let x = domain.denormalize(tau);
        assert_abs_diff_eq!(
            fast.evaluate(x).unwrap(),
            generic.evaluate(x).unwrap(),
            epsilon = 1e-12
        );
    }
}

#[test]
fn interpolation_lobatto_constructor_reproduces_samples() {
    use cheby::core::{nodes_mapped, Domain, NodeKind};
    const N: usize = 9;
    let domain = Domain::new(0.0_f64, 2.0);
    let xs = nodes_mapped::<f64, N>(domain, NodeKind::Lobatto);
    let ys: [f64; N] = std::array::from_fn(|k| xs[k].sin());
    let interp = BarycentricInterpolator::on_lobatto_nodes(domain, ys).unwrap();
    for k in 0..N {
        assert_abs_diff_eq!(interp.evaluate(xs[k]).unwrap(), ys[k], epsilon = 1e-12);
    }
}

#[test]
fn derivative_of_integral_recovers_series_approximately() {
    let series = ChebySeries::new([1.0, -0.5, 0.25, 0.0, 0.0]);
    let recovered = series.integral(0.0).derivative();
    for tau in [-0.75, -0.25, 0.25, 0.75] {
        assert_abs_diff_eq!(
            recovered.evaluate(tau),
            series.evaluate(tau),
            epsilon = 1e-12
        );
    }
}

#[test]
fn piecewise_lookup_chooses_correct_segment() {
    let table: ChebySegmentTable<f64, f64, 5> = ChebySegmentTable::from_fn(|x| x, 0.0, 4.0, 1.0);
    assert_abs_diff_eq!(table.get_segment(0.25).unwrap().domain().start(), 0.0);
    assert_abs_diff_eq!(table.get_segment(1.25).unwrap().domain().start(), 1.0);
    assert_abs_diff_eq!(table.get_segment(3.25).unwrap().domain().start(), 3.0);
}

#[test]
#[cfg(feature = "adaptive")]
fn adaptive_fitting_converges_for_smooth_function() {
    let result =
        cheby::approx::adaptive::AdaptiveFit::new(Domain::new(-1.0, 1.0), |x: f64| x.sin())
            .max_degree(64)
            .absolute_tolerance(1e-9)
            .relative_tolerance(1e-9)
            .build::<f64>()
            .unwrap();
    assert!(result.report.converged);
}

#[test]
#[cfg(feature = "adaptive")]
fn adaptive_fitting_honors_non_power_of_two_limit() {
    let result =
        cheby::approx::adaptive::AdaptiveFit::new(Domain::new(-1.0, 1.0), |x: f64| (3.0 * x).sin())
            .max_degree(20)
            .absolute_tolerance(0.0)
            .relative_tolerance(0.0)
            .build::<f64>()
            .unwrap();
    assert_eq!(result.report.degree, 20);
}

#[test]
#[cfg(feature = "adaptive")]
fn adaptive_fitting_allows_small_degree_limit() {
    let result =
        cheby::approx::adaptive::AdaptiveFit::new(Domain::new(-1.0, 1.0), |x: f64| x.sin())
            .max_degree(7)
            .absolute_tolerance(0.0)
            .relative_tolerance(0.0)
            .build::<f64>()
            .unwrap();
    assert_eq!(result.report.degree, 7);
}