ndarray/numeric/impl_numeric.rs
1// Copyright 2014-2016 bluss and ndarray developers.
2//
3// Licensed under the Apache License, Version 2.0 <LICENSE-APACHE or
4// http://www.apache.org/licenses/LICENSE-2.0> or the MIT license
5// <LICENSE-MIT or http://opensource.org/licenses/MIT>, at your
6// option. This file may not be copied, modified, or distributed
7// except according to those terms.
8
9#[cfg(feature = "std")]
10use num_traits::Float;
11use num_traits::{self, FromPrimitive, Zero};
12use std::ops::{Add, Div, Mul};
13
14use crate::imp_prelude::*;
15use crate::numeric_util;
16
17/// # Numerical Methods for Arrays
18impl<A, S, D> ArrayBase<S, D>
19where
20 S: Data<Elem = A>,
21 D: Dimension,
22{
23 /// Return the sum of all elements in the array.
24 ///
25 /// ```
26 /// use ndarray::arr2;
27 ///
28 /// let a = arr2(&[[1., 2.],
29 /// [3., 4.]]);
30 /// assert_eq!(a.sum(), 10.);
31 /// ```
32 pub fn sum(&self) -> A
33 where
34 A: Clone + Add<Output = A> + num_traits::Zero,
35 {
36 if let Some(slc) = self.as_slice_memory_order() {
37 return numeric_util::unrolled_fold(slc, A::zero, A::add);
38 }
39 let mut sum = A::zero();
40 for row in self.rows() {
41 if let Some(slc) = row.as_slice() {
42 sum = sum + numeric_util::unrolled_fold(slc, A::zero, A::add);
43 } else {
44 sum = sum + row.iter().fold(A::zero(), |acc, elt| acc + elt.clone());
45 }
46 }
47 sum
48 }
49
50 /// Return the sum of all elements in the array.
51 ///
52 /// *This method has been renamed to `.sum()`*
53 #[deprecated(note="renamed to `sum`", since="0.15.0")]
54 pub fn scalar_sum(&self) -> A
55 where
56 A: Clone + Add<Output = A> + num_traits::Zero,
57 {
58 self.sum()
59 }
60
61 /// Returns the [arithmetic mean] x̅ of all elements in the array:
62 ///
63 /// ```text
64 /// 1 n
65 /// x̅ = ― ∑ xᵢ
66 /// n i=1
67 /// ```
68 ///
69 /// If the array is empty, `None` is returned.
70 ///
71 /// **Panics** if `A::from_usize()` fails to convert the number of elements in the array.
72 ///
73 /// [arithmetic mean]: https://en.wikipedia.org/wiki/Arithmetic_mean
74 pub fn mean(&self) -> Option<A>
75 where
76 A: Clone + FromPrimitive + Add<Output = A> + Div<Output = A> + Zero,
77 {
78 let n_elements = self.len();
79 if n_elements == 0 {
80 None
81 } else {
82 let n_elements = A::from_usize(n_elements)
83 .expect("Converting number of elements to `A` must not fail.");
84 Some(self.sum() / n_elements)
85 }
86 }
87
88 /// Return the product of all elements in the array.
89 ///
90 /// ```
91 /// use ndarray::arr2;
92 ///
93 /// let a = arr2(&[[1., 2.],
94 /// [3., 4.]]);
95 /// assert_eq!(a.product(), 24.);
96 /// ```
97 pub fn product(&self) -> A
98 where
99 A: Clone + Mul<Output = A> + num_traits::One,
100 {
101 if let Some(slc) = self.as_slice_memory_order() {
102 return numeric_util::unrolled_fold(slc, A::one, A::mul);
103 }
104 let mut sum = A::one();
105 for row in self.rows() {
106 if let Some(slc) = row.as_slice() {
107 sum = sum * numeric_util::unrolled_fold(slc, A::one, A::mul);
108 } else {
109 sum = sum * row.iter().fold(A::one(), |acc, elt| acc * elt.clone());
110 }
111 }
112 sum
113 }
114
115 /// Return variance of elements in the array.
116 ///
117 /// The variance is computed using the [Welford one-pass
118 /// algorithm](https://www.jstor.org/stable/1266577).
119 ///
120 /// The parameter `ddof` specifies the "delta degrees of freedom". For
121 /// example, to calculate the population variance, use `ddof = 0`, or to
122 /// calculate the sample variance, use `ddof = 1`.
123 ///
124 /// The variance is defined as:
125 ///
126 /// ```text
127 /// 1 n
128 /// variance = ―――――――― ∑ (xᵢ - x̅)²
129 /// n - ddof i=1
130 /// ```
131 ///
132 /// where
133 ///
134 /// ```text
135 /// 1 n
136 /// x̅ = ― ∑ xᵢ
137 /// n i=1
138 /// ```
139 ///
140 /// and `n` is the length of the array.
141 ///
142 /// **Panics** if `ddof` is less than zero or greater than `n`
143 ///
144 /// # Example
145 ///
146 /// ```
147 /// use ndarray::array;
148 /// use approx::assert_abs_diff_eq;
149 ///
150 /// let a = array![1., -4.32, 1.14, 0.32];
151 /// let var = a.var(1.);
152 /// assert_abs_diff_eq!(var, 6.7331, epsilon = 1e-4);
153 /// ```
154 #[cfg(feature = "std")]
155 pub fn var(&self, ddof: A) -> A
156 where
157 A: Float + FromPrimitive,
158 {
159 let zero = A::from_usize(0).expect("Converting 0 to `A` must not fail.");
160 let n = A::from_usize(self.len()).expect("Converting length to `A` must not fail.");
161 assert!(
162 !(ddof < zero || ddof > n),
163 "`ddof` must not be less than zero or greater than the length of \
164 the axis",
165 );
166 let dof = n - ddof;
167 let mut mean = A::zero();
168 let mut sum_sq = A::zero();
169 let mut i = 0;
170 self.for_each(|&x| {
171 let count = A::from_usize(i + 1).expect("Converting index to `A` must not fail.");
172 let delta = x - mean;
173 mean = mean + delta / count;
174 sum_sq = (x - mean).mul_add(delta, sum_sq);
175 i += 1;
176 });
177 sum_sq / dof
178 }
179
180 /// Return standard deviation of elements in the array.
181 ///
182 /// The standard deviation is computed from the variance using
183 /// the [Welford one-pass algorithm](https://www.jstor.org/stable/1266577).
184 ///
185 /// The parameter `ddof` specifies the "delta degrees of freedom". For
186 /// example, to calculate the population standard deviation, use `ddof = 0`,
187 /// or to calculate the sample standard deviation, use `ddof = 1`.
188 ///
189 /// The standard deviation is defined as:
190 ///
191 /// ```text
192 /// ⎛ 1 n ⎞
193 /// stddev = sqrt ⎜ ―――――――― ∑ (xᵢ - x̅)²⎟
194 /// ⎝ n - ddof i=1 ⎠
195 /// ```
196 ///
197 /// where
198 ///
199 /// ```text
200 /// 1 n
201 /// x̅ = ― ∑ xᵢ
202 /// n i=1
203 /// ```
204 ///
205 /// and `n` is the length of the array.
206 ///
207 /// **Panics** if `ddof` is less than zero or greater than `n`
208 ///
209 /// # Example
210 ///
211 /// ```
212 /// use ndarray::array;
213 /// use approx::assert_abs_diff_eq;
214 ///
215 /// let a = array![1., -4.32, 1.14, 0.32];
216 /// let stddev = a.std(1.);
217 /// assert_abs_diff_eq!(stddev, 2.59483, epsilon = 1e-4);
218 /// ```
219 #[cfg(feature = "std")]
220 pub fn std(&self, ddof: A) -> A
221 where
222 A: Float + FromPrimitive,
223 {
224 self.var(ddof).sqrt()
225 }
226
227 /// Return sum along `axis`.
228 ///
229 /// ```
230 /// use ndarray::{aview0, aview1, arr2, Axis};
231 ///
232 /// let a = arr2(&[[1., 2., 3.],
233 /// [4., 5., 6.]]);
234 /// assert!(
235 /// a.sum_axis(Axis(0)) == aview1(&[5., 7., 9.]) &&
236 /// a.sum_axis(Axis(1)) == aview1(&[6., 15.]) &&
237 ///
238 /// a.sum_axis(Axis(0)).sum_axis(Axis(0)) == aview0(&21.)
239 /// );
240 /// ```
241 ///
242 /// **Panics** if `axis` is out of bounds.
243 pub fn sum_axis(&self, axis: Axis) -> Array<A, D::Smaller>
244 where
245 A: Clone + Zero + Add<Output = A>,
246 D: RemoveAxis,
247 {
248 let min_stride_axis = self.dim.min_stride_axis(&self.strides);
249 if axis == min_stride_axis {
250 crate::Zip::from(self.lanes(axis)).map_collect(|lane| lane.sum())
251 } else {
252 let mut res = Array::zeros(self.raw_dim().remove_axis(axis));
253 for subview in self.axis_iter(axis) {
254 res = res + &subview;
255 }
256 res
257 }
258 }
259
260 /// Return mean along `axis`.
261 ///
262 /// Return `None` if the length of the axis is zero.
263 ///
264 /// **Panics** if `axis` is out of bounds or if `A::from_usize()`
265 /// fails for the axis length.
266 ///
267 /// ```
268 /// use ndarray::{aview0, aview1, arr2, Axis};
269 ///
270 /// let a = arr2(&[[1., 2., 3.],
271 /// [4., 5., 6.]]);
272 /// assert!(
273 /// a.mean_axis(Axis(0)).unwrap() == aview1(&[2.5, 3.5, 4.5]) &&
274 /// a.mean_axis(Axis(1)).unwrap() == aview1(&[2., 5.]) &&
275 ///
276 /// a.mean_axis(Axis(0)).unwrap().mean_axis(Axis(0)).unwrap() == aview0(&3.5)
277 /// );
278 /// ```
279 pub fn mean_axis(&self, axis: Axis) -> Option<Array<A, D::Smaller>>
280 where
281 A: Clone + Zero + FromPrimitive + Add<Output = A> + Div<Output = A>,
282 D: RemoveAxis,
283 {
284 let axis_length = self.len_of(axis);
285 if axis_length == 0 {
286 None
287 } else {
288 let axis_length =
289 A::from_usize(axis_length).expect("Converting axis length to `A` must not fail.");
290 let sum = self.sum_axis(axis);
291 Some(sum / aview0(&axis_length))
292 }
293 }
294
295 /// Return variance along `axis`.
296 ///
297 /// The variance is computed using the [Welford one-pass
298 /// algorithm](https://www.jstor.org/stable/1266577).
299 ///
300 /// The parameter `ddof` specifies the "delta degrees of freedom". For
301 /// example, to calculate the population variance, use `ddof = 0`, or to
302 /// calculate the sample variance, use `ddof = 1`.
303 ///
304 /// The variance is defined as:
305 ///
306 /// ```text
307 /// 1 n
308 /// variance = ―――――――― ∑ (xᵢ - x̅)²
309 /// n - ddof i=1
310 /// ```
311 ///
312 /// where
313 ///
314 /// ```text
315 /// 1 n
316 /// x̅ = ― ∑ xᵢ
317 /// n i=1
318 /// ```
319 ///
320 /// and `n` is the length of the axis.
321 ///
322 /// **Panics** if `ddof` is less than zero or greater than `n`, if `axis`
323 /// is out of bounds, or if `A::from_usize()` fails for any any of the
324 /// numbers in the range `0..=n`.
325 ///
326 /// # Example
327 ///
328 /// ```
329 /// use ndarray::{aview1, arr2, Axis};
330 ///
331 /// let a = arr2(&[[1., 2.],
332 /// [3., 4.],
333 /// [5., 6.]]);
334 /// let var = a.var_axis(Axis(0), 1.);
335 /// assert_eq!(var, aview1(&[4., 4.]));
336 /// ```
337 #[cfg(feature = "std")]
338 pub fn var_axis(&self, axis: Axis, ddof: A) -> Array<A, D::Smaller>
339 where
340 A: Float + FromPrimitive,
341 D: RemoveAxis,
342 {
343 let zero = A::from_usize(0).expect("Converting 0 to `A` must not fail.");
344 let n = A::from_usize(self.len_of(axis)).expect("Converting length to `A` must not fail.");
345 assert!(
346 !(ddof < zero || ddof > n),
347 "`ddof` must not be less than zero or greater than the length of \
348 the axis",
349 );
350 let dof = n - ddof;
351 let mut mean = Array::<A, _>::zeros(self.dim.remove_axis(axis));
352 let mut sum_sq = Array::<A, _>::zeros(self.dim.remove_axis(axis));
353 for (i, subview) in self.axis_iter(axis).enumerate() {
354 let count = A::from_usize(i + 1).expect("Converting index to `A` must not fail.");
355 azip!((mean in &mut mean, sum_sq in &mut sum_sq, &x in &subview) {
356 let delta = x - *mean;
357 *mean = *mean + delta / count;
358 *sum_sq = (x - *mean).mul_add(delta, *sum_sq);
359 });
360 }
361 sum_sq.mapv_into(|s| s / dof)
362 }
363
364 /// Return standard deviation along `axis`.
365 ///
366 /// The standard deviation is computed from the variance using
367 /// the [Welford one-pass algorithm](https://www.jstor.org/stable/1266577).
368 ///
369 /// The parameter `ddof` specifies the "delta degrees of freedom". For
370 /// example, to calculate the population standard deviation, use `ddof = 0`,
371 /// or to calculate the sample standard deviation, use `ddof = 1`.
372 ///
373 /// The standard deviation is defined as:
374 ///
375 /// ```text
376 /// ⎛ 1 n ⎞
377 /// stddev = sqrt ⎜ ―――――――― ∑ (xᵢ - x̅)²⎟
378 /// ⎝ n - ddof i=1 ⎠
379 /// ```
380 ///
381 /// where
382 ///
383 /// ```text
384 /// 1 n
385 /// x̅ = ― ∑ xᵢ
386 /// n i=1
387 /// ```
388 ///
389 /// and `n` is the length of the axis.
390 ///
391 /// **Panics** if `ddof` is less than zero or greater than `n`, if `axis`
392 /// is out of bounds, or if `A::from_usize()` fails for any any of the
393 /// numbers in the range `0..=n`.
394 ///
395 /// # Example
396 ///
397 /// ```
398 /// use ndarray::{aview1, arr2, Axis};
399 ///
400 /// let a = arr2(&[[1., 2.],
401 /// [3., 4.],
402 /// [5., 6.]]);
403 /// let stddev = a.std_axis(Axis(0), 1.);
404 /// assert_eq!(stddev, aview1(&[2., 2.]));
405 /// ```
406 #[cfg(feature = "std")]
407 pub fn std_axis(&self, axis: Axis, ddof: A) -> Array<A, D::Smaller>
408 where
409 A: Float + FromPrimitive,
410 D: RemoveAxis,
411 {
412 self.var_axis(axis, ddof).mapv_into(|x| x.sqrt())
413 }
414}