fusion_core/
jacobian.rs

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

fusion_core/
jacobian.rs
