solver/numeric/finite_difference/discretization.rs
1//! Shared finite-difference discretization primitives.
2//!
3//! `DimensionKind`, `Transport`, and the boundary-condition types are shared
4//! by the time-dependent [`super::pde::PdeProblem`] path and the stationary
5//! elliptic solver. They live here so both solution strategies use the same
6//! spatial stencils without reimplementing them.
7
8/// Declarative discretization kind for one state coordinate.
9///
10/// The finite-difference solvers use this to select the stencil without
11/// inspecting which derivative components a model happens to populate.
12#[derive(Clone, Copy, Debug, PartialEq, Eq)]
13pub enum DimensionKind {
14 /// Integer state with unit jumps under a point-process intensity. The
15 /// forward and backward transition rates are used directly, not divided
16 /// by a spacing.
17 DiscreteJump,
18 /// Integer state with arbitrary-amplitude jumps under a point-process
19 /// intensity. The jump amplitudes and rates are supplied by the
20 /// [`JumpKernel`] associated with the dimension.
21 Jump,
22 /// Continuous state driven by Brownian diffusion. Uses central second
23 /// differences plus upwind first differences.
24 Diffusion,
25 /// Continuous state with deterministic drift and no diffusion. Uses
26 /// upwind first differences only.
27 DeterministicDrift,
28}
29
30/// Per-dimension transport coefficients for a finite-difference operator.
31///
32/// The operator is decomposed as
33///
34/// ```text
35/// L V = sum_i [ plus_i (V_{i,+} - V) + minus_i (V_{i,-} - V) ] + source
36/// ```
37///
38/// `plus_i` and `minus_i` are the forward and backward transport coefficients
39/// documented in `docs/src/solver/solver_pde.md`. `source` collects the purely
40/// local terms.
41///
42/// The transport coefficients and the source have clean semantics: transport
43/// goes into the operator and the source goes into the right-hand side. There
44/// is no hidden cancellation term.
45#[derive(Clone, Copy, Debug)]
46pub struct Transport<const N: usize> {
47 /// Forward transport coefficient in each dimension.
48 pub plus: [f64; N],
49 /// Backward transport coefficient in each dimension.
50 pub minus: [f64; N],
51 /// Purely local source term.
52 pub source: f64,
53}
54
55impl<const N: usize> Transport<N> {
56 /// Builds a transport bundle from raw coefficients.
57 ///
58 /// # Examples
59 ///
60 /// ```
61 /// use solver::numeric::finite_difference::discretization::{Transport, DimensionKind};
62 /// let t = Transport::<2>::new([0.5, 0.1], [0.25, 0.0], 1.0);
63 /// assert_eq!(t.plus, [0.5, 0.1]);
64 /// assert_eq!(t.minus, [0.25, 0.0]);
65 /// assert_eq!(t.source, 1.0);
66 /// let _kind = DimensionKind::Diffusion;
67 /// ```
68 pub fn new(plus: [f64; N], minus: [f64; N], source: f64) -> Self {
69 Self {
70 plus,
71 minus,
72 source,
73 }
74 }
75
76 /// An empty transport bundle with a zero source.
77 pub fn zero() -> Self {
78 Self {
79 plus: [0.0; N],
80 minus: [0.0; N],
81 source: 0.0,
82 }
83 }
84}
85
86/// A single jump transition in a jump-kernel dimension.
87///
88/// `amplitude` is the integer change in the state coordinate (positive for an
89/// upward jump, negative for a downward jump); `rate` is the non-negative
90/// point-process intensity at which the jump occurs. The state is interpreted
91/// on an integer lattice, so the amplitude is applied directly without a grid
92/// spacing.
93#[derive(Clone, Copy, Debug, PartialEq)]
94pub struct JumpTransition {
95 /// Integer change in the state coordinate.
96 pub amplitude: i32,
97 /// Transition intensity (arrival rate).
98 pub rate: f64,
99}
100
101impl JumpTransition {
102 /// Builds a jump transition.
103 ///
104 /// # Examples
105 ///
106 /// ```
107 /// use solver::numeric::finite_difference::discretization::JumpTransition;
108 /// let jump = JumpTransition::new(2, 1.5);
109 /// assert_eq!(jump.amplitude, 2);
110 /// assert_eq!(jump.rate, 1.5);
111 /// ```
112 pub fn new(amplitude: i32, rate: f64) -> Self {
113 Self { amplitude, rate }
114 }
115}
116
117/// Per-dimension jump kernel for arbitrary-amplitude jump processes.
118///
119/// Each dimension carries an ordered list of [`JumpTransition`] pairs. For a
120/// `DiscreteJump` dimension the kernel is exactly the unit up/down pair
121/// `[(+1, plus), (-1, minus)]`; a `Jump` dimension may carry any finite set of
122/// amplitudes and rates. The finite-difference operator and the BSDE forward
123/// step both sum or sample over this kernel.
124#[derive(Clone, Debug)]
125pub struct JumpKernel<const N: usize> {
126 /// Jump transitions for each dimension.
127 pub dims: [Vec<JumpTransition>; N],
128}
129
130impl<const N: usize> Default for JumpKernel<N> {
131 fn default() -> Self {
132 Self::empty()
133 }
134}
135
136impl<const N: usize> JumpKernel<N> {
137 /// Builds an empty jump kernel (no jumps in any dimension).
138 ///
139 /// # Examples
140 ///
141 /// ```
142 /// use solver::numeric::finite_difference::discretization::JumpKernel;
143 /// let kernel = JumpKernel::<2>::empty();
144 /// assert!(kernel.is_empty(0));
145 /// assert!(kernel.is_empty(1));
146 /// ```
147 pub fn empty() -> Self {
148 Self {
149 dims: std::array::from_fn(|_| Vec::new()),
150 }
151 }
152
153 /// Builds a jump kernel from a transport bundle, folding each dimension's
154 /// `plus`/`minus` rates into the unit `+1`/`-1` transitions used by the
155 /// `DiscreteJump` convention.
156 ///
157 /// # Examples
158 ///
159 /// ```
160 /// use solver::numeric::finite_difference::discretization::{JumpKernel, Transport};
161 /// let t = Transport::<2>::new([1.5, 0.0], [0.0, 2.0], 0.0);
162 /// let kernel = JumpKernel::from_transport(&t);
163 /// assert_eq!(kernel.dims[0][0].amplitude, 1);
164 /// assert_eq!(kernel.dims[1][0].amplitude, -1);
165 /// ```
166 pub fn from_transport(transport: &Transport<N>) -> Self {
167 let mut kernel = Self::empty();
168 for (dim, jumps) in kernel.dims.iter_mut().enumerate() {
169 if transport.plus[dim] > 0.0 {
170 jumps.push(JumpTransition::new(1, transport.plus[dim]));
171 }
172 if transport.minus[dim] > 0.0 {
173 jumps.push(JumpTransition::new(-1, transport.minus[dim]));
174 }
175 }
176 kernel
177 }
178
179 /// Number of jump transitions in `dim`.
180 ///
181 /// # Examples
182 ///
183 /// ```
184 /// use solver::numeric::finite_difference::discretization::{JumpKernel, Transport};
185 /// let t = Transport::<2>::new([1.5, 0.0], [0.0, 2.0], 0.0);
186 /// let kernel = JumpKernel::from_transport(&t);
187 /// assert_eq!(kernel.len(0), 1);
188 /// assert_eq!(kernel.len(1), 1);
189 /// ```
190 pub fn len(&self, dim: usize) -> usize {
191 self.dims[dim].len()
192 }
193
194 /// Whether a dimension has no jump transitions.
195 ///
196 /// # Examples
197 ///
198 /// ```
199 /// use solver::numeric::finite_difference::discretization::{JumpKernel, Transport};
200 /// let t = Transport::<2>::new([0.0, 0.0], [2.0, 0.0], 0.0);
201 /// let kernel = JumpKernel::from_transport(&t);
202 /// assert!(!kernel.is_empty(0));
203 /// assert!(kernel.is_empty(1));
204 /// ```
205 pub fn is_empty(&self, dim: usize) -> bool {
206 self.dims[dim].is_empty()
207 }
208
209 /// Total intensity of all transitions in `dim`.
210 ///
211 /// # Examples
212 ///
213 /// ```
214 /// use solver::numeric::finite_difference::discretization::{JumpKernel, Transport};
215 /// let t = Transport::<2>::new([1.5, 0.0], [0.0, 2.0], 0.0);
216 /// let kernel = JumpKernel::from_transport(&t);
217 /// assert_eq!(kernel.total_rate(0), 1.5);
218 /// assert_eq!(kernel.total_rate(1), 2.0);
219 /// ```
220 pub fn total_rate(&self, dim: usize) -> f64 {
221 self.dims[dim].iter().map(|j| j.rate).sum()
222 }
223}
224
225/// Samples a single jump from a set of competing Poisson transitions.
226///
227/// Each transition fires with intensity `rate`; the first to fire wins. Under
228/// competing exponentials the winner is transition `i` with probability
229/// `rate_i / sum_j rate_j`, so the amplitude is drawn by mapping a uniform
230/// `u` (via the standard-normal CDF of `noise`) to the cumulative rate. When
231/// the total rate is zero or `u` falls outside `[0, 1)` no jump occurs.
232/// # Examples
233///
234/// ```
235/// use solver::numeric::finite_difference::discretization::{
236/// JumpTransition, sample_jump_kernel,
237/// };
238/// let jumps = [JumpTransition::new(2, 1.0), JumpTransition::new(-3, 1.0)];
239/// // Noise far negative maps the CDF to zero, selecting the first transition.
240/// assert_eq!(sample_jump_kernel(&jumps, -10.0), 2);
241/// // Noise far positive maps the CDF to one, selecting the last transition.
242/// assert_eq!(sample_jump_kernel(&jumps, 10.0), -3);
243/// assert_eq!(sample_jump_kernel(&[], 0.0), 0);
244/// ```
245pub fn sample_jump_kernel(jumps: &[JumpTransition], noise: f64) -> i32 {
246 let total: f64 = jumps.iter().map(|j| j.rate).sum();
247 if total <= 0.0 {
248 return 0;
249 }
250 // Map the noise through the standard-normal CDF to a uniform draw, then
251 // clamp into the half-open interval `[0, 1)` so `u = 1` (a CDF that rounds
252 // to exactly one for large positive noise) falls to the last transition
253 // rather than being treated as "no jump".
254 let u = crate::models::normal_cdf(noise);
255 let u = if u >= 1.0 {
256 1.0 - f64::EPSILON
257 } else if u < 0.0 {
258 0.0
259 } else {
260 u
261 };
262 let mut cumulative = 0.0;
263 for jump in jumps {
264 cumulative += jump.rate;
265 if u * total < cumulative {
266 return jump.amplitude;
267 }
268 }
269 jumps.last().map_or(0, |j| j.amplitude)
270}
271
272/// A boundary condition for one side of one coordinate axis.
273///
274/// The derivative in `Neumann` and `Robin` is taken with respect to the
275/// increasing coordinate direction `x`, not an outward normal. This keeps the
276/// Cartesian interpretation unambiguous.
277#[derive(Clone, Copy, Debug, PartialEq)]
278pub enum BoundaryCondition {
279 /// `u = value` at the boundary.
280 Dirichlet(f64),
281 /// `du/dx = value` at the boundary.
282 Neumann(f64),
283 /// `alpha * u + beta * du/dx = value` at the boundary.
284 Robin {
285 /// Right-hand side of the Robin condition.
286 value: f64,
287 /// Coefficient of the value `u`.
288 alpha: f64,
289 /// Coefficient of the derivative `du/dx`.
290 beta: f64,
291 },
292}
293
294/// Per-dimension lower and upper boundary conditions.
295#[derive(Clone, Debug)]
296pub struct BoundaryConditions<const N: usize> {
297 /// Condition at the lower (`x_i = min_i`) boundary of each dimension.
298 pub lower: [BoundaryCondition; N],
299 /// Condition at the upper (`x_i = max_i`) boundary of each dimension.
300 pub upper: [BoundaryCondition; N],
301}
302
303impl<const N: usize> BoundaryConditions<N> {
304 /// Builds explicit lower and upper boundary condition arrays.
305 pub fn new(lower: [BoundaryCondition; N], upper: [BoundaryCondition; N]) -> Self {
306 Self { lower, upper }
307 }
308}
309
310impl<const N: usize> Default for BoundaryConditions<N> {
311 /// Defaults every side to `Neumann(0.0)` (zero derivative).
312 fn default() -> Self {
313 Self {
314 lower: std::array::from_fn(|_| BoundaryCondition::Neumann(0.0)),
315 upper: std::array::from_fn(|_| BoundaryCondition::Neumann(0.0)),
316 }
317 }
318}
319
320#[cfg(test)]
321mod tests {
322 use super::*;
323
324 #[test]
325 fn default_boundary_conditions_are_zero_neumann() {
326 let b = BoundaryConditions::<2>::default();
327 for cond in b.lower.iter().chain(b.upper.iter()) {
328 assert_eq!(*cond, BoundaryCondition::Neumann(0.0));
329 }
330 }
331
332 #[test]
333 fn transport_zero_is_empty() {
334 let t = Transport::<3>::zero();
335 assert_eq!(t.plus, [0.0; 3]);
336 assert_eq!(t.minus, [0.0; 3]);
337 assert_eq!(t.source, 0.0);
338 }
339
340 #[test]
341 fn jump_kernel_from_transport_folds_unit_jumps() {
342 let t = Transport::<2>::new([1.5, 0.0], [0.0, 2.0], 0.0);
343 let k = JumpKernel::from_transport(&t);
344 assert_eq!(k.dims[0], vec![JumpTransition::new(1, 1.5)]);
345 assert_eq!(k.dims[1], vec![JumpTransition::new(-1, 2.0)]);
346 }
347
348 #[test]
349 fn jump_kernel_empty_by_default() {
350 let k = JumpKernel::<2>::empty();
351 assert!(k.is_empty(0));
352 assert!(k.is_empty(1));
353 }
354
355 #[test]
356 fn sample_jump_kernel_selects_by_cumulative_rate() {
357 let jumps = vec![JumpTransition::new(2, 1.0), JumpTransition::new(-3, 1.0)];
358 // noise negative -> cdf ~ 0 -> first transition (cumulative reaches 1.0).
359 assert_eq!(sample_jump_kernel(&jumps, -10.0), 2);
360 // noise positive -> cdf ~ 1 -> second transition.
361 assert_eq!(sample_jump_kernel(&jumps, 10.0), -3);
362 // No transitions -> no jump.
363 assert_eq!(sample_jump_kernel(&[], 0.0), 0);
364 // Zero total rate -> no jump.
365 let zero = vec![JumpTransition::new(5, 0.0)];
366 assert_eq!(sample_jump_kernel(&zero, 0.0), 0);
367 }
368}