hstats 0.3.0

use core::fmt::{Error, Formatter};
use core::num::NonZeroU64;
use core::{
    fmt::{Debug, Display},
    ops::AddAssign,
};

use alloc::format;
use alloc::string::{String, ToString};
use alloc::vec;
use alloc::vec::Vec;

use num_traits::{Float, FromPrimitive};
use rolling_stats::Stats;

const DEFAULT_BAR_CHAR: &str = "░";
const DEFAULT_PRECISION: usize = 2;
const DEFAULT_DISPLAY_PERCENTILES: &[u64] = &[25, 50, 75, 90, 99];

/// `Hstats` is a struct for creating and managing histograms of data.
///
/// The generic `T` refers to the type of the data this histogram manages.
/// `T` must be a float type that implements Float, AddAssign, FromPrimitive,
///  Debug , Display
///
/// The struct includes fields for managing the histogram bins, underflow,
/// overflow, and other statistics.
#[derive(Debug, Clone)]
pub struct Hstats<T>
where
    T: Float + AddAssign + FromPrimitive + Debug + Display,
{
    start: T,
    end: T,
    bin_count: usize,
    bin_width: T,
    bins: Vec<u64>,
    underflow: u64,
    overflow: u64,
    stats: Stats<T>,
    precision: usize,
    bar_char: String,
    display_percentiles: Vec<u64>,
}

impl<T> Hstats<T>
where
    T: Float + AddAssign + FromPrimitive + Debug + Display,
{
    /// Constructs a new `Hstats` instance with specified start and end points and bin count.
    ///
    /// # Arguments
    ///
    /// * `start`: Lower bound of the range for the histogram bins.
    /// * `end`: Upper bound of the range for the histogram bins.
    /// * `bin_count`: Number of bins in the histogram.
    ///
    /// # Panics
    ///
    /// Panics if `start` >= `end` or if `bin_count` <= 0.
    pub fn new(start: T, end: T, bin_count: usize) -> Self {
        assert!(start < end, "start ({start}) must be less than end ({end})");
        assert!(
            bin_count > 0,
            "bin_count ({bin_count}) must be greater than 0"
        );

        let bin_width = (end - start) / T::from(bin_count).unwrap();

        Self {
            start,
            end,
            bin_count,
            bin_width,
            bins: vec![0; bin_count],
            underflow: 0,
            overflow: 0,
            stats: Stats::new(),
            precision: DEFAULT_PRECISION,
            bar_char: DEFAULT_BAR_CHAR.to_string(),
            display_percentiles: DEFAULT_DISPLAY_PERCENTILES.to_vec(),
        }
    }

    /// Adds a value to the histogram and updates the statistics.
    ///
    /// # Arguments
    ///
    /// * `value`: Value to be added to the histogram.
    pub fn add(&mut self, value: T) {
        if value.is_nan() {
            return;
        }

        self.stats.update(value);

        if value < self.start {
            self.underflow += 1;
        } else if value >= self.end {
            self.overflow += 1;
        } else {
            let index = ((value - self.start) / self.bin_width)
                .floor()
                .to_usize()
                .unwrap();
            self.bins[index] += 1;
        }
    }

    /// Merges this histogram with another histogram.
    ///
    /// # Arguments
    ///
    /// * `other`: Another `Hstats` instance to merge with this histogram.
    ///
    /// # Returns
    ///
    /// A new `Hstats` instance resulting from the merge.
    ///
    /// # Panics
    ///
    /// Panics if the `start`, `end`, and `bin_count` of the two histograms aren't equal.
    pub fn merge(&self, other: &Self) -> Self {
        assert_eq!(self.start, other.start, "Starts must be equal");
        assert_eq!(self.end, other.end, "Ends must be equal");
        assert_eq!(self.bin_count, other.bin_count, "Bin counts must be equal");

        let mut merged = Hstats::new(self.start, self.end, self.bin_count);
        merged.precision = self.precision;
        merged.bar_char = self.bar_char.clone();
        merged.display_percentiles = self.display_percentiles.clone();

        merged.underflow = self.underflow + other.underflow;
        merged.overflow = self.overflow + other.overflow;

        for (i, (left, right)) in (self.bins.iter().zip(other.bins.iter())).enumerate() {
            merged.bins[i] = *left + *right;
        }

        merged.stats = self.stats.merge(&other.stats);

        merged
    }

    /// Returns the number of bins in the histogram.
    /// Same as the `bin_count` argument passed to `new()`.
    pub fn bin_count(&self) -> usize {
        self.bin_count
    }

    /// Returns the width of each bin in the histogram.
    /// Same as `(end - start) / bin_count`.
    pub fn bin_width(&self) -> T {
        self.bin_width
    }

    /// Returns the start of the range for the histogram bins.
    ///
    /// Values < `start` are counted in the underflow. Values >= `start` are counted in the bins.
    pub fn start(&self) -> T {
        self.start
    }

    /// Returns the end of the range for the histogram bins.
    ///
    /// Values < `end` are counted in the bins. Values >= `end` are counted in the overflow.
    pub fn end(&self) -> T {
        self.end
    }

    /// Returns the ranges and counts for the histogram bins.
    ///
    /// # Returns
    ///
    /// A vector of tuples. Each tuple has lower bound, upper bound, and count for each bin.
    pub fn bins(&self) -> Vec<(T, T, u64)> {
        let mut bins = Vec::with_capacity(self.bin_count + 2);

        // From negative infinity to the start of the first bin
        bins.push((T::neg_infinity(), self.start, self.underflow));

        // From the start of the first bin to the end of the last bin
        let mut lower = self.start;
        let mut upper = self.start + self.bin_width;

        for count in &self.bins {
            bins.push((lower, upper, *count));
            lower = upper;
            upper += self.bin_width;
        }

        // From the end of the last bin to positive infinity
        bins.push((self.end, T::infinity(), self.overflow));

        bins
    }

    /// Returns the same bins as `bins`, but with cumulative counts for each bin.
    ///
    /// # Returns
    ///
    /// A vector of tuples. Each tuple has lower bound, upper bound and cumulative count for each bin.
    pub fn bins_cumulative(&self) -> Vec<(T, T, u64)> {
        let bins = self.bins();
        let mut cumul_bins = Vec::<(T, T, u64)>::with_capacity(bins.len());

        let mut aggregate = 0;
        for (lower, upper, count) in bins {
            aggregate += count;
            cumul_bins.push((lower, upper, aggregate));
        }

        cumul_bins
    }

    /// Shortcut to `bins_at_quantiles`, scaled to a 10,000 for myriatiles.
    ///
    /// Example:
    /// ```
    /// use hstats::Hstats;
    ///
    /// let mut hstats = Hstats::new(0.0, 10.0, 10);
    /// hstats.add(5.0);
    ///
    /// assert_eq!(
    ///     hstats.bins_at_myriatiles(&[0,5_000, 10_000]),
    ///     hstats.bins_at_centiles(&[0, 50, 100])
    /// );
    ///
    /// // Bin at 99.99% ( 9999 / 10_000 ):
    /// assert_eq!(
    ///     hstats.bins_at_myriatiles(&[9_999])[0],
    ///     hstats.bins_cumulative()[6]
    /// );
    /// ```
    pub fn bins_at_myriatiles(&self, myriatiles: &[u64]) -> Vec<(T, T, u64)> {
        self.bins_at_quantiles(myriatiles, NonZeroU64::new(10_000).unwrap())
    }

    /// Shortcut to `bins_at_quantiles`, scaled to a 100 for percentiles.
    ///
    /// Example:
    /// ```
    /// use hstats::Hstats;
    ///
    /// let mut hstats = Hstats::new(0.0, 10.0, 10);
    /// hstats.add(5.0);
    ///
    /// assert_eq!(hstats.bins_at_centiles(&[0])[0] , hstats.bins_cumulative()[0]);
    /// assert_eq!(hstats.bins_at_centiles(&[50])[0] , hstats.bins_cumulative()[6]);
    /// assert_eq!(hstats.bins_at_centiles(&[100])[0] , hstats.bins_cumulative()[6]);
    /// ```
    pub fn bins_at_centiles(&self, percentiles: &[u64]) -> Vec<(T, T, u64)> {
        self.bins_at_quantiles(percentiles, NonZeroU64::new(100).unwrap())
    }

    /// Shortcut to `bins_at_quantiles`, scaled to a 4 for quartiles.
    ///
    /// Example:
    /// ```
    /// use hstats::Hstats;
    ///
    /// let mut hstats = Hstats::new(0.0, 10.0, 10);
    /// hstats.add(5.0);
    ///
    /// assert_eq!(hstats.bins_at_centiles(&[0])[0] , hstats.bins_at_quartiles(&[0])[0]);
    /// assert_eq!(hstats.bins_at_centiles(&[50])[0] , hstats.bins_at_quartiles(&[2])[0]);
    /// assert_eq!(hstats.bins_at_centiles(&[100])[0] , hstats.bins_at_quartiles(&[4])[0]);
    /// ```
    pub fn bins_at_quartiles(&self, quartiles: &[u64]) -> Vec<(T, T, u64)> {
        self.bins_at_quantiles(quartiles, NonZeroU64::new(4).unwrap())
    }

    /// Returns the bins (potentially duplicated) covering the given quantile and scale.
    ///
    /// Quantile values are clamped to be in [0,scale] (boundaries included)
    ///
    /// # Returns
    ///
    /// A vector of tuples. Each tuple has lower bound, upper bound and cumulative count for each bin and contains the given quantile/scale.
    ///
    /// Example:
    /// ```
    /// use hstats::Hstats;
    /// use core::num::NonZeroU64;
    ///
    /// let mut hstats = Hstats::new(0.0, 10.0, 10);
    /// hstats.add(5.0);
    ///
    /// // Get the 99.999% quantile bin:
    /// assert_eq!(
    ///     hstats.bins_at_quantiles(&[99_999], NonZeroU64::new(100_000).unwrap())[0],
    ///     hstats.bins_cumulative()[6]
    /// );
    ///
    /// ```
    pub fn bins_at_quantiles(&self, quantiles: &[u64], scale: NonZeroU64) -> Vec<(T, T, u64)> {
        // Compute the count boundaries and save the original indices/order for quantiles.
        let mut quantiles_counts = {
            let count = self.count();
            quantiles
                .iter()
                .map(|&pc| {
                    pc.clamp(0, scale.get())
                        .saturating_mul(count as u64)
                        .div_ceil(scale.get())
                })
                .enumerate()
                .collect::<Vec<_>>()
        };

        // Sort by increasing value for the algorithm below.
        quantiles_counts.sort_unstable_by_key(|&(_, v)| v);

        let mut res = Vec::with_capacity(quantiles_counts.len());

        let mut cumul_bins = self.bins_cumulative().into_iter().peekable();

        for quantile_count in quantiles_counts {
            // Is the stable next value still good for the percentile?
            while let Some(&bin) = cumul_bins.peek() {
                // Advance if the bin is not good.
                if bin.2 < quantile_count.1 {
                    cumul_bins.next();
                } else {
                    // Bin is good. save and stay put for the next count.
                    res.push((quantile_count.0, bin));
                    break;
                }
            }
        }

        // Reorder the bins by original quantiles order
        res.sort_unstable_by_key(|&(idx, _)| idx);

        debug_assert_eq!(res.len(), quantiles.len());

        res.into_iter().map(|(_, bin)| bin).collect()
    }

    /// Maximum value seen so far.
    pub fn max(&self) -> T {
        self.stats.max
    }

    /// Minimum value seen so far.
    pub fn min(&self) -> T {
        self.stats.min
    }

    /// Mean value calculated so far.
    pub fn mean(&self) -> T {
        self.stats.mean
    }

    /// Standard deviation calculated so far.
    pub fn std_dev(&self) -> T {
        self.stats.std_dev
    }

    /// Number of values seen so far.
    pub fn count(&self) -> usize {
        self.stats.count
    }

    /// Modifies the decimal precision used in display output.
    pub fn with_precision(mut self, precision: usize) -> Self {
        self.precision = precision;
        self
    }

    /// Modifies the bar character used in display output.
    pub fn with_bar_char(mut self, bar_char: &str) -> Self {
        self.bar_char = bar_char.to_string();
        self
    }

    /// Modifies the percentiles shown in display output.
    /// Pass an empty slice to hide percentiles.
    pub fn with_display_percentiles(mut self, percentiles: &[u64]) -> Self {
        self.display_percentiles = percentiles.to_vec();
        self
    }
}

impl<T> Default for Hstats<T>
where
    T: Float + AddAssign + FromPrimitive + Debug + Display,
{
    fn default() -> Self {
        Self::new(T::zero(), T::one(), 10)
    }
}

/// Display the histogram as a text-based histogram.
impl<T> Display for Hstats<T>
where
    T: Float + AddAssign + FromPrimitive + Debug + Display,
{
    fn fmt(&self, f: &mut Formatter<'_>) -> Result<(), Error> {
        const MAX_BAR_SIZE: usize = 60;

        let max_count = *self.bins.iter().max().unwrap_or(&0);
        let total_count = self.count();
        let bins = self.bins();

        let col1 = bins
            .iter()
            .map(|(start, _, _)| format!("{:.*}", self.precision, start).len())
            .max()
            .unwrap_or(5);

        let col2 = bins
            .iter()
            .map(|(_, end, _)| format!("{:.*}", self.precision, end).len())
            .max()
            .unwrap_or(5);

        let precision = self.precision;

        writeln!(f, "{:^col1$} | {:^col2$}", "Start", "End")?;
        writeln!(f, "{:-^col1$}-|-{:-^col2$}-", "", "")?;
        for (range_start, range_end, count) in &bins {
            let bar_length = if max_count > 0 {
                ((*count as f64 / max_count as f64) * MAX_BAR_SIZE as f64) as usize
            } else {
                0
            };

            let percent = if total_count > 0 {
                *count as f64 / total_count as f64 * 100.0
            } else {
                0.0
            };

            let bar = self.bar_char.repeat(bar_length);

            writeln!(
                f,
                "{range_start:>col1$.precision$} | {range_end:>col2$.precision$} | {bar} {count} ({percent:.2}%)",
            )?;
        }
        writeln!(f)?;
        write!(f, "Total Count: {}", total_count)?;
        write!(f, " Min: {:.*}", self.precision, self.min())?;
        write!(f, " Max: {:.*}", self.precision, self.max())?;
        write!(f, " Mean: {:.*}", self.precision, self.mean())?;
        writeln!(f, " Std Dev: {:.*}", self.precision, self.std_dev())?;
        writeln!(f)?;

        if total_count > 0 && !self.display_percentiles.is_empty() {
            let pcts = self.bins_at_centiles(&self.display_percentiles);
            writeln!(f, "Percentiles:")?;
            let two = T::from(2.0).unwrap();
            for (p, (lower, upper, _)) in self.display_percentiles.iter().zip(pcts.iter()) {
                let midpoint = (*lower + *upper) / two;
                writeln!(f, "  p{p}: ~{:.*}", self.precision, midpoint)?;
            }
        }

        Ok(())
    }
}

#[cfg(test)]
mod tests {
    use float_cmp::ApproxEq;
    use rand::SeedableRng;
    use rand_distr::{Distribution, Normal};
    use rayon::{prelude::ParallelIterator, slice::ParallelSlice};

    use super::*;

    // Tests for Hstats::new
    #[test]
    fn test_new() {
        let hstats = Hstats::new(0.0, 10.0, 10);

        assert_eq!(hstats.start, 0.0);
        assert_eq!(hstats.end, 10.0);
        assert_eq!(hstats.bin_count, 10);
        assert_eq!(hstats.bin_width, 1.0);
        assert_eq!(hstats.bins.len(), 10);
    }

    #[test]
    #[should_panic(expected = "start (0) must be less than end (0)")]
    fn test_new_start_equal_end() {
        let _ = Hstats::new(0.0, 0.0, 10);
    }

    #[test]
    #[should_panic(expected = "bin_count (0) must be greater than 0")]
    fn test_new_bin_count_zero() {
        let _ = Hstats::new(0.0, 10.0, 0);
    }

    // Tests for Hstats::add
    #[test]
    fn test_add() {
        let mut hstats = Hstats::new(0.0, 10.0, 10);
        hstats.add(5.0);

        assert_eq!(hstats.bins[5], 1);
        assert_eq!(hstats.count(), 1);
    }

    #[test]
    fn test_add_underflow() {
        let mut hstats = Hstats::new(0.0, 10.0, 10);
        hstats.add(-1.0);
        hstats.add(0.0);

        assert_eq!(hstats.underflow, 1);
        assert_eq!(hstats.count(), 2);
    }

    #[test]
    fn test_add_overflow() {
        let mut hstats = Hstats::new(0.0, 10.0, 10);
        hstats.add(11.0);
        hstats.add(10.0);

        assert_eq!(hstats.overflow, 2);
        assert_eq!(hstats.count(), 2);
    }

    // Test for Hstats::merge
    #[test]
    fn test_merge() {
        let mut hstats1 = Hstats::new(0.0, 10.0, 10);
        hstats1.add(5.0);
        let mut hstats2 = Hstats::new(0.0, 10.0, 10);
        hstats2.add(6.0);

        let merged = hstats1.merge(&hstats2);
        assert_eq!(merged.bins[5], 1);
        assert_eq!(merged.bins[6], 1);
        assert_eq!(merged.count(), 2);
    }

    #[test]
    #[should_panic(expected = "Starts must be equal")]
    fn test_merge_different_start() {
        let hstats1 = Hstats::new(0.0, 10.0, 10);
        let hstats2 = Hstats::new(1.0, 10.0, 10);

        let _ = hstats1.merge(&hstats2);
    }

    #[test]
    #[should_panic(expected = "Ends must be equal")]
    fn test_merge_different_end() {
        let hstats1 = Hstats::new(0.0, 10.0, 10);
        let hstats2 = Hstats::new(0.0, 11.0, 10);

        let _ = hstats1.merge(&hstats2);
    }

    #[test]
    #[should_panic(expected = "Bin counts must be equal")]
    fn test_merge_different_bin_count() {
        let hstats1 = Hstats::new(0.0, 10.0, 10);
        let hstats2 = Hstats::new(0.0, 10.0, 11);

        let _ = hstats1.merge(&hstats2);
    }

    #[test]
    fn test_default() {
        let hstats: Hstats<f64> = Hstats::default();
        assert_eq!(hstats.start(), 0.0);
        assert_eq!(hstats.end(), 1.0);
        assert_eq!(hstats.bin_count(), 10);
        assert_eq!(hstats.count(), 0);
    }

    #[test]
    fn test_add_nan() {
        let mut hstats = Hstats::new(0.0, 10.0, 10);
        hstats.add(5.0);
        hstats.add(f64::NAN);
        hstats.add(3.0);

        assert_eq!(hstats.count(), 2);
        assert_eq!(hstats.bins[5], 1);
        assert_eq!(hstats.bins[3], 1);

        // 33% percentile gives the first bin with non zero count.
        let first_non_zero = hstats
            .bins_cumulative()
            .into_iter()
            .find(|&bin| bin.2 >= 1)
            .unwrap();
        assert_eq!(hstats.bins_at_centiles(&[33])[0], first_non_zero);
    }

    #[test]
    fn test_quantiles_empty_histogram() {
        let hstats = Hstats::new(0.0, 10.0, 10);
        let cumulative = hstats.bins_cumulative();

        assert_eq!(hstats.bins_at_centiles(&[0])[0], cumulative[0]);
        assert_eq!(
            hstats.bins_at_centiles(&[50, 100]),
            vec![cumulative[0], cumulative[0]]
        );
    }

    #[test]
    fn test_quantiles_single_element() {
        let mut hstats = Hstats::new(0.0, 10.0, 10);
        hstats.add(5.0);
        let cumulative = hstats.bins_cumulative();

        // Bin 6 in cumulative (underflow + bins 0..=5) contains the element
        let bin_with_element = cumulative[6];
        assert_eq!(bin_with_element.2, 1);

        // All non-zero percentiles should resolve to the bin containing the element
        assert_eq!(hstats.bins_at_centiles(&[1])[0], bin_with_element);
        assert_eq!(hstats.bins_at_centiles(&[50])[0], bin_with_element);
        assert_eq!(hstats.bins_at_centiles(&[100])[0], bin_with_element);

        // 0th percentile returns the first bin (underflow, count 0)
        assert_eq!(hstats.bins_at_centiles(&[0])[0], cumulative[0]);
    }

    #[test]
    fn test_merge_preserves_settings() {
        let mut h1 = Hstats::new(0.0, 10.0, 10)
            .with_precision(4)
            .with_bar_char("#")
            .with_display_percentiles(&[50, 99]);
        h1.add(5.0);
        let mut h2 = Hstats::new(0.0, 10.0, 10);
        h2.add(6.0);

        let merged = h1.merge(&h2);
        assert_eq!(merged.precision, 4);
        assert_eq!(merged.bar_char, "#");
        assert_eq!(merged.display_percentiles, vec![50, 99]);
    }

    #[test]
    fn test_display_empty_histogram() {
        let hstats = Hstats::new(0.0, 10.0, 5);
        let output = format!("{}", hstats);
        assert!(output.contains("Total Count: 0"));
    }

    #[test]
    fn stats_for_large_random_data() {
        type T = f64;

        // Define some constants
        // Note that the random data mean is not centered on the middle of the START..END interval.
        const MEAN: T = 2.0;
        const STD_DEV: T = 3.0;
        const SEED: u64 = 42;
        const NUM_SAMPLES: usize = 10_000;
        const NUM_BINS: usize = 100;
        const START: T = -10.0;
        const END: T = 10.0;

        let mut hstats = Hstats::new(START, END, NUM_BINS);

        let mut rng = rand::rngs::StdRng::seed_from_u64(SEED);

        let normal = Normal::<T>::new(MEAN, STD_DEV).unwrap();

        // Generate some random data
        let random_data: Vec<T> = (0..NUM_SAMPLES).map(|_x| normal.sample(&mut rng)).collect();

        // Update the stats
        random_data.iter().for_each(|v| hstats.add(*v));

        // Check the standard deviation against the stats' standard deviation
        assert!(hstats.std_dev().approx_eq(STD_DEV, (1.0e-2, 2)));

        // Check the mean against the stats' mean
        assert!(hstats.mean().approx_eq(MEAN, (1.0e-1, 2)));

        // Check the counts
        assert_eq!(hstats.count(), random_data.len());

        let count_from_bins: u64 =
            hstats.bins.iter().copied().sum::<u64>() + hstats.underflow + hstats.overflow;

        assert_eq!(count_from_bins as usize, random_data.len());

        let cumulative = hstats.bins_cumulative();

        assert_eq!(hstats.bins().len(), cumulative.len());

        // Check cumulative bins counts increase monotonically
        let mut prev_cumul = u64::MIN;
        for cbin in cumulative.iter() {
            assert!(cbin.2 >= prev_cumul);
            prev_cumul = cbin.2;
        }

        // Check last cumulative bin count is the actual count
        assert_eq!(cumulative.last().unwrap().2 as usize, hstats.count());

        // bins_at_percentiles
        // 0 should give the first bin.
        assert_eq!(hstats.bins_at_centiles(&[0])[0], cumulative[0]);
        // 100 should give the last bin
        assert_eq!(
            &hstats.bins_at_centiles(&[100])[0],
            cumulative.last().unwrap()
        );

        // Several times the same percentage should give several time
        // the same bin, and the order of the requested percentiles is preserved.
        assert_eq!(
            &hstats.bins_at_centiles(&[100, 0, 100])[0],
            cumulative.last().unwrap()
        );
        assert_eq!(hstats.bins_at_centiles(&[100, 0, 100])[1], cumulative[0]);
        assert_eq!(
            &hstats.bins_at_centiles(&[100, 0, 100])[2],
            cumulative.last().unwrap()
        );

        // The bin 61 is at the 50th percentile.
        // It is not the bin 50, because the mean of the random data
        // is 2.0 , and not zero (bins range from -10.0 to +10.0)
        assert_eq!(hstats.bins_at_centiles(&[50, 50])[0], cumulative[61]);
        // The second value refers to the same bin
        assert_eq!(hstats.bins_at_centiles(&[50, 50])[1], cumulative[61]);

        // Check equivalence between centiles and quartiles
        assert_eq!(
            hstats.bins_at_centiles(&[0, 25, 50, 75, 100]),
            hstats.bins_at_quartiles(&[0, 1, 2, 3, 4])
        );

        // Min, Max, Mean, and StdDev are tested by rolling-stats tests, so we don't need to test them here
    }

    #[test]
    fn stats_parallel() {
        type T = f64;

        // Define some constants
        const MEAN: T = 2.0;
        const STD_DEV: T = 3.0;
        const SEED: u64 = 42;
        const NUM_SAMPLES: usize = 10_000;
        const NUM_BINS: usize = 100;
        const START: T = -10.0;
        const END: T = 10.0;

        // let mut hstats = Hstats::new(START, END, NUM_BINS);

        let mut rng = rand::rngs::StdRng::seed_from_u64(SEED);

        let normal = Normal::<T>::new(MEAN, STD_DEV).unwrap();

        // Generate some random data
        let random_data: Vec<T> = (0..NUM_SAMPLES).map(|_x| normal.sample(&mut rng)).collect();

        // Update the stats

        let thread_count = rayon::current_num_threads() * 2;
        let chunk_size = random_data.len() / thread_count;

        // Update the stats for each chunk in parallel
        let hstats_list: Vec<Hstats<T>> = random_data
            .par_chunks(chunk_size)
            .map(|chunk| {
                let mut hstats = Hstats::new(START, END, NUM_BINS);
                chunk.iter().for_each(|v| hstats.add(*v));
                hstats
            })
            .collect();

        // Merge the stats from each chunk
        let merged = hstats_list
            .into_iter()
            .reduce(|hstats1, hstats2| hstats1.merge(&hstats2))
            .unwrap();

        // Check the standard deviation against the stats' standard deviation
        assert!(merged.std_dev().approx_eq(STD_DEV, (1.0e-2, 2)));

        // Check the mean against the stats' mean
        assert!(merged.mean().approx_eq(MEAN, (1.0e-1, 2)));

        // Check the count
        assert_eq!(merged.count(), random_data.len());

        // Quantiles on merged histogram should match a single-thread histogram
        let mut single = Hstats::new(START, END, NUM_BINS);
        random_data.iter().for_each(|v| single.add(*v));

        assert_eq!(
            merged.bins_at_centiles(&[0, 25, 50, 75, 99, 100]),
            single.bins_at_centiles(&[0, 25, 50, 75, 99, 100])
        );
    }
}