Skip to content

Commit

Permalink
refactor: remove Integraal::compute_error (#28)
Browse files Browse the repository at this point in the history 
* delete `Integraal::compute_error`

* update tests

* rewrite `almost_equal!` macro

* delete `is_within_tolerance`
  • Loading branch information
@imrn99
imrn99 authored Aug 1, 2024
1 parent bdf219b commit 940f180
Show file tree
Hide file tree
Showing 2 changed files with 30 additions and 113 deletions.
82 changes: 0 additions & 82 deletions integraal/src/structure/implementations.rs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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<X, IntegraalError> {
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<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!(),
})
.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",
))
}
}
}

// ---
Expand Down
61 changes: 30 additions & 31 deletions integraal/src/structure/tests.rs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,9 +2,7 @@
// ------ IMPORTS

use crate::{
ComputeMethod, DomainDescriptor, FunctionDescriptor, Integraal, IntegraalError, Scalar,
};
use crate::{ComputeMethod, DomainDescriptor, FunctionDescriptor, Integraal, IntegraalError};

// ------ CONTENT

Expand Down Expand Up @@ -101,23 +99,14 @@ fn inconsistent_parameters() {

// correct usages

fn is_within_tolerance<T: Scalar>(
mut integraal: Integraal<T>,
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
};
}

Expand Down Expand Up @@ -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
);
}
};
}
Expand Down Expand Up @@ -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!(
Expand All @@ -217,7 +211,8 @@ mod a_rectangleleft {
step: STEP,
n_step: (1000. * std::f64::consts::PI) as usize,
},
ComputeMethod::RectangleLeft
ComputeMethod::RectangleLeft,
RECTANGLE_TOLERANCE
);
}

Expand All @@ -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!(
Expand All @@ -270,7 +266,8 @@ mod a_rectangleright {
step: STEP,
n_step: (1000. * std::f64::consts::PI) as usize,
},
ComputeMethod::RectangleRight
ComputeMethod::RectangleRight,
RECTANGLE_TOLERANCE
);
}

Expand All @@ -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!(
Expand All @@ -323,7 +321,8 @@ mod a_trapezoid {
step: STEP,
n_step: (1000. * std::f64::consts::PI) as usize,
},
ComputeMethod::Trapezoid
ComputeMethod::Trapezoid,
TRAPEZOID_TOLERANCE
);
}

Expand Down Expand Up @@ -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()
);
Expand Down

0 comments on commit 940f180

Please sign in to comment.