From 940f180f620daad20de728f45e910c95a8300458 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Isa=C3=AFe?= Date: Thu, 1 Aug 2024 07:07:00 +0200 Subject: [PATCH] refactor: remove `Integraal::compute_error` (#28) * delete `Integraal::compute_error` * update tests * rewrite `almost_equal!` macro * delete `is_within_tolerance` --- integraal/src/structure/implementations.rs | 82 ---------------------- integraal/src/structure/tests.rs | 61 ++++++++-------- 2 files changed, 30 insertions(+), 113 deletions(-) diff --git a/integraal/src/structure/implementations.rs b/integraal/src/structure/implementations.rs index 24b7231..f1d1fc3 100644 --- a/integraal/src/structure/implementations.rs +++ b/integraal/src/structure/implementations.rs @@ -87,88 +87,6 @@ impl<'a, X: Scalar> Integraal<'a, X> { self.function = None; // is this really useful? we could wire returns directly using `?` if this wasn't here Ok(res) } - - #[allow( - clippy::missing_errors_doc, - clippy::missing_panics_doc, - clippy::too_many_lines - )] - /// This method attempts to compute the numerical error one can expect from the approximation. - /// - /// Error formulae were taken from the respective methods' Wikipedia pages. - /// - /// # Return / Errors - /// - /// This method returns a `Result` taking the following values: - /// - `Ok(X: Scalar)` -- The computation was successfuly done - /// - `Err(IntegraalError)` -- The computation failed for the reason specified by the enum. - pub fn compute_error(&self) -> Result { - if self.domain.is_none() | self.function.is_none() | self.method.is_none() { - return Err(IntegraalError::MissingParameters( - "one or more parameter is missing", - )); - } - if let Some(DomainDescriptor::Uniform { - start, - step, - n_step, - }) = self.domain - { - let res: X = match self.method { - // ref: https://en.wikipedia.org/wiki/Riemann_sum#Riemann_summation_methods - Some(ComputeMethod::RectangleLeft | ComputeMethod::RectangleRight) => { - let m1: X = (1..n_step) - .map(|step_id| match &self.function { - Some(FunctionDescriptor::Closure(f)) => { - (f(start + step * X::from_usize(step_id).unwrap()) - - f(start + step * X::from_usize(step_id - 1).unwrap())) - / step - } - Some(FunctionDescriptor::Values(v)) => { - (v[step_id] - v[step_id - 1]) / step - } - None => unreachable!(), - }) - .max_by(|t1, t2| t1.abs().partial_cmp(&t2.abs()).unwrap()) - .unwrap(); - let end = start + step * X::from_usize(n_step).unwrap(); - m1 * (end - start).powi(2) / X::from_usize(2 * n_step).unwrap() - } - // ref: https://en.wikipedia.org/wiki/Trapezoidal_rule#Error_analysis - Some(ComputeMethod::Trapezoid) => { - let d1: Vec = (1..n_step) - .map(|step_id| match &self.function { - Some(FunctionDescriptor::Closure(f)) => { - (f(start + step * X::from_usize(step_id).unwrap()) - - f(start + step * X::from_usize(step_id - 1).unwrap())) - / step - } - Some(FunctionDescriptor::Values(v)) => { - (v[step_id] - v[step_id - 1]) / step - } - None => unreachable!(), - }) - .collect(); - let m2: X = (1..n_step - 2) - .map(|step_id| d1[step_id] - d1[step_id - 1] / step) - .max_by(|t1, t2| t1.abs().partial_cmp(&t2.abs()).unwrap()) - .unwrap(); - let end = start + step * X::from_usize(n_step).unwrap(); - -m2 * (end - start).powi(3) / X::from_usize(24 * n_step.pow(2)).unwrap() - } - #[cfg(feature = "montecarlo")] - Some(ComputeMethod::MonteCarlo { .. }) => { - todo!() - } - None => unreachable!(), - }; - Ok(res) - } else { - Err(IntegraalError::InconsistentParameters( - "numerical error computation in not supported for non-uniform domains", - )) - } - } } // --- diff --git a/integraal/src/structure/tests.rs b/integraal/src/structure/tests.rs index d17d7f1..03cf1ee 100644 --- a/integraal/src/structure/tests.rs +++ b/integraal/src/structure/tests.rs @@ -2,9 +2,7 @@ // ------ IMPORTS -use crate::{ - ComputeMethod, DomainDescriptor, FunctionDescriptor, Integraal, IntegraalError, Scalar, -}; +use crate::{ComputeMethod, DomainDescriptor, FunctionDescriptor, Integraal, IntegraalError}; // ------ CONTENT @@ -101,23 +99,14 @@ fn inconsistent_parameters() { // correct usages -fn is_within_tolerance( - mut integraal: Integraal, - expected_result: T, -) -> (bool, String) { - let tolerance = integraal.compute_error().unwrap(); - let computed_result = integraal.compute().unwrap(); - let sum = expected_result.abs() + computed_result.abs(); - let delta = (computed_result - expected_result).abs(); - ( - delta < tolerance + T::epsilon() * sum.min(T::max_value()), - format!("computed value: {computed_result:?}\nexpected value {expected_result:?}\ncomputed tolerance: {tolerance:?}"), - ) -} - +// works for F: Float macro_rules! almost_equal { - ($v1: expr, $v2: expr, $tol: ident) => { - ($v1 - $v2).abs() < $tol + ($ft: ty, $v1: expr, $v2: expr) => { + ($v1 - $v2).abs() < ($v1.abs() + $v2.abs()).min(<$ft as num_traits::Float>::max_value()) + }; + ($ft: ty,$v1: expr, $v2: expr, $tol: ident) => { + ($v1 - $v2).abs() + < ($v1.abs() + $v2.abs()).min(<$ft as num_traits::Float>::max_value()) + $tol }; } @@ -145,25 +134,29 @@ macro_rules! generate_test { let res = integraal.compute(); assert!(res.is_ok()); assert!( - almost_equal!(res.unwrap(), 2.0, $tol), + almost_equal!(f64, res.unwrap(), 2.0_f64, $tol), "left: {} \nright: 2.0", res.unwrap() ); } }; - ($name: ident, $fnd: expr, $dmd: expr, $met: expr) => { + ($name: ident, $fnd: expr, $dmd: expr, $met: expr, $tol: ident) => { #[allow(non_snake_case)] #[test] fn $name() { let functiond = $fnd; let domaind = $dmd; let computem = $met; - let integraal = Integraal::default() + let mut integraal = Integraal::default() .function(functiond) .domain(domaind) .method(computem); - let (res, msg) = is_within_tolerance(integraal, 2.0); - assert!(res, "{msg}"); + let res = integraal.compute().unwrap(); + assert!( + almost_equal!(f64, res, 2.0_f64, $tol), + "computed value: {res:?}\nexpected value: 2.0\ntolerance: {:?}", + $tol + ); } }; } @@ -191,7 +184,8 @@ mod a_rectangleleft { step: STEP, n_step: (1000. * std::f64::consts::PI) as usize, }, - ComputeMethod::RectangleLeft + ComputeMethod::RectangleLeft, + RECTANGLE_TOLERANCE ); generate_test!( @@ -217,7 +211,8 @@ mod a_rectangleleft { step: STEP, n_step: (1000. * std::f64::consts::PI) as usize, }, - ComputeMethod::RectangleLeft + ComputeMethod::RectangleLeft, + RECTANGLE_TOLERANCE ); } @@ -244,7 +239,8 @@ mod a_rectangleright { step: STEP, n_step: (1000. * std::f64::consts::PI) as usize, }, - ComputeMethod::RectangleRight + ComputeMethod::RectangleRight, + RECTANGLE_TOLERANCE ); generate_test!( @@ -270,7 +266,8 @@ mod a_rectangleright { step: STEP, n_step: (1000. * std::f64::consts::PI) as usize, }, - ComputeMethod::RectangleRight + ComputeMethod::RectangleRight, + RECTANGLE_TOLERANCE ); } @@ -297,7 +294,8 @@ mod a_trapezoid { step: STEP, n_step: (1000. * std::f64::consts::PI) as usize, }, - ComputeMethod::Trapezoid + ComputeMethod::Trapezoid, + TRAPEZOID_TOLERANCE ); generate_test!( @@ -323,7 +321,8 @@ mod a_trapezoid { step: STEP, n_step: (1000. * std::f64::consts::PI) as usize, }, - ComputeMethod::Trapezoid + ComputeMethod::Trapezoid, + TRAPEZOID_TOLERANCE ); } @@ -351,7 +350,7 @@ fn B_Closure_Explicit_Rectangle() { .compute(); assert!(res.is_ok()); assert!( - almost_equal!(res.unwrap(), 0.0, RECTANGLE_TOLERANCE), + almost_equal!(f64, res.unwrap(), 0.0_f64, RECTANGLE_TOLERANCE), "left: {} \nright: 0.0", res.unwrap() );