diff --git a/.github/workflows/benchmarks.yml b/.github/workflows/benchmarks.yml index 4d9b619..1e6720c 100644 --- a/.github/workflows/benchmarks.yml +++ b/.github/workflows/benchmarks.yml @@ -14,15 +14,18 @@ permissions: env: CARGO_TERM_COLOR: always - # Ten percent is a useful default for noisy shared CI runners. Adjust this - # repository-level value if the project needs a stricter or looser boundary. - BENCHER_THRESHOLD: "0.10" + # Forty percent keeps normal shared-runner noise from failing the report. + # Adjust this repository-level value if a stricter or looser boundary is needed. + BENCHER_THRESHOLD: "0.30" + # Keep every benchmark target in CI, but omit the largest input size. Local + # runs use all sizes unless this variable is set explicitly. + BENCHMARK_MAX_INPUT_SIZE: "1000" # A relative path is valid from the checkout and is available at workflow # parse time, unlike the runner context used by the previous configuration. CARGO_TARGET_DIR: target jobs: - benchmark: + benchmark_tests: name: Test and benchmark runs-on: ubuntu-latest steps: @@ -42,6 +45,33 @@ jobs: - name: Run functionality tests run: cargo test --locked + benchmark: + name: Benchmark (${{ matrix.target }}) + needs: benchmark_tests + runs-on: ubuntu-latest + strategy: + fail-fast: false + max-parallel: 4 + matrix: + target: + - sort + - summary_statistics + - deviation + - core_operations + steps: + - name: Check out source + uses: actions/checkout@v4 + with: + fetch-depth: 0 + + - name: Install Rust + uses: dtolnay/rust-toolchain@master + with: + toolchain: stable + + - name: Cache Rust build artifacts + uses: Swatinem/rust-cache@v2 + # Bencher Cloud requires a project and API key. Keeping the check in a # step makes forked pull requests safely fall back to artifacts because # GitHub does not expose repository secrets to them. @@ -68,6 +98,7 @@ jobs: env: BENCHER_API_KEY: ${{ secrets.BENCHER_API_KEY }} BENCHER_PROJECT: ${{ vars.BENCHER_PROJECT }} + BENCHER_CI_ID: ${{ matrix.target }} GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} run: | set -euo pipefail @@ -89,6 +120,7 @@ jobs: --threshold-upper-boundary "$BENCHER_THRESHOLD" --thresholds-reset --adapter rust_criterion + --ci-id "$BENCHER_CI_ID" --github-actions "$GITHUB_TOKEN" ) @@ -104,7 +136,7 @@ jobs: # Do not expose the Bencher key to the benchmarked crate. env -u BENCHER_API_KEY bencher "${bencher_args[@]}" \ - 'for bench in sort summary_statistics deviation core_operations; do cargo bench --locked --bench "$bench"; done' + 'cargo bench --locked --bench "$BENCHER_CI_ID"' - name: Fetch master for local Criterion comparison if: steps.bencher.outputs.enabled != 'true' && (github.event_name != 'push' || github.ref_name != 'master') @@ -117,18 +149,17 @@ jobs: set -euo pipefail master_worktree="$RUNNER_TEMP/ndarray-stats-master" - baseline_targets="$RUNNER_TEMP/criterion-baseline-targets" - : > "$baseline_targets" + baseline_target="$RUNNER_TEMP/criterion-baseline-target" + benchmark_target="${{ matrix.target }}" + : > "$baseline_target" git worktree add --detach "$master_worktree" origin/master - for bench in sort summary_statistics deviation core_operations; do - if grep -q "name = \"$bench\"" "$master_worktree/Cargo.toml"; then - cargo bench --manifest-path "$master_worktree/Cargo.toml" \ - --locked --bench "$bench" -- --save-baseline master - echo "$bench" >> "$baseline_targets" - else - echo "Skipping $bench: it is not present on master yet." - fi - done + if grep -q "name = \"$benchmark_target\"" "$master_worktree/Cargo.toml"; then + cargo bench --manifest-path "$master_worktree/Cargo.toml" \ + --locked --bench "$benchmark_target" -- --save-baseline master + echo "$benchmark_target" >> "$baseline_target" + else + echo "Skipping $benchmark_target: it is not present on master yet." + fi git worktree remove --force "$master_worktree" - name: Run local Criterion comparison @@ -137,23 +168,20 @@ jobs: run: | set -euo pipefail + benchmark_target="${{ matrix.target }}" if [[ "$GITHUB_EVENT_NAME" != "push" || "$GITHUB_REF_NAME" != "master" ]]; then - for bench in sort summary_statistics deviation core_operations; do - if grep -Fxq "$bench" "$RUNNER_TEMP/criterion-baseline-targets"; then - cargo bench --locked --bench "$bench" -- --baseline master 2>&1 | \ - tee -a "$RUNNER_TEMP/criterion-summary.txt" - else - echo "No master baseline for $bench; measuring without comparison." | \ - tee -a "$RUNNER_TEMP/criterion-summary.txt" - cargo bench --locked --bench "$bench" 2>&1 | \ - tee -a "$RUNNER_TEMP/criterion-summary.txt" - fi - done - else - for bench in sort summary_statistics deviation core_operations; do - cargo bench --locked --bench "$bench" 2>&1 | \ + if grep -Fxq "$benchmark_target" "$RUNNER_TEMP/criterion-baseline-target"; then + cargo bench --locked --bench "$benchmark_target" -- --baseline master 2>&1 | \ + tee -a "$RUNNER_TEMP/criterion-summary.txt" + else + echo "No master baseline for $benchmark_target; measuring without comparison." | \ tee -a "$RUNNER_TEMP/criterion-summary.txt" - done + cargo bench --locked --bench "$benchmark_target" 2>&1 | \ + tee -a "$RUNNER_TEMP/criterion-summary.txt" + fi + else + cargo bench --locked --bench "$benchmark_target" 2>&1 | \ + tee -a "$RUNNER_TEMP/criterion-summary.txt" fi - name: Create custom benchmark summary @@ -162,6 +190,7 @@ jobs: env: BENCHER_ENABLED: ${{ steps.bencher.outputs.enabled }} BENCHER_PROJECT: ${{ vars.BENCHER_PROJECT }} + BENCHER_CI_ID: ${{ matrix.target }} run: | set -euo pipefail @@ -297,13 +326,14 @@ jobs: context_row "Runner" "$RUNNER_OS / ubuntu-latest" context_row "Rust" "$rust_version" context_row "Features" "default" - context_row "Benchmark targets" "sort, summary_statistics, deviation, core_operations" + context_row "Benchmark target" "$BENCHER_CI_ID" context_row "Input shapes" "1-D n; correlation 3 x n; histogram n x 2" + context_row "Input size cap" "n <= $BENCHMARK_MAX_INPUT_SIZE" context_row "Backend" "$backend" context_row "Comparison" "$baseline_info" context_row "Bencher project" "${BENCHER_PROJECT:-not configured}" if [[ "$BENCHER_ENABLED" == "true" ]]; then - context_row "Threshold" "10% upper boundary" + context_row "Threshold" "40% upper boundary" else context_row "Threshold" "not enforced by artifact fallback" fi @@ -315,10 +345,55 @@ jobs: if: always() && steps.bencher.outputs.enabled != 'true' uses: actions/upload-artifact@v4 with: - name: criterion-${{ github.run_id }} + name: criterion-${{ github.run_id }}-${{ matrix.target }} path: | ${{ runner.temp }}/criterion-summary.txt ${{ runner.temp }}/custom-benchmark-summary.md target/criterion if-no-files-found: warn retention-days: 14 + + bencher_report: + name: Bencher Report + if: always() && github.event_name == 'pull_request' + needs: benchmark + runs-on: ubuntu-latest + steps: + # Bencher intentionally marks its GitHub Check as failed when it finds an + # alert. Keep those reports visible, but make alerts informational for PRs. + - name: Keep Bencher reports non-blocking + uses: actions/github-script@v7 + with: + github-token: ${{ secrets.GITHUB_TOKEN }} + script: | + try { + const { data } = await github.rest.checks.listForRef({ + owner: context.repo.owner, + repo: context.repo.repo, + ref: context.sha, + per_page: 100, + }); + const bencherReports = data.check_runs.filter(({ name }) => + name.startsWith('Bencher Report (') + ); + + for (const report of bencherReports) { + try { + await github.rest.checks.update({ + owner: context.repo.owner, + repo: context.repo.repo, + check_run_id: report.id, + conclusion: 'neutral', + }); + core.info(`${report.name} is informational for this pull request.`); + } catch (error) { + core.warning(`Could not update ${report.name}: ${error.message}`); + } + } + + if (bencherReports.length === 0) { + core.info(`No external Bencher report checks found on ${context.sha}.`); + } + } catch (error) { + core.warning(`Could not update Bencher report checks: ${error.message}`); + } diff --git a/benches/common/mod.rs b/benches/common/mod.rs new file mode 100644 index 0000000..4fbc6d6 --- /dev/null +++ b/benches/common/mod.rs @@ -0,0 +1,22 @@ +use std::env; + +/// Return the benchmark input sizes, optionally applying the CI size cap. +pub(crate) fn benchmark_lengths() -> Vec { + const DEFAULT_LENGTHS: [usize; 4] = [10, 100, 1_000, 10_000]; + + if let Ok(max_length) = env::var("BENCHMARK_MAX_INPUT_SIZE") { + if let Ok(max_length) = max_length.parse::() { + let lengths: Vec<_> = DEFAULT_LENGTHS + .iter() + .copied() + .filter(|length| *length <= max_length) + .collect(); + + if !lengths.is_empty() { + return lengths; + } + } + } + + DEFAULT_LENGTHS.to_vec() +} diff --git a/benches/core_operations.rs b/benches/core_operations.rs index 66be7a2..54b86bc 100644 --- a/benches/core_operations.rs +++ b/benches/core_operations.rs @@ -8,8 +8,10 @@ use ndarray_stats::histogram::{strategies::Auto, GridBuilder, HistogramExt}; use ndarray_stats::{interpolate::Linear, CorrelationExt, DeviationExt, EntropyExt, Quantile1dExt}; use noisy_float::types::n64; +mod common; + fn mean(c: &mut Criterion) { - let lens = vec![10, 100, 1000, 10000]; + let lens = common::benchmark_lengths(); let mut group = c.benchmark_group("mean"); group.plot_config(PlotConfiguration::default().summary_scale(AxisScale::Logarithmic)); for len in &lens { @@ -23,7 +25,7 @@ fn mean(c: &mut Criterion) { } fn quantiles_mut(c: &mut Criterion) { - let lens = vec![10, 100, 1000, 10000]; + let lens = common::benchmark_lengths(); let quantile_indexes = Array1::from_vec(vec![n64(0.25), n64(0.5), n64(0.75)]); let mut group = c.benchmark_group("quantiles_mut"); group.plot_config(PlotConfiguration::default().summary_scale(AxisScale::Logarithmic)); @@ -46,7 +48,7 @@ fn quantiles_mut(c: &mut Criterion) { } fn pearson_correlation(c: &mut Criterion) { - let lens = vec![10, 100, 1000, 10000]; + let lens = common::benchmark_lengths(); let mut group = c.benchmark_group("pearson_correlation"); group.plot_config(PlotConfiguration::default().summary_scale(AxisScale::Logarithmic)); for len in &lens { @@ -59,7 +61,7 @@ fn pearson_correlation(c: &mut Criterion) { } fn spearman_correlation(c: &mut Criterion) { - let lens = vec![10, 100, 1000, 10000]; + let lens = common::benchmark_lengths(); let mut group = c.benchmark_group("spearman_correlation"); group.plot_config(PlotConfiguration::default().summary_scale(AxisScale::Logarithmic)); for len in &lens { @@ -74,7 +76,7 @@ fn spearman_correlation(c: &mut Criterion) { } fn kendall_tau(c: &mut Criterion) { - let lens = vec![10, 100, 1000, 10000]; + let lens = common::benchmark_lengths(); let mut group = c.benchmark_group("kendall_tau"); group.plot_config(PlotConfiguration::default().summary_scale(AxisScale::Logarithmic)); for len in &lens { @@ -89,7 +91,7 @@ fn kendall_tau(c: &mut Criterion) { } fn entropy(c: &mut Criterion) { - let lens = vec![10, 100, 1000, 10000]; + let lens = common::benchmark_lengths(); let mut group = c.benchmark_group("entropy"); group.plot_config(PlotConfiguration::default().summary_scale(AxisScale::Logarithmic)); for len in &lens { @@ -103,7 +105,7 @@ fn entropy(c: &mut Criterion) { } fn histogram(c: &mut Criterion) { - let lens = vec![10, 100, 1000, 10000]; + let lens = common::benchmark_lengths(); let mut group = c.benchmark_group("histogram"); group.plot_config(PlotConfiguration::default().summary_scale(AxisScale::Logarithmic)); for len in &lens { @@ -123,7 +125,7 @@ fn histogram(c: &mut Criterion) { } fn l1_dist(c: &mut Criterion) { - let lens = vec![10, 100, 1000, 10000]; + let lens = common::benchmark_lengths(); let mut group = c.benchmark_group("l1_dist"); group.plot_config(PlotConfiguration::default().summary_scale(AxisScale::Logarithmic)); for len in &lens { diff --git a/benches/deviation.rs b/benches/deviation.rs index 2cd9b91..3eee04c 100644 --- a/benches/deviation.rs +++ b/benches/deviation.rs @@ -6,8 +6,10 @@ use ndarray_rand::rand_distr::Uniform; use ndarray_rand::RandomExt; use ndarray_stats::DeviationExt; +mod common; + fn sq_l2_dist(c: &mut Criterion) { - let lens = vec![10, 100, 1000, 10000]; + let lens = common::benchmark_lengths(); let mut group = c.benchmark_group("sq_l2_dist"); group.plot_config(PlotConfiguration::default().summary_scale(AxisScale::Logarithmic)); for len in &lens { diff --git a/benches/sort.rs b/benches/sort.rs index 1a2f442..ba7e9ee 100644 --- a/benches/sort.rs +++ b/benches/sort.rs @@ -5,8 +5,10 @@ use ndarray::prelude::*; use ndarray_stats::Sort1dExt; use rand::prelude::*; +mod common; + fn get_from_sorted_mut(c: &mut Criterion) { - let lens = vec![10, 100, 1000, 10000]; + let lens = common::benchmark_lengths(); let mut group = c.benchmark_group("get_from_sorted_mut"); group.plot_config(PlotConfiguration::default().summary_scale(AxisScale::Logarithmic)); for len in &lens { @@ -30,7 +32,7 @@ fn get_from_sorted_mut(c: &mut Criterion) { } fn get_many_from_sorted_mut(c: &mut Criterion) { - let lens = vec![10, 100, 1000, 10000]; + let lens = common::benchmark_lengths(); let mut group = c.benchmark_group("get_many_from_sorted_mut"); group.plot_config(PlotConfiguration::default().summary_scale(AxisScale::Logarithmic)); for len in &lens { diff --git a/benches/summary_statistics.rs b/benches/summary_statistics.rs index 64f22be..b1b607b 100644 --- a/benches/summary_statistics.rs +++ b/benches/summary_statistics.rs @@ -4,10 +4,49 @@ use criterion::{ use ndarray::prelude::*; use ndarray_rand::rand_distr::Uniform; use ndarray_rand::RandomExt; -use ndarray_stats::SummaryStatisticsExt; +use ndarray_stats::{DescriptiveStatistics, QuantileExt, SummaryStatisticsExt}; + +mod common; + +fn score_f64(summary: &DescriptiveStatistics) -> f64 { + summary.count() as f64 + + summary.mean() + + summary.min() + + summary.max() + + summary.population_variance() + + summary.sample_variance() + + summary.population_std() + + summary.sample_std() +} + +fn score_f32(summary: &DescriptiveStatistics) -> f32 { + summary.count() as f32 + + summary.mean() + + summary.min() + + summary.max() + + summary.population_variance() + + summary.sample_variance() + + summary.population_std() + + summary.sample_std() +} + +fn score_repeated_axis(data: &Array2, axis: Axis, weights: &Array1) -> f64 { + data.lanes(axis) + .into_iter() + .map(|lane| { + lane.mean().unwrap() + + *lane.min().unwrap() + + *lane.max().unwrap() + + lane.weighted_var(weights, 0.0).unwrap() + + lane.weighted_var(weights, 1.0).unwrap() + + lane.weighted_std(weights, 0.0).unwrap() + + lane.weighted_std(weights, 1.0).unwrap() + }) + .sum() +} fn weighted_std(c: &mut Criterion) { - let lens = vec![10, 100, 1000, 10000]; + let lens = common::benchmark_lengths(); let mut group = c.benchmark_group("weighted_std"); group.plot_config(PlotConfiguration::default().summary_scale(AxisScale::Logarithmic)); for len in &lens { @@ -27,9 +66,105 @@ fn weighted_std(c: &mut Criterion) { group.finish(); } +fn descriptive_statistics(c: &mut Criterion) { + let lens = common::benchmark_lengths(); + let mut group = c.benchmark_group("descriptive_statistics"); + group.plot_config(PlotConfiguration::default().summary_scale(AxisScale::Logarithmic)); + + for &len in &lens { + let data = Array::random(len, Uniform::new(0.0, 1.0).unwrap()); + let weights = Array1::ones(len); + + group.bench_with_input(format!("fused/{len}"), &data, |b, data| { + b.iter(|| { + let summary = black_box(data.descriptive_statistics().unwrap()); + black_box(score_f64(&summary)); + }) + }); + + group.bench_with_input(format!("repeated/{len}"), &data, |b, data| { + b.iter(|| { + let result = ( + data.mean().unwrap(), + *data.min().unwrap(), + *data.max().unwrap(), + data.weighted_var(&weights, 0.0).unwrap(), + data.weighted_var(&weights, 1.0).unwrap(), + data.weighted_std(&weights, 0.0).unwrap(), + data.weighted_std(&weights, 1.0).unwrap(), + ); + black_box(result); + }) + }); + } + + group.finish(); +} + +fn descriptive_statistics_axis(c: &mut Criterion) { + let lens = common::benchmark_lengths(); + let mut group = c.benchmark_group("descriptive_statistics_axis"); + group.plot_config(PlotConfiguration::default().summary_scale(AxisScale::Logarithmic)); + + for &len in &lens { + let data = Array::random((len, 8), Uniform::new(0.0, 1.0).unwrap()); + let axis_zero_weights = Array1::ones(len); + let axis_one_weights = Array1::ones(8); + + group.bench_with_input(format!("fused_axis0/{len}"), &data, |b, data| { + b.iter(|| { + let summaries = black_box(data.descriptive_statistics_axis(Axis(0)).unwrap()); + let score = summaries.iter().map(score_f64).sum::(); + black_box(score); + }) + }); + + group.bench_with_input(format!("repeated_axis0/{len}"), &data, |b, data| { + b.iter(|| { + black_box(score_repeated_axis(data, Axis(0), &axis_zero_weights)); + }) + }); + + group.bench_with_input(format!("fused_axis1/{len}"), &data, |b, data| { + b.iter(|| { + let summaries = black_box(data.descriptive_statistics_axis(Axis(1)).unwrap()); + let score = summaries.iter().map(score_f64).sum::(); + black_box(score); + }) + }); + + group.bench_with_input(format!("repeated_axis1/{len}"), &data, |b, data| { + b.iter(|| { + black_box(score_repeated_axis(data, Axis(1), &axis_one_weights)); + }) + }); + } + + group.finish(); +} + +fn descriptive_statistics_f32(c: &mut Criterion) { + let lens = common::benchmark_lengths(); + let mut group = c.benchmark_group("descriptive_statistics_f32"); + group.plot_config(PlotConfiguration::default().summary_scale(AxisScale::Logarithmic)); + + for &len in &lens { + group.bench_function(format!("fused/{len}"), |b| { + let data: Array1 = Array::random(len, Uniform::new(0.0, 1.0).unwrap()); + b.iter(|| { + let summary = black_box(data.descriptive_statistics().unwrap()); + black_box(score_f32(&summary)); + }) + }); + } + + group.finish(); +} + criterion_group! { name = benches; config = Criterion::default(); - targets = weighted_std + targets = weighted_std, descriptive_statistics, descriptive_statistics_axis, + descriptive_statistics_f32 } criterion_main!(benches); diff --git a/src/errors.rs b/src/errors.rs index e2617f3..1ad91d6 100644 --- a/src/errors.rs +++ b/src/errors.rs @@ -15,6 +15,34 @@ impl fmt::Display for EmptyInput { impl Error for EmptyInput {} +/// An error returned when computing a descriptive-statistics summary. +#[derive(Clone, Copy, Debug, Eq, PartialEq)] +pub enum SummaryStatisticsError { + /// The input array, or one of its summary lanes, was empty. + EmptyInput, + /// A pairwise ordering required for the minimum or maximum was undefined. + UndefinedOrder, +} + +impl fmt::Display for SummaryStatisticsError { + fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { + match self { + SummaryStatisticsError::EmptyInput => write!(f, "Empty input."), + SummaryStatisticsError::UndefinedOrder => { + write!(f, "Undefined ordering between a tested pair of values.") + } + } + } +} + +impl Error for SummaryStatisticsError {} + +impl From for SummaryStatisticsError { + fn from(_: EmptyInput) -> Self { + SummaryStatisticsError::EmptyInput + } +} + /// An error computing a minimum/maximum value. #[derive(Clone, Debug, Eq, PartialEq)] pub enum MinMaxError { diff --git a/src/lib.rs b/src/lib.rs index b02ddad..21369a3 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -36,7 +36,7 @@ pub use crate::histogram::HistogramExt; pub use crate::maybe_nan::{MaybeNan, MaybeNanExt}; pub use crate::quantile::{interpolate, Quantile1dExt, QuantileExt}; pub use crate::sort::Sort1dExt; -pub use crate::summary_statistics::SummaryStatisticsExt; +pub use crate::summary_statistics::{DescriptiveStatistics, SummaryStatisticsExt}; #[cfg(test)] #[macro_use] diff --git a/src/summary_statistics/descriptive.rs b/src/summary_statistics/descriptive.rs new file mode 100644 index 0000000..c0de50f --- /dev/null +++ b/src/summary_statistics/descriptive.rs @@ -0,0 +1,141 @@ +use crate::errors::SummaryStatisticsError; +use num_traits::{Float, FromPrimitive}; +use std::cmp::Ordering; + +/// A fused descriptive-statistics summary for a non-empty floating-point input. +/// +/// The value stores the observation count, mean, minimum, maximum, and the +/// second central-moment accumulator used to derive population and sample +/// variance. Construct values with +/// [`crate::SummaryStatisticsExt::descriptive_statistics`] or +/// [`crate::SummaryStatisticsExt::descriptive_statistics_axis`]. +#[derive(Clone, Copy, Debug, PartialEq)] +pub struct DescriptiveStatistics { + count: usize, + mean: A, + m2: A, + min: A, + max: A, +} + +impl DescriptiveStatistics +where + A: Float + FromPrimitive, +{ + /// Returns the number of observations represented by this summary. + pub fn count(&self) -> usize { + self.count + } + + /// Returns the arithmetic mean. + pub fn mean(&self) -> A { + self.mean + } + + /// Returns the minimum observation. + pub fn min(&self) -> A { + self.min + } + + /// Returns the maximum observation. + pub fn max(&self) -> A { + self.max + } + + /// Returns the population variance, dividing the second-moment + /// accumulator by the number of observations. + pub fn population_variance(&self) -> A { + self.m2 / Self::from_usize(self.count) + } + + /// Returns the sample variance, dividing the second-moment accumulator by + /// one fewer than the number of observations. + /// + /// For a one-observation summary, this follows floating-point division + /// semantics and is therefore `NaN`. + pub fn sample_variance(&self) -> A { + self.m2 / Self::from_usize(self.count - 1) + } + + /// Returns the population standard deviation. + pub fn population_std(&self) -> A { + self.population_variance().sqrt() + } + + /// Returns the sample standard deviation. + /// + /// For a one-observation summary, this follows floating-point division + /// semantics and is therefore `NaN`. + pub fn sample_std(&self) -> A { + self.sample_variance().sqrt() + } + + pub(super) fn from_iter(values: I) -> Result + where + I: IntoIterator, + { + let mut values = values.into_iter(); + let first = values.next().ok_or(SummaryStatisticsError::EmptyInput)?; + let mut accumulator = Accumulator { + count: 1, + mean: first, + m2: A::zero(), + min: first, + max: first, + }; + + for value in values { + accumulator.update(value)?; + } + + Ok(accumulator.finish()) + } + + fn from_usize(value: usize) -> A { + A::from_usize(value).expect("Converting an observation count to `A` must not fail.") + } +} + +struct Accumulator { + count: usize, + mean: A, + m2: A, + min: A, + max: A, +} + +impl Accumulator +where + A: Float + FromPrimitive, +{ + fn update(&mut self, value: A) -> Result<(), SummaryStatisticsError> { + self.count += 1; + let count = DescriptiveStatistics::::from_usize(self.count); + let delta = value - self.mean; + self.mean = self.mean + delta / count; + self.m2 = self.m2 + delta * (value - self.mean); + + match value.partial_cmp(&self.min) { + Some(Ordering::Less) => self.min = value, + Some(_) => {} + None => return Err(SummaryStatisticsError::UndefinedOrder), + } + match value.partial_cmp(&self.max) { + Some(Ordering::Greater) => self.max = value, + Some(_) => {} + None => return Err(SummaryStatisticsError::UndefinedOrder), + } + + Ok(()) + } + + fn finish(self) -> DescriptiveStatistics { + DescriptiveStatistics { + count: self.count, + mean: self.mean, + m2: self.m2, + min: self.min, + max: self.max, + } + } +} diff --git a/src/summary_statistics/means.rs b/src/summary_statistics/means.rs index a2e8ab0..c7ed653 100644 --- a/src/summary_statistics/means.rs +++ b/src/summary_statistics/means.rs @@ -1,5 +1,6 @@ +use super::DescriptiveStatistics; use super::SummaryStatisticsExt; -use crate::errors::{EmptyInput, MultiInputError, ShapeMismatch}; +use crate::errors::{EmptyInput, MultiInputError, ShapeMismatch, SummaryStatisticsError}; use ndarray::{Array, ArrayBase, ArrayRef, Axis, Data, Dimension, Ix1, RemoveAxis}; use num_integer::IterBinomial; use num_traits::{Float, FromPrimitive, Zero}; @@ -9,6 +10,36 @@ impl SummaryStatisticsExt for ArrayRef where D: Dimension, { + fn descriptive_statistics(&self) -> Result, SummaryStatisticsError> + where + A: Float + FromPrimitive, + { + DescriptiveStatistics::from_iter(self.iter().copied()) + } + + fn descriptive_statistics_axis( + &self, + axis: Axis, + ) -> Result, D::Smaller>, SummaryStatisticsError> + where + A: Float + FromPrimitive, + D: RemoveAxis, + { + if self.is_empty() { + return Err(SummaryStatisticsError::EmptyInput); + } + + let shape = self.raw_dim().remove_axis(axis); + let summaries = self + .lanes(axis) + .into_iter() + .map(|lane| DescriptiveStatistics::from_iter(lane.iter().copied())) + .collect::, _>>()?; + + Ok(Array::from_shape_vec(shape, summaries) + .expect("descriptive-statistics lanes must match the output shape")) + } + fn mean(&self) -> Result where A: Clone + FromPrimitive + Add + Div + Zero, diff --git a/src/summary_statistics/mod.rs b/src/summary_statistics/mod.rs index 239440e..ca570da 100644 --- a/src/summary_statistics/mod.rs +++ b/src/summary_statistics/mod.rs @@ -1,15 +1,41 @@ //! Summary statistics (e.g. mean, variance, etc.). -use crate::errors::{EmptyInput, MultiInputError}; +use crate::errors::{EmptyInput, MultiInputError, SummaryStatisticsError}; use ndarray::{Array, ArrayRef, Axis, Dimension, Ix1, RemoveAxis}; use num_traits::{Float, FromPrimitive, Zero}; use std::ops::{Add, AddAssign, Div, Mul}; +mod descriptive; + +pub use self::descriptive::DescriptiveStatistics; + /// Extension trait for `ArrayRef` providing methods /// to compute several summary statistics (e.g. mean, variance, etc.). pub trait SummaryStatisticsExt where D: Dimension, { + /// Returns a fused descriptive-statistics summary of all elements in the array. + /// + /// If the array is empty, `SummaryStatisticsError::EmptyInput` is returned. + /// If a required minimum or maximum comparison has undefined ordering, + /// `SummaryStatisticsError::UndefinedOrder` is returned. + fn descriptive_statistics(&self) -> Result, SummaryStatisticsError> + where + A: Float + FromPrimitive; + + /// Returns a descriptive-statistics summary for every lane along `axis`. + /// + /// The returned array has the input shape with `axis` removed. The method + /// panics if `axis` is out of bounds and returns the first summary error + /// encountered while processing the lanes. + fn descriptive_statistics_axis( + &self, + axis: Axis, + ) -> Result, D::Smaller>, SummaryStatisticsError> + where + A: Float + FromPrimitive, + D: RemoveAxis; + /// Returns the [`arithmetic mean`] x̅ of all elements in the array: /// /// ```text diff --git a/tests/summary_statistics.rs b/tests/summary_statistics.rs index af6256a..7c5daeb 100644 --- a/tests/summary_statistics.rs +++ b/tests/summary_statistics.rs @@ -3,13 +3,153 @@ use ndarray::{arr0, array, Array, Array1, Array2, Axis}; use ndarray_rand::rand_distr::Uniform; use ndarray_rand::RandomExt; use ndarray_stats::{ - errors::{EmptyInput, MultiInputError, ShapeMismatch}, + errors::{EmptyInput, MultiInputError, ShapeMismatch, SummaryStatisticsError}, SummaryStatisticsExt, }; use noisy_float::types::N64; use quickcheck::{quickcheck, TestResult}; use std::f64; +#[test] +fn descriptive_statistics_known_values() { + let a = array![1.0, 2.0, 2.0, 3.0]; + let summary = a.descriptive_statistics().unwrap(); + + assert_eq!(summary.count(), 4); + assert_eq!(summary.min(), 1.0); + assert_eq!(summary.max(), 3.0); + assert_abs_diff_eq!(summary.mean(), 2.0, epsilon = 1e-12); + assert_abs_diff_eq!(summary.population_variance(), 0.5, epsilon = 1e-12); + assert_abs_diff_eq!(summary.sample_variance(), 2.0 / 3.0, epsilon = 1e-12); + assert_abs_diff_eq!(summary.population_std(), 0.5_f64.sqrt(), epsilon = 1e-12); + assert_abs_diff_eq!( + summary.sample_std(), + (2.0_f64 / 3.0).sqrt(), + epsilon = 1e-12 + ); +} + +#[test] +fn descriptive_statistics_supports_f32() { + let summary = array![1.0_f32, 2.0, 3.0].descriptive_statistics().unwrap(); + + assert_eq!(summary.count(), 3); + assert_abs_diff_eq!(summary.mean(), 2.0_f32, epsilon = 1e-6); + assert_abs_diff_eq!(summary.population_variance(), 2.0_f32 / 3.0, epsilon = 1e-6); + assert_abs_diff_eq!(summary.sample_variance(), 1.0_f32, epsilon = 1e-6); +} + +#[test] +fn descriptive_statistics_reports_empty_input() { + let a: Array1 = array![]; + + assert_eq!( + a.descriptive_statistics(), + Err(SummaryStatisticsError::EmptyInput) + ); + assert_eq!( + a.descriptive_statistics_axis(Axis(0)), + Err(SummaryStatisticsError::EmptyInput) + ); +} + +#[test] +fn descriptive_statistics_preserves_undefined_order_behavior() { + let a = array![1.0, f64::NAN]; + assert_eq!( + a.descriptive_statistics(), + Err(SummaryStatisticsError::UndefinedOrder) + ); + assert_eq!( + a.descriptive_statistics_axis(Axis(0)), + Err(SummaryStatisticsError::UndefinedOrder) + ); + + let singleton = array![f64::NAN]; + let summary = singleton.descriptive_statistics().unwrap(); + assert!(summary.mean().is_nan()); + assert!(summary.min().is_nan()); + assert!(summary.max().is_nan()); +} + +#[test] +fn descriptive_statistics_sample_one_observation_uses_float_semantics() { + let summary = array![42.0_f64].descriptive_statistics().unwrap(); + + assert_eq!(summary.population_variance(), 0.0); + assert_eq!(summary.population_std(), 0.0); + assert!(summary.sample_variance().is_nan()); + assert!(summary.sample_std().is_nan()); +} + +#[test] +fn descriptive_statistics_axis_preserves_shape_and_values() { + let a = array![[1.0, 2.0, 3.0], [4.0, 5.0, 6.0]]; + + let columns = a.descriptive_statistics_axis(Axis(0)).unwrap(); + assert_eq!(columns.shape(), &[3]); + for (summary, (mean, min, max)) in + columns + .iter() + .zip([(2.5, 1.0, 4.0), (3.5, 2.0, 5.0), (4.5, 3.0, 6.0)]) + { + assert_eq!(summary.count(), 2); + assert_abs_diff_eq!(summary.mean(), mean, epsilon = 1e-12); + assert_abs_diff_eq!(summary.min(), min, epsilon = 1e-12); + assert_abs_diff_eq!(summary.max(), max, epsilon = 1e-12); + assert_abs_diff_eq!(summary.population_variance(), 2.25, epsilon = 1e-12); + assert_abs_diff_eq!(summary.sample_variance(), 4.5, epsilon = 1e-12); + } + + let rows = a.descriptive_statistics_axis(Axis(1)).unwrap(); + assert_eq!(rows.shape(), &[2]); + assert_abs_diff_eq!(rows[0].mean(), 2.0, epsilon = 1e-12); + assert_abs_diff_eq!(rows[1].mean(), 5.0, epsilon = 1e-12); + assert_eq!(rows[0].min(), 1.0); + assert_eq!(rows[1].max(), 6.0); +} + +#[test] +fn descriptive_statistics_matches_reference_for_finite_values() { + fn prop(values: Vec) -> TestResult { + if values.is_empty() { + return TestResult::discard(); + } + + let values: Vec = values.into_iter().map(|value| value % 100.0).collect(); + let array = Array1::from(values.clone()); + let summary = array.descriptive_statistics().unwrap(); + let count = values.len() as f64; + let mean = values.iter().sum::() / count; + let sum_squared = values + .iter() + .map(|value| (value - mean).powi(2)) + .sum::(); + let min = values.iter().copied().fold(f64::INFINITY, f64::min); + let max = values.iter().copied().fold(f64::NEG_INFINITY, f64::max); + + TestResult::from_bool( + summary.count() == values.len() + && abs_diff_eq!(summary.mean(), mean, epsilon = 1e-10) + && abs_diff_eq!( + summary.population_variance(), + sum_squared / count, + epsilon = 1e-8 + ) + && (summary.sample_variance().is_nan() + || abs_diff_eq!( + summary.sample_variance(), + sum_squared / (count - 1.0), + epsilon = 1e-8 + )) + && summary.min() == min + && summary.max() == max, + ) + } + + quickcheck(prop as fn(Vec) -> TestResult); +} + #[test] fn test_with_nan_values() { let a = array![f64::NAN, 1.];