diff --git a/Cargo.lock b/Cargo.lock index e91fe100e..8ea8bb5ca 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -4709,6 +4709,7 @@ dependencies = [ "digital-muon-common", "digital-muon-streaming-types", "git-version", + "libm", "miette", "num", "rand 0.9.2", diff --git a/simulator/Cargo.toml b/simulator/Cargo.toml index 9ca2e2738..a3837f1b0 100644 --- a/simulator/Cargo.toml +++ b/simulator/Cargo.toml @@ -8,6 +8,7 @@ edition.workspace = true chrono.workspace = true clap.workspace = true git-version.workspace = true +libm.workspace = true miette = { workspace = true, features = ["fancy"] } num.workspace = true rand.workspace = true diff --git a/simulator/src/integrated/simulation.rs b/simulator/src/integrated/simulation.rs index 6128cbf9c..fffced854 100644 --- a/simulator/src/integrated/simulation.rs +++ b/simulator/src/integrated/simulation.rs @@ -149,11 +149,12 @@ mod tests { } }, "pulses": [{ - "pulse-type": "biexp", - "height": { "random-type": "uniform-float", "min": { "const": 30 }, "max": { "const": 70 } }, - "start": { "random-type": "exponential", "lifetime": { "const": 2200 } }, - "rise": { "random-type": "uniform-float", "min": { "const": 20 }, "max": { "const": 30 } }, - "decay": { "random-type": "uniform-float", "min": { "const": 5 }, "max": { "const": 10 } } + "pulse-type": "back-to-back-exp", + "spread": { "random-type": "constant-float", "value": { "const": 5.5 } }, + "rising": { "random-type": "constant-float", "value": { "const": 3.5 } }, + "falling": { "random-type": "constant-float", "value": { "const": 2.4 } }, + "peak_time": { "random-type": "exponential", "lifetime": { "const": 2200 } }, + "peak_height": { "random-type": "uniform-float", "min": { "const": 250 }, "max": { "const": 1100 } } }, { "pulse-type": "flat", diff --git a/simulator/src/integrated/simulation_elements/pulses.rs b/simulator/src/integrated/simulation_elements/pulses.rs index 7d3357c13..ae11dc755 100644 --- a/simulator/src/integrated/simulation_elements/pulses.rs +++ b/simulator/src/integrated/simulation_elements/pulses.rs @@ -1,3 +1,5 @@ +use core::f64; + use super::{FloatRandomDistribution, utils::JsonValueError}; use digital_muon_common::{Intensity, Time}; use serde::Deserialize; @@ -21,11 +23,12 @@ pub(crate) enum PulseTemplate { peak_time: FloatRandomDistribution, sd: FloatRandomDistribution, }, - Biexp { - start: FloatRandomDistribution, - decay: FloatRandomDistribution, - rise: FloatRandomDistribution, - height: FloatRandomDistribution, + BackToBackExp { + peak_height: FloatRandomDistribution, + peak_time: FloatRandomDistribution, + spread: FloatRandomDistribution, + falling: FloatRandomDistribution, + rising: FloatRandomDistribution, }, } @@ -49,14 +52,16 @@ pub(crate) enum PulseEvent { sd: f64, peak_amplitude: f64, }, - Biexp { + BackToBackExp { start: f64, stop: f64, - decay: f64, - rise: f64, - peak_height: f64, - coef: f64, peak_time: f64, + falling: f64, + rising: f64, + normalising_factor: f64, + rising_spread: f64, + falling_spread: f64, + frac_1_sqrt_2_spread: f64, }, } @@ -97,37 +102,63 @@ impl PulseEvent { } => { let mean = peak_time.sample(frame)?; let sd = sd.sample(frame)?; + let peak_amplitude = height.sample(frame)?; + let distance_to_value_of_one = 2.0 * sd * peak_amplitude.ln().sqrt(); Ok(Self::Gaussian { - start: mean - 4.0 * sd, - stop: mean + 4.0 * sd, + start: mean - distance_to_value_of_one, + stop: mean + distance_to_value_of_one, mean, sd, - peak_amplitude: height.sample(frame)?, + peak_amplitude, }) } - PulseTemplate::Biexp { - start, - decay, - rise, - height, + PulseTemplate::BackToBackExp { + peak_height, + peak_time, + spread, + falling, + rising, } => { - let start = start.sample(frame)?; - let decay = decay.sample(frame)?; - let rise = rise.sample(frame)?; - let peak_height = height.sample(frame)?; - let ratio = decay / rise; - let coef = peak_height - / (f64::powf(ratio, 1.0 / ratio - 1.0) - f64::powf(ratio, 1.0 - ratio)); - let peak_time = f64::ln(f64::powf(ratio, decay * rise / (decay - rise))); - let stop = f64::MAX; // This needs to be changed - Ok(Self::Biexp { + let rising = rising.sample(frame)?; + let falling = falling.sample(frame)?; + let peak_height = peak_height.sample(frame)?; + let spread = spread.sample(frame)?; + let peak_time = peak_time.sample(frame)?; + + let rising_spread = rising * spread.powi(2); + let falling_spread = falling * spread.powi(2); + let frac_1_sqrt_2_spread = std::f64::consts::FRAC_1_SQRT_2 / spread; + + let normalising_factor = { + let rising_erfc = libm::erfc(rising_spread * frac_1_sqrt_2_spread); + let rising_exp = if rising_erfc == 0.0 { + 0.0 + } else { + f64::exp(0.5 * rising * rising_spread) + }; + let falling_erfc = libm::erfc(falling_spread * frac_1_sqrt_2_spread); + let falling_exp = if falling_erfc == 0.0 { + 0.0 + } else { + f64::exp(0.5 * falling * falling_spread) + }; + + peak_height / (rising_exp * rising_erfc + falling_exp * falling_erfc) + }; + + let start = peak_time - 0.5 * rising_spread - normalising_factor.ln() / rising; + let stop = peak_time + 0.5 * falling_spread + normalising_factor.ln() / falling; + + Ok(Self::BackToBackExp { start, stop, - decay, - rise, - peak_height, - coef, peak_time, + falling, + rising, + normalising_factor, + rising_spread, + falling_spread, + frac_1_sqrt_2_spread, }) } } @@ -138,7 +169,7 @@ impl PulseEvent { Self::Flat { start, .. } => *start, Self::Triangular { start, .. } => *start, Self::Gaussian { start, .. } => *start, - Self::Biexp { start, .. } => *start, + Self::BackToBackExp { start, .. } => *start, }) as Time } @@ -147,7 +178,7 @@ impl PulseEvent { Self::Flat { stop, .. } => *stop, Self::Triangular { stop, .. } => *stop, Self::Gaussian { stop, .. } => *stop, - Self::Biexp { stop, .. } => *stop, + Self::BackToBackExp { stop, .. } => *stop, }) as Time } @@ -156,46 +187,51 @@ impl PulseEvent { Self::Flat { start, .. } => *start, Self::Triangular { peak_time, .. } => *peak_time, Self::Gaussian { mean, .. } => *mean, - Self::Biexp { - start, peak_time, .. - } => *start + *peak_time / 2.0, + Self::BackToBackExp { peak_time, .. } => *peak_time, }) as Time } pub(crate) fn intensity(&self) -> Intensity { - *match self { - Self::Flat { amplitude, .. } => amplitude, - Self::Triangular { amplitude, .. } => amplitude, - Self::Gaussian { peak_amplitude, .. } => peak_amplitude, - Self::Biexp { peak_height, .. } => peak_height, - } as Intensity + (match self { + Self::Flat { amplitude, .. } => *amplitude, + Self::Triangular { amplitude, .. } => *amplitude, + Self::Gaussian { peak_amplitude, .. } => *peak_amplitude, + Self::BackToBackExp { + falling, + rising, + normalising_factor, + rising_spread, + falling_spread, + frac_1_sqrt_2_spread, + .. + } => { + let rising_exp = f64::exp(rising * (0.5 * rising_spread)); + let rising_erfc = libm::erfc(rising_spread * frac_1_sqrt_2_spread); + let falling_exp = f64::exp(falling * (0.5 * falling_spread)); + let falling_erfc = libm::erfc(falling_spread * frac_1_sqrt_2_spread); + + normalising_factor * (rising_exp * rising_erfc + falling_exp * falling_erfc) + } + }) as Intensity } pub(crate) fn get_value_at(&self, time: f64) -> f64 { + if self.get_start() > time as Time || time as Time > self.get_end() { + return Default::default(); + } + match *self { - Self::Flat { - start, - stop, - amplitude, - } => { - if start <= time && time < stop { - amplitude - } else { - f64::default() - } - } + Self::Flat { amplitude, .. } => amplitude, Self::Triangular { start, peak_time, stop, amplitude, } => { - if start <= time && time < peak_time { + if time < peak_time { amplitude * (time - start) / (peak_time - start) - } else if peak_time <= time && time < stop { - amplitude * (stop - time) / (stop - peak_time) } else { - f64::default() + amplitude * (stop - time) / (stop - peak_time) } } Self::Gaussian { @@ -203,27 +239,152 @@ impl PulseEvent { sd, peak_amplitude, .. - } => { - if mean - 6.0 * sd > time || time > mean + 6.0 * sd { - f64::default() - } else { - peak_amplitude * f64::exp(-f64::powi(0.5 * (time - mean) / sd, 2)) - } - } - Self::Biexp { - start, - decay, - rise, - coef, + } => peak_amplitude * f64::exp(-f64::powi(0.5 * (time - mean) / sd, 2)), + Self::BackToBackExp { + peak_time, + falling, + rising, + normalising_factor, + rising_spread, + falling_spread, + frac_1_sqrt_2_spread, .. } => { - if time < start { - f64::default() + let time_shift = time - peak_time; // -2.9539084192897 + + let rising_erfc = libm::erfc((rising_spread + time_shift) * frac_1_sqrt_2_spread); + // Guard against rising_exp being NaN or Inf + let rising_exp = if rising_erfc >= f64::EPSILON { + f64::exp(rising * (0.5 * rising_spread + time_shift)) } else { - let time = time - start; - coef * (f64::exp(-time / decay) - f64::exp(-time / rise)) - } + Default::default() + }; + + let falling_erfc = libm::erfc((falling_spread - time_shift) * frac_1_sqrt_2_spread); + // Guard against falling_exp being NaN or Inf + let falling_exp = if falling_erfc >= f64::EPSILON { + f64::exp(falling * (0.5 * falling_spread - time_shift)) + } else { + Default::default() + }; + + normalising_factor * (rising_exp * rising_erfc + falling_exp * falling_erfc) } } } } + +#[cfg(test)] +mod tests { + use crate::integrated::simulation_elements::NumExpression; + + use super::*; + + const TEMPLATE: PulseTemplate = PulseTemplate::BackToBackExp { + peak_height: FloatRandomDistribution::ConstantFloat { + value: NumExpression::Const(2100.0), + }, + peak_time: FloatRandomDistribution::ConstantFloat { + value: NumExpression::Const(2200.0), + }, + spread: FloatRandomDistribution::ConstantFloat { + value: NumExpression::Const(3.0), + }, + falling: FloatRandomDistribution::ConstantFloat { + value: NumExpression::Const(2.5), + }, + rising: FloatRandomDistribution::ConstantFloat { + value: NumExpression::Const(1.5), + }, + }; + + #[test] + fn back_to_back_exp_template() { + let pulse = PulseEvent::sample(&TEMPLATE, 0); + assert!(pulse.is_ok()); + let pulse = pulse.unwrap(); + assert_eq!(pulse.get_start(), 2187); + assert_eq!(pulse.get_end(), 2214); + assert_eq!(pulse.intensity(), 2100); + assert_eq!(pulse.time(), 2200); + } + + #[test] + fn back_to_back_exp_values() { + let pulse = PulseEvent::sample(&TEMPLATE, 0).unwrap(); + const VALUES: [Intensity; 27] = [ + 0, 1, 5, 14, 35, 78, 159, 292, 487, 730, 988, 1793, 2044, 2100, 1942, 1616, 1211, 816, + 495, 270, 132, 58, 23, 8, 2, 0, 0, + ]; + for (t, &v) in VALUES.iter().enumerate() { + assert_eq!( + pulse.get_value_at((pulse.get_start() + t as Time) as f64) as Intensity, + v + ); + } + } + + #[test] + fn back_to_back_exp_extreme_rising() { + let template = PulseTemplate::BackToBackExp { + peak_height: FloatRandomDistribution::ConstantFloat { + value: NumExpression::Const(2100.0), + }, + peak_time: FloatRandomDistribution::ConstantFloat { + value: NumExpression::Const(2200.0), + }, + spread: FloatRandomDistribution::ConstantFloat { + value: NumExpression::Const(3.0), + }, + falling: FloatRandomDistribution::ConstantFloat { + value: NumExpression::Const(12.5), + }, + rising: FloatRandomDistribution::ConstantFloat { + value: NumExpression::Const(1.5), + }, + }; + let event = PulseEvent::sample(&template, 0).unwrap(); + const VALUES1: [Intensity; 27] = [ + 0, 2, 7, 20, 50, 112, 228, 420, 699, 1049, 1419, 1728, 1893, 1866, 1652, 1314, 938, + 601, 346, 178, 82, 34, 12, 4, 1, 0, 0, + ]; + for (t, &v) in VALUES1.iter().enumerate() { + assert_eq!( + event.get_value_at(event.get_start() as f64 + t as f64) as Intensity, + v + ); + } + } + + #[test] + fn back_to_back_exp_extreme_falling() { + let template = PulseTemplate::BackToBackExp { + peak_height: FloatRandomDistribution::ConstantFloat { + value: NumExpression::Const(2100.0), + }, + peak_time: FloatRandomDistribution::ConstantFloat { + value: NumExpression::Const(2200.0), + }, + spread: FloatRandomDistribution::ConstantFloat { + value: NumExpression::Const(4.0), + }, + falling: FloatRandomDistribution::ConstantFloat { + value: NumExpression::Const(12.5), + }, + rising: FloatRandomDistribution::ConstantFloat { + value: NumExpression::Const(1.5), + }, + }; + let event = PulseEvent::sample(&template, 0).unwrap(); + const VALUES2: [Intensity; 27] = [ + 0, 0, 0, 1, 4, 10, 21, 43, 83, 151, 258, 414, 626, 892, 1197, 1512, 1798, 2012, 2119, + 2100, 1957, 1716, 1416, 1098, 801, 550, 355, + ]; + for (t, &v) in VALUES2.iter().enumerate() { + assert_eq!( + event.get_value_at(event.get_start() as f64 + t as f64) as Intensity, + v + ); + } + } +}