xsf 0.5.0+0.2.0

Rust implementations and bindings for scipy.special functions.
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
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139

//! <https://github.com/numpy/numpy/blob/4992eaf/numpy/_core/src/npymath/npy_math_internal.h.src>

use core::f64::consts::{LN_2, LOG2_E};

mod sealed {
    pub trait Sealed {}
    impl Sealed for f32 {}
    impl Sealed for f64 {}
}

pub trait LogAddExpArg: sealed::Sealed {
    const ZERO: Self;
    fn npy_logaddexp(self, other: Self) -> Self;
    fn npy_logaddexp2(self, other: Self) -> Self;
}

impl LogAddExpArg for f32 {
    const ZERO: Self = 0.0;

    #[inline]
    #[allow(clippy::cast_possible_truncation)]
    fn npy_logaddexp(self, other: Self) -> Self {
        f64::from(self).npy_logaddexp(f64::from(other)) as f32
    }

    #[inline]
    #[allow(clippy::cast_possible_truncation)]
    fn npy_logaddexp2(self, other: Self) -> Self {
        f64::from(self).npy_logaddexp2(f64::from(other)) as f32
    }
}

impl LogAddExpArg for f64 {
    const ZERO: Self = 0.0;

    #[inline]
    fn npy_logaddexp(self, other: Self) -> Self {
        #[allow(clippy::float_cmp)]
        if self == other {
            // Handles infinities of the same sign without warnings
            self + LN_2
        } else {
            let tmp = self - other;
            if tmp > 0.0 {
                self + (-tmp).exp().ln_1p()
            } else if tmp <= 0.0 {
                other + tmp.exp().ln_1p()
            } else {
                // NaNs
                tmp
            }
        }
    }

    #[inline]
    fn npy_logaddexp2(self, other: Self) -> Self {
        #[allow(clippy::float_cmp)]
        if self == other {
            // Handles infinities of the same sign without warnings
            self + 1.0
        } else {
            let tmp = self - other;
            if tmp > 0.0 {
                self + (-tmp).exp2().ln_1p() * LOG2_E
            } else if tmp <= 0.0 {
                other + tmp.exp2().ln_1p() * LOG2_E
            } else {
                // NaNs
                tmp
            }
        }
    }
}

/// ln(e<sup>x</sup> + e<sup>y</sup>)
///
/// This function is useful in statistics where the calculated probabilities of events may be so
/// small as to exceed the range of normal floating point numbers. In such cases the logarithm of
/// the calculated probability is stored. This function allows adding probabilities stored in such
/// a fashion.
///
/// # See also
/// - [`logaddexp2`]
/// - [NumPy documentation](https://numpy.org/doc/stable/reference/generated/numpy.logaddexp.html)
///
#[inline]
pub fn logaddexp<T: LogAddExpArg>(x: T, y: T) -> T {
    x.npy_logaddexp(y)
}

/// log<sub>2</sub>(2<sup>x</sup> + 2<sup>y</sup>)
///
/// This function is useful in machine learning when the calculated probabilities of events may be
/// so small as to exceed the range of normal floating point numbers. In such cases the base-2
/// logarithm of the calculated probability can be used instead. This function allows adding
/// probabilities stored in such a fashion.
///
/// # See also
/// - [`logaddexp`]
/// - [NumPy documentation](https://numpy.org/doc/stable/reference/generated/numpy.logaddexp2.html)
///
#[inline]
pub fn logaddexp2<T: LogAddExpArg>(x: T, y: T) -> T {
    x.npy_logaddexp2(y)
}

#[cfg(test)]
mod tests {
    use crate::np_assert_allclose;

    #[test]
    fn test_logaddexp_f32() {
        let x = [1.0, 2.0, 3.0, 4.0, 5.0];
        let y = [5.0, 4.0, 3.0, 2.0, 1.0];
        let z = [6.0, 6.0, 6.0, 6.0, 6.0];

        let xf = x.map(f64::log2);
        let yf = y.map(f64::log2);
        let zf = z.map(f64::log2);

        #[allow(clippy::cast_possible_truncation)]
        let zr: [f32; 5] = std::array::from_fn(|i| crate::logaddexp2(xf[i] as f32, yf[i] as f32));
        np_assert_allclose!(zr.map(f64::from), zf, atol = 1.5e-7);
    }

    #[test]
    fn test_logaddexp_f64() {
        let x = [1.0, 2.0, 3.0, 4.0, 5.0];
        let y = [5.0, 4.0, 3.0, 2.0, 1.0];
        let z = [6.0, 6.0, 6.0, 6.0, 6.0];

        let xf = x.map(f64::log2);
        let yf = y.map(f64::log2);
        let zf = z.map(f64::log2);

        let zr: [f64; 5] = std::array::from_fn(|i| crate::logaddexp2(xf[i], yf[i]));
        np_assert_allclose!(zr, zf, atol = 1e-15);
    }
}