Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

refactor: remove Integraal::compute_error #28

Merged
merged 4 commits into from
Aug 1, 2024
Merged

refactor: remove Integraal::compute_error #28

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
8a1b4c9
delete `Integraal::compute_error`
imrn99 Aug 1, 2024
301e1a9
update tests
imrn99 Aug 1, 2024
84ed73f
rewrite `almost_equal!` macro
imrn99 Aug 1, 2024
e5532f2
delete `is_within_tolerance`
imrn99 Aug 1, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
82 changes: 0 additions & 82 deletions integraal/src/structure/implementations.rs
View file
Open in desktop
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
Open in desktop
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