axiolid_numeric/error.rs
1//! Refusals and convergence status shared by every routine.
2
3use core::fmt;
4
5/// Why a numeric routine refused its input or could not produce a value.
6///
7/// Every variant names the offending argument or condition. A routine never
8/// substitutes a default, clamps silently, or returns NaN in place of one of
9/// these.
10#[derive(Debug, Clone, PartialEq)]
11#[non_exhaustive]
12pub enum NumericError {
13 /// An input value was NaN or infinite. `name` is the argument.
14 NonFiniteInput {
15 /// The argument that was not finite, such as `"lower bound"`.
16 name: &'static str,
17 },
18 /// A caller-supplied function returned NaN or infinity.
19 NonFiniteEvaluation {
20 /// Which function, such as `"integrand"` or `"residuals"`.
21 name: &'static str,
22 /// The scalar argument it was evaluated at, when it has one.
23 at: Option<f64>,
24 },
25 /// An argument was finite but outside its valid range.
26 InvalidArgument {
27 /// The argument, such as `"x tolerance"`.
28 name: &'static str,
29 /// What it must satisfy.
30 reason: &'static str,
31 },
32 /// Sizes of two arguments disagree.
33 DimensionMismatch {
34 /// The argument whose size is wrong.
35 name: &'static str,
36 /// The size it must have.
37 expected: usize,
38 /// The size it has.
39 found: usize,
40 },
41 /// The function has the same sign at both ends of the bracket, so a
42 /// root is not guaranteed between them.
43 NoSignChange {
44 /// Value at the lower end.
45 f_lower: f64,
46 /// Value at the upper end.
47 f_upper: f64,
48 },
49 /// Every coefficient of the polynomial is zero: every point is a root.
50 ZeroPolynomial,
51 /// A square matrix is singular to working precision.
52 Singular {
53 /// Which matrix, such as `"system matrix"`.
54 name: &'static str,
55 },
56 /// A least-squares matrix has fewer independent columns than unknowns,
57 /// so the solution is not unique.
58 RankDeficient {
59 /// Which matrix.
60 name: &'static str,
61 /// Numerical rank found.
62 rank: usize,
63 /// Rank required (the number of unknowns or constraints).
64 required: usize,
65 },
66 /// A matrix given to a symmetric solver is not symmetric.
67 NotSymmetric {
68 /// Row of the first asymmetric pair found.
69 row: usize,
70 /// Column of the first asymmetric pair found.
71 column: usize,
72 },
73 /// A symmetric matrix is not positive definite to working precision.
74 NotPositiveDefinite {
75 /// Index of the first non-positive pivot.
76 pivot: usize,
77 },
78 /// An iterative routine stopped without meeting its tolerance. Returned
79 /// only by the `converged()` accessors; the routines themselves report
80 /// this case in their result's status with the estimate they reached.
81 NotConverged {
82 /// Which routine.
83 name: &'static str,
84 /// The error estimate it reached.
85 error_estimate: f64,
86 },
87}
88
89impl fmt::Display for NumericError {
90 fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
91 match self {
92 Self::NonFiniteInput { name } => write!(f, "{name} is not finite"),
93 Self::NonFiniteEvaluation { name, at: Some(x) } => {
94 write!(f, "{name} returned a non-finite value at {x}")
95 }
96 Self::NonFiniteEvaluation { name, at: None } => {
97 write!(f, "{name} returned a non-finite value")
98 }
99 Self::InvalidArgument { name, reason } => write!(f, "{name} {reason}"),
100 Self::DimensionMismatch {
101 name,
102 expected,
103 found,
104 } => write!(f, "{name} has size {found}, expected {expected}"),
105 Self::NoSignChange { f_lower, f_upper } => write!(
106 f,
107 "function values {f_lower} and {f_upper} at the bracket ends do not change sign"
108 ),
109 Self::ZeroPolynomial => f.write_str("polynomial is identically zero"),
110 Self::Singular { name } => write!(f, "{name} is singular to working precision"),
111 Self::RankDeficient {
112 name,
113 rank,
114 required,
115 } => write!(f, "{name} has numerical rank {rank}, {required} required"),
116 Self::NotSymmetric { row, column } => {
117 write!(f, "matrix is not symmetric at ({row}, {column})")
118 }
119 Self::NotPositiveDefinite { pivot } => {
120 write!(f, "matrix is not positive definite (pivot {pivot})")
121 }
122 Self::NotConverged {
123 name,
124 error_estimate,
125 } => write!(
126 f,
127 "{name} did not converge (error estimate {error_estimate})"
128 ),
129 }
130 }
131}
132
133impl std::error::Error for NumericError {}
134
135/// Result of a numeric routine.
136pub type NumericResult<T> = Result<T, NumericError>;
137
138/// How an iterative routine stopped.
139#[derive(Debug, Clone, Copy, PartialEq, Eq)]
140pub enum Status {
141 /// The requested tolerance was met.
142 Converged,
143 /// The iteration, evaluation or subdivision budget ran out first. The
144 /// result still carries the best value and its error estimate.
145 BudgetExhausted,
146 /// Rounding error stopped further progress before the tolerance was
147 /// met; the error estimate is the best `f64` can certify here.
148 RoundoffLimited,
149}
150
151/// Refuse a non-finite scalar by name.
152pub(crate) fn finite(value: f64, name: &'static str) -> NumericResult<f64> {
153 if value.is_finite() {
154 Ok(value)
155 } else {
156 Err(NumericError::NonFiniteInput { name })
157 }
158}
159
160/// Refuse a slice containing a non-finite value by name.
161pub(crate) fn all_finite(values: &[f64], name: &'static str) -> NumericResult<()> {
162 if values.iter().all(|v| v.is_finite()) {
163 Ok(())
164 } else {
165 Err(NumericError::NonFiniteInput { name })
166 }
167}
168
169/// Refuse a tolerance that is negative or not finite.
170pub(crate) fn tolerance(value: f64, name: &'static str) -> NumericResult<f64> {
171 finite(value, name)?;
172 if value < 0.0 {
173 return Err(NumericError::InvalidArgument {
174 name,
175 reason: "must be non-negative",
176 });
177 }
178 Ok(value)
179}