comnum_math/
log_domain.rs

1//! The log domain of probabilities and LLRs, without overflow: [`ln_1p_exp`],
2//! [`max_star`] and its choice [`Jacobian`] (of `comnum-core`), [`log_sum_exp`] and [`log_sum_exp_iter`]
3//! (re-exported at the root).
4
5/// `ln(1 + e^x)` (the softplus) without overflow.
6///
7/// # Examples
8/// ```
9/// use comnum_math::ln_1p_exp;
10/// assert!((ln_1p_exp(0.0) - 2f64.ln()).abs() < 1e-15);
11/// assert_eq!(ln_1p_exp(1000.0), 1000.0);
12/// ```
13#[inline]
14pub fn ln_1p_exp(x: f64) -> f64 {
15    x.max(0.0) + (-x.abs()).exp().ln_1p()
16}
17
18/// The Jacobian logarithm, exact or max-log: defined in `comnum-core`, since the soft
19/// demapper trait [`comnum_core::traits::Demap`] takes it, and re-exported here for
20/// [`max_star`], the BCJR and turbo decoders of `comnum-fec` and the demappers of
21/// `comnum-mapping`.
22pub use comnum_core::traits::Jacobian;
23
24/// The `max*` operation `ln(e^a + e^b)` of the log-MAP algorithms; `max(a, b)` with
25/// [`Jacobian::MaxLog`] or when the maximum is not finite.
26///
27/// # Examples
28/// ```
29/// use comnum_math::{Jacobian, max_star};
30/// assert!((max_star(0.0, 0.0, Jacobian::Exact) - 2f64.ln()).abs() < 1e-15);
31/// assert_eq!(max_star(0.0, -1.0, Jacobian::MaxLog), 0.0);
32/// ```
33#[inline]
34pub fn max_star(a: f64, b: f64, jacobian: Jacobian) -> f64 {
35    let m = a.max(b);
36    if !jacobian.is_exact() || !m.is_finite() { m } else { m + (-(a - b).abs()).exp().ln_1p() }
37}
38
39/// `ln Σ_i e^{v_i}` without overflow (the log domain of probabilities and LLRs).
40///
41/// # Examples
42/// ```
43/// use comnum_math::log_sum_exp;
44/// assert!((log_sum_exp(&[1000.0, 1000.0]) - (1000.0 + 2f64.ln())).abs() < 1e-9);
45/// assert_eq!(log_sum_exp(&[]), f64::NEG_INFINITY);
46/// ```
47pub fn log_sum_exp(v: &[f64]) -> f64 {
48    log_sum_exp_iter(v.iter().copied())
49}
50
51/// [`log_sum_exp`] of the values of an iterator, traversed twice (no allocation).
52///
53/// # Examples
54/// ```
55/// use comnum_math::{log_sum_exp, log_sum_exp_iter};
56/// let v = [0.5, -2.0, 3.0, 1.0];
57/// // the terms of even index only
58/// let even = log_sum_exp_iter(v.iter().step_by(2).copied());
59/// assert_eq!(even, log_sum_exp(&[0.5, 3.0]));
60/// ```
61pub fn log_sum_exp_iter<I: Iterator<Item = f64> + Clone>(v: I) -> f64 {
62    let m = v.clone().fold(f64::NEG_INFINITY, f64::max);
63    if !m.is_finite() {
64        return m;
65    }
66    m + v.map(|x| (x - m).exp()).sum::<f64>().ln()
67}
comnum_math/log_domain.rs

comnum_math/
log_domain.rs
