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
65
66
67
68
69
70
71

use std::convert::Into;
use std::os::raw::c_int;

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

use super::{make_gsl_function, GSLIntegrationError, GSLIntegrationWorkspace};

/// Quadrature Adaptive General integration for Infinite intervals. It applies
/// the QAGS algorithm to a transformation of the input integral, such that
/// integrating the transformed function from 0 to 1 yields the infinite
/// integral of the original function.
///
/// See GSL docs
/// [here](https://www.gnu.org/software/gsl/doc/html/integration.html#c.gsl_integration_qagi).
#[derive(Debug, Clone)]
pub struct QAGI {
    wkspc: GSLIntegrationWorkspace,
}

impl QAGI {
    /// Creates a new QAGI with enough memory for `nintervals` subintervals.
    pub fn new(nintervals: usize) -> Self {
        QAGI {
            wkspc: GSLIntegrationWorkspace::new(nintervals)
        }
    }

    /// Discards the old workspace and allocates a new one with enough memory
    /// for `nintervals` subintervals.
    pub fn with_nintervals(self, nintervals: usize) -> Self {
        QAGI {
            wkspc: GSLIntegrationWorkspace::new(nintervals),
            ..self
        }
    }
}

impl Integrator for QAGI {
    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 lp = LandingPad::new(fun);
        let retcode = unsafe {
            let mut gslfn = make_gsl_function(&mut lp, self.range_low, self.range_high)?;
            bindings::gsl_integration_qagi(&mut gslfn.function,
                                           epsabs, epsrel,
                                           self.wkspc.nintervals,
                                           self.wkspc.wkspc,
                                           &mut value,
                                           &mut error)
        };
        lp.maybe_resume_unwind();

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