diff --git a/docs/build.sh b/docs/build.sh
index 63d5e20..3f59413 100644
--- a/docs/build.sh
+++ b/docs/build.sh
@@ -1,2 +1,3 @@
cp docs/index.html target/doc
-cp docs/index.css target/doc
\ No newline at end of file
+cp docs/index.css target/doc
+cp docs/compute_methods.svg target/doc
diff --git a/docs/compute_methods.svg b/docs/compute_methods.svg
new file mode 100644
index 0000000..1436e97
--- /dev/null
+++ b/docs/compute_methods.svg
@@ -0,0 +1,423 @@
+
+
+
+
diff --git a/examples/examples/rectangle/integraal.rs b/examples/examples/rectangle/integraal.rs
index 8a4bcdf..fa22bee 100644
--- a/examples/examples/rectangle/integraal.rs
+++ b/examples/examples/rectangle/integraal.rs
@@ -8,7 +8,7 @@ fn main() {
n_step: 100_001,
};
let function = FunctionDescriptor::Closure(Box::new(|x: f64| 2.0 * x));
- let method = ComputeMethod::Rectangle;
+ let method = ComputeMethod::RectangleLeft;
// build the integral
let mut integral = Integraal::default();
diff --git a/integraal/src/parameters.rs b/integraal/src/parameters.rs
index f9eddbe..8353574 100644
--- a/integraal/src/parameters.rs
+++ b/integraal/src/parameters.rs
@@ -36,10 +36,30 @@ pub enum FunctionDescriptor {
}
/// Numerical integration method enum
+///
+/// # Note on computations
+///
+/// For the considered integral to be consistent across compute methods for a given description,
+/// the left-rectangle (resp. right-rectangle) has to ignore the last (resp. first) value given
+/// in the descriptors. This can be visualized in the following example:
+///
+/// 
+///
+/// Out of 11 samples, both methods compute the area of 10 polygons. In the case where the domain
+/// is uniform & described using a step, the eleventh sample value is useless (for a left-rectangle
+/// method).
+///
+/// The crate assumes that the first and last samples making up your domain corresponds to the
+/// limits of the integral. Therefore, these values will be ignored when computing the integral
+/// using rectangles.
#[derive(Debug, Clone, Copy)]
pub enum ComputeMethod {
- /// Rectangle method -- [reference](https://en.wikipedia.org/wiki/Riemann_sum)
- Rectangle,
+ /// Rectangle method, using the left rule --
+ /// [reference](https://en.wikipedia.org/wiki/Riemann_sum#Left_rule)
+ RectangleLeft,
+ /// Rectangle method, using the right rule --
+ /// [reference](https://en.wikipedia.org/wiki/Riemann_sum#Right_rule)
+ RectangleRight,
/// Trapezoid method [reference](https://en.wikipedia.org/wiki/Trapezoidal_rule)
Trapezoid,
#[cfg(feature = "montecarlo")]
diff --git a/integraal/src/structure/implementations.rs b/integraal/src/structure/implementations.rs
index c923b33..05edb36 100644
--- a/integraal/src/structure/implementations.rs
+++ b/integraal/src/structure/implementations.rs
@@ -51,12 +51,18 @@ impl<'a> Integraal<'a> {
// because the domain may be not uniform, we have to compute step values
match &self.method {
- Some(ComputeMethod::Rectangle) => (1..n_sample)
+ Some(ComputeMethod::RectangleLeft) => (1..n_sample)
.map(|idx| {
let step = args[idx] - args[idx - 1];
step * vals[idx - 1]
})
.sum(),
+ Some(ComputeMethod::RectangleRight) => (1..n_sample)
+ .map(|idx| {
+ let step = args[idx] - args[idx - 1];
+ step * vals[idx]
+ })
+ .sum(),
Some(ComputeMethod::Trapezoid) => (1..n_sample)
.map(|idx| {
let step = args[idx] - args[idx - 1];
@@ -88,8 +94,13 @@ impl<'a> Integraal<'a> {
// we can use the uniform domain's step & number of step to compute areas
match &self.method {
- Some(ComputeMethod::Rectangle) => {
- (0..*n_step).map(|step_id| vals[step_id] * step).sum()
+ Some(ComputeMethod::RectangleLeft) => {
+ // ignore the last value since its a left rule
+ (0..*n_step - 1).map(|step_id| vals[step_id] * step).sum()
+ }
+ Some(ComputeMethod::RectangleRight) => {
+ // ignore the last value since its a left rule
+ (1..*n_step).map(|step_id| vals[step_id] * step).sum()
}
Some(ComputeMethod::Trapezoid) => (1..*n_step)
.map(|step_id| {
@@ -109,12 +120,18 @@ impl<'a> Integraal<'a> {
Some(FunctionDescriptor::Closure(closure)),
Some(DomainDescriptor::Explicit(args)),
) => match &self.method {
- Some(ComputeMethod::Rectangle) => (1..args.len())
+ Some(ComputeMethod::RectangleLeft) => (1..args.len())
.map(|idx| {
let step = args[idx] - args[idx - 1];
step * closure(args[idx - 1])
})
.sum(),
+ Some(ComputeMethod::RectangleRight) => (1..args.len())
+ .map(|idx| {
+ let step = args[idx] - args[idx - 1];
+ step * closure(args[idx])
+ })
+ .sum(),
Some(ComputeMethod::Trapezoid) => (1..args.len())
.map(|idx| {
let step = args[idx] - args[idx - 1];
@@ -139,7 +156,13 @@ impl<'a> Integraal<'a> {
) => {
// compute args
match &self.method {
- Some(ComputeMethod::Rectangle) => (0..*n_step)
+ Some(ComputeMethod::RectangleLeft) => (0..*n_step - 1)
+ .map(|step_id| {
+ let x = start + step * step_id as f64;
+ closure(x) * step
+ })
+ .sum(),
+ Some(ComputeMethod::RectangleRight) => (1..*n_step)
.map(|step_id| {
let x = start + step * step_id as f64;
closure(x) * step
diff --git a/integraal/src/structure/tests.rs b/integraal/src/structure/tests.rs
index bb41ee6..fe6b2cf 100644
--- a/integraal/src/structure/tests.rs
+++ b/integraal/src/structure/tests.rs
@@ -22,7 +22,7 @@ macro_rules! generate_sample_descriptors {
($f: ident, $d: ident, $c: ident) => {
let $f = FunctionDescriptor::Closure(Box::new(|x| x));
let $d = DomainDescriptor::Explicit(&[]);
- let $c = ComputeMethod::Rectangle;
+ let $c = ComputeMethod::RectangleLeft;
};
}
@@ -69,7 +69,7 @@ fn missing_parameters() {
#[test]
fn inconsistent_parameters() {
- let method = ComputeMethod::Rectangle;
+ let method = ComputeMethod::RectangleLeft;
let function = FunctionDescriptor::Values(vec![1., 1., 1., 1., 1., 1.]);
let domain = vec![0.0, 0.1, 0.2, 0.3, 0.4]; // missing the last x value
let domain = DomainDescriptor::Explicit(&domain);
@@ -158,7 +158,7 @@ macro_rules! generate_test {
}
#[allow(clippy::cast_sign_loss, clippy::cast_possible_truncation)]
-mod a_rectangle {
+mod a_rectangleleft {
use super::*;
generate_test!(
@@ -168,7 +168,7 @@ mod a_rectangle {
.collect(),
FunctionDescriptor::Closure(Box::new(f64::sin)),
DomainDescriptor::Explicit(&domain),
- ComputeMethod::Rectangle,
+ ComputeMethod::RectangleLeft,
RECTANGLE_TOLERANCE
);
@@ -180,7 +180,7 @@ mod a_rectangle {
step: STEP,
n_step: (1000. * std::f64::consts::PI) as usize,
},
- ComputeMethod::Rectangle,
+ ComputeMethod::RectangleLeft,
RECTANGLE_TOLERANCE
);
@@ -191,7 +191,7 @@ mod a_rectangle {
.collect(),
FunctionDescriptor::Values(domain.iter().copied().map(f64::sin).collect()),
DomainDescriptor::Explicit(&domain),
- ComputeMethod::Rectangle,
+ ComputeMethod::RectangleLeft,
RECTANGLE_TOLERANCE
);
@@ -207,7 +207,62 @@ mod a_rectangle {
step: STEP,
n_step: (1000. * std::f64::consts::PI) as usize,
},
- ComputeMethod::Rectangle,
+ ComputeMethod::RectangleLeft,
+ RECTANGLE_TOLERANCE
+ );
+}
+
+#[allow(clippy::cast_sign_loss, clippy::cast_possible_truncation)]
+mod a_rectangleright {
+ use super::*;
+
+ generate_test!(
+ ClosureExplicit,
+ let domain: Vec = (0..(std::f64::consts::PI * 1000.) as usize)
+ .map(|step_id| step_id as f64 * STEP)
+ .collect(),
+ FunctionDescriptor::Closure(Box::new(f64::sin)),
+ DomainDescriptor::Explicit(&domain),
+ ComputeMethod::RectangleRight,
+ RECTANGLE_TOLERANCE
+ );
+
+ generate_test!(
+ ClosureUniform,
+ FunctionDescriptor::Closure(Box::new(f64::sin)),
+ DomainDescriptor::Uniform {
+ start: 0.,
+ step: STEP,
+ n_step: (1000. * std::f64::consts::PI) as usize,
+ },
+ ComputeMethod::RectangleRight,
+ RECTANGLE_TOLERANCE
+ );
+
+ generate_test!(
+ ValuesExplicit,
+ let domain: Vec = (0..(std::f64::consts::PI * 1000.) as usize)
+ .map(|step_id| step_id as f64 * STEP)
+ .collect(),
+ FunctionDescriptor::Values(domain.iter().copied().map(f64::sin).collect()),
+ DomainDescriptor::Explicit(&domain),
+ ComputeMethod::RectangleRight,
+ RECTANGLE_TOLERANCE
+ );
+
+ generate_test!(
+ ValuesUniform,
+ FunctionDescriptor::Values(
+ (0..(1000. * std::f64::consts::PI) as usize)
+ .map(|step_id| (step_id as f64 * STEP).sin())
+ .collect()
+ ),
+ DomainDescriptor::Uniform {
+ start: 0.,
+ step: STEP,
+ n_step: (1000. * std::f64::consts::PI) as usize,
+ },
+ ComputeMethod::RectangleRight,
RECTANGLE_TOLERANCE
);
}
@@ -282,7 +337,7 @@ fn B_Closure_Explicit_Rectangle() {
.map(|step_id| -1. + f64::from(step_id) * STEP)
.collect();
let domaind = DomainDescriptor::Explicit(&domain);
- let computem = ComputeMethod::Rectangle;
+ let computem = ComputeMethod::RectangleLeft;
let mut integraal = Integraal::default();
let res = integraal
.function(functiond)