From 8a1b4c90ad1da87dd09010ea1021d55df729f06f Mon Sep 17 00:00:00 2001 From: imrn99 <95699343+imrn99@users.noreply.github.com> Date: Thu, 1 Aug 2024 06:34:24 +0200 Subject: [PATCH 1/4] delete `Integraal::compute_error` --- integraal/src/structure/implementations.rs | 82 ---------------------- 1 file changed, 82 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", - )) - } - } } // --- From 301e1a90eb704411a9ef436912f8d5f1246bd233 Mon Sep 17 00:00:00 2001 From: imrn99 <95699343+imrn99@users.noreply.github.com> Date: Thu, 1 Aug 2024 06:42:52 +0200 Subject: [PATCH 2/4] update tests --- integraal/src/structure/tests.rs | 32 +++++++++++++++++++------------- 1 file changed, 19 insertions(+), 13 deletions(-) diff --git a/integraal/src/structure/tests.rs b/integraal/src/structure/tests.rs index d17d7f1..0d8eab4 100644 --- a/integraal/src/structure/tests.rs +++ b/integraal/src/structure/tests.rs @@ -102,16 +102,15 @@ fn inconsistent_parameters() { // correct usages fn is_within_tolerance( - mut integraal: Integraal, + computed_result: T, expected_result: T, + tolerance: 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:?}"), + format!("computed value: {computed_result:?}\nexpected value {expected_result:?}\ntolerance: {tolerance:?}"), ) } @@ -151,18 +150,19 @@ macro_rules! generate_test { ); } }; - ($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); + let res = integraal.compute().unwrap(); + let (res, msg) = is_within_tolerance(res, 2.0, $tol); assert!(res, "{msg}"); } }; @@ -191,7 +191,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 +218,8 @@ mod a_rectangleleft { step: STEP, n_step: (1000. * std::f64::consts::PI) as usize, }, - ComputeMethod::RectangleLeft + ComputeMethod::RectangleLeft, + RECTANGLE_TOLERANCE ); } @@ -244,7 +246,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 +273,8 @@ mod a_rectangleright { step: STEP, n_step: (1000. * std::f64::consts::PI) as usize, }, - ComputeMethod::RectangleRight + ComputeMethod::RectangleRight, + RECTANGLE_TOLERANCE ); } @@ -297,7 +301,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 +328,8 @@ mod a_trapezoid { step: STEP, n_step: (1000. * std::f64::consts::PI) as usize, }, - ComputeMethod::Trapezoid + ComputeMethod::Trapezoid, + TRAPEZOID_TOLERANCE ); } From 84ed73f11bd20cc8ad70178bf363a642160fa587 Mon Sep 17 00:00:00 2001 From: imrn99 <95699343+imrn99@users.noreply.github.com> Date: Thu, 1 Aug 2024 06:52:03 +0200 Subject: [PATCH 3/4] rewrite `almost_equal!` macro --- integraal/src/structure/tests.rs | 13 +++++++++---- 1 file changed, 9 insertions(+), 4 deletions(-) diff --git a/integraal/src/structure/tests.rs b/integraal/src/structure/tests.rs index 0d8eab4..d036245 100644 --- a/integraal/src/structure/tests.rs +++ b/integraal/src/structure/tests.rs @@ -114,9 +114,14 @@ fn is_within_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 }; } @@ -144,7 +149,7 @@ 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() ); @@ -357,7 +362,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() ); From e5532f22429077ffb152d5941fd37af28595f174 Mon Sep 17 00:00:00 2001 From: imrn99 <95699343+imrn99@users.noreply.github.com> Date: Thu, 1 Aug 2024 06:54:43 +0200 Subject: [PATCH 4/4] delete `is_within_tolerance` --- integraal/src/structure/tests.rs | 24 ++++++------------------ 1 file changed, 6 insertions(+), 18 deletions(-) diff --git a/integraal/src/structure/tests.rs b/integraal/src/structure/tests.rs index d036245..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,19 +99,6 @@ fn inconsistent_parameters() { // correct usages -fn is_within_tolerance( - computed_result: T, - expected_result: T, - tolerance: T, -) -> (bool, String) { - 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:?}\ntolerance: {tolerance:?}"), - ) -} - // works for F: Float macro_rules! almost_equal { ($ft: ty, $v1: expr, $v2: expr) => { @@ -167,8 +152,11 @@ macro_rules! generate_test { .domain(domaind) .method(computem); let res = integraal.compute().unwrap(); - let (res, msg) = is_within_tolerance(res, 2.0, $tol); - assert!(res, "{msg}"); + assert!( + almost_equal!(f64, res, 2.0_f64, $tol), + "computed value: {res:?}\nexpected value: 2.0\ntolerance: {:?}", + $tol + ); } }; }