quad_rs/
solve.rs

1use nalgebra::{ComplexField, RealField};
2use num_traits::FromPrimitive;
3
4use std::ops::Range;
5use trellis_runner::{Calculation, Problem};
6
7use crate::{
8    contour::split_range_around_singularities, AccumulateError, Contour, GaussKronrod,
9    GaussKronrodCore, Integrable, IntegrableFloat, IntegrationError, IntegrationOutput,
10    IntegrationResult, IntegrationState, RescaleError,
11};
12
13#[derive(Debug)]
14pub struct AdaptiveIntegrator<F: ComplexField> {
15    limits: Range<F>,
16    max_elements: usize,
17    minimum_element_size: <F as ComplexField>::RealField,
18    singularities: Vec<F>,
19    integrator: GaussKronrod<<F as ComplexField>::RealField>,
20    relative_tolerance: <F as ComplexField>::RealField,
21    absolute_tolerance: <F as ComplexField>::RealField,
22}
23
24impl<F: ComplexField> AdaptiveIntegrator<F> {
25    pub fn new_real(
26        Range { start, end }: Range<<F as ComplexField>::RealField>,
27        max_elements: usize,
28        minimum_element_size: <F as ComplexField>::RealField,
29        singularities: Vec<F>,
30        relative_tolerance: <F as ComplexField>::RealField,
31        absolute_tolerance: <F as ComplexField>::RealField,
32    ) -> Self {
33        Self {
34            limits: Range {
35                start: F::from_real(start),
36                end: F::from_real(end),
37            },
38            max_elements,
39            minimum_element_size,
40            singularities,
41            integrator: GaussKronrod::default(),
42            relative_tolerance,
43            absolute_tolerance,
44        }
45    }
46
47    pub fn new(
48        limits: Range<F>,
49        max_elements: usize,
50        minimum_element_size: <F as ComplexField>::RealField,
51        singularities: Vec<F>,
52        relative_tolerance: <F as ComplexField>::RealField,
53        absolute_tolerance: <F as ComplexField>::RealField,
54    ) -> Self {
55        Self {
56            limits,
57            max_elements,
58            minimum_element_size,
59            singularities,
60            integrator: GaussKronrod::default(),
61            relative_tolerance,
62            absolute_tolerance,
63        }
64    }
65}
66
67impl<P> Integrable for Problem<P>
68where
69    P: Integrable,
70{
71    type Input = P::Input;
72    type Output = P::Output;
73
74    fn integrand(
75        &self,
76        input: &Self::Input,
77    ) -> Result<Self::Output, crate::EvaluationError<Self::Input>> {
78        self.as_ref().integrand(input)
79    }
80}
81
82impl<I, O, F, P> Calculation<P, IntegrationState<I, O, F>> for AdaptiveIntegrator<I>
83where
84    P: Integrable<Input = I, Output = O> + Send + Sync,
85    I: ComplexField<RealField = F> + FromPrimitive + Copy,
86    O: IntegrationOutput<Float = F, Scalar = I>,
87    F: IntegrableFloat + RealField + FromPrimitive + AccumulateError<F> + RescaleError,
88{
89    const NAME: &'static str = "Adaptive Integrator";
90    type Error = IntegrationError<I>;
91    type Output = IntegrationResult<I, O>;
92
93    fn initialise(
94        &mut self,
95        problem: &mut Problem<P>,
96        state: IntegrationState<I, O, F>,
97    ) -> Result<IntegrationState<I, O, F>, Self::Error> {
98        let initial_segments = match self.singularities.len() {
99            0 => self
100                .integrator
101                // .gauss_kronrod(|x| problem.integrand(&x).unwrap(), self.limits.clone())
102                .gauss_kronrod(&*problem, self.limits.clone())?,
103            _ => split_range_around_singularities(self.limits.clone(), self.singularities.clone())
104                .into_iter()
105                .map(|range| {
106                    self.integrator
107                        // .gauss_kronrod(|x| problem.integrand(&x).unwrap(), range)
108                        .gauss_kronrod(&*problem, range)
109                })
110                .collect::<Result<Vec<_>, _>>()?
111                .into_iter()
112                .flatten()
113                .collect(),
114        };
115
116        Ok(state.segments(initial_segments))
117    }
118
119    fn next(
120        &mut self,
121        problem: &mut Problem<P>,
122        mut state: IntegrationState<I, O, F>,
123    ) -> Result<IntegrationState<I, O, F>, Self::Error> {
124        let worst_segment = state.pop_worst_segment().unwrap();
125        let new_segments = self.integrator.split_segment(&*problem, worst_segment)?;
126        Ok(state.segments(new_segments))
127    }
128
129    fn finalise(
130        &mut self,
131        _problem: &mut Problem<P>,
132        state: IntegrationState<I, O, F>,
133    ) -> Result<Self::Output, Self::Error> {
134        Ok(state.into())
135    }
136}
137
138#[derive(Debug)]
139pub struct AdaptiveRectangularContourIntegrator<F: ComplexField> {
140    pub(crate) contour: Contour<F>,
141    pub(crate) max_elements: usize,
142    pub(crate) minimum_element_size: <F as ComplexField>::RealField,
143    pub(crate) integrator: GaussKronrod<<F as ComplexField>::RealField>,
144    pub(crate) relative_tolerance: <F as ComplexField>::RealField,
145    pub(crate) absolute_tolerance: <F as ComplexField>::RealField,
146}
147
148impl<F> AdaptiveRectangularContourIntegrator<F>
149where
150    F: ComplexField + Copy,
151    <F as ComplexField>::RealField: Copy,
152{
153    pub fn new(
154        contour: Contour<F>,
155        max_elements: usize,
156        minimum_element_size: <F as ComplexField>::RealField,
157        relative_tolerance: <F as ComplexField>::RealField,
158        absolute_tolerance: <F as ComplexField>::RealField,
159    ) -> Self {
160        Self {
161            contour,
162            max_elements,
163            minimum_element_size,
164            integrator: GaussKronrod::default(),
165            relative_tolerance,
166            absolute_tolerance,
167        }
168    }
169    pub fn initialise<I: Into<Contour<F>>>(
170        input: I,
171        max_elements: usize,
172        minimum_element_size: <F as ComplexField>::RealField,
173        relative_tolerance: <F as ComplexField>::RealField,
174        absolute_tolerance: <F as ComplexField>::RealField,
175    ) -> Self {
176        Self {
177            contour: input.into(),
178            max_elements,
179            minimum_element_size,
180            integrator: GaussKronrod::default(),
181            relative_tolerance,
182            absolute_tolerance,
183        }
184    }
185}
186
187impl<I, O, F, P> Calculation<P, IntegrationState<I, O, F>>
188    for AdaptiveRectangularContourIntegrator<I>
189where
190    P: Integrable<Input = I, Output = O> + Send + Sync,
191    I: ComplexField<RealField = F> + FromPrimitive + Copy,
192    O: IntegrationOutput<Float = F, Scalar = I>,
193    F: IntegrableFloat + RealField + FromPrimitive + AccumulateError<F> + RescaleError,
194{
195    const NAME: &'static str = "Adaptive Contour Integrator";
196    type Error = IntegrationError<I>;
197    type Output = IntegrationResult<I, O>;
198
199    fn initialise(
200        &mut self,
201        problem: &mut Problem<P>,
202        state: IntegrationState<I, O, F>,
203    ) -> Result<IntegrationState<I, O, F>, Self::Error> {
204        let initial_segments = self
205            .contour
206            .range
207            .iter()
208            .map(|range| {
209                self.integrator
210                    // .gauss_kronrod(|x| problem.integrand(&x).unwrap(), range.clone())
211                    .gauss_kronrod(&*problem, range.clone())
212            })
213            .collect::<Result<Vec<_>, _>>()?
214            .into_iter()
215            .flatten()
216            .collect();
217        Ok(state.segments(initial_segments))
218    }
219
220    fn next(
221        &mut self,
222        problem: &mut Problem<P>,
223        mut state: IntegrationState<I, O, F>,
224    ) -> Result<IntegrationState<I, O, F>, Self::Error> {
225        let worst_segment = state.pop_worst_segment().unwrap();
226        let new_segments = self
227            .integrator
228            // .split_segment(|x| problem.integrand(&x).unwrap(), worst_segment)
229            .split_segment(&*problem, worst_segment)?;
230        Ok(state.segments(new_segments))
231    }
232
233    fn finalise(
234        &mut self,
235        _problem: &mut Problem<P>,
236        state: IntegrationState<I, O, F>,
237    ) -> Result<Self::Output, Self::Error> {
238        Ok(state.into())
239    }
240}
quad_rs/solve.rs

quad_rs/
solve.rs
