diff --git a/Cargo.toml b/Cargo.toml index 9df854b..6fe7da5 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -26,3 +26,4 @@ integraal-examples = { version = "0.0.2", path = "./examples" } # external rand = "0.9.0-alpha.1" rustversion = "1.0.15" +num = "0.4.3" diff --git a/integraal/Cargo.toml b/integraal/Cargo.toml index fe63f68..164ee97 100644 --- a/integraal/Cargo.toml +++ b/integraal/Cargo.toml @@ -20,6 +20,7 @@ montecarlo = ["dep:rand"] # DEPS [dependencies] +num.workspace = true rand = { workspace = true, features = ["small_rng"], optional = true } [build-dependencies] diff --git a/integraal/src/lib.rs b/integraal/src/lib.rs index a031fd6..279643f 100644 --- a/integraal/src/lib.rs +++ b/integraal/src/lib.rs @@ -22,8 +22,10 @@ mod parameters; mod structure; +mod traits; // --- RE-EXPORTS pub use parameters::{ComputeMethod, DomainDescriptor, FunctionDescriptor}; pub use structure::{Integraal, IntegraalError}; +pub use traits::Scalar; diff --git a/integraal/src/parameters.rs b/integraal/src/parameters.rs index 8353574..94f8a79 100644 --- a/integraal/src/parameters.rs +++ b/integraal/src/parameters.rs @@ -1,5 +1,7 @@ //! integral parameterization code +use crate::Scalar; + /// Domain description enum /// /// This is essentially a discretization of the integrated space. @@ -8,15 +10,15 @@ /// `f64` values (i.e. the type used for further computations). In the future, adding support /// for higher dimension & generic value type can be considered. #[derive(Debug, Clone)] -pub enum DomainDescriptor<'a> { +pub enum DomainDescriptor<'a, T: Scalar> { /// List of values taken by the variable on which we integrate. - Explicit(&'a [f64]), + Explicit(&'a [T]), /// Description of a uniform discretization over a certain range of values. Uniform { /// First value of the range - start: f64, + start: T, /// Step between each value of the range - step: f64, + step: T, /// Total number of values n_step: usize, }, @@ -26,13 +28,16 @@ pub enum DomainDescriptor<'a> { /// /// This enum is used to provide either the values taken by the function or describe of to compute /// those. -pub enum FunctionDescriptor { +pub enum FunctionDescriptor +where + X: Scalar, +{ /// Direct expression of the function, taking a value of the domain as input & returning the /// image of that value through the function. - Closure(Box f64>), + Closure(Box X>), /// List of values taken by the function. The coherence with the domain description must /// be ensured by the user in this case. - Values(Vec), + Values(Vec), } /// Numerical integration method enum diff --git a/integraal/src/structure/definitions.rs b/integraal/src/structure/definitions.rs index b83ab7f..841e0ce 100644 --- a/integraal/src/structure/definitions.rs +++ b/integraal/src/structure/definitions.rs @@ -2,7 +2,7 @@ // ------ IMPORTS -use crate::{ComputeMethod, DomainDescriptor, FunctionDescriptor}; +use crate::{ComputeMethod, DomainDescriptor, FunctionDescriptor, Scalar}; // ------ CONTENT @@ -56,8 +56,8 @@ pub enum IntegraalError { /// # } /// ``` #[derive(Default)] -pub struct Integraal<'a> { - pub(crate) domain: Option>, - pub(crate) function: Option, +pub struct Integraal<'a, X: Scalar> { + pub(crate) domain: Option>, + pub(crate) function: Option>, pub(crate) method: Option, } diff --git a/integraal/src/structure/implementations.rs b/integraal/src/structure/implementations.rs index 05edb36..f3ccc58 100644 --- a/integraal/src/structure/implementations.rs +++ b/integraal/src/structure/implementations.rs @@ -2,19 +2,23 @@ // ------ IMPORTS -use crate::{ComputeMethod, DomainDescriptor, FunctionDescriptor, Integraal, IntegraalError}; +use crate::{ + ComputeMethod, DomainDescriptor, FunctionDescriptor, Integraal, IntegraalError, Scalar, +}; +use num::abs; +use std::ops::Deref; // ------ CONTENT -impl<'a> Integraal<'a> { +impl<'a, X: Scalar> Integraal<'a, X> { /// Set the domain descriptor. - pub fn domain(&mut self, domain_descriptor: DomainDescriptor<'a>) -> &mut Self { + pub fn domain(&mut self, domain_descriptor: DomainDescriptor<'a, X>) -> &mut Self { self.domain = Some(domain_descriptor); self } /// Set the function descriptor. - pub fn function(&mut self, function_descriptor: FunctionDescriptor) -> &mut Self { + pub fn function(&mut self, function_descriptor: FunctionDescriptor) -> &mut Self { self.function = Some(function_descriptor); self } @@ -25,7 +29,11 @@ impl<'a> Integraal<'a> { self } - #[allow(clippy::missing_errors_doc, clippy::too_many_lines)] + #[allow( + clippy::missing_errors_doc, + clippy::missing_panics_doc, + clippy::too_many_lines + )] /// This method attempts to compute the integral. If it is successful, it will clear the /// internal [`FunctionDescriptor`] object before returning the result. /// @@ -34,7 +42,7 @@ impl<'a> Integraal<'a> { /// This method returns a `Result` taking the following values: /// - `Ok(f64)` -- The computation was successfuly done /// - `Err(IntegraalError)` -- The computation failed for the reason specified by the enum. - pub fn compute(&mut self) -> Result { + pub fn compute(&mut 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", @@ -54,13 +62,13 @@ impl<'a> Integraal<'a> { Some(ComputeMethod::RectangleLeft) => (1..n_sample) .map(|idx| { let step = args[idx] - args[idx - 1]; - step * vals[idx - 1] + vals[idx - 1] * step }) .sum(), Some(ComputeMethod::RectangleRight) => (1..n_sample) .map(|idx| { let step = args[idx] - args[idx - 1]; - step * vals[idx] + vals[idx] * step }) .sum(), Some(ComputeMethod::Trapezoid) => (1..n_sample) @@ -68,7 +76,7 @@ impl<'a> Integraal<'a> { let step = args[idx] - args[idx - 1]; let y1 = vals[idx - 1]; let y2 = vals[idx]; - step * (y1.min(y2) + (y1 - y2).abs() / 2.0) + (y1.min(y2) + abs(y1 - y2) / X::from_f32(2.0).unwrap()) * step }) .sum(), #[cfg(feature = "montecarlo")] @@ -96,17 +104,17 @@ impl<'a> Integraal<'a> { match &self.method { Some(ComputeMethod::RectangleLeft) => { // ignore the last value since its a left rule - (0..*n_step - 1).map(|step_id| vals[step_id] * step).sum() + (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() + (1..*n_step).map(|step_id| vals[step_id] * *step).sum() } Some(ComputeMethod::Trapezoid) => (1..*n_step) .map(|step_id| { let y1 = vals[step_id - 1]; let y2 = vals[step_id]; - step * (y1.min(y2) + (y1 - y2).abs() / 2.0) + (y1.min(y2) + (y1 - y2).abs() / X::from_f32(2.0).unwrap()) * *step }) .sum(), #[cfg(feature = "montecarlo")] @@ -123,21 +131,21 @@ impl<'a> Integraal<'a> { Some(ComputeMethod::RectangleLeft) => (1..args.len()) .map(|idx| { let step = args[idx] - args[idx - 1]; - step * closure(args[idx - 1]) + closure(args[idx - 1]) * step }) .sum(), Some(ComputeMethod::RectangleRight) => (1..args.len()) .map(|idx| { let step = args[idx] - args[idx - 1]; - step * closure(args[idx]) + closure(args[idx]) * step }) .sum(), Some(ComputeMethod::Trapezoid) => (1..args.len()) .map(|idx| { let step = args[idx] - args[idx - 1]; - let y1 = closure(args[idx - 1]); + let y1 = closure.deref()(args[idx - 1]); let y2 = closure(args[idx]); - step * (y1.min(y2) + (y1 - y2).abs() / 2.0) + (y1.min(y2) + (y1 - y2).abs() / X::from_f32(2.0).unwrap()) * step }) .sum(), #[cfg(feature = "montecarlo")] @@ -158,23 +166,23 @@ impl<'a> Integraal<'a> { match &self.method { Some(ComputeMethod::RectangleLeft) => (0..*n_step - 1) .map(|step_id| { - let x = start + step * step_id as f64; - closure(x) * step + let x = *start + *step * X::from_usize(step_id).unwrap(); + closure(x) * *step }) .sum(), Some(ComputeMethod::RectangleRight) => (1..*n_step) .map(|step_id| { - let x = start + step * step_id as f64; - closure(x) * step + let x = *start + *step * X::from_usize(step_id).unwrap(); + closure(x) * *step }) .sum(), Some(ComputeMethod::Trapezoid) => (1..*n_step) .map(|step_id| { - let x1 = start + step * (step_id - 1) as f64; - let x2 = start + step * step_id as f64; - let y1 = closure(x1); + let x1 = *start + *step * X::from_usize(step_id - 1).unwrap(); + let x2 = *start + *step * X::from_usize(step_id).unwrap(); + let y1 = closure.deref()(x1); let y2 = closure(x2); - step * (y1.min(y2) + (y1 - y2).abs() / 2.0) + (y1.min(y2) + (y1 - y2).abs() / X::from_f32(2.0).unwrap()) * *step }) .sum(), #[cfg(feature = "montecarlo")] diff --git a/integraal/src/structure/tests.rs b/integraal/src/structure/tests.rs index fe6b2cf..3a18e59 100644 --- a/integraal/src/structure/tests.rs +++ b/integraal/src/structure/tests.rs @@ -20,15 +20,15 @@ macro_rules! almost_equal { 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::RectangleLeft; + let $f: FunctionDescriptor = FunctionDescriptor::Closure(Box::new(|x| x)); + let $d: DomainDescriptor<'_, f64> = DomainDescriptor::Explicit(&[]); + let $c: ComputeMethod = ComputeMethod::RectangleLeft; }; } macro_rules! generate_missing { ($a: ident, $b: ident) => { - let mut integral = Integraal::default(); + let mut integral: Integraal<'_, f64> = Integraal::default(); integral.$a($a).$b($b); assert_eq!( integral.compute(), @@ -57,7 +57,7 @@ fn missing_parameters() { generate_missing!(function, domain); // missing all but one - let mut integral = Integraal::default(); + let mut integral: Integraal<'_, f64> = Integraal::default(); integral.method(method); assert_eq!( integral.compute(), diff --git a/integraal/src/traits.rs b/integraal/src/traits.rs new file mode 100644 index 0000000..0d0ecfd --- /dev/null +++ b/integraal/src/traits.rs @@ -0,0 +1,29 @@ +//! module doc + +/// Scalar value trait. +/// +/// This trait is automatically implemented for all types implementing its requirements. +pub trait Scalar: + Clone + + Copy + + num::Float + + num::Signed + + num::FromPrimitive + + std::ops::Sub + + std::ops::Mul + + std::iter::Sum +{ +} + +impl< + X: Clone + + Copy + + num::Float + + num::Signed + + num::FromPrimitive + + std::ops::Sub + + std::ops::Mul + + std::iter::Sum, + > Scalar for X +{ +}