integrators 0.0.3

Generic traits to allow easy plugging of different numerical integration algorithms, and wrappers for the Cuba and GSL integrators.
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64

use ::bindings;
use ::ffi::LandingPad;
use ::traits::{IntegrandInput, IntegrandOutput};
use ::{IntegrationResult, Integrator, Real};
use super::{make_gsl_function, GSLIntegrationError};

/// Quadrature Non-adaptive General-use integrator. Iteratively
/// applies Gauss-Kronrod quadrature rules of successively higher
/// order until a desired precision is reached. If such precision
/// is not reached by the 87-point quadrature rule, the integration
/// fails.
///
/// See GSL Docs [here](https://www.gnu.org/software/gsl/doc/html/integration.html#c.gsl_integration_qng).
#[derive(Debug, Copy, Clone)]
pub struct QNG {
    range_low: Real,
    range_high: Real
}

impl QNG {
    /// Creates a new QNG integrator which will integrate a 1-dimensional
    /// function from `range_low` to `range_high`.
    pub fn new(range_low: Real, range_high: Real) -> Self {
        QNG { range_low, range_high }
    }

    /// Builder pattern to change the integration range of the QNG.
    pub fn with_range(self, range_low: Real, range_high: Real) -> Self {
        QNG { range_low, range_high }
    }
}

impl Integrator for QNG {
    type Success = IntegrationResult;
    type Failure = GSLIntegrationError;
    fn integrate<A, B, F: FnMut(A) -> B>(&mut self, fun: F, epsrel: Real, epsabs: Real) -> Result<Self::Success, Self::Failure>
        where A: IntegrandInput,
              B: IntegrandOutput
    {
        let mut value: Real = 0.0;
        let mut error: Real = 0.0;
        let mut neval: usize = 0;

        let mut lp = LandingPad::new(fun);
        let retcode = unsafe {
            let mut gslfn = make_gsl_function(&mut lp, self.range_low, self.range_high)?;
            bindings::gsl_integration_qng(&mut gslfn.function,
                                          self.range_low, self.range_high,
                                          epsabs, epsrel,
                                          &mut value,
                                          &mut error,
                                          &mut neval)
        };
        lp.maybe_resume_unwind();

        if retcode != bindings::GSL_SUCCESS {
            Err(GSLIntegrationError::GSLError(retcode.into()))
        } else {
            Ok(IntegrationResult {
                value, error
            })
        }
    }
}