quad_rs/
output.rs

1//! Output abstractions for integration results.
2//!
3//! This module defines [`IntegrationOutput`], the trait implemented by values
4//! that can be returned from an integrand.
5//!
6//! An integration output may be:
7//!
8//! - a real scalar, such as `f64`,
9//! - a complex scalar, such as `Complex<f64>`,
10//! - a vector, matrix, or higher-dimensional array.
11//!
12//! The integrator needs to add, subtract, scale, and divide these values while
13//! forming Gauss–Kronrod estimates. It also needs to reduce output-valued error
14//! estimates to a scalar error used by the adaptive controller.
15//!
16//! [`ErrorNorm`] controls how componentwise errors are reduced for vector-like
17//! outputs.
18
19use nalgebra::ComplexField;
20use num_complex::Complex;
21
22use crate::IntegrableFloat;
23
24/// Strategy used to reduce componentwise output errors to one scalar.
25///
26/// For scalar outputs, both variants are equivalent.
27///
28/// For vector, matrix, or array outputs, this determines how local integration
29/// error is converted to the scalar value used by the adaptive controller.
30#[derive(Copy, Clone, Debug, PartialEq, Default)]
31pub enum ErrorNorm {
32    /// Use the arithmetic mean of component magnitudes.
33    Mean,
34
35    /// Use the largest component magnitude.
36    #[default]
37    Max,
38}
39
40/// Type that can be used as an integrand output.
41///
42/// The integrator forms linear combinations of integrand outputs using scalar
43/// weights from the integration domain. Therefore an output must support
44/// addition, subtraction, multiplication by its scalar type, and division by its
45/// scalar type.
46///
47/// `IntegrationOutput` also provides scalar diagnostics:
48///
49/// - [`modulus`](Self::modulus), used for absolute integral estimates,
50/// - [`is_finite`](Self::is_finite), used to detect invalid integrand values,
51/// - [`reduce_error`](Self::reduce_error), used to convert an output-valued
52///   local error estimate into a scalar adaptive error.
53///
54/// # Associated types
55///
56/// - [`Float`](Self::Float): underlying real floating-point type.
57pub trait IntegrationOutput<S>: Clone + Default {
58    /// Underlying real floating-point type.
59    type Float: IntegrableFloat;
60
61    fn add(&self, other: &Self) -> Self;
62    fn sub(&self, other: &Self) -> Self;
63    fn mul_scalar(&self, scalar: &S) -> Self;
64
65    /// Returns a scalar magnitude for this output.
66    ///
67    /// For scalar outputs this is the absolute value or complex modulus. For
68    /// array-like outputs this should return a norm-like aggregate magnitude.
69    fn modulus(&self) -> Self::Float;
70
71    /// Returns `false` if this output contains `NaN` or infinity.
72    fn is_finite(&self) -> bool;
73
74    /// Returns the largest component magnitude.
75    fn max_component(&self) -> Self::Float;
76
77    /// Returns the mean component magnitude.
78    fn mean_component(&self) -> Self::Float;
79
80    /// Reduces this output to a scalar error according to `mode`.
81    fn reduce_error(&self, mode: ErrorNorm) -> Self::Float {
82        match mode {
83            ErrorNorm::Mean => self.mean_component(),
84            ErrorNorm::Max => self.max_component(),
85        }
86    }
87}
88
89impl IntegrationOutput<Complex<f32>> for Complex<f32> {
90    type Float = f32;
91
92    fn add(&self, other: &Self) -> Self {
93        self + other
94    }
95    fn sub(&self, other: &Self) -> Self {
96        self - other
97    }
98    fn mul_scalar(&self, scalar: &Complex<f32>) -> Self {
99        self * scalar
100    }
101
102    fn modulus(&self) -> Self::Float {
103        <Self as ComplexField>::modulus(*self)
104    }
105
106    fn is_finite(&self) -> bool {
107        ComplexField::is_finite(self)
108    }
109
110    fn max_component(&self) -> Self::Float {
111        <Self as IntegrationOutput<Complex<f32>>>::modulus(self)
112    }
113
114    fn mean_component(&self) -> Self::Float {
115        <Self as IntegrationOutput<Complex<f32>>>::modulus(self)
116    }
117}
118
119impl IntegrationOutput<f32> for Complex<f32> {
120    type Float = f32;
121
122    fn add(&self, other: &Self) -> Self {
123        self + other
124    }
125    fn sub(&self, other: &Self) -> Self {
126        self - other
127    }
128    fn mul_scalar(&self, scalar: &f32) -> Self {
129        self * scalar
130    }
131
132    fn modulus(&self) -> Self::Float {
133        <Self as ComplexField>::modulus(*self)
134    }
135
136    fn is_finite(&self) -> bool {
137        ComplexField::is_finite(self)
138    }
139
140    fn max_component(&self) -> Self::Float {
141        <Self as IntegrationOutput<f32>>::modulus(self)
142    }
143
144    fn mean_component(&self) -> Self::Float {
145        <Self as IntegrationOutput<f32>>::modulus(self)
146    }
147}
148
149impl IntegrationOutput<f32> for f32 {
150    type Float = Self;
151
152    fn add(&self, other: &Self) -> Self {
153        self + other
154    }
155    fn sub(&self, other: &Self) -> Self {
156        self - other
157    }
158    fn mul_scalar(&self, scalar: &f32) -> Self {
159        self * scalar
160    }
161
162    fn modulus(&self) -> Self::Float {
163        <Self as ComplexField>::modulus(*self)
164    }
165
166    fn is_finite(&self) -> bool {
167        ComplexField::is_finite(self)
168    }
169
170    fn max_component(&self) -> Self::Float {
171        self.modulus()
172    }
173
174    fn mean_component(&self) -> Self::Float {
175        self.modulus()
176    }
177}
178
179impl IntegrationOutput<Complex<f64>> for Complex<f64> {
180    type Float = f64;
181
182    fn add(&self, other: &Self) -> Self {
183        self + other
184    }
185    fn sub(&self, other: &Self) -> Self {
186        self - other
187    }
188    fn mul_scalar(&self, scalar: &Complex<f64>) -> Self {
189        self * scalar
190    }
191
192    fn modulus(&self) -> Self::Float {
193        <Self as ComplexField>::modulus(*self)
194    }
195
196    fn is_finite(&self) -> bool {
197        ComplexField::is_finite(self)
198    }
199
200    fn max_component(&self) -> Self::Float {
201        <Self as IntegrationOutput<Complex<f64>>>::modulus(self)
202    }
203
204    fn mean_component(&self) -> Self::Float {
205        <Self as IntegrationOutput<Complex<f64>>>::modulus(self)
206    }
207}
208
209impl IntegrationOutput<f64> for Complex<f64> {
210    type Float = f64;
211
212    fn add(&self, other: &Self) -> Self {
213        self + other
214    }
215    fn sub(&self, other: &Self) -> Self {
216        self - other
217    }
218    fn mul_scalar(&self, scalar: &f64) -> Self {
219        self * scalar
220    }
221
222    fn modulus(&self) -> Self::Float {
223        <Self as ComplexField>::modulus(*self)
224    }
225
226    fn is_finite(&self) -> bool {
227        ComplexField::is_finite(self)
228    }
229
230    fn max_component(&self) -> Self::Float {
231        <Self as IntegrationOutput<f64>>::modulus(self)
232    }
233
234    fn mean_component(&self) -> Self::Float {
235        <Self as IntegrationOutput<f64>>::modulus(self)
236    }
237}
238
239impl IntegrationOutput<f64> for f64 {
240    type Float = Self;
241
242    fn add(&self, other: &Self) -> Self {
243        self + other
244    }
245    fn sub(&self, other: &Self) -> Self {
246        self - other
247    }
248    fn mul_scalar(&self, scalar: &f64) -> Self {
249        self * scalar
250    }
251
252    fn modulus(&self) -> Self::Float {
253        <Self as ComplexField>::modulus(*self)
254    }
255
256    fn is_finite(&self) -> bool {
257        ComplexField::is_finite(self)
258    }
259
260    fn max_component(&self) -> Self::Float {
261        <Self as IntegrationOutput<f64>>::modulus(self)
262    }
263
264    fn mean_component(&self) -> Self::Float {
265        <Self as IntegrationOutput<f64>>::modulus(self)
266    }
267}
268
269#[cfg(feature = "ndarray")]
270use ndarray::{Array, Dimension};
271
272#[cfg(feature = "ndarray")]
273impl<D> IntegrationOutput<f32> for Array<f32, D>
274where
275    D: Dimension,
276{
277    type Float = f32;
278
279    fn add(&self, other: &Self) -> Self {
280        self + other
281    }
282
283    fn sub(&self, other: &Self) -> Self {
284        self - other
285    }
286
287    fn mul_scalar(&self, scalar: &f32) -> Self {
288        self.mapv(|value| value * *scalar)
289    }
290
291    fn modulus(&self) -> Self::Float {
292        self.iter()
293            .map(|value| {
294                let x = value.modulus();
295                x * x
296            })
297            .sum::<Self::Float>()
298            .sqrt()
299    }
300
301    fn is_finite(&self) -> bool {
302        self.iter().all(|value| f32::is_finite(*value))
303    }
304
305    fn max_component(&self) -> Self::Float {
306        self.iter().map(|value| value.modulus()).fold(
307            <Self::Float as num_traits::Float>::neg_zero(),
308            |acc, value| {
309                if value > acc { value } else { acc }
310            },
311        )
312    }
313
314    fn mean_component(&self) -> Self::Float {
315        if self.is_empty() {
316            return <Self::Float as num_traits::Float>::neg_zero();
317        }
318
319        let sum = self
320            .iter()
321            .map(|value| value.modulus())
322            .sum::<Self::Float>();
323        sum / <Self::Float as num_traits::FromPrimitive>::from_usize(self.len()).unwrap()
324    }
325}
326
327#[cfg(feature = "ndarray")]
328impl<D> IntegrationOutput<f32> for Array<Complex<f32>, D>
329where
330    D: Dimension,
331{
332    type Float = f32;
333
334    fn add(&self, other: &Self) -> Self {
335        self + other
336    }
337
338    fn sub(&self, other: &Self) -> Self {
339        self - other
340    }
341
342    fn mul_scalar(&self, scalar: &f32) -> Self {
343        self.mapv(|value| value * *scalar)
344    }
345
346    fn modulus(&self) -> Self::Float {
347        self.iter()
348            .map(|value| {
349                let x = value.abs();
350                x * x
351            })
352            .sum::<Self::Float>()
353            .sqrt()
354    }
355
356    fn is_finite(&self) -> bool {
357        self.iter().all(|value| Complex::is_finite(*value))
358    }
359
360    fn max_component(&self) -> Self::Float {
361        self.iter().map(|value| value.abs()).fold(
362            <Self::Float as num_traits::Float>::neg_zero(),
363            |acc, value| {
364                if value > acc { value } else { acc }
365            },
366        )
367    }
368
369    fn mean_component(&self) -> Self::Float {
370        if self.is_empty() {
371            return <Self::Float as num_traits::Float>::neg_zero();
372        }
373
374        let sum = self.iter().map(|value| value.abs()).sum::<Self::Float>();
375        sum / <Self::Float as num_traits::FromPrimitive>::from_usize(self.len()).unwrap()
376    }
377}
378
379#[cfg(feature = "ndarray")]
380impl<D> IntegrationOutput<Complex<f32>> for Array<Complex<f32>, D>
381where
382    D: Dimension,
383{
384    type Float = f32;
385
386    fn add(&self, other: &Self) -> Self {
387        self + other
388    }
389
390    fn sub(&self, other: &Self) -> Self {
391        self - other
392    }
393
394    fn mul_scalar(&self, scalar: &Complex<f32>) -> Self {
395        self.mapv(|value| value * *scalar)
396    }
397
398    fn modulus(&self) -> Self::Float {
399        self.iter()
400            .map(|value| {
401                let x = value.abs();
402                x * x
403            })
404            .sum::<Self::Float>()
405            .sqrt()
406    }
407
408    fn is_finite(&self) -> bool {
409        self.iter().all(|value| Complex::is_finite(*value))
410    }
411
412    fn max_component(&self) -> Self::Float {
413        self.iter().map(|value| value.abs()).fold(
414            <Self::Float as num_traits::Float>::neg_zero(),
415            |acc, value| {
416                if value > acc { value } else { acc }
417            },
418        )
419    }
420
421    fn mean_component(&self) -> Self::Float {
422        if self.is_empty() {
423            return <Self::Float as num_traits::Float>::neg_zero();
424        }
425
426        let sum = self.iter().map(|value| value.abs()).sum::<Self::Float>();
427        sum / <Self::Float as num_traits::FromPrimitive>::from_usize(self.len()).unwrap()
428    }
429}
430
431#[cfg(feature = "ndarray")]
432impl<D> IntegrationOutput<f64> for Array<f64, D>
433where
434    D: Dimension,
435{
436    type Float = f64;
437
438    fn add(&self, other: &Self) -> Self {
439        self + other
440    }
441
442    fn sub(&self, other: &Self) -> Self {
443        self - other
444    }
445
446    fn mul_scalar(&self, scalar: &f64) -> Self {
447        self.mapv(|value| value * *scalar)
448    }
449
450    fn modulus(&self) -> Self::Float {
451        self.iter()
452            .map(|value| {
453                let x = value.abs();
454                x * x
455            })
456            .sum::<Self::Float>()
457            .sqrt()
458    }
459
460    fn is_finite(&self) -> bool {
461        self.iter().all(|value| f64::is_finite(*value))
462    }
463
464    fn max_component(&self) -> Self::Float {
465        self.iter().map(|value| value.abs()).fold(
466            <Self::Float as num_traits::Float>::neg_zero(),
467            |acc, value| {
468                if value > acc { value } else { acc }
469            },
470        )
471    }
472
473    fn mean_component(&self) -> Self::Float {
474        if self.is_empty() {
475            return <Self::Float as num_traits::Float>::neg_zero();
476        }
477
478        let sum = self.iter().map(|value| value.abs()).sum::<Self::Float>();
479        sum / <Self::Float as num_traits::FromPrimitive>::from_usize(self.len()).unwrap()
480    }
481}
482
483#[cfg(feature = "ndarray")]
484impl<D> IntegrationOutput<f64> for Array<Complex<f64>, D>
485where
486    D: Dimension,
487{
488    type Float = f64;
489
490    fn add(&self, other: &Self) -> Self {
491        self + other
492    }
493
494    fn sub(&self, other: &Self) -> Self {
495        self - other
496    }
497
498    fn mul_scalar(&self, scalar: &f64) -> Self {
499        self.mapv(|value| value * *scalar)
500    }
501
502    fn modulus(&self) -> Self::Float {
503        self.iter()
504            .map(|value| {
505                let x = value.abs();
506                x * x
507            })
508            .sum::<Self::Float>()
509            .sqrt()
510    }
511
512    fn is_finite(&self) -> bool {
513        self.iter().all(|value| Complex::is_finite(*value))
514    }
515
516    fn max_component(&self) -> Self::Float {
517        self.iter().map(|value| value.abs()).fold(
518            <Self::Float as num_traits::Float>::neg_zero(),
519            |acc, value| {
520                if value > acc { value } else { acc }
521            },
522        )
523    }
524
525    fn mean_component(&self) -> Self::Float {
526        if self.is_empty() {
527            return <Self::Float as num_traits::Float>::neg_zero();
528        }
529
530        let sum = self.iter().map(|value| value.abs()).sum::<Self::Float>();
531        sum / <Self::Float as num_traits::FromPrimitive>::from_usize(self.len()).unwrap()
532    }
533}
534
535#[cfg(feature = "ndarray")]
536impl<D> IntegrationOutput<Complex<f64>> for Array<Complex<f64>, D>
537where
538    D: Dimension,
539{
540    type Float = f64;
541
542    fn add(&self, other: &Self) -> Self {
543        self + other
544    }
545
546    fn sub(&self, other: &Self) -> Self {
547        self - other
548    }
549
550    fn mul_scalar(&self, scalar: &Complex<f64>) -> Self {
551        self.mapv(|value| value * *scalar)
552    }
553
554    fn modulus(&self) -> Self::Float {
555        self.iter()
556            .map(|value| {
557                let x = value.abs();
558                x * x
559            })
560            .sum::<Self::Float>()
561            .sqrt()
562    }
563
564    fn is_finite(&self) -> bool {
565        self.iter().all(|value| Complex::is_finite(*value))
566    }
567
568    fn max_component(&self) -> Self::Float {
569        self.iter().map(|value| value.abs()).fold(
570            <Self::Float as num_traits::Float>::neg_zero(),
571            |acc, value| {
572                if value > acc { value } else { acc }
573            },
574        )
575    }
576
577    fn mean_component(&self) -> Self::Float {
578        if self.is_empty() {
579            return <Self::Float as num_traits::Float>::neg_zero();
580        }
581
582        let sum = self.iter().map(|value| value.abs()).sum::<Self::Float>();
583        sum / <Self::Float as num_traits::FromPrimitive>::from_usize(self.len()).unwrap()
584    }
585}
586
587#[cfg(test)]
588mod tests {
589    use super::*;
590    use num_complex::Complex;
591    use num_traits::{Float, FromPrimitive};
592
593    fn assert_close<F: Float + FromPrimitive + std::fmt::Display>(a: F, b: F) {
594        assert!(
595            (a - b).abs() < F::from_f64(1e-12).unwrap(),
596            "expected {b}, got {a}, diff = {}",
597            (a - b).abs()
598        );
599    }
600
601    #[test]
602    fn f64_output_reductions_are_absolute_value() {
603        let x = -3.5_f64;
604
605        assert_close(x.modulus(), 3.5);
606        assert_close(x.max_component(), 3.5);
607        assert_close(x.mean_component(), 3.5);
608        assert_close(x.reduce_error(ErrorNorm::Max), 3.5);
609        assert_close(x.reduce_error(ErrorNorm::Mean), 3.5);
610        assert!(x.is_finite());
611    }
612
613    #[test]
614    fn complex_output_reductions_use_modulus() {
615        let z: Complex<f64> = Complex::new(3.0_f64, 4.0);
616
617        assert_close(<Complex<f64> as IntegrationOutput<f64>>::modulus(&z), 5.0);
618        assert_close(
619            <Complex<f64> as IntegrationOutput<f64>>::max_component(&z),
620            5.0,
621        );
622        assert_close(
623            <Complex<f64> as IntegrationOutput<f64>>::mean_component(&z),
624            5.0,
625        );
626        assert_close(
627            <Complex<f64> as IntegrationOutput<f64>>::reduce_error(&z, ErrorNorm::Max),
628            5.0,
629        );
630        assert_close(
631            <Complex<f64> as IntegrationOutput<f64>>::reduce_error(&z, ErrorNorm::Mean),
632            5.0,
633        );
634        assert!(z.is_finite());
635    }
636
637    #[test]
638    fn non_finite_scalars_are_detected() {
639        assert!(!f64::NAN.is_finite());
640        assert!(!f64::INFINITY.is_finite());
641
642        let z = Complex::new(1.0_f64, f64::NAN);
643        assert!(!<Complex<f64> as IntegrationOutput<f64>>::is_finite(&z));
644    }
645
646    #[cfg(feature = "ndarray")]
647    mod ndarray_tests {
648        use super::*;
649        use ndarray::{Array1, Array2, array};
650
651        #[test]
652        fn array1_f64_real_output_reductions_work() {
653            let x: Array1<f64> = array![-1.0, 2.0, -3.0];
654
655            assert_close(
656                <Array1<f64> as IntegrationOutput<f64>>::modulus(&x),
657                Float::sqrt(14.0),
658            );
659            assert_close(
660                <Array1<f64> as IntegrationOutput<f64>>::max_component(&x),
661                3.0,
662            );
663            assert_close(
664                <Array1<f64> as IntegrationOutput<f64>>::mean_component(&x),
665                2.0,
666            );
667            assert_close(
668                <Array1<f64> as IntegrationOutput<f64>>::reduce_error(&x, ErrorNorm::Max),
669                3.0,
670            );
671            assert_close(
672                <Array1<f64> as IntegrationOutput<f64>>::reduce_error(&x, ErrorNorm::Mean),
673                2.0,
674            );
675            assert!(x.is_finite());
676        }
677
678        #[test]
679        fn array1_f32_real_output_reductions_work() {
680            let x: Array1<f32> = array![-1.0, 2.0, -3.0];
681
682            assert_close(
683                <Array1<f32> as IntegrationOutput<f32>>::modulus(&x),
684                Float::sqrt(14.0),
685            );
686            assert_close(
687                <Array1<f32> as IntegrationOutput<f32>>::max_component(&x),
688                3.0,
689            );
690            assert_close(
691                <Array1<f32> as IntegrationOutput<f32>>::mean_component(&x),
692                2.0,
693            );
694            assert_close(
695                <Array1<f32> as IntegrationOutput<f32>>::reduce_error(&x, ErrorNorm::Max),
696                3.0,
697            );
698            assert_close(
699                <Array1<f32> as IntegrationOutput<f32>>::reduce_error(&x, ErrorNorm::Mean),
700                2.0,
701            );
702            assert!(x.is_finite());
703        }
704
705        #[test]
706        fn array2_f64_real_output_reductions_work() {
707            let x: Array2<f64> = array![[1.0, -2.0], [-3.0, 4.0],];
708
709            assert_close(
710                <Array2<f64> as IntegrationOutput<f64>>::modulus(&x),
711                Float::sqrt(30.0),
712            );
713            assert_close(
714                <Array2<f64> as IntegrationOutput<f64>>::max_component(&x),
715                4.0,
716            );
717            assert_close(
718                <Array2<f64> as IntegrationOutput<f64>>::mean_component(&x),
719                2.5,
720            );
721            assert_close(
722                <Array2<f64> as IntegrationOutput<f64>>::reduce_error(&x, ErrorNorm::Max),
723                4.0,
724            );
725            assert_close(
726                <Array2<f64> as IntegrationOutput<f64>>::reduce_error(&x, ErrorNorm::Mean),
727                2.5,
728            );
729            assert!(x.is_finite());
730        }
731
732        #[test]
733        fn array2_f32_real_output_reductions_work() {
734            let x: Array2<f32> = array![[1.0, -2.0], [-3.0, 4.0],];
735
736            assert_close(
737                <Array2<f32> as IntegrationOutput<f32>>::modulus(&x),
738                Float::sqrt(30.0),
739            );
740            assert_close(
741                <Array2<f32> as IntegrationOutput<f32>>::max_component(&x),
742                4.0,
743            );
744            assert_close(
745                <Array2<f32> as IntegrationOutput<f32>>::mean_component(&x),
746                2.5,
747            );
748            assert_close(
749                <Array2<f32> as IntegrationOutput<f32>>::reduce_error(&x, ErrorNorm::Max),
750                4.0,
751            );
752            assert_close(
753                <Array2<f32> as IntegrationOutput<f32>>::reduce_error(&x, ErrorNorm::Mean),
754                2.5,
755            );
756            assert!(x.is_finite());
757        }
758
759        #[test]
760        fn array1_complex_f64_output_reductions_with_real_scalar_work() {
761            let x: Array1<Complex<f64>> = array![
762                Complex::new(3.0, 4.0),
763                Complex::new(5.0, 12.0),
764                Complex::new(8.0, 15.0),
765            ];
766
767            assert_close(
768                <Array1<Complex<f64>> as IntegrationOutput<f64>>::modulus(&x),
769                (x[0] * x[0].conj() + x[1] * x[1].conj() + x[2] * x[2].conj())
770                    .sqrt()
771                    .re,
772            );
773            assert_close(
774                <Array1<Complex<f64>> as IntegrationOutput<f64>>::max_component(&x),
775                17.0,
776            );
777            assert_close(
778                <Array1<Complex<f64>> as IntegrationOutput<f64>>::mean_component(&x),
779                (5.0 + 13.0 + 17.0) / 3.0,
780            );
781            assert!(<Array1<Complex<f64>> as IntegrationOutput<f64>>::is_finite(
782                &x
783            ));
784        }
785
786        #[test]
787        fn array1_complex_f64_output_reductions_with_complex_scalar_work() {
788            let x: Array1<Complex<f64>> = array![
789                Complex::new(3.0, 4.0),
790                Complex::new(5.0, 12.0),
791                Complex::new(8.0, 15.0),
792            ];
793
794            assert_close(
795                <Array1<Complex<f64>> as IntegrationOutput<Complex<f64>>>::modulus(&x),
796                (x[0] * x[0].conj() + x[1] * x[1].conj() + x[2] * x[2].conj())
797                    .sqrt()
798                    .re,
799            );
800            assert_close(
801                <Array1<Complex<f64>> as IntegrationOutput<Complex<f64>>>::max_component(&x),
802                17.0,
803            );
804            assert_close(
805                <Array1<Complex<f64>> as IntegrationOutput<Complex<f64>>>::mean_component(&x),
806                (5.0 + 13.0 + 17.0) / 3.0,
807            );
808            assert!(<Array1<Complex<f64>> as IntegrationOutput<Complex<f64>>>::is_finite(&x));
809        }
810
811        #[test]
812        fn array1_complex_f32_output_reductions_with_real_scalar_work() {
813            let x: Array1<Complex<f32>> = array![
814                Complex::new(3.0, 4.0),
815                Complex::new(5.0, 12.0),
816                Complex::new(8.0, 15.0),
817            ];
818
819            assert_close(
820                <Array1<Complex<f32>> as IntegrationOutput<f32>>::modulus(&x),
821                (x[0] * x[0].conj() + x[1] * x[1].conj() + x[2] * x[2].conj())
822                    .sqrt()
823                    .re,
824            );
825            assert_close(
826                <Array1<Complex<f32>> as IntegrationOutput<f32>>::max_component(&x),
827                17.0,
828            );
829            assert_close(
830                <Array1<Complex<f32>> as IntegrationOutput<f32>>::mean_component(&x),
831                (5.0 + 13.0 + 17.0) / 3.0,
832            );
833            assert!(<Array1<Complex<f32>> as IntegrationOutput<f32>>::is_finite(
834                &x
835            ));
836        }
837
838        #[test]
839        fn array1_complex_f32_output_reductions_with_complex_scalar_work() {
840            let x: Array1<Complex<f32>> = array![
841                Complex::new(3.0, 4.0),
842                Complex::new(5.0, 12.0),
843                Complex::new(8.0, 15.0),
844            ];
845
846            assert_close(
847                <Array1<Complex<f32>> as IntegrationOutput<Complex<f32>>>::modulus(&x),
848                (x[0] * x[0].conj() + x[1] * x[1].conj() + x[2] * x[2].conj())
849                    .sqrt()
850                    .re,
851            );
852            assert_close(
853                <Array1<Complex<f32>> as IntegrationOutput<Complex<f32>>>::max_component(&x),
854                17.0,
855            );
856            assert_close(
857                <Array1<Complex<f32>> as IntegrationOutput<Complex<f32>>>::mean_component(&x),
858                (5.0 + 13.0 + 17.0) / 3.0,
859            );
860            assert!(<Array1<Complex<f32>> as IntegrationOutput<Complex<f32>>>::is_finite(&x));
861        }
862
863        #[test]
864        fn array2_complex_f64_output_reductions_with_real_scalar_work() {
865            let x: Array2<Complex<f64>> = array![
866                [Complex::new(3.0, 4.0), Complex::new(0.0, 2.0)],
867                [Complex::new(5.0, 12.0), Complex::new(1.0, 0.0)],
868            ];
869
870            let modulus = x
871                .iter()
872                .map(|each| each * each.conj())
873                .map(|each| each.re)
874                .sum::<f64>()
875                .sqrt();
876
877            assert_close(
878                <Array2<Complex<f64>> as IntegrationOutput<f64>>::modulus(&x),
879                modulus,
880            );
881            assert_close(
882                <Array2<Complex<f64>> as IntegrationOutput<f64>>::max_component(&x),
883                13.0,
884            );
885            assert_close(
886                <Array2<Complex<f64>> as IntegrationOutput<f64>>::mean_component(&x),
887                (5.0 + 2.0 + 13.0 + 1.0) / 4.0,
888            );
889            assert!(<Array2<Complex<f64>> as IntegrationOutput<f64>>::is_finite(
890                &x
891            ));
892        }
893
894        #[test]
895        fn array2_complex_f64_output_reductions_with_complex_scalar_work() {
896            let x: Array2<Complex<f64>> = array![
897                [Complex::new(3.0, 4.0), Complex::new(0.0, 2.0)],
898                [Complex::new(5.0, 12.0), Complex::new(1.0, 0.0)],
899            ];
900
901            let modulus = x
902                .iter()
903                .map(|each| each * each.conj())
904                .map(|each| each.re)
905                .sum::<f64>()
906                .sqrt();
907
908            assert_close(
909                <Array2<Complex<f64>> as IntegrationOutput<Complex<f64>>>::modulus(&x),
910                modulus,
911            );
912            assert_close(
913                <Array2<Complex<f64>> as IntegrationOutput<Complex<f64>>>::max_component(&x),
914                13.0,
915            );
916            assert_close(
917                <Array2<Complex<f64>> as IntegrationOutput<Complex<f64>>>::mean_component(&x),
918                (5.0 + 2.0 + 13.0 + 1.0) / 4.0,
919            );
920            assert!(<Array2<Complex<f64>> as IntegrationOutput<Complex<f64>>>::is_finite(&x));
921        }
922
923        #[test]
924        fn array2_complex_f32_output_reductions_with_real_scalar_work() {
925            let x: Array2<Complex<f32>> = array![
926                [Complex::new(3.0, 4.0), Complex::new(0.0, 2.0)],
927                [Complex::new(5.0, 12.0), Complex::new(1.0, 0.0)],
928            ];
929
930            let modulus = x
931                .iter()
932                .map(|each| each * each.conj())
933                .map(|each| each.re)
934                .sum::<f32>()
935                .sqrt();
936
937            assert_close(
938                <Array2<Complex<f32>> as IntegrationOutput<f32>>::modulus(&x),
939                modulus,
940            );
941            assert_close(
942                <Array2<Complex<f32>> as IntegrationOutput<f32>>::max_component(&x),
943                13.0,
944            );
945            assert_close(
946                <Array2<Complex<f32>> as IntegrationOutput<f32>>::mean_component(&x),
947                (5.0 + 2.0 + 13.0 + 1.0) / 4.0,
948            );
949            assert!(<Array2<Complex<f32>> as IntegrationOutput<f32>>::is_finite(
950                &x
951            ));
952        }
953
954        #[test]
955        fn array2_complex_f32_output_reductions_with_complex_scalar_work() {
956            let x: Array2<Complex<f32>> = array![
957                [Complex::new(3.0, 4.0), Complex::new(0.0, 2.0)],
958                [Complex::new(5.0, 12.0), Complex::new(1.0, 0.0)],
959            ];
960
961            let modulus = x
962                .iter()
963                .map(|each| each * each.conj())
964                .map(|each| each.re)
965                .sum::<f32>()
966                .sqrt();
967
968            assert_close(
969                <Array2<Complex<f32>> as IntegrationOutput<Complex<f32>>>::modulus(&x),
970                modulus,
971            );
972            assert_close(
973                <Array2<Complex<f32>> as IntegrationOutput<Complex<f32>>>::max_component(&x),
974                13.0,
975            );
976            assert_close(
977                <Array2<Complex<f32>> as IntegrationOutput<Complex<f32>>>::mean_component(&x),
978                (5.0 + 2.0 + 13.0 + 1.0) / 4.0,
979            );
980            assert!(<Array2<Complex<f32>> as IntegrationOutput<Complex<f32>>>::is_finite(&x));
981        }
982
983        #[test]
984        fn array_non_finite_values_are_detected() {
985            let x: Array1<f64> = array![1.0, f64::NAN, 3.0];
986            assert!(!x.is_finite());
987
988            let z: Array1<Complex<f64>> =
989                array![Complex::new(1.0, 0.0), Complex::new(f64::INFINITY, 0.0),];
990            assert!(!<Array1<Complex<f64>> as IntegrationOutput<f64>>::is_finite(&z));
991            assert!(!<Array1<Complex<f64>> as IntegrationOutput<Complex<f64>>>::is_finite(&z));
992        }
993
994        #[test]
995        fn empty_array_reductions_are_well_defined() {
996            let x: Array1<f64> = Array1::from_vec(vec![]);
997
998            assert_close(<Array1<f64> as IntegrationOutput<f64>>::modulus(&x), 0.0);
999            assert_close(
1000                <Array1<f64> as IntegrationOutput<f64>>::max_component(&x),
1001                0.0,
1002            );
1003            assert_close(
1004                <Array1<f64> as IntegrationOutput<f64>>::mean_component(&x),
1005                0.0,
1006            );
1007            assert!(<Array1<f64> as IntegrationOutput<f64>>::is_finite(&x));
1008        }
1009    }
1010}
quad_rs/output.rs

quad_rs/
output.rs
