tensor4all_tensorci/
integration.rs

1//! Numerical integration using TCI and Gauss-Kronrod quadrature.
2//!
3//! [`integrate`] approximates a multi-dimensional integral over a hypercube
4//! by building a tensor-train representation of the integrand on
5//! Gauss-Kronrod nodes and then summing the tensor train.
6//!
7//! Supported quadrature orders: 15 and 31 (standard Gauss-Kronrod rules).
8
9use crate::error::{Result, TCIError};
10use crate::tensorci2::{crossinterpolate2, TCI2Options};
11use tensor4all_simplett::{AbstractTensorTrain, TTScalar};
12use tensor4all_tcicore::{
13    DenseFaerLuKernel, LazyBlockRookKernel, MatrixLuciScalar, MultiIndex, PivotKernel, Scalar,
14};
15
16/// Gauss-Kronrod 15-point rule: nodes on [-1, 1]
17const GK15_NODES: [f64; 15] = [
18    -0.991_455_371_120_812_6,
19    -0.949_107_912_342_758_5,
20    -0.864_864_423_359_769_1,
21    -0.741_531_185_599_394_5,
22    -0.586_087_235_467_691_1,
23    -0.405_845_151_377_397_2,
24    -0.207_784_955_007_898_48,
25    0.000000000000000000000000000000000,
26    0.207_784_955_007_898_48,
27    0.405_845_151_377_397_2,
28    0.586_087_235_467_691_1,
29    0.741_531_185_599_394_5,
30    0.864_864_423_359_769_1,
31    0.949_107_912_342_758_5,
32    0.991_455_371_120_812_6,
33];
34
35/// Gauss-Kronrod 15-point rule: weights on [-1, 1]
36const GK15_WEIGHTS: [f64; 15] = [
37    0.022_935_322_010_529_224,
38    0.063_092_092_629_978_56,
39    0.104_790_010_322_250_19,
40    0.140_653_259_715_525_92,
41    0.169_004_726_639_267_9,
42    0.190_350_578_064_785_42,
43    0.204_432_940_075_298_89,
44    0.209_482_141_084_727_82,
45    0.204_432_940_075_298_89,
46    0.190_350_578_064_785_42,
47    0.169_004_726_639_267_9,
48    0.140_653_259_715_525_92,
49    0.104_790_010_322_250_19,
50    0.063_092_092_629_978_56,
51    0.022_935_322_010_529_224,
52];
53
54/// Gauss-Kronrod 31-point rule: nodes on [-1, 1]
55/// Source: QUADPACK (Piessens et al., 1983)
56const GK31_NODES: [f64; 31] = [
57    -0.998_002_298_693_397_1,
58    -0.987_992_518_020_485_4,
59    -0.967_739_075_679_139_1,
60    -0.937_273_392_400_706,
61    -0.897_264_532_344_081_9,
62    -0.848_206_583_410_427_2,
63    -0.790_418_501_442_466,
64    -0.724_417_731_360_170_1,
65    -0.650_996_741_297_416_9,
66    -0.570_972_172_608_539,
67    -0.485_081_863_640_239_7,
68    -0.394_151_347_077_563_4,
69    -0.299_180_007_153_168_8,
70    -0.201_194_093_997_434_5,
71    -0.101_142_066_918_717_5,
72    0.0,
73    0.101_142_066_918_717_5,
74    0.201_194_093_997_434_5,
75    0.299_180_007_153_168_8,
76    0.394_151_347_077_563_4,
77    0.485_081_863_640_239_7,
78    0.570_972_172_608_539,
79    0.650_996_741_297_416_9,
80    0.724_417_731_360_170_1,
81    0.790_418_501_442_466,
82    0.848_206_583_410_427_2,
83    0.897_264_532_344_081_9,
84    0.937_273_392_400_706,
85    0.967_739_075_679_139_1,
86    0.987_992_518_020_485_4,
87    0.998_002_298_693_397_1,
88];
89
90/// Gauss-Kronrod 31-point rule: weights on [-1, 1]
91const GK31_WEIGHTS: [f64; 31] = [
92    0.005_377_479_872_923_349,
93    0.015_007_947_329_316_122,
94    0.025_460_847_326_715_32,
95    0.035_346_360_791_375_85,
96    0.044_589_751_324_764_88,
97    0.053_481_524_690_928_09,
98    0.062_009_567_800_670_64,
99    0.069_854_121_318_728_26,
100    0.076_849_680_757_720_38,
101    0.083_080_502_823_133_17,
102    0.088_564_443_056_211_77,
103    0.093_126_598_170_825_32,
104    0.096_642_726_983_623_68,
105    0.099_173_598_721_791_96,
106    0.100_769_845_523_875_6,
107    0.101_330_007_014_791_55,
108    0.100_769_845_523_875_6,
109    0.099_173_598_721_791_96,
110    0.096_642_726_983_623_68,
111    0.093_126_598_170_825_32,
112    0.088_564_443_056_211_77,
113    0.083_080_502_823_133_17,
114    0.076_849_680_757_720_38,
115    0.069_854_121_318_728_26,
116    0.062_009_567_800_670_64,
117    0.053_481_524_690_928_09,
118    0.044_589_751_324_764_88,
119    0.035_346_360_791_375_85,
120    0.025_460_847_326_715_32,
121    0.015_007_947_329_316_122,
122    0.005_377_479_872_923_349,
123];
124
125/// Get Gauss-Kronrod nodes and weights for order `gk_order`.
126///
127/// Supported orders: 15, 31.
128fn gk_nodes_weights(gk_order: usize) -> Result<(&'static [f64], &'static [f64])> {
129    match gk_order {
130        15 => Ok((&GK15_NODES, &GK15_WEIGHTS)),
131        31 => Ok((&GK31_NODES, &GK31_WEIGHTS)),
132        _ => Err(TCIError::InvalidOperation {
133            message: format!(
134                "GK order {} not supported. Supported orders: 15, 31.",
135                gk_order
136            ),
137        }),
138    }
139}
140
141/// Integrate a function over a hypercube `[a, b]` using TCI and
142/// Gauss-Kronrod quadrature.
143///
144/// The integrand is sampled at Gauss-Kronrod nodes in each dimension,
145/// and a tensor-train approximation is built via [`crossinterpolate2`].
146/// The integral is then computed as a weighted sum of the tensor train.
147///
148/// # Arguments
149///
150/// * `f` -- function to integrate, takes a coordinate slice `&[f64]`
151/// * `a` -- lower bounds for each dimension
152/// * `b` -- upper bounds for each dimension
153/// * `gk_order` -- Gauss-Kronrod quadrature order (15 or 31)
154/// * `tci_options` -- options for the TCI approximation
155///
156/// # Returns
157///
158/// The approximate integral value.
159///
160/// # Errors
161///
162/// Returns an error if `a` and `b` have different lengths, or if
163/// `gk_order` is not 15 or 31.
164///
165/// # Examples
166///
167/// ```
168/// use tensor4all_tensorci::integration::integrate;
169/// use tensor4all_tensorci::TCI2Options;
170///
171/// // Integrate (x^2 + y^2) over [0,1]^2 = 2/3
172/// let f = |x: &[f64]| x.iter().map(|&xi| xi * xi).sum::<f64>();
173/// let a = vec![0.0, 0.0];
174/// let b = vec![1.0, 1.0];
175/// let options = TCI2Options { tolerance: 1e-10, ..TCI2Options::default() };
176///
177/// let result: f64 = integrate(&f, &a, &b, 15, options).unwrap();
178/// assert!((result - 2.0 / 3.0).abs() < 1e-8);
179/// ```
180pub fn integrate<T, F>(
181    f: &F,
182    a: &[f64],
183    b: &[f64],
184    gk_order: usize,
185    tci_options: TCI2Options,
186) -> Result<T>
187where
188    T: Scalar + TTScalar + Default + MatrixLuciScalar,
189    DenseFaerLuKernel: PivotKernel<T>,
190    LazyBlockRookKernel: PivotKernel<T>,
191    F: Fn(&[f64]) -> T,
192{
193    if a.len() != b.len() {
194        return Err(TCIError::DimensionMismatch {
195            message: format!(
196                "Lower bounds ({}) and upper bounds ({}) must have same length",
197                a.len(),
198                b.len()
199            ),
200        });
201    }
202
203    let ndims = a.len();
204    let (nodes_ref, weights_ref) = gk_nodes_weights(gk_order)?;
205    let n_nodes = nodes_ref.len();
206
207    // Transform nodes from [-1, 1] to [a_i, b_i] and weights
208    let mut nodes: Vec<Vec<f64>> = Vec::with_capacity(ndims);
209    let mut weights: Vec<Vec<f64>> = Vec::with_capacity(ndims);
210    for d in 0..ndims {
211        let mut dim_nodes = Vec::with_capacity(n_nodes);
212        let mut dim_weights = Vec::with_capacity(n_nodes);
213        for k in 0..n_nodes {
214            // x = (b - a) * (node + 1) / 2 + a
215            dim_nodes.push((b[d] - a[d]) * (nodes_ref[k] + 1.0) / 2.0 + a[d]);
216            // w = (b - a) * weight / 2
217            dim_weights.push((b[d] - a[d]) * weights_ref[k] / 2.0);
218        }
219        nodes.push(dim_nodes);
220        weights.push(dim_weights);
221    }
222
223    let normalization = (gk_order as f64).powi(ndims as i32);
224
225    // Create wrapper function: F(indices) = w_prod * f(x) * normalization
226    let wrapper = |idx: &MultiIndex| -> T {
227        let x: Vec<f64> = (0..ndims).map(|d| nodes[d][idx[d]]).collect();
228        let w: f64 = (0..ndims).map(|d| weights[d][idx[d]]).product();
229        let fval = f(&x);
230        // fval * (w * normalization)
231        let scale = <T as Scalar>::from_f64(w * normalization);
232        fval * scale
233    };
234
235    let local_dims = vec![n_nodes; ndims];
236
237    let (tci, _ranks, _errors) = crossinterpolate2::<T, _, fn(&[MultiIndex]) -> Vec<T>>(
238        wrapper,
239        None,
240        local_dims,
241        vec![],
242        tci_options,
243    )?;
244
245    let tt = tci.to_tensor_train()?;
246    let sum = tt.sum();
247
248    // Return sum / normalization
249    let inv_norm = <T as Scalar>::from_f64(1.0 / normalization);
250    Ok(sum * inv_norm)
251}
252
253#[cfg(test)]
254mod tests {
255    use super::*;
256
257    #[test]
258    fn test_integrate_constant() {
259        // Integral of 1.0 over [0, 1]^2 = 1.0
260        let f = |_x: &[f64]| 1.0f64;
261        let a = vec![0.0, 0.0];
262        let b = vec![1.0, 1.0];
263        let options = TCI2Options {
264            tolerance: 1e-10,
265            ..TCI2Options::default()
266        };
267
268        let result: f64 = integrate(&f, &a, &b, 15, options).unwrap();
269        assert!((result - 1.0).abs() < 1e-8, "Expected 1.0, got {}", result);
270    }
271
272    #[test]
273    fn test_integrate_polynomial() {
274        // Integral of x*y over [0, 1]^2 = 1/4
275        let f = |x: &[f64]| x[0] * x[1];
276        let a = vec![0.0, 0.0];
277        let b = vec![1.0, 1.0];
278        let options = TCI2Options {
279            tolerance: 1e-10,
280            ..TCI2Options::default()
281        };
282
283        let result: f64 = integrate(&f, &a, &b, 15, options).unwrap();
284        assert!(
285            (result - 0.25).abs() < 1e-8,
286            "Expected 0.25, got {}",
287            result
288        );
289    }
290
291    #[test]
292    fn test_integrate_5d_polynomial() {
293        // Integral of product of polynomials over [0, 1]^5
294        // p(x) = 1 + x + x^2, integral = 1 + 1/2 + 1/3 = 11/6
295        // 5D integral = (11/6)^5
296        let polynomial = |x: f64| 1.0 + x + x * x;
297        let f = |coords: &[f64]| coords.iter().map(|&x| polynomial(x)).product::<f64>();
298
299        let n = 5;
300        let a = vec![0.0; n];
301        let b = vec![1.0; n];
302        let options = TCI2Options {
303            tolerance: 1e-10,
304            ..TCI2Options::default()
305        };
306
307        let result: f64 = integrate(&f, &a, &b, 15, options).unwrap();
308        let expected = (11.0 / 6.0_f64).powi(5);
309        assert!(
310            (result - expected).abs() < 1e-6,
311            "Expected {}, got {}",
312            expected,
313            result
314        );
315    }
316
317    #[test]
318    fn test_integrate_with_gk31() {
319        // Same polynomial test with GK31 — should be more accurate
320        let polynomial = |x: f64| 1.0 + x + x * x;
321        let f = |coords: &[f64]| coords.iter().map(|&x| polynomial(x)).product::<f64>();
322
323        let n = 3;
324        let a = vec![0.0; n];
325        let b = vec![1.0; n];
326        let options = TCI2Options {
327            tolerance: 1e-12,
328            ..TCI2Options::default()
329        };
330
331        let result: f64 = integrate(&f, &a, &b, 31, options).unwrap();
332        let expected = (11.0 / 6.0_f64).powi(3);
333        assert!(
334            (result - expected).abs() < 1e-10,
335            "GK31: Expected {}, got {}, diff={}",
336            expected,
337            result,
338            (result - expected).abs()
339        );
340    }
341
342    #[test]
343    fn test_gk_order_unsupported() {
344        let f = |_x: &[f64]| 1.0f64;
345        let a = vec![0.0];
346        let b = vec![1.0];
347        let options = TCI2Options::default();
348        let result: std::result::Result<f64, _> = integrate(&f, &a, &b, 7, options);
349        assert!(result.is_err());
350    }
351
352    /// Port of Julia test_integration.jl: "Integrate real polynomials" with arbitrary bounds
353    #[test]
354    fn test_integrate_arbitrary_bounds() {
355        // p(x) = 1 + 0.5*x + 0.3*x^2
356        // integral of p from a to b = [x + 0.25*x^2 + 0.1*x^3]_a^b
357        let polynomial = |x: f64| 1.0 + 0.5 * x + 0.3 * x * x;
358        let poly_integral = |x: f64| x + 0.25 * x * x + 0.1 * x * x * x;
359        let f = |coords: &[f64]| coords.iter().map(|&x| polynomial(x)).product::<f64>();
360
361        let n = 3;
362        let a = vec![0.2, 0.5, 0.1];
363        let b = vec![0.8, 0.9, 0.7];
364        let expected: f64 = (0..n)
365            .map(|i| poly_integral(b[i]) - poly_integral(a[i]))
366            .product();
367
368        let options = TCI2Options {
369            tolerance: 1e-10,
370            ..TCI2Options::default()
371        };
372        let result: f64 = integrate(&f, &a, &b, 15, options).unwrap();
373        assert!(
374            (result - expected).abs() < 1e-8,
375            "Arbitrary bounds: expected {expected}, got {result}"
376        );
377    }
378
379    /// Port of Julia test_integration.jl: "Integrate complex polynomials"
380    #[test]
381    fn test_integrate_complex_polynomial() {
382        use num_complex::Complex64;
383
384        let polynomial = |x: f64| Complex64::new(1.0 + x, 0.5 * x * x);
385        let poly_integral = |x: f64| Complex64::new(x + 0.5 * x * x, 0.5 * x * x * x / 3.0);
386        let f = |coords: &[f64]| coords.iter().map(|&x| polynomial(x)).product::<Complex64>();
387
388        let n = 3;
389        let a = vec![0.0; n];
390        let b = vec![1.0; n];
391        let expected: Complex64 = (0..n)
392            .map(|_| poly_integral(1.0) - poly_integral(0.0))
393            .product();
394
395        let options = TCI2Options {
396            tolerance: 1e-10,
397            ..TCI2Options::default()
398        };
399        let result: Complex64 = integrate(&f, &a, &b, 15, options).unwrap();
400        assert!(
401            (result - expected).norm() < 1e-8,
402            "Complex: expected {expected}, got {result}"
403        );
404    }
405}
tensor4all_tensorci/integration.rs

tensor4all_tensorci/
integration.rs
