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
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132

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 rule to apply for QAG integration. Rules are supported for 15,
/// 21, 31, 41, 51, 61 points.
#[derive(Debug, Hash, Copy, Clone, PartialEq, Eq)]
pub enum QAGRule {
    Gauss15,
    Gauss21,
    Gauss31,
    Gauss41,
    Gauss51,
    Gauss61
}

impl Into<c_int> for QAGRule {
    fn into(self) -> c_int {
        use self::QAGRule::*;
        match self {
            Gauss15 => bindings::GSL_INTEG_GAUSS15 as c_int,
            Gauss21 => bindings::GSL_INTEG_GAUSS21 as c_int,
            Gauss31 => bindings::GSL_INTEG_GAUSS31 as c_int,
            Gauss41 => bindings::GSL_INTEG_GAUSS41 as c_int,
            Gauss51 => bindings::GSL_INTEG_GAUSS51 as c_int,
            Gauss61 => bindings::GSL_INTEG_GAUSS61 as c_int,
        }
    }
}

/// Quadrature Adaptive General integration. Iteratively applies quadrature
/// rules to sub-regions of a function until it converges to the requested
/// uncertainty. This method will work reasonably well for most well-behaved
/// functions.
///
/// ```
/// use integrators::{Integrator, Real};
/// let mut qag = integrators::QAG::new(1000);
///
/// let res1 = qag.integrate(|a: Real| a * a, 1e-6, 1e-10)
///               .unwrap();
///
/// assert!((res1.value - 3f64.recip()).abs() < res1.error);
///
/// let res2 = qag.with_range(3.0, 10.0)
///               .integrate(|a: Real| a * a, 1e-6, 1e-10)
///               .unwrap();
/// assert!((res2.value - (1000f64 - 27f64) / 3.).abs() < res2.error);
/// ```
#[derive(Debug, Clone)]
pub struct QAG {
    range_low: Real,
    range_high: Real,
    wkspc: GSLIntegrationWorkspace,
    rule: QAGRule
}

impl QAG {
    /// Creates a new QAG with enough memory for `nintervals` subintervals.
    /// This will create a QAG to integrate the range [0, 1], and using the
    /// 61-point (highest-order) Gauss-Kronrod rule. To change the integration
    /// bounds, see `with_range`, and to change the quadrature rule, see
    /// `with_rule`.
    pub fn new(nintervals: usize) -> Self {
        QAG {
            range_low: 0.0,
            range_high: 1.0,
            rule: QAGRule::Gauss61,
            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 {
        QAG {
            wkspc: GSLIntegrationWorkspace::new(nintervals),
            ..self
        }
    }

    /// Use a different integration range. (Default = [0, 1])
    pub fn with_range(self, range_low: Real, range_high: Real) -> Self {
        QAG { range_low, range_high, ..self }
    }


    /// Use a specific quadrature rule. (Default = `QAGRule::Gauss61`)
    pub fn with_rule(self, rule: QAGRule) -> Self {
        QAG { rule, ..self }
    }
}

impl Integrator for QAG {
    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_qag(&mut gslfn.function,
                                          self.range_low, self.range_high,
                                          epsabs, epsrel,
                                          self.wkspc.nintervals,
                                          self.rule.into(),
                                          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
            })
        }
    }
}