Repository: rust-ndarray/ndarray
Pull request #1491 (2025-03-18 21:41:42)

Add cumprod
# Add cumprod() to match NumPy behavior

Hello ndarray community!

This is implementation of the `cumprod()` function to match NumPy's behavior regarding axis handling.

## Related links
- #1354 
- [numpy doc](https://numpy.org/doc/stable/reference/generated/numpy.cumulative_prod.html)

## Details
1. ~~check `if self.ndim() > 1 && axis.is_none() { panic!(...) }`~~- removed since we now use `Axis` directly
2. ~~`Some(axis)` branch to use input array values instead of result array values~~ - removed since we now use `Axis` directly
3. tests to match NumPy's behavior

## Example Usage
```rust
use ndarray::{array, Axis};

// 1D array - Axis(0) works
let a = array![1, 2, 3, 4];
assert_eq!(a.cumprod(Axis(0)), array![1, 2, 6, 24]);

// 2D array - must specify axis
let b = array![[1, 2], [3, 4]];
assert_eq!(b.cumprod(Axis(0)), array![[1, 2], [3, 8]]);
assert_eq!(b.cumprod(Axis(1)), array![[1, 2], [3, 12]]);
```

## Testing
All existing tests pass,
* 1D array behavior with `Axis(0)`
* N-D array behavior requiring axis specification
* Error cases for invalid axis usage

The implementation maintains consistency with NumPy's behavior, but there might be more elegant way to implement this(especially for the nested loop). please give me any comment/feedback! 

@akern40 what's up! I love the automated CI pipeline, thanks 👍

NewBornRustacean: @nilgoyette thanks for your quick review! I was seeking for what would be a better approach, your comment was very helpful! I tried to vectorize it a bit more with `fold_axis()`
NewBornRustacean: @nilgoyette @akern40  sorry to belated reply! thanks to suggestions 🥇  most of the suggestions are included and passes all the tests. 
NewBornRustacean: > The code is really nice with `accumulate_axis_inplace`! Sorry for the detour. I guess we all learned something ;)

yes indeed! thank you so much 😃 
akern40: Thanks for your persistence! Last thing: I don't think `Copy` should be needed. I am assuming it prompted you for it when you wrote `*curr = *curr * prev`? If so, I think you can either a) use `*curr = curr.clone() * prev`, in which case Rust should elide the clone; or b) call `*curr *= prev`, which would require a `MulAssign` trait bound. I think I'm partial to option (a), since Rust _should_ still just do a multiply-assign operation when possible, and you have fewer trait bounds. But I'm open to ideas.
NewBornRustacean: > Thanks for your persistence! Last thing: I don't think `Copy` should be needed. I am assuming it prompted you for it when you wrote `*curr = *curr * prev`? If so, I think you can either a) use `*curr = curr.clone() * prev`, in which case Rust should elide the clone; or b) call `*curr *= prev`, which would require a `MulAssign` trait bound. I think I'm partial to option (a), since Rust _should_ still just do a multiply-assign operation when possible, and you have fewer trait bounds. But I'm open to ideas.

thanks for the suggestion! 😺 how about this one?  
```rust
    #[track_caller]
    pub fn cumprod(&self, axis: Axis) -> Array<A, D>
    where
        A: Clone + Mul<Output = A> + MulAssign,
        D: Dimension + RemoveAxis,
    {
        if axis.0 >= self.ndim() {
            panic!("axis is out of bounds for array of dimension");
        }

        let mut result = self.to_owned();
        result.accumulate_axis_inplace(axis, |prev, curr| *curr *= prev.clone());
        result
    }
```


akern40: Looks good to me! Thanks again for your work 😄

Patch:
diff --git a/src/numeric/impl_numeric.rs b/src/numeric/impl_numeric.rs
index 27c5687e..ae82a482 100644
--- a/src/numeric/impl_numeric.rs
+++ b/src/numeric/impl_numeric.rs
@@ -10,7 +10,7 @@
 use num_traits::Float;
 use num_traits::One;
 use num_traits::{FromPrimitive, Zero};
-use std::ops::{Add, Div, Mul, Sub};
+use std::ops::{Add, Div, Mul, MulAssign, Sub};
 
 use crate::imp_prelude::*;
 use crate::numeric_util;
@@ -97,6 +97,45 @@ where D: Dimension
         sum
     }
 
+    /// Return the cumulative product of elements along a given axis.
+    ///
+    /// ```
+    /// use ndarray::{arr2, Axis};
+    ///
+    /// let a = arr2(&[[1., 2., 3.],
+    ///                [4., 5., 6.]]);
+    ///
+    /// // Cumulative product along rows (axis 0)
+    /// assert_eq!(
+    ///     a.cumprod(Axis(0)),
+    ///     arr2(&[[1., 2., 3.],
+    ///           [4., 10., 18.]])
+    /// );
+    ///
+    /// // Cumulative product along columns (axis 1)
+    /// assert_eq!(
+    ///     a.cumprod(Axis(1)),
+    ///     arr2(&[[1., 2., 6.],
+    ///           [4., 20., 120.]])
+    /// );
+    /// ```
+    ///
+    /// **Panics** if `axis` is out of bounds.
+    #[track_caller]
+    pub fn cumprod(&self, axis: Axis) -> Array<A, D>
+    where
+        A: Clone + Mul<Output = A> + MulAssign,
+        D: Dimension + RemoveAxis,
+    {
+        if axis.0 >= self.ndim() {
+            panic!("axis is out of bounds for array of dimension");
+        }
+
+        let mut result = self.to_owned();
+        result.accumulate_axis_inplace(axis, |prev, curr| *curr *= prev.clone());
+        result
+    }
+
     /// Return variance of elements in the array.
     ///
     /// The variance is computed using the [Welford one-pass


Test patch:
diff --git a/tests/numeric.rs b/tests/numeric.rs
index 839aba58..7e696481 100644
--- a/tests/numeric.rs
+++ b/tests/numeric.rs
@@ -75,6 +75,76 @@ fn sum_mean_prod_empty()
     assert_eq!(a, None);
 }
 
+#[test]
+fn test_cumprod_1d()
+{
+    let a = array![1, 2, 3, 4];
+    let result = a.cumprod(Axis(0));
+    assert_eq!(result, array![1, 2, 6, 24]);
+}
+
+#[test]
+fn test_cumprod_2d()
+{
+    let a = array![[1, 2], [3, 4]];
+
+    let result_axis0 = a.cumprod(Axis(0));
+    assert_eq!(result_axis0, array![[1, 2], [3, 8]]);
+
+    let result_axis1 = a.cumprod(Axis(1));
+    assert_eq!(result_axis1, array![[1, 2], [3, 12]]);
+}
+
+#[test]
+fn test_cumprod_3d()
+{
+    let a = array![[[1, 2], [3, 4]], [[5, 6], [7, 8]]];
+
+    let result_axis0 = a.cumprod(Axis(0));
+    assert_eq!(result_axis0, array![[[1, 2], [3, 4]], [[5, 12], [21, 32]]]);
+
+    let result_axis1 = a.cumprod(Axis(1));
+    assert_eq!(result_axis1, array![[[1, 2], [3, 8]], [[5, 6], [35, 48]]]);
+
+    let result_axis2 = a.cumprod(Axis(2));
+    assert_eq!(result_axis2, array![[[1, 2], [3, 12]], [[5, 30], [7, 56]]]);
+}
+
+#[test]
+fn test_cumprod_empty()
+{
+    // For 2D empty array
+    let b: Array2<i32> = Array2::zeros((0, 0));
+    let result_axis0 = b.cumprod(Axis(0));
+    assert_eq!(result_axis0, Array2::zeros((0, 0)));
+    let result_axis1 = b.cumprod(Axis(1));
+    assert_eq!(result_axis1, Array2::zeros((0, 0)));
+}
+
+#[test]
+fn test_cumprod_1_element()
+{
+    // For 1D array with one element
+    let a = array![5];
+    let result = a.cumprod(Axis(0));
+    assert_eq!(result, array![5]);
+
+    // For 2D array with one element
+    let b = array![[5]];
+    let result_axis0 = b.cumprod(Axis(0));
+    let result_axis1 = b.cumprod(Axis(1));
+    assert_eq!(result_axis0, array![[5]]);
+    assert_eq!(result_axis1, array![[5]]);
+}
+
+#[test]
+#[should_panic(expected = "axis is out of bounds for array of dimension")]
+fn test_cumprod_axis_out_of_bounds()
+{
+    let a = array![[1, 2], [3, 4]];
+    let _result = a.cumprod(Axis(2));
+}
+
 #[test]
 #[cfg(feature = "std")]
 fn var()