noobase 0.0.6

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
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174

//! Cross-module integration tests for the convolve surface.
//!
//! These exercise the *public* API only (`noobase::convolve::*`,
//! `noobase::image::convolve_*`, `noobase::spectroscopy::Spectrum::
//! convolve_lsf`) to pin that the layers compose across module
//! boundaries as documented: the bare kernels feed the NaN wrappers,
//! which feed the image domain entries; and the LSF entry composes the
//! `bins` Grid with the `convolve` Gaussian. Per-module behavior is
//! already covered by the in-crate `mod tests`; this file is the
//! outside-the-crate contract.

use ndarray::{Array1, Array2, Axis};
use noobase::axis::{Grid, GridKind};
use noobase::convolve::{Boundary, GaussianSampling, Normalization, conv2d, gaussian1d};
use noobase::image::{convolve_gaussian_axis, convolve_psf};
use noobase::spectroscopy::{LsfSpec, Spectrum};

fn approx_eq(a: f64, b: f64, tol: f64) -> bool {
    (a - b).abs() <= tol * a.abs().max(b.abs()).max(1.0)
}

#[test]
fn convolve_lsf_conserves_constant_via_public_api_f64() {
    // Spectrum (spectroscopy) + Grid (bins) + Gaussian (convolve) all
    // composed through the public surface; constant template stays
    // constant including the edges.
    let wavelength = Grid::<f64>::logspace(4000.0, 7000.0, 200, GridKind::Centers);
    let flux = Array1::<f64>::from_elem(200, 2.5);
    let spectrum = Spectrum::new(wavelength, flux, None, None).unwrap();
    let broadened = spectrum.convolve_lsf(LsfSpec::ConstantR(2500.0)).unwrap();
    for value in broadened.flux().iter() {
        assert!(approx_eq(*value, 2.5, 1e-9));
    }
    assert!(broadened.error().is_none());
    assert!(broadened.mask().is_none());
}

#[test]
fn convolve_lsf_rejects_non_template_via_public_api() {
    let wavelength = Grid::<f64>::logspace(4000.0, 7000.0, 32, GridKind::Centers);
    let flux = Array1::<f64>::from_elem(32, 1.0);
    let error = Array1::<f64>::from_elem(32, 0.1);
    let observed = Spectrum::new(wavelength, flux, Some(error), None).unwrap();
    assert!(observed.convolve_lsf(LsfSpec::ConstantR(2000.0)).is_err());
}

#[test]
fn convolve_psf_interior_composes_bare_conv2d_f64() {
    // image::convolve_psf flips the PSF and renormalizes. With a
    // Sum-normalized PSF and a NaN-free image, the interior weight is 1,
    // so the interior must equal the bare conv2d of the flipped PSF —
    // pinning the documented cross-module composition.
    let mut image = Array2::<f64>::zeros((8, 9));
    for ((i, j), v) in image.indexed_iter_mut() {
        *v = (i * 4 + j * 3) as f64 * 0.1 - 1.5;
    }
    let mut psf = ndarray::array![
        [0.05_f64, 0.10, 0.02],
        [0.15, 0.30, 0.08],
        [0.04, 0.18, 0.08]
    ];
    let sum: f64 = psf.iter().sum();
    psf.mapv_inplace(|v| v / sum);

    let via_entry = convolve_psf(image.view(), psf.view());
    let flipped = psf.slice(ndarray::s![..;-1, ..;-1]).to_owned();
    let via_bare = conv2d(image.view(), flipped.view(), Boundary::Zero);
    for i in 1..7 {
        for j in 1..8 {
            assert!(
                approx_eq(via_entry[(i, j)], via_bare[(i, j)], 1e-9),
                "interior mismatch at ({i},{j})"
            );
        }
    }
}

#[test]
fn convolve_gaussian_axis_matched_filter_peaks_at_injected_sigma_f64() {
    let rows = 121usize;
    let cols = 4usize;
    let line_row = 60.0_f64;
    let injected_sigma = 5.0_f64;
    let mut image = Array2::<f64>::zeros((rows, cols));
    for i in 0..rows {
        let z = (i as f64 - line_row) / injected_sigma;
        let amplitude = (-0.5 * z * z).exp();
        for j in 0..cols {
            image[(i, j)] = amplitude;
        }
    }
    let peak = |probe: f64| {
        convolve_gaussian_axis(
            image.view(),
            probe,
            Axis(0),
            Normalization::L2,
            Boundary::Zero,
            false,
        )[(60, 0)]
    };
    let matched = peak(injected_sigma);
    assert!(matched > peak(injected_sigma * 0.5));
    assert!(matched > peak(injected_sigma * 2.0));
}

#[test]
fn gaussian1d_public_normalizations() {
    // Confirms the convolve re-exports are reachable at the crate path
    // and the two kernel normalizations behave.
    let sum_kernel = gaussian1d::<f64>(
        1.7,
        4.0,
        GaussianSampling::ErfIntegrated,
        Normalization::Sum,
    );
    let l2_kernel = gaussian1d::<f64>(1.7, 4.0, GaussianSampling::ErfIntegrated, Normalization::L2);
    let sum: f64 = sum_kernel.iter().sum();
    let sum_sq: f64 = l2_kernel.iter().map(|v| v * v).sum();
    assert!(approx_eq(sum, 1.0, 1e-12));
    assert!(approx_eq(sum_sq, 1.0, 1e-12));
}

#[test]
fn f32_and_f64_agree_within_f32_precision() {
    // The convolve entries are generic over the float channel; the same
    // data as f32 or f64 must agree to within f32 precision.
    let mut image_f64 = Array2::<f64>::zeros((10, 11));
    for ((i, j), v) in image_f64.indexed_iter_mut() {
        *v = (i as f64) * 0.3 + (j as f64) * 0.2 + 1.0;
    }
    let image_f32 = image_f64.mapv(|v| v as f32);
    let psf_f64 = Array2::<f64>::from_elem((3, 3), 1.0 / 9.0);
    let psf_f32 = psf_f64.mapv(|v| v as f32);

    let out_f64 = convolve_psf(image_f64.view(), psf_f64.view());
    let out_f32 = convolve_psf(image_f32.view(), psf_f32.view());
    let tolerance = 1e-5;
    for i in 0..10 {
        for j in 0..11 {
            let a = out_f64[(i, j)];
            let b = out_f32[(i, j)] as f64;
            let scale = a.abs().max(b.abs()).max(1.0);
            assert!(
                (a - b).abs() <= tolerance * scale,
                "f32/f64 mismatch at ({i},{j}): {a} vs {b}"
            );
        }
    }

    // Same for the LSF entry.
    let wavelength_f64 = Grid::<f64>::logspace(4000.0, 7000.0, 128, GridKind::Centers);
    let wavelength_f32 = Grid::<f32>::logspace(4000.0, 7000.0, 128, GridKind::Centers);
    let flux_f64: Array1<f64> = (0..128)
        .map(|i| 1.0 + 0.3 * (i as f64 / 7.0).sin())
        .collect();
    let flux_f32 = flux_f64.mapv(|v| v as f32);
    let lsf_f64 = Spectrum::new(wavelength_f64, flux_f64, None, None)
        .unwrap()
        .convolve_lsf(LsfSpec::ConstantR(3000.0))
        .unwrap();
    let lsf_f32 = Spectrum::new(wavelength_f32, flux_f32, None, None)
        .unwrap()
        .convolve_lsf(LsfSpec::ConstantR(3000.0))
        .unwrap();
    for (a, b) in lsf_f64.flux().iter().zip(lsf_f32.flux().iter()) {
        let b = *b as f64;
        let scale = a.abs().max(b.abs()).max(1.0);
        assert!(
            (a - b).abs() <= 1e-4 * scale,
            "lsf f32/f64 mismatch: {a} vs {b}"
        );
    }
}