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}