scirs2_core/ndarray_ext/stats/descriptive.rs
1//! Descriptive statistics for ndarray arrays
2//!
3//! This module provides descriptive statistical functions such as mean, median,
4//! standard deviation, variance, min, max, etc. for ndarray arrays.
5
6use ndarray::{Array, ArrayView, Axis, Dimension, Ix1, Ix2};
7use num_traits::{Float, FromPrimitive};
8
9/// Calculate the mean of array elements (2D arrays)
10///
11/// # Errors
12///
13/// Returns an error if the array is empty or if conversion fails.
14///
15/// # Panics
16///
17/// Panics if type conversion from usize fails.
18///
19/// # Arguments
20///
21/// * `array` - The input 2D array
22/// * `axis` - Optional axis along which to compute the mean (None for global mean)
23///
24/// # Returns
25///
26/// The mean of the array elements
27///
28/// # Examples
29///
30/// ```
31/// use ndarray::{array, Axis};
32/// use scirs2_core::ndarray_ext::stats::mean_2d;
33///
34/// let a = array![[1.0, 2.0, 3.0], [4.0, 5.0, 6.0]];
35///
36/// // Global mean
37/// let global_mean = mean_2d(&a.view(), None).unwrap();
38/// assert_eq!(global_mean[0], 3.5);
39///
40/// // Mean along axis 0 (columns)
41/// let col_means = mean_2d(&a.view(), Some(Axis(0))).unwrap();
42/// assert_eq!(col_means.len(), 3);
43/// assert_eq!(col_means[0], 2.5);
44/// assert_eq!(col_means[1], 3.5);
45/// assert_eq!(col_means[2], 4.5);
46/// ```
47#[allow(dead_code)]
48pub fn mean_2d<T>(
49 array: &ArrayView<T, Ix2>,
50 axis: Option<Axis>,
51) -> Result<Array<T, Ix1>, &'static str>
52where
53 T: Clone + Float + FromPrimitive,
54{
55 if array.is_empty() {
56 return Err("Cannot compute mean of an empty array");
57 }
58
59 if let Some(ax) = axis {
60 let (rows, cols) = (array.shape()[0], array.shape()[1]);
61
62 match ax.index() {
63 0 => {
64 // Mean along axis 0 (columns)
65 let mut result = Array::<T, Ix1>::zeros(cols);
66 let n = T::from_usize(rows).unwrap();
67
68 for j in 0..cols {
69 let mut sum = T::zero();
70 for i in 0..rows {
71 sum = sum + array[[i, j]];
72 }
73 result[j] = sum / n;
74 }
75
76 Ok(result)
77 }
78 1 => {
79 // Mean along axis 1 (rows)
80 let mut result = Array::<T, Ix1>::zeros(rows);
81 let n = T::from_usize(cols).unwrap();
82
83 for i in 0..rows {
84 let mut sum = T::zero();
85 for j in 0..cols {
86 sum = sum + array[[i, j]];
87 }
88 result[0] = sum / n;
89 }
90
91 Ok(result)
92 }
93 _ => Err("Axis index out of bounds for 2D array"),
94 }
95 } else {
96 // Global mean
97 let total_elements = array.len();
98 let mut sum = T::zero();
99
100 for &val in array {
101 sum = sum + val;
102 }
103
104 let count = T::from_usize(total_elements).ok_or("Cannot convert array length to T")?;
105 Ok(Array::from_elem(1, sum / count))
106 }
107}
108
109/// Calculate the median of array elements (2D arrays)
110///
111/// # Errors
112///
113/// Returns an error if the array is empty.
114///
115/// # Panics
116///
117/// Panics if partial comparison fails or type conversion fails.
118///
119/// # Arguments
120///
121/// * `array` - The input 2D array
122/// * `axis` - Optional axis along which to compute the median (None for global median)
123///
124/// # Returns
125///
126/// The median of the array elements
127///
128/// # Examples
129///
130/// ```
131/// use ndarray::{array, Axis};
132/// use scirs2_core::ndarray_ext::stats::median_2d;
133///
134/// let a = array![[1.0, 3.0, 5.0], [2.0, 4.0, 6.0]];
135///
136/// // Global median
137/// let global_median = median_2d(&a.view(), None).unwrap();
138/// assert_eq!(global_median[0], 3.5);
139///
140/// // Median along axis 0 (columns)
141/// let col_medians = median_2d(&a.view(), Some(Axis(0))).unwrap();
142/// assert_eq!(col_medians.len(), 3);
143/// assert_eq!(col_medians[0], 1.5);
144/// assert_eq!(col_medians[1], 3.5);
145/// assert_eq!(col_medians[2], 5.5);
146/// ```
147#[allow(dead_code)]
148pub fn median_2d<T>(
149 array: &ArrayView<T, Ix2>,
150 axis: Option<Axis>,
151) -> Result<Array<T, Ix1>, &'static str>
152where
153 T: Clone + Float + FromPrimitive,
154{
155 if array.is_empty() {
156 return Err("Cannot compute median of an empty array");
157 }
158
159 if let Some(ax) = axis {
160 let (rows, cols) = (array.shape()[0], array.shape()[1]);
161
162 match ax.index() {
163 0 => {
164 // Median along axis 0 (columns)
165 let mut result = Array::<T, Ix1>::zeros(cols);
166
167 for j in 0..cols {
168 let mut column_values = Vec::with_capacity(rows);
169 for i in 0..rows {
170 column_values.push(array[[i, j]]);
171 }
172
173 column_values
174 .sort_by(|a, b| a.partial_cmp(b).unwrap_or(std::cmp::Ordering::Equal));
175
176 let median_value = if column_values.len() % 2 == 0 {
177 let mid = column_values.len() / 2;
178 (column_values[mid - 1] + column_values[mid]) / T::from_f64(2.0).unwrap()
179 } else {
180 column_values[column_values.len() / 2]
181 };
182
183 result[j] = median_value;
184 }
185
186 Ok(result)
187 }
188 1 => {
189 // Median along axis 1 (rows)
190 let mut result = Array::<T, Ix1>::zeros(rows);
191
192 for i in 0..rows {
193 let mut row_values = Vec::with_capacity(cols);
194 for j in 0..cols {
195 row_values.push(array[[i, j]]);
196 }
197
198 row_values
199 .sort_by(|a, b| a.partial_cmp(b).unwrap_or(std::cmp::Ordering::Equal));
200
201 let median_value = if row_values.len() % 2 == 0 {
202 let mid = row_values.len() / 2;
203 (row_values[mid - 1] + row_values[mid]) / T::from_f64(2.0).unwrap()
204 } else {
205 row_values[row_values.len() / 2]
206 };
207
208 result[0] = median_value;
209 }
210
211 Ok(result)
212 }
213 _ => Err("Axis index out of bounds for 2D array"),
214 }
215 } else {
216 // Global median
217 let mut values: Vec<_> = array.iter().copied().collect();
218 values.sort_by(|a, b| a.partial_cmp(b).unwrap_or(std::cmp::Ordering::Equal));
219
220 let median_value = if values.len() % 2 == 0 {
221 let mid = values.len() / 2;
222 (values[mid - 1] + values[mid]) / T::from_f64(2.0).unwrap()
223 } else {
224 values[values.len() / 2]
225 };
226
227 Ok(Array::from_elem(1, median_value))
228 }
229}
230
231/// Calculate the standard deviation of array elements (2D arrays)
232///
233/// # Errors
234///
235/// Returns an error if the array is empty or variance calculation fails.
236///
237/// # Panics
238///
239/// Panics if variance calculation panics.
240///
241/// # Arguments
242///
243/// * `array` - The input 2D array
244/// * `axis` - Optional axis along which to compute the std dev (None for global std dev)
245/// * `ddof` - Delta degrees of freedom (default 0)
246///
247/// # Returns
248///
249/// The standard deviation of the array elements
250///
251/// # Examples
252///
253/// ```
254/// use ndarray::{array, Axis};
255/// use scirs2_core::ndarray_ext::stats::std_dev_2d;
256///
257/// let a = array![[1.0, 2.0, 3.0], [4.0, 5.0, 6.0]];
258///
259/// // Global standard deviation
260/// let global_std = std_dev_2d(&a.view(), None, 1).unwrap();
261/// assert!((global_std[0] - 1.87082869339_f64).abs() < 1e-10);
262///
263/// // Standard deviation along axis 0 (columns)
264/// let col_stds = std_dev_2d(&a.view(), Some(Axis(0)), 1).unwrap();
265/// assert_eq!(col_stds.len(), 3);
266/// ```
267#[allow(dead_code)]
268pub fn std_dev_2d<T>(
269 array: &ArrayView<T, Ix2>,
270 axis: Option<Axis>,
271 ddof: usize,
272) -> Result<Array<T, Ix1>, &'static str>
273where
274 T: Clone + Float + FromPrimitive,
275{
276 let var_result = variance_2d(array, axis, ddof)?;
277 Ok(var_result.mapv(|x| x.sqrt()))
278}
279
280/// Calculate the variance of array elements (2D arrays)
281///
282/// # Errors
283///
284/// Returns an error if the array is empty or if conversion fails.
285///
286/// # Panics
287///
288/// Panics if type conversion from usize fails.
289///
290/// # Arguments
291///
292/// * `array` - The input 2D array
293/// * `axis` - Optional axis along which to compute the variance (None for global variance)
294/// * `ddof` - Delta degrees of freedom (default 0)
295///
296/// # Returns
297///
298/// The variance of the array elements
299///
300/// # Examples
301///
302/// ```
303/// use ndarray::{array, Axis};
304/// use scirs2_core::ndarray_ext::stats::variance_2d;
305///
306/// let a = array![[1.0, 2.0, 3.0], [4.0, 5.0, 6.0]];
307///
308/// // Global variance
309/// let global_var = variance_2d(&a.view(), None, 1).unwrap();
310/// assert!((global_var[0] - 3.5_f64).abs() < 1e-10);
311///
312/// // Variance along axis 0 (columns)
313/// let col_vars = variance_2d(&a.view(), Some(Axis(0)), 1).unwrap();
314/// assert_eq!(col_vars.len(), 3);
315/// ```
316#[allow(dead_code)]
317pub fn variance_2d<T>(
318 array: &ArrayView<T, Ix2>,
319 axis: Option<Axis>,
320 ddof: usize,
321) -> Result<Array<T, Ix1>, &'static str>
322where
323 T: Clone + Float + FromPrimitive,
324{
325 if array.is_empty() {
326 return Err("Cannot compute variance of an empty array");
327 }
328
329 if let Some(ax) = axis {
330 let (rows, cols) = (array.shape()[0], array.shape()[1]);
331
332 match ax.index() {
333 0 => {
334 // Variance along axis 0 (columns)
335 let means = mean_2d(array, Some(ax))?;
336
337 if rows <= ddof {
338 return Err("Not enough data points for variance calculation with given ddof");
339 }
340
341 let mut result = Array::<T, Ix1>::zeros(cols);
342
343 for j in 0..cols {
344 let mut sum_sq_diff = T::zero();
345 for i in 0..rows {
346 let diff = array[[i, j]] - means[j];
347 sum_sq_diff = sum_sq_diff + (diff * diff);
348 }
349
350 let divisor = T::from_usize(rows - ddof).unwrap();
351 result[j] = sum_sq_diff / divisor;
352 }
353
354 Ok(result)
355 }
356 1 => {
357 // Variance along axis 1 (rows)
358 let means = mean_2d(array, Some(ax))?;
359
360 if cols <= ddof {
361 return Err("Not enough data points for variance calculation with given ddof");
362 }
363
364 let mut result = Array::<T, Ix1>::zeros(rows);
365
366 for i in 0..rows {
367 let mut sum_sq_diff = T::zero();
368 for j in 0..cols {
369 let diff = array[[i, j]] - means[i];
370 sum_sq_diff = sum_sq_diff + (diff * diff);
371 }
372
373 let divisor = T::from_usize(cols - ddof).unwrap();
374 result[0] = sum_sq_diff / divisor;
375 }
376
377 Ok(result)
378 }
379 _ => Err("Axis index out of bounds for 2D array"),
380 }
381 } else {
382 // Global variance
383 let total_elements = array.len();
384
385 if total_elements <= ddof {
386 return Err("Not enough data points for variance calculation with given ddof");
387 }
388
389 // Calculate global mean
390 let global_mean = mean_2d(array, None)?[0];
391
392 // Calculate sum of squared differences from the mean
393 let mut sum_sq_diff = T::zero();
394 for &val in array {
395 let diff = val - global_mean;
396 sum_sq_diff = sum_sq_diff + (diff * diff);
397 }
398
399 let divisor = T::from_usize(total_elements - ddof).unwrap();
400
401 Ok(Array::from_elem(1, sum_sq_diff / divisor))
402 }
403}
404
405// Need to implement the following functions:
406// min_2d, max_2d, sum_2d, product_2d, percentile_2d, mean, median, variance, std_dev, min, max, percentile
407
408// Let's continue with min_2d and max_2d
409
410/// Calculate the minimum value(s) of array elements (2D arrays)
411///
412/// # Errors
413///
414/// Returns an error if the array is empty.
415///
416/// # Arguments
417///
418/// * `array` - The input 2D array
419/// * `axis` - Optional axis along which to compute the minimum (None for global minimum)
420///
421/// # Returns
422///
423/// The minimum value(s) of the array elements
424#[allow(dead_code)]
425pub fn min_2d<T>(
426 array: &ArrayView<T, Ix2>,
427 axis: Option<Axis>,
428) -> Result<Array<T, Ix1>, &'static str>
429where
430 T: Clone + Float,
431{
432 if array.is_empty() {
433 return Err("Cannot compute minimum of an empty array");
434 }
435
436 match axis {
437 Some(ax) => {
438 let (rows, cols) = (array.shape()[0], array.shape()[1]);
439
440 match ax.index() {
441 0 => {
442 // Min along axis 0 (columns)
443 let mut result = Array::<T, Ix1>::zeros(cols);
444
445 for j in 0..cols {
446 let mut min_val = array[[0, j]];
447 for i in 1..rows {
448 if array[[i, j]] < min_val {
449 min_val = array[[i, j]];
450 }
451 }
452 result[j] = min_val;
453 }
454
455 Ok(result)
456 }
457 1 => {
458 // Min along axis 1 (rows)
459 let mut result = Array::<T, Ix1>::zeros(rows);
460
461 for i in 0..rows {
462 let mut min_val = array[[i, 0]];
463 for j in 1..cols {
464 if array[[i, j]] < min_val {
465 min_val = array[[i, j]];
466 }
467 }
468 result[i] = min_val;
469 }
470
471 Ok(result)
472 }
473 _ => Err("Axis index out of bounds for 2D array"),
474 }
475 }
476 None => {
477 // Global min
478 let mut min_val = array[[0, 0]];
479
480 for &val in array {
481 if val < min_val {
482 min_val = val;
483 }
484 }
485
486 Ok(Array::from_elem(1, min_val))
487 }
488 }
489}
490
491/// Calculate the maximum value(s) of array elements (2D arrays)
492///
493/// # Errors
494///
495/// Returns an error if the array is empty.
496///
497/// # Arguments
498///
499/// * `array` - The input 2D array
500/// * `axis` - Optional axis along which to compute the maximum (None for global maximum)
501///
502/// # Returns
503///
504/// The maximum value(s) of the array elements
505#[allow(dead_code)]
506pub fn max_2d<T>(
507 array: &ArrayView<T, Ix2>,
508 axis: Option<Axis>,
509) -> Result<Array<T, Ix1>, &'static str>
510where
511 T: Clone + Float,
512{
513 if array.is_empty() {
514 return Err("Cannot compute maximum of an empty array");
515 }
516
517 match axis {
518 Some(ax) => {
519 let (rows, cols) = (array.shape()[0], array.shape()[1]);
520
521 match ax.index() {
522 0 => {
523 // Max along axis 0 (columns)
524 let mut result = Array::<T, Ix1>::zeros(cols);
525
526 for j in 0..cols {
527 let mut max_val = array[[0, j]];
528 for i in 1..rows {
529 if array[[i, j]] > max_val {
530 max_val = array[[i, j]];
531 }
532 }
533 result[j] = max_val;
534 }
535
536 Ok(result)
537 }
538 1 => {
539 // Max along axis 1 (rows)
540 let mut result = Array::<T, Ix1>::zeros(rows);
541
542 for i in 0..rows {
543 let mut max_val = array[[i, 0]];
544 for j in 1..cols {
545 if array[[i, j]] > max_val {
546 max_val = array[[i, j]];
547 }
548 }
549 result[i] = max_val;
550 }
551
552 Ok(result)
553 }
554 _ => Err("Axis index out of bounds for 2D array"),
555 }
556 }
557 None => {
558 // Global max
559 let mut max_val = array[[0, 0]];
560
561 for &val in array {
562 if val > max_val {
563 max_val = val;
564 }
565 }
566
567 Ok(Array::from_elem(1, max_val))
568 }
569 }
570}
571
572/// Calculate the sum of array elements (2D arrays)
573///
574/// # Errors
575///
576/// Returns an error if the array is empty.
577///
578/// # Arguments
579///
580/// * `array` - The input 2D array
581/// * `axis` - Optional axis along which to compute the sum (None for global sum)
582///
583/// # Returns
584///
585/// The sum of the array elements
586#[allow(dead_code)]
587pub fn sum_2d<T>(
588 array: &ArrayView<T, Ix2>,
589 axis: Option<Axis>,
590) -> Result<Array<T, Ix1>, &'static str>
591where
592 T: Clone + Float,
593{
594 if array.is_empty() {
595 return Err("Cannot compute sum of an empty array");
596 }
597
598 match axis {
599 Some(ax) => {
600 let (rows, cols) = (array.shape()[0], array.shape()[1]);
601
602 match ax.index() {
603 0 => {
604 // Sum along axis 0 (columns)
605 let mut result = Array::<T, Ix1>::zeros(cols);
606
607 for j in 0..cols {
608 let mut sum = T::zero();
609 for i in 0..rows {
610 sum = sum + array[[i, j]];
611 }
612 result[j] = sum;
613 }
614
615 Ok(result)
616 }
617 1 => {
618 // Sum along axis 1 (rows)
619 let mut result = Array::<T, Ix1>::zeros(rows);
620
621 for i in 0..rows {
622 let mut sum = T::zero();
623 for j in 0..cols {
624 sum = sum + array[[i, j]];
625 }
626 result[0] = sum;
627 }
628
629 Ok(result)
630 }
631 _ => Err("Axis index out of bounds for 2D array"),
632 }
633 }
634 None => {
635 // Global sum
636 let mut sum = T::zero();
637
638 for &val in array {
639 sum = sum + val;
640 }
641
642 Ok(Array::from_elem(1, sum))
643 }
644 }
645}
646
647/// Calculate the product of array elements (2D arrays)
648///
649/// # Errors
650///
651/// Returns an error if the array is empty.
652///
653/// # Arguments
654///
655/// * `array` - The input 2D array
656/// * `axis` - Optional axis along which to compute the product (None for global product)
657///
658/// # Returns
659///
660/// The product of the array elements
661#[allow(dead_code)]
662pub fn product_2d<T>(
663 array: &ArrayView<T, Ix2>,
664 axis: Option<Axis>,
665) -> Result<Array<T, Ix1>, &'static str>
666where
667 T: Clone + Float,
668{
669 if array.is_empty() {
670 return Err("Cannot compute product of an empty array");
671 }
672
673 match axis {
674 Some(ax) => {
675 let (rows, cols) = (array.shape()[0], array.shape()[1]);
676
677 match ax.index() {
678 0 => {
679 // Product along axis 0 (columns)
680 let mut result = Array::<T, Ix1>::from_elem(cols, T::one());
681
682 for j in 0..cols {
683 for i in 0..rows {
684 result[j] = result[j] * array[[i, j]];
685 }
686 }
687
688 Ok(result)
689 }
690 1 => {
691 // Product along axis 1 (rows)
692 let mut result = Array::<T, Ix1>::from_elem(rows, T::one());
693
694 for i in 0..rows {
695 for j in 0..cols {
696 result[i] = result[i] * array[[i, j]];
697 }
698 }
699
700 Ok(result)
701 }
702 _ => Err("Axis index out of bounds for 2D array"),
703 }
704 }
705 None => {
706 // Global product
707 let mut product = T::one();
708
709 for &val in array {
710 product = product * val;
711 }
712
713 Ok(Array::from_elem(1, product))
714 }
715 }
716}
717
718/// Calculate the percentile of array elements (2D arrays)
719///
720/// # Errors
721///
722/// Returns an error if the array is empty or percentile is invalid.
723///
724/// # Panics
725///
726/// Panics if type conversion or partial comparison fails.
727///
728/// # Arguments
729///
730/// * `array` - The input 2D array
731/// * `q` - The percentile to compute (0 to 100)
732/// * `axis` - Optional axis along which to compute the percentile (None for global percentile)
733///
734/// # Returns
735///
736/// The percentile value(s) of the array elements
737#[allow(dead_code)]
738pub fn percentile_2d<T>(
739 array: &ArrayView<T, Ix2>,
740 q: f64,
741 axis: Option<Axis>,
742) -> Result<Array<T, Ix1>, &'static str>
743where
744 T: Clone + Float + FromPrimitive,
745{
746 if array.is_empty() {
747 return Err("Cannot compute percentile of an empty array");
748 }
749
750 if !(0.0..=100.0).contains(&q) {
751 return Err("Percentile must be between 0 and 100");
752 }
753
754 match axis {
755 Some(ax) => {
756 let (rows, cols) = (array.shape()[0], array.shape()[1]);
757
758 match ax.index() {
759 0 => {
760 // Percentile along axis 0 (columns)
761 let mut result = Array::<T, Ix1>::zeros(cols);
762
763 for j in 0..cols {
764 let mut column_values = Vec::with_capacity(rows);
765 for i in 0..rows {
766 column_values.push(array[[i, j]]);
767 }
768
769 column_values
770 .sort_by(|a, b| a.partial_cmp(b).unwrap_or(std::cmp::Ordering::Equal));
771
772 // Linear interpolation
773 let pos = (q / 100.0) * (column_values.len() as f64 - 1.0);
774 let idx_low = pos.floor() as usize;
775 let idx_high = pos.ceil() as usize;
776
777 if idx_low == idx_high {
778 result[j] = column_values[idx_low];
779 } else {
780 let weight_high = pos - (idx_low as f64);
781 let weight_low = 1.0 - weight_high;
782
783 result[j] = column_values[idx_low] * T::from_f64(weight_low).unwrap()
784 + column_values[idx_high] * T::from_f64(weight_high).unwrap();
785 }
786 }
787
788 Ok(result)
789 }
790 1 => {
791 // Percentile along axis 1 (rows)
792 let mut result = Array::<T, Ix1>::zeros(rows);
793
794 for i in 0..rows {
795 let mut row_values = Vec::with_capacity(cols);
796 for j in 0..cols {
797 row_values.push(array[[i, j]]);
798 }
799
800 row_values
801 .sort_by(|a, b| a.partial_cmp(b).unwrap_or(std::cmp::Ordering::Equal));
802
803 // Linear interpolation
804 let pos = (q / 100.0) * (row_values.len() as f64 - 1.0);
805 let idx_low = pos.floor() as usize;
806 let idx_high = pos.ceil() as usize;
807
808 if idx_low == idx_high {
809 result[0] = row_values[idx_low];
810 } else {
811 let weight_high = pos - (idx_low as f64);
812 let weight_low = 1.0 - weight_high;
813
814 result[0] = row_values[idx_low] * T::from_f64(weight_low).unwrap()
815 + row_values[idx_high] * T::from_f64(weight_high).unwrap();
816 }
817 }
818
819 Ok(result)
820 }
821 _ => Err("Axis index out of bounds for 2D array"),
822 }
823 }
824 None => {
825 // Global percentile
826 let mut values: Vec<_> = array.iter().copied().collect();
827 values.sort_by(|a, b| a.partial_cmp(b).unwrap_or(std::cmp::Ordering::Equal));
828
829 // Linear interpolation
830 let pos = (q / 100.0) * (values.len() as f64 - 1.0);
831 let idx_low = pos.floor() as usize;
832 let idx_high = pos.ceil() as usize;
833
834 let result = if idx_low == idx_high {
835 values[idx_low]
836 } else {
837 let weight_high = pos - (idx_low as f64);
838 let weight_low = 1.0 - weight_high;
839
840 values[idx_low] * T::from_f64(weight_low).unwrap()
841 + values[idx_high] * T::from_f64(weight_high).unwrap()
842 };
843
844 Ok(Array::from_elem(1, result))
845 }
846 }
847}
848
849// Generic implementations for n-dimensional arrays
850
851/// Calculate the mean of array elements
852///
853/// # Errors
854///
855/// Returns an error if the array is empty or if conversion fails.
856///
857/// # Panics
858///
859/// Panics if type conversion from usize fails.
860///
861/// # Arguments
862///
863/// * `array` - The input array
864/// * `axis` - Optional axis along which to compute the mean (None for global mean)
865///
866/// # Returns
867///
868/// The mean of the array elements
869#[allow(dead_code)]
870pub fn mean<T, D>(
871 array: &ArrayView<T, D>,
872 axis: Option<Axis>,
873) -> Result<Array<T, Ix1>, &'static str>
874where
875 T: Clone + Float + FromPrimitive,
876 D: Dimension + ndarray::RemoveAxis,
877{
878 if array.is_empty() {
879 return Err("Cannot compute mean of an empty array");
880 }
881
882 match axis {
883 Some(ax) => {
884 // Axis-specific mean implementation for arbitrary dimensions
885 let ndim = array.ndim();
886 if ax.index() >= ndim {
887 return Err("Axis index out of bounds");
888 }
889
890 // Create output shape by removing the specified axis
891 let mut outputshape = array.shape().to_vec();
892 outputshape.remove(ax.index());
893
894 // Handle case where removing axis results in scalar
895 if outputshape.is_empty() {
896 outputshape.push(1);
897 }
898
899 // Calculate mean along specified axis using ndarray's mean_axis
900 let result = array
901 .mean_axis(ax)
902 .ok_or("Failed to compute mean along axis")?;
903
904 // Convert to 1D array as expected by function signature
905 let flat_result = result
906 .to_shape((result.len(),))
907 .map_err(|_| "Failed to reshape result to 1D")?;
908
909 Ok(flat_result.into_owned())
910 }
911 None => {
912 // Global mean
913 let total_elements = array.len();
914 let mut sum = T::zero();
915
916 for &val in array {
917 sum = sum + val;
918 }
919
920 let count = T::from_usize(total_elements).ok_or("Cannot convert array length to T")?;
921 Ok(Array::from_elem(1, sum / count))
922 }
923 }
924}
925
926/// Calculate the median of array elements
927///
928/// # Arguments
929///
930/// * `array` - The input array
931/// * `axis` - Optional axis along which to compute the median (None for global median)
932///
933/// # Returns
934///
935/// The median of the array elements
936///
937/// # Errors
938///
939/// Returns an error if the array is empty.
940///
941/// # Panics
942///
943/// Panics if type conversion or partial comparison fails.
944#[allow(dead_code)]
945pub fn median<T, D>(
946 array: &ArrayView<T, D>,
947 axis: Option<Axis>,
948) -> Result<Array<T, Ix1>, &'static str>
949where
950 T: Clone + Float + FromPrimitive,
951 D: Dimension + ndarray::RemoveAxis,
952{
953 if array.is_empty() {
954 return Err("Cannot compute median of an empty array");
955 }
956
957 match axis {
958 Some(_) => {
959 // For higher dimensional arrays, we need to implement axis-specific logic
960 // This is a placeholder for now
961 Err("Axis-specific median for arbitrary dimension arrays not yet implemented")
962 }
963 None => {
964 // Global median
965 let mut values: Vec<_> = array.iter().copied().collect();
966 values.sort_by(|a, b| a.partial_cmp(b).unwrap_or(std::cmp::Ordering::Equal));
967
968 let median_value = if values.len() % 2 == 0 {
969 let mid = values.len() / 2;
970 (values[mid - 1] + values[mid]) / T::from_f64(2.0).unwrap()
971 } else {
972 values[values.len() / 2]
973 };
974
975 Ok(Array::from_elem(1, median_value))
976 }
977 }
978}
979
980/// Calculate the variance of array elements
981///
982/// # Arguments
983///
984/// * `array` - The input array
985/// * `axis` - Optional axis along which to compute the variance (None for global variance)
986/// * `ddof` - Delta degrees of freedom (default 0)
987///
988/// # Returns
989///
990/// The variance of the array elements
991///
992/// # Errors
993///
994/// Returns an error if the array is empty or if conversion fails.
995///
996/// # Panics
997///
998/// Panics if type conversion from usize fails.
999#[allow(dead_code)]
1000pub fn variance<T, D>(
1001 array: &ArrayView<T, D>,
1002 axis: Option<Axis>,
1003 ddof: usize,
1004) -> Result<Array<T, Ix1>, &'static str>
1005where
1006 T: Clone + Float + FromPrimitive,
1007 D: Dimension + ndarray::RemoveAxis,
1008{
1009 if array.is_empty() {
1010 return Err("Cannot compute variance of an empty array");
1011 }
1012
1013 match axis {
1014 Some(ax) => {
1015 // Axis-specific variance implementation for arbitrary dimensions
1016 let ndim = array.ndim();
1017 if ax.index() >= ndim {
1018 return Err("Axis index out of bounds");
1019 }
1020
1021 // Create output shape by removing the specified axis
1022 let mut outputshape = array.shape().to_vec();
1023 outputshape.remove(ax.index());
1024
1025 // Handle case where removing axis results in scalar
1026 if outputshape.is_empty() {
1027 outputshape.push(1);
1028 }
1029
1030 // Calculate variance along specified axis using ndarray's var_axis
1031 let result = array.var_axis(ax, T::from_usize(ddof).unwrap());
1032
1033 // Convert to 1D array as expected by function signature
1034 let flat_result = result
1035 .to_shape((result.len(),))
1036 .map_err(|_| "Failed to reshape variance result to 1D")?;
1037
1038 Ok(flat_result.into_owned())
1039 }
1040 None => {
1041 // Global variance
1042 let total_elements = array.len();
1043
1044 if total_elements <= ddof {
1045 return Err("Not enough data points for variance calculation with given ddof");
1046 }
1047
1048 // Calculate global mean
1049 let global_mean = mean(array, None)?[0];
1050
1051 // Calculate sum of squared differences from the mean
1052 let mut sum_sq_diff = T::zero();
1053 for &val in array {
1054 let diff = val - global_mean;
1055 sum_sq_diff = sum_sq_diff + (diff * diff);
1056 }
1057
1058 let divisor = T::from_usize(total_elements - ddof).unwrap();
1059
1060 Ok(Array::from_elem(1, sum_sq_diff / divisor))
1061 }
1062 }
1063}
1064
1065/// Calculate the standard deviation of array elements
1066///
1067/// # Arguments
1068///
1069/// * `array` - The input array
1070/// * `axis` - Optional axis along which to compute the std dev (None for global std dev)
1071/// * `ddof` - Delta degrees of freedom (default 0)
1072///
1073/// # Returns
1074///
1075/// The standard deviation of the array elements
1076///
1077/// # Errors
1078///
1079/// Returns an error if the array is empty or variance calculation fails.
1080///
1081/// # Panics
1082///
1083/// Panics if variance calculation panics.
1084#[allow(dead_code)]
1085pub fn std_dev<T, D>(
1086 array: &ArrayView<T, D>,
1087 axis: Option<Axis>,
1088 ddof: usize,
1089) -> Result<Array<T, Ix1>, &'static str>
1090where
1091 T: Clone + Float + FromPrimitive,
1092 D: Dimension + ndarray::RemoveAxis,
1093{
1094 let var_result = variance(array, axis, ddof)?;
1095 Ok(var_result.mapv(|x| x.sqrt()))
1096}
1097
1098/// Calculate the minimum value(s) of array elements
1099///
1100/// # Arguments
1101///
1102/// * `array` - The input array
1103/// * `axis` - Optional axis along which to compute the minimum (None for global minimum)
1104///
1105/// # Returns
1106///
1107/// The minimum value(s) of the array elements
1108///
1109/// # Errors
1110///
1111/// Returns an error if the array is empty.
1112#[allow(dead_code)]
1113pub fn min<T, D>(array: &ArrayView<T, D>, axis: Option<Axis>) -> Result<Array<T, Ix1>, &'static str>
1114where
1115 T: Clone + Float,
1116 D: Dimension + ndarray::RemoveAxis,
1117{
1118 if array.is_empty() {
1119 return Err("Cannot compute minimum of an empty array");
1120 }
1121
1122 match axis {
1123 Some(ax) => {
1124 // Axis-specific minimum implementation for arbitrary dimensions
1125 let ndim = array.ndim();
1126 if ax.index() >= ndim {
1127 return Err("Axis index out of bounds");
1128 }
1129
1130 // Create output shape by removing the specified axis
1131 let mut outputshape = array.shape().to_vec();
1132 outputshape.remove(ax.index());
1133
1134 // Handle case where removing axis results in scalar
1135 if outputshape.is_empty() {
1136 outputshape.push(1);
1137 }
1138
1139 // Use ndarray's fold_axis to compute minimum along specified axis
1140 let result = array.fold_axis(ax, T::infinity(), |&a, &b| if a < b { a } else { b });
1141
1142 // Convert to 1D array as expected by function signature
1143 let flat_result = result
1144 .to_shape((result.len(),))
1145 .map_err(|_| "Failed to reshape minimum result to 1D")?;
1146
1147 Ok(flat_result.into_owned())
1148 }
1149 None => {
1150 // Global minimum
1151 let mut min_val = *array.iter().next().unwrap();
1152
1153 for &val in array {
1154 if val < min_val {
1155 min_val = val;
1156 }
1157 }
1158
1159 Ok(Array::from_elem(1, min_val))
1160 }
1161 }
1162}
1163
1164/// Calculate the maximum value(s) of array elements
1165///
1166/// # Arguments
1167///
1168/// * `array` - The input array
1169/// * `axis` - Optional axis along which to compute the maximum (None for global maximum)
1170///
1171/// # Returns
1172///
1173/// The maximum value(s) of the array elements
1174///
1175/// # Errors
1176///
1177/// Returns an error if the array is empty.
1178#[allow(dead_code)]
1179pub fn max<T, D>(array: &ArrayView<T, D>, axis: Option<Axis>) -> Result<Array<T, Ix1>, &'static str>
1180where
1181 T: Clone + Float,
1182 D: Dimension + ndarray::RemoveAxis,
1183{
1184 if array.is_empty() {
1185 return Err("Cannot compute maximum of an empty array");
1186 }
1187
1188 match axis {
1189 Some(ax) => {
1190 // Axis-specific maximum implementation for arbitrary dimensions
1191 let ndim = array.ndim();
1192 if ax.index() >= ndim {
1193 return Err("Axis index out of bounds");
1194 }
1195
1196 // Create output shape by removing the specified axis
1197 let mut outputshape = array.shape().to_vec();
1198 outputshape.remove(ax.index());
1199
1200 // Handle case where removing axis results in scalar
1201 if outputshape.is_empty() {
1202 outputshape.push(1);
1203 }
1204
1205 // Use ndarray's fold_axis to compute maximum along specified axis
1206 let result = array.fold_axis(ax, T::neg_infinity(), |&a, &b| if a > b { a } else { b });
1207
1208 // Convert to 1D array as expected by function signature
1209 let flat_result = result
1210 .to_shape((result.len(),))
1211 .map_err(|_| "Failed to reshape maximum result to 1D")?;
1212
1213 Ok(flat_result.into_owned())
1214 }
1215 None => {
1216 // Global maximum
1217 let mut max_val = *array.iter().next().unwrap();
1218
1219 for &val in array {
1220 if val > max_val {
1221 max_val = val;
1222 }
1223 }
1224
1225 Ok(Array::from_elem(1, max_val))
1226 }
1227 }
1228}
1229
1230/// Calculate the percentile of array elements
1231///
1232/// # Arguments
1233///
1234/// * `array` - The input array
1235/// * `q` - The percentile to compute (0 to 100)
1236/// * `axis` - Optional axis along which to compute the percentile (None for global percentile)
1237///
1238/// # Returns
1239///
1240/// The percentile value(s) of the array elements
1241///
1242/// # Errors
1243///
1244/// Returns an error if the array is empty or percentile is invalid.
1245///
1246/// # Panics
1247///
1248/// Panics if type conversion or partial comparison fails.
1249#[allow(dead_code)]
1250pub fn percentile<T, D>(
1251 array: &ArrayView<T, D>,
1252 q: f64,
1253 axis: Option<Axis>,
1254) -> Result<Array<T, Ix1>, &'static str>
1255where
1256 T: Clone + Float + FromPrimitive,
1257 D: Dimension + ndarray::RemoveAxis,
1258{
1259 if array.is_empty() {
1260 return Err("Cannot compute percentile of an empty array");
1261 }
1262
1263 if !(0.0..=100.0).contains(&q) {
1264 return Err("Percentile must be between 0 and 100");
1265 }
1266
1267 match axis {
1268 Some(_) => {
1269 // For higher dimensional arrays, we need to implement axis-specific logic
1270 // This is a placeholder for now
1271 Err("Axis-specific percentile for arbitrary dimension arrays not yet implemented")
1272 }
1273 None => {
1274 // Global percentile
1275 let mut values: Vec<_> = array.iter().copied().collect();
1276 values.sort_by(|a, b| a.partial_cmp(b).unwrap_or(std::cmp::Ordering::Equal));
1277
1278 // Linear interpolation
1279 let pos = (q / 100.0) * (values.len() as f64 - 1.0);
1280 let idx_low = pos.floor() as usize;
1281 let idx_high = pos.ceil() as usize;
1282
1283 let result = if idx_low == idx_high {
1284 values[idx_low]
1285 } else {
1286 let weight_high = pos - (idx_low as f64);
1287 let weight_low = 1.0 - weight_high;
1288
1289 values[idx_low] * T::from_f64(weight_low).unwrap()
1290 + values[idx_high] * T::from_f64(weight_high).unwrap()
1291 };
1292
1293 Ok(Array::from_elem(1, result))
1294 }
1295 }
1296}