fix logsumexp and add a test (#11)

* fix logsumexp and add a test

* drop printlns from test

* bump version

* format

* changelog

* add link

Author: cscherrer

CHANGELOG.md

@@ -1,5 +1,10 @@
 # Changelog
 
+## [0.18.1] - 2025-01-07
+
+### Fixed
+- Fixed logsumexp function for -inf values
+
 ## [0.18.0] - 2024-06-24
 
 ### Added
@@ -215,6 +220,7 @@
 - Remove dependency on `quadrature` crate in favor of hand-rolled adaptive
     Simpson's rule, which handles multimodal distributions better.
 
+[0.18.1]: https://github.com/promise-ai/rv/compare/v0.18.0...v0.18.1
 [0.18.0]: https://github.com/promise-ai/rv/compare/v0.17.0...v0.18.0
 [0.17.0]: https://github.com/promise-ai/rv/compare/v0.16.5...v0.17.0
 [0.16.5]: https://github.com/promise-ai/rv/compare/v0.16.4...v0.16.5

Cargo.lock

@@ -1,6 +1,6 @@
 [package]
 name = "rv"
-version = "0.18.0"
+version = "0.18.1"
 authors = ["Baxter Eaves", "Michael Schmidt", "Chad Scherrer"]
 description = "Random variables"
 repository = "https://github.com/promised-ai/rv"

Cargo.toml

@@ -127,7 +127,10 @@ where
         let (alpha, r) =
             self.fold((f64::NEG_INFINITY, 0.0), |(alpha, r), x| {
                 let x = *x.borrow();
-                if x <= alpha {
+
+                if x == f64::NEG_INFINITY {
+                    return (alpha, r);
+                } else if x <= alpha {
                     (alpha, r + (x - alpha).exp())
                 } else {
                     (x, (alpha - x).exp().mul_add(r, 1.0))
@@ -864,6 +867,27 @@ mod tests {
         assert_eq!(argmax(&xs), vec![4, 6]);
     }
 
+    #[test]
+    fn logsumexp_nan_handling() {
+        let a: f64 = -3.0;
+        let b: f64 = -7.0;
+        let target: f64 = logaddexp(a, b);
+        let xs = [
+            -f64::INFINITY,
+            a,
+            -f64::INFINITY,
+            b,
+            -f64::INFINITY,
+            -f64::INFINITY,
+            -f64::INFINITY,
+            -f64::INFINITY,
+            -f64::INFINITY,
+            -f64::INFINITY,
+        ];
+        let result = xs.iter().logsumexp();
+        assert!((result - target).abs() < 1e-12);
+    }
+
     proptest! {
         #[test]
         fn proptest_logsumexp(xs in prop::collection::vec(-1e10_f64..1e10_f64, 0..100)) {