fusion_core/
jacobian.rs

1// ─────────────────────────────────────────────────────────────────────
2// SCPN Fusion Core — Analytical Jacobian Utilities
3// © 1998–2026 Miroslav Šotek. All rights reserved.
4// Contact: www.anulum.li | protoscience@anulum.li
5// ORCID: https://orcid.org/0009-0009-3560-0851
6// License: GNU AGPL v3 | Commercial licensing available
7// ─────────────────────────────────────────────────────────────────────
8//! Jacobian builders for inverse profile reconstruction.
9//!
10//! The forward model used in this crate-level inverse module maps probe-space
11//! normalized flux samples to synthetic measurements:
12//!   m_i = p(psi_i) + ff(psi_i)
13//! where p and ff are mTanh profiles with independent parameter sets.
14
15use crate::source::{
16    mtanh_profile, mtanh_profile_derivatives as source_mtanh_profile_derivatives, ProfileParams,
17};
18use fusion_types::error::{FusionError, FusionResult};
19
20fn validate_profile_params(params: &ProfileParams, label: &str) -> FusionResult<()> {
21    if !params.ped_top.is_finite() || params.ped_top <= 0.0 {
22        return Err(FusionError::ConfigError(format!(
23            "{label}.ped_top must be finite and > 0, got {}",
24            params.ped_top
25        )));
26    }
27    if !params.ped_width.is_finite() || params.ped_width <= 0.0 {
28        return Err(FusionError::ConfigError(format!(
29            "{label}.ped_width must be finite and > 0, got {}",
30            params.ped_width
31        )));
32    }
33    if !params.ped_height.is_finite() || !params.core_alpha.is_finite() {
34        return Err(FusionError::ConfigError(format!(
35            "{label}.ped_height/core_alpha must be finite"
36        )));
37    }
38    Ok(())
39}
40
41fn validate_probe_psi_norm(probe_psi_norm: &[f64]) -> FusionResult<()> {
42    if probe_psi_norm.is_empty() {
43        return Err(FusionError::ConfigError(
44            "probe_psi_norm must be non-empty".to_string(),
45        ));
46    }
47    if probe_psi_norm
48        .iter()
49        .any(|psi| !psi.is_finite() || *psi < 0.0 || *psi > 1.0)
50    {
51        return Err(FusionError::ConfigError(
52            "probe_psi_norm values must be finite and within [0, 1]".to_string(),
53        ));
54    }
55    Ok(())
56}
57
58/// Re-export profile derivatives for callers expecting this symbol in `jacobian`.
59pub fn mtanh_profile_derivatives(psi_norm: f64, params: &ProfileParams) -> [f64; 4] {
60    source_mtanh_profile_derivatives(psi_norm, params)
61}
62
63/// Forward model used by inverse reconstruction.
64pub fn forward_model_response(
65    probe_psi_norm: &[f64],
66    params_p: &ProfileParams,
67    params_ff: &ProfileParams,
68) -> FusionResult<Vec<f64>> {
69    validate_probe_psi_norm(probe_psi_norm)?;
70    validate_profile_params(params_p, "params_p")?;
71    validate_profile_params(params_ff, "params_ff")?;
72
73    let out: Vec<f64> = probe_psi_norm
74        .iter()
75        .map(|&psi| mtanh_profile(psi, params_p) + mtanh_profile(psi, params_ff))
76        .collect();
77    if out.iter().any(|v| !v.is_finite()) {
78        return Err(FusionError::ConfigError(
79            "forward model response contains non-finite values".to_string(),
80        ));
81    }
82    Ok(out)
83}
84
85/// Build analytical Jacobian with 8 columns:
86/// [p_ped_height, p_ped_top, p_ped_width, p_core_alpha,
87///  ff_ped_height, ff_ped_top, ff_ped_width, ff_core_alpha].
88pub fn compute_analytical_jacobian(
89    params_p: &ProfileParams,
90    params_ff: &ProfileParams,
91    probe_psi_norm: &[f64],
92) -> FusionResult<Vec<Vec<f64>>> {
93    validate_probe_psi_norm(probe_psi_norm)?;
94    validate_profile_params(params_p, "params_p")?;
95    validate_profile_params(params_ff, "params_ff")?;
96
97    let mut jac = Vec::with_capacity(probe_psi_norm.len());
98    for &psi in probe_psi_norm {
99        let dp = mtanh_profile_derivatives(psi, params_p);
100        let df = mtanh_profile_derivatives(psi, params_ff);
101        if dp.iter().any(|v| !v.is_finite()) || df.iter().any(|v| !v.is_finite()) {
102            return Err(FusionError::ConfigError(
103                "analytical jacobian derivative evaluation became non-finite".to_string(),
104            ));
105        }
106        jac.push(vec![dp[0], dp[1], dp[2], dp[3], df[0], df[1], df[2], df[3]]);
107    }
108    if jac.iter().flatten().any(|v| !v.is_finite()) {
109        return Err(FusionError::ConfigError(
110            "analytical jacobian contains non-finite values".to_string(),
111        ));
112    }
113    Ok(jac)
114}
115
116fn perturb_param(
117    params_p: &ProfileParams,
118    params_ff: &ProfileParams,
119    idx: usize,
120    delta: f64,
121) -> (ProfileParams, ProfileParams) {
122    let mut p = *params_p;
123    let mut ff = *params_ff;
124    match idx {
125        0 => p.ped_height += delta,
126        1 => p.ped_top += delta,
127        2 => p.ped_width += delta,
128        3 => p.core_alpha += delta,
129        4 => ff.ped_height += delta,
130        5 => ff.ped_top += delta,
131        6 => ff.ped_width += delta,
132        7 => ff.core_alpha += delta,
133        _ => {}
134    }
135    (p, ff)
136}
137
138/// Finite-difference Jacobian using forward difference.
139pub fn compute_fd_jacobian(
140    params_p: &ProfileParams,
141    params_ff: &ProfileParams,
142    probe_psi_norm: &[f64],
143    fd_step: f64,
144) -> FusionResult<Vec<Vec<f64>>> {
145    if !fd_step.is_finite() || fd_step <= 0.0 {
146        return Err(FusionError::ConfigError(
147            "jacobian fd_step must be finite and > 0".to_string(),
148        ));
149    }
150    validate_probe_psi_norm(probe_psi_norm)?;
151    validate_profile_params(params_p, "params_p")?;
152    validate_profile_params(params_ff, "params_ff")?;
153
154    let base = forward_model_response(probe_psi_norm, params_p, params_ff)?;
155    let h = fd_step;
156    const PARAM_INDICES: [usize; 8] = [0, 1, 2, 3, 4, 5, 6, 7];
157
158    let mut jac = vec![vec![0.0; 8]; probe_psi_norm.len()];
159    for col in PARAM_INDICES {
160        let (p_pert, ff_pert) = perturb_param(params_p, params_ff, col, h);
161        let f_pert = forward_model_response(probe_psi_norm, &p_pert, &ff_pert)?;
162        for (row, (&fp, &b)) in jac.iter_mut().zip(f_pert.iter().zip(base.iter())) {
163            row[col] = (fp - b) / h;
164        }
165    }
166    if jac.iter().flatten().any(|v| !v.is_finite()) {
167        return Err(FusionError::ConfigError(
168            "finite-difference jacobian contains non-finite values".to_string(),
169        ));
170    }
171    Ok(jac)
172}
173
174#[cfg(test)]
175mod tests {
176    use super::*;
177
178    #[test]
179    fn test_analytical_vs_fd_jacobian() {
180        let p = ProfileParams {
181            ped_top: 0.9,
182            ped_width: 0.08,
183            ped_height: 1.1,
184            core_alpha: 0.25,
185        };
186        let ff = ProfileParams {
187            ped_top: 0.85,
188            ped_width: 0.06,
189            ped_height: 0.95,
190            core_alpha: 0.15,
191        };
192        let probes: Vec<f64> = (0..32).map(|i| i as f64 / 31.0).collect();
193
194        let ja = compute_analytical_jacobian(&p, &ff, &probes).expect("valid analytical jacobian");
195        let jf = compute_fd_jacobian(&p, &ff, &probes, 1e-6).expect("valid fd_step");
196        for i in 0..probes.len() {
197            for j in 0..8 {
198                let abs = (ja[i][j] - jf[i][j]).abs();
199                let denom = jf[i][j].abs().max(1e-8);
200                let rel = abs / denom;
201                assert!(
202                    abs < 1e-8 || rel < 5e-2,
203                    "Mismatch at row {i}, col {j}: analytical={}, fd={}, abs={}, rel={}",
204                    ja[i][j],
205                    jf[i][j],
206                    abs,
207                    rel
208                );
209            }
210        }
211    }
212
213    #[test]
214    fn test_analytical_jacobian_symmetry() {
215        let params = ProfileParams {
216            ped_top: 0.92,
217            ped_width: 0.07,
218            ped_height: 1.0,
219            core_alpha: 0.2,
220        };
221        let probes: Vec<f64> = (0..16).map(|i| i as f64 / 15.0).collect();
222        let jac = compute_analytical_jacobian(&params, &params, &probes)
223            .expect("valid analytical jacobian");
224
225        // If p and ff parameters are identical, corresponding Jacobian columns should match.
226        for row in &jac {
227            for k in 0..4 {
228                assert!(
229                    (row[k] - row[k + 4]).abs() < 1e-12,
230                    "Expected symmetry in Jacobian columns {k} and {}",
231                    k + 4
232                );
233            }
234        }
235    }
236
237    #[test]
238    fn test_fd_jacobian_rejects_invalid_fd_step() {
239        let p = ProfileParams {
240            ped_top: 0.9,
241            ped_width: 0.08,
242            ped_height: 1.1,
243            core_alpha: 0.25,
244        };
245        let ff = p;
246        let probes: Vec<f64> = (0..8).map(|i| i as f64 / 7.0).collect();
247
248        assert!(compute_fd_jacobian(&p, &ff, &probes, 0.0).is_err());
249        assert!(compute_fd_jacobian(&p, &ff, &probes, -1e-6).is_err());
250        assert!(compute_fd_jacobian(&p, &ff, &probes, f64::NAN).is_err());
251    }
252
253    #[test]
254    fn test_forward_model_rejects_invalid_runtime_inputs() {
255        let p = ProfileParams {
256            ped_top: 0.9,
257            ped_width: 0.08,
258            ped_height: 1.1,
259            core_alpha: 0.25,
260        };
261        let ff = p;
262        assert!(forward_model_response(&[], &p, &ff).is_err());
263        assert!(forward_model_response(&[0.1, f64::INFINITY], &p, &ff).is_err());
264        let bad = ProfileParams {
265            ped_top: 0.9,
266            ped_width: 0.0,
267            ped_height: 1.0,
268            core_alpha: 0.2,
269        };
270        assert!(forward_model_response(&[0.2, 0.4], &bad, &ff).is_err());
271        assert!(compute_analytical_jacobian(&p, &ff, &[0.2, -0.1]).is_err());
272    }
273}
fusion_core/jacobian.rs

fusion_core/
jacobian.rs
