rustnum 0.1.0

NumPy-style n-dimensional arrays for Rust: one array type, NumPy names, pure Rust.
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104

//! Shared assertions for the NumPy-parity tests.
//!
//! Every test creates a [`Report`] first; failed comparisons are collected and printed together
//! when the report is dropped, so one run shows all disagreements with NumPy.

use std::cell::RefCell;

use super::complex::Complex;
use crate::Float;

thread_local! {
    static FAILS: RefCell<Vec<String>> = const { RefCell::new(Vec::new()) };
}

fn fail(msg: String) {
    FAILS.with(|f| f.borrow_mut().push(msg));
}

/// Panics on drop if any comparison of this test failed.
pub struct Report;

impl Drop for Report {
    fn drop(&mut self) {
        if std::thread::panicking() {
            return;
        }
        let fails = FAILS.with(|f| core::mem::take(&mut *f.borrow_mut()));
        if !fails.is_empty() {
            let shown = fails
                .iter()
                .take(40)
                .cloned()
                .collect::<Vec<_>>()
                .join("\n");
            panic!(
                "{} comparison(s) disagree with NumPy:\n{shown}",
                fails.len()
            );
        }
    }
}

/// Compare `got` with NumPy's `want`.
///
/// NaN matches NaN (payload ignored). Infinities and zeros must match exactly, including the
/// sign of zero. Other values must agree within `atol + rtol * |want|`; with `rtol = atol = 0`
/// the bit patterns must be identical.
pub fn check<F: Float>(what: &str, input: &str, got: F, want: F, rtol: f64, atol: f64) {
    let ok = if want.is_nan() {
        got.is_nan()
    } else if want.is_infinite() || want == F::ZERO || (rtol == 0.0 && atol == 0.0) {
        got.to_bits64() == want.to_bits64()
    } else {
        let (g, w) = (got.to_f64(), want.to_f64());
        (g - w).abs() <= atol + rtol * w.abs()
    };
    if !ok {
        fail(format!("{what}({input}): got {got:?}, NumPy says {want:?}"));
    }
}

/// Like [`check`] for the two parts of a complex result.
///
/// Finite parts are compared against the modulus of the expected value, so cancellation in one
/// part does not cause false alarms. When `loose_signs` is set (used for non-finite inputs,
/// where C99 Annex G leaves the sign of a zero or infinity unspecified and libms disagree),
/// zeros and infinities only have to match in magnitude.
pub fn check_c<F: Float>(
    what: &str,
    input: &str,
    got: Complex<F>,
    want: (F, F),
    rtol: f64,
    loose_signs: bool,
) {
    let norm = if want.0.is_finite() && want.1.is_finite() {
        want.0.to_f64().hypot(want.1.to_f64())
    } else {
        0.0
    };
    for (part, g, w) in [("re", got.re, want.0), ("im", got.im, want.1)] {
        let ok = if w.is_nan() {
            g.is_nan()
        } else if w.is_infinite() || w == F::ZERO || rtol == 0.0 {
            if loose_signs {
                g.abs().to_bits64() == w.abs().to_bits64()
            } else {
                g.to_bits64() == w.to_bits64()
            }
        } else {
            (g.to_f64() - w.to_f64()).abs() <= rtol * norm.max(w.to_f64().abs())
        };
        if !ok {
            fail(format!(
                "{what}({input}).{part}: got {g:?}, NumPy says {w:?} (full got {got:?}, want {want:?})"
            ));
        }
    }
}

/// True when both parts of `z` are finite (the strict-sign regime).
pub fn finite_input<F: Float>(z: Complex<F>) -> bool {
    z.re.is_finite() && z.im.is_finite()
}