semiflow 0.15.0-beta

no_std Rust PDE solver: evolution equations via Chernoff approximation of operator semigroups (heat, diffusion, Schrodinger)
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

//! Unit tests for `graph_sensitivity.rs`.
//!
//! Moved here to keep `graph_sensitivity.rs` within the 500-line suckless limit.

#[cfg(test)]
mod tests {
    use alloc::{sync::Arc, vec};

    use crate::{
        graph::{Graph, Laplacian},
        graph_sensitivity::{
            adjoint_state_gradient, apply_edge_weight_deriv, magnus_step_jvp_into,
            EdgeWeightSensitivity, GeneratorSensitivity, NodeTimescaleSensitivity,
        },
        graph_signal::GraphSignal,
        scratch::ScratchPool,
        state::State,
    };

    #[test]
    fn rank1_stencil_edge() {
        let n = 3;
        let v = vec![1.0_f64, 2.0, 3.0];
        let mut out = vec![0.0_f64; n];
        apply_edge_weight_deriv(0, 1, &v, &mut out).unwrap();
        let diff = v[0] - v[1]; // = -1.0
        assert!((out[0] - diff).abs() < 1e-14, "out[0]={}", out[0]);
        assert!((out[1] + diff).abs() < 1e-14, "out[1]={}", out[1]);
        assert!(out[2].abs() < 1e-14, "out[2]={}", out[2]);
        assert!(apply_edge_weight_deriv(1, 1, &v, &mut out).is_err());
        assert!(apply_edge_weight_deriv(0, 5, &v, &mut out).is_err());
    }

    #[test]
    fn node_timescale_finite() {
        let n = 3;
        let g = Arc::new(Graph::<f64>::path(n));
        let bare_lap = Laplacian::assemble_combinatorial(&g);
        let sqrt_a = vec![1.0_f64, 2.0_f64.sqrt(), 0.5_f64.sqrt()];
        let sens = NodeTimescaleSensitivity { sqrt_a, bare_lap };
        let v = vec![1.0_f64, -0.5, 0.3];
        let mut out = vec![0.0_f64; n];
        for k in 0..n {
            sens.apply_param_deriv(k, 0.0, &v, &mut out).unwrap();
            assert!(out.iter().all(|x| x.is_finite()), "node {k}: {out:?}");
        }
    }

    /// `accumulate_bilinear` validates `a.len() == b.len() == n_nodes` (ADR-0203).
    #[test]
    fn node_timescale_bilinear_rejects_length_mismatch() {
        let g = Arc::new(Graph::<f64>::path(3));
        let bare_lap = Laplacian::assemble_combinatorial(&g);
        let sens = NodeTimescaleSensitivity {
            sqrt_a: vec![1.0_f64; 3],
            bare_lap,
        };
        let mut scratch = ScratchPool::new();
        let mut grad = vec![0.0_f64; 3];
        let ok = [1.0_f64, 2.0, 3.0];
        let short = [1.0_f64, 2.0];
        assert!(sens
            .accumulate_bilinear(0.1, 1.0, &ok, &ok, &mut grad, &mut scratch)
            .is_ok());
        assert!(sens
            .accumulate_bilinear(0.1, 1.0, &short, &short, &mut grad, &mut scratch)
            .is_err());
        assert!(sens
            .accumulate_bilinear(0.1, 1.0, &ok, &short, &mut grad, &mut scratch)
            .is_err());
        let mut bad_grad = vec![0.0_f64; 2];
        assert!(sens
            .accumulate_bilinear(0.1, 1.0, &ok, &ok, &mut bad_grad, &mut scratch)
            .is_err());
    }

    #[test]
    fn adjoint_grad_finite() {
        use crate::magnus_graph::{LaplacianAtTime, MagnusGraphHeatChernoff};
        let n = 3;
        let g = Arc::new(Graph::<f64>::path(n));
        let g2 = Arc::clone(&g);
        let lap_fn: LaplacianAtTime<f64> =
            alloc::boxed::Box::new(move |_t| Arc::new(Laplacian::assemble_combinatorial(&g2)));
        let mc = MagnusGraphHeatChernoff::new(Arc::clone(&g), lap_fn, 3.0, true).unwrap();
        let u0 = GraphSignal::from_fn(Arc::clone(&g), |i| f64::from(i) * 0.3 + 0.1);
        let tgt = GraphSignal::from_fn(Arc::clone(&g), |i| f64::from(i) * 0.1 - 0.2);
        let n_steps = 2;
        let tau = 0.02_f64;
        let mut scratch = ScratchPool::new();
        let mut cur = u0.clone();
        let mut nxt = GraphSignal::zeros(Arc::clone(&g));
        for k in 0..n_steps {
            #[allow(clippy::cast_precision_loss)]
            mc.apply_into_at(k as f64 * tau, tau, &cur, &mut nxt, &mut scratch)
                .unwrap();
            core::mem::swap(&mut cur, &mut nxt);
        }
        let mut dj = cur;
        State::<f64>::axpy_into(&mut dj, -1.0, &tgt);
        let sens = EdgeWeightSensitivity {
            params: vec![(0, 1), (1, 2)],
            n_nodes: n,
        };
        let mut grad = vec![0.0_f64; 2];
        adjoint_state_gradient(&mc, &u0, n_steps, tau, &dj, &sens, &mut grad, &mut scratch)
            .unwrap();
        assert!(grad.iter().all(|x| x.is_finite()), "gradient: {grad:?}");
    }

    #[test]
    fn jvp_into_finite() {
        let n = 3;
        let g = Arc::new(Graph::<f64>::path(n));
        let lap = Laplacian::assemble_combinatorial(&g);
        let u = vec![1.0_f64, 0.5, -0.3];
        let mut out = vec![0.0_f64; n];
        let mut scratch = ScratchPool::new();
        magnus_step_jvp_into(&lap, &lap, &lap, &lap, 0.01, &u, &mut out, &mut scratch).unwrap();
        assert!(out.iter().all(|x| x.is_finite()), "jvp: {out:?}");
    }
}