noobase 0.0.3

Foundational pure-function utilities for astronomy analysis
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
144
145
146
147
148
149
150
151
152

//! Shared separable bicubic Catmull-Rom interpolation kernel.
//!
//! This is the single source of the interpolation weights used by both
//! the forward operator ([`super::render`]) and its adjoint
//! ([`super::accumulate`]). Sharing one weight function is the structural
//! guarantee that `accumulate` is the *exact* transpose of `render`:
//! `render` gathers `value += epsf[ru, rv] * w_u * w_v` while
//! `accumulate` scatters `epsf_acc[ru, rv] += coeff * w_u * w_v` over the
//! very same per-tap weights, the same `floor`-derived tap base, and the
//! same zero-padding rule. Were the two to compute weights independently,
//! any drift between them would silently break the inner-product
//! transpose identity that the solver relies on.
//!
//! Keys cubic convolution with `a = -1/2` (the Catmull-Rom member),
//! four taps per axis. It is a pure-Rust, psf-module-internal helper that
//! never crosses the FFI boundary and is never visible to the caller.

use ndarray::ArrayView2;

/// Keys cubic-convolution (Catmull-Rom, `a = -1/2`) weights for the four
/// taps surrounding a fractional position `frac` in `[0, 1)`. The
/// returned array is ordered `[w(-1), w(0), w(+1), w(+2)]` relative to
/// `floor(coordinate)`.
///
/// The weights sum to one (partition of unity) and at `frac = 0` collapse
/// to `[0, 1, 0, 0]`, which is what makes an integer-aligned sample an
/// exact, interpolation-free grid pick (and, on the adjoint side, an
/// exact grid scatter).
#[inline]
pub(super) fn catmull_rom_weights(frac: f64) -> [f64; 4] {
    // Catmull-Rom basis on the unit interval, parameter t = frac:
    //   w(-1) = -0.5 t^3 +      t^2 - 0.5 t
    //   w( 0) =  1.5 t^3 -  2.5 t^2       + 1
    //   w(+1) = -1.5 t^3 +  2.0 t^2 + 0.5 t
    //   w(+2) =  0.5 t^3 -  0.5 t^2
    let t = frac;
    [
        ((-0.5 * t + 1.0) * t - 0.5) * t,
        ((1.5 * t - 2.5) * t) * t + 1.0,
        ((-1.5 * t + 2.0) * t + 0.5) * t,
        ((0.5 * t - 0.5) * t) * t,
    ]
}

/// Derivatives of the four Keys cubic-convolution (Catmull-Rom,
/// `a = -1/2`) tap weights with respect to the fractional position
/// `frac`. The returned array is ordered `[w'(-1), w'(0), w'(+1),
/// w'(+2)]`, tap-for-tap aligned with [`catmull_rom_weights`].
///
/// These are the exact term-by-term derivatives of the polynomials in
/// [`catmull_rom_weights`]. Differentiating the partition-of-unity
/// identity `sum w = 1` gives `sum w' = 0`, so the four returned values
/// always sum to zero. They are the structural source of the analytic
/// `dS/ddelta` Jacobian used by [`super::nuisance`]: with the render
/// sampling coordinate `k = K_c + os * ((i - c_det) - delta)` we have
/// `dk/ddelta = -os` and `dfrac/dk = 1`, so the centroid Jacobian is
/// these weight derivatives scaled by `-os` (pinned against a finite
/// difference of `render` in the nuisance tests). Sharing this single
/// derivative source with the forward kernel is the same structural
/// guarantee that keeps `accumulate` an exact transpose of `render`.
#[inline]
pub(super) fn catmull_rom_weight_derivatives(frac: f64) -> [f64; 4] {
    // d/dt of the basis written out in `catmull_rom_weights`:
    //   w(-1)' = -1.5 t^2 + 2.0 t - 0.5
    //   w( 0)' =  4.5 t^2 - 5.0 t
    //   w(+1)' = -4.5 t^2 + 4.0 t + 0.5
    //   w(+2)' =  1.5 t^2 - 1.0 t
    let t = frac;
    [
        (-1.5 * t + 2.0) * t - 0.5,
        (4.5 * t - 5.0) * t,
        (-4.5 * t + 4.0) * t + 0.5,
        (1.5 * t - 1.0) * t,
    ]
}

/// Separable bicubic Catmull-Rom point sample of `epsf` at the continuous
/// model coordinate `(k_u, k_v)`. Taps that fall outside the model grid
/// are zero (zero padding).
///
/// The adjoint scatters the identical per-tap weights returned by
/// [`catmull_rom_weights`] at the identical tap base derived here, so this
/// function and the scatter stay exact transposes by construction.
pub(super) fn catmull_rom_sample(epsf: &ArrayView2<f64>, k_u: f64, k_v: f64) -> f64 {
    let side = epsf.shape()[0] as i64; // square: shape()[0] == shape()[1]

    let u_floor = k_u.floor();
    let v_floor = k_v.floor();
    let weights_u = catmull_rom_weights(k_u - u_floor);
    let weights_v = catmull_rom_weights(k_v - v_floor);

    // The four taps of a Catmull-Rom segment sit at integer offsets
    // -1, 0, +1, +2 from floor(k); `weights_*[t]` is the weight of tap
    // `base + t`. The adjoint must scatter these very same weights.
    let base_u = u_floor as i64 - 1;
    let base_v = v_floor as i64 - 1;

    let mut value = 0.0_f64;
    for (tap_u, &weight_u) in weights_u.iter().enumerate() {
        let row = base_u + tap_u as i64;
        if row < 0 || row >= side {
            continue; // zero-pad outside the model grid
        }
        for (tap_v, &weight_v) in weights_v.iter().enumerate() {
            let column = base_v + tap_v as i64;
            if column < 0 || column >= side {
                continue; // zero-pad outside the model grid
            }
            value += epsf[(row as usize, column as usize)] * weight_u * weight_v;
        }
    }
    value
}

#[cfg(test)]
mod tests {
    use super::{catmull_rom_weight_derivatives, catmull_rom_weights};

    #[test]
    fn weight_derivatives_match_central_difference() {
        // Pin the analytic 4-tap weight derivatives against a central
        // finite difference of `catmull_rom_weights`. A sign or
        // coefficient slip in the polynomial derivative (which would
        // flip the GN-delta step) is gross relative to this bound.
        let step = 1e-6;
        for &frac in &[0.05, 0.2, 0.37, 0.5, 0.63, 0.8, 0.95] {
            let analytic = catmull_rom_weight_derivatives(frac);
            let forward = catmull_rom_weights(frac + step);
            let backward = catmull_rom_weights(frac - step);
            for tap in 0..4 {
                let finite_difference = (forward[tap] - backward[tap]) / (2.0 * step);
                assert!(
                    (analytic[tap] - finite_difference).abs() < 1e-6,
                    "tap {tap} at frac {frac}: analytic {} vs finite difference {}",
                    analytic[tap],
                    finite_difference
                );
            }
        }
    }

    #[test]
    fn weight_derivatives_sum_to_zero() {
        // Differentiating the partition-of-unity identity (sum w = 1)
        // gives sum w' = 0 for every fractional position.
        for &frac in &[0.0, 0.1, 0.25, 0.5, 0.75, 0.999] {
            let derivatives = catmull_rom_weight_derivatives(frac);
            let sum: f64 = derivatives.iter().sum();
            assert!(sum.abs() < 1e-12, "sum {sum} at frac {frac}");
        }
    }
}