diffusers 0.3.1

//! The double exponential algorithm is naturally adaptive, it stops calling the integrand when the error is reduced to below the desired threshold.
//! It also does not allocate. No box, no vec, etc.
//! It has a hard coded maximum of approximately 350 function evaluations. This guarantees that the algorithm will return.
//! The error in the algorithm decreases exponentially in the number of function evaluations, specifically O(exp(-cN/log(N))).
//! So if 350 function evaluations is not giving the desired accuracy than the programmer probably needs to give some guidance by splitting up the range at singularities or [other preparation techniques](http://www.johndcook.com/blog/2012/02/21/care-and-treatment-of-singularities/).
//!
//! This code is mostly copied and adapted from `quadrature` https://github.com/Eh2406/quadrature
use std::f64;

#[derive(Clone, Copy, Debug)]
pub struct Output {
    pub num_function_evaluations: u32,
    pub error_estimate: f64,
    pub integral: f64,
}

impl Output {
    fn scale(self, c: f64) -> Self {
        Output {
            num_function_evaluations: self.num_function_evaluations,
            error_estimate: c * self.error_estimate,
            integral: c * self.integral,
        }
    }
}

pub fn integrate<F>(f: F, a: f64, b: f64, target_absolute_error: f64) -> Output
where
    F: Fn(f64) -> f64,
{
    // Apply the linear change of variables x = ct + d
    // $$\int_a^b f(x) dx = c \int_{-1}^1 f( ct + d ) dt$$
    // c = (b-a)/2, d = (a+b)/2
    let c = 0.5 * (b - a);
    let d = 0.5 * (a + b);
    integrate_core(
        |x| {
            let out = f(c * x + d);
            if out.is_finite() {
                out
            } else {
                0.0
            }
        },
        0.25 * target_absolute_error / c,
    )
    .scale(c)
}

/// Integrate f(x) from [-1.0, 1.0]
fn integrate_core<F>(f: F, target_absolute_error: f64) -> Output
where
    F: Fn(f64) -> f64,
{
    let mut error_estimate = f64::MAX;
    let mut num_function_evaluations = 1;
    let mut current_delta = f64::MAX;

    let mut integral = 2.0 * f64::consts::FRAC_PI_2 * f(0.0);

    for &weight in &WEIGHTS {
        let new_contribution =
            weight.iter().map(|&(w, x)| w * (f(x) + f(-x))).fold(0.0, |sum, x| sum + x);
        num_function_evaluations += 2 * weight.len();

        // difference in consecutive integral estimates
        let previous_delta_ln = current_delta.ln();
        current_delta = (0.5 * integral - new_contribution).abs();
        integral = 0.5 * integral + new_contribution;

        // Once convergence kicks in, error is approximately squared at each step.
        // Determine whether we're in the convergent region by looking at the trend in the error.
        if num_function_evaluations <= 13 {
            // level <= 1
            continue; // previousDelta meaningless, so cannot check convergence.
        }

        // Exact comparison with zero is harmless here.  Could possibly be replaced with
        // a small positive upper limit on the size of currentDelta, but determining
        // that upper limit would be difficult.  At worse, the loop is executed more
        // times than necessary.  But no infinite loop can result since there is
        // an upper bound on the loop variable.
        if current_delta == 0.0 {
            error_estimate = 0.0;
            break;
        }
        // previousDelta != 0 or would have been kicked out previously
        let r = current_delta.ln() / previous_delta_ln;

        if r > 1.9 && r < 2.1 {
            // If convergence theory applied perfectly, r would be 2 in the convergence region.
            // r close to 2 is good enough. We expect the difference between this integral estimate
            // and the next one to be roughly delta^2.
            error_estimate = current_delta * current_delta;
        } else {
            // Not in the convergence region.  Assume only that error is decreasing.
            error_estimate = current_delta;
        }

        if error_estimate < target_absolute_error {
            break;
        }
    }

    Output { num_function_evaluations: num_function_evaluations as u32, error_estimate, integral }
}

pub const WEIGHTS: [&[(f64, f64)]; 7] = [
    &[
        // First layer weights
        (0.230_022_394_514_788_68, 0.951_367_964_072_746_9),
        (0.000_266_200_513_752_716_93, 0.999_977_477_192_461_6),
        (1.358_178_427_453_909_1e-12, 0.999_999_999_999_957),
    ],
    &[
        // 2nd layer weights and abscissas: transformed 1/2, 3/2, 5/2
        (0.5 * 0.965_976_579_412_301_2, 0.674_271_492_248_435_9),
        (0.5 * 0.018_343_166_989_927_842, 0.997_514_856_457_224_4),
        (0.5 * 2.143_120_455_694_304e-7, 0.999_999_988_875_664_9),
    ],
    &[
        // 3rd layer weights and abscissas: transformed 1/4, 3/4, ...
        (0.25 * 1.389_614_759_247_256_3, 0.377_209_738_164_034_2),
        (0.25 * 0.531_078_275_428_054, 0.859_569_058_689_896_6),
        (0.25 * 0.076_385_743_570_832_3, 0.987_040_560_507_376_9),
        (0.25 * 0.002_902_517_747_901_313_7, 0.999_688_264_028_353_2),
        (0.25 * 0.000_011_983_701_363_170_72, 0.999_999_204_737_114_7),
        (0.25 * 1.163_116_581_425_578_2e-9, 0.999_999_999_952_856_5),
    ],
    &[
        // 4th layer weights and abscissas: transformed 1/8, 3/8, ...
        (0.125 * 1.523_283_718_634_705_2, 0.194_357_003_324_935_44),
        (0.125 * 1.193_463_025_849_157, 0.539_146_705_387_967_7),
        (0.125 * 0.737_437_848_361_547_8, 0.780_607_438_983_200_3),
        (0.125 * 0.360_461_418_469_343_65, 0.914_879_263_264_574_6),
        (0.125 * 0.137_422_107_733_167_74, 0.973_966_868_195_677_5),
        (0.125 * 0.039_175_005_493_600_78, 0.994_055_506_631_402_2),
        (0.125 * 0.007_742_601_026_064_241, 0.999_065_196_455_785_8),
        (0.125 * 0.000_949_946_804_283_468_7, 0.999_909_384_695_144),
        (0.125 * 0.000_062_482_559_240_744_09, 0.999_995_316_041_220_5),
        (0.125 * 1.826_332_059_371_066e-6, 0.999_999_892_781_612_4),
        (0.125 * 1.868_728_226_873_641e-8, 0.999_999_999_142_705_1),
        (0.125 * 4.937_853_877_663_192_6e-11, 0.999_999_999_998_232_2),
    ],
    &[
        //  5th layer weights and abscissa: transformed 1/16, 3/16, ...
        (0.0625 * 1.558_773_355_533_33, 0.097_923_885_287_832_33),
        (0.0625 * 1.466_014_426_716_965_7, 0.287_879_932_742_715_9),
        (0.0625 * 1.297_475_750_424_978, 0.461_253_543_939_585_73),
        (0.0625 * 1.081_634_985_490_070_4, 0.610_273_657_500_639),
        (0.0625 * 0.850_172_856_456_620_1, 0.731_018_034_792_561_6),
        (0.0625 * 0.630_405_135_164_743_7, 0.823_317_005_506_402_4),
        (0.0625 * 0.440_833_236_273_858_23, 0.889_891_402_784_260_2),
        (0.0625 * 0.290_240_679_312_454_2, 0.935_160_857_521_984_6),
        (0.0625 * 0.179_324_412_110_728_3, 0.964_112_164_223_547_3),
        (0.0625 * 0.103_432_154_223_332_9, 0.981_454_826_677_335_2),
        (0.0625 * 0.055_289_683_742_240_58, 0.991_126_992_441_698_8),
        (0.0625 * 0.027_133_510_013_712_003, 0.996_108_665_437_508_5),
        (0.0625 * 0.012_083_543_599_157_953, 0.998_454_208_767_697_7),
        (0.0625 * 0.004_816_298_143_928_463, 0.999_451_434_435_274_6),
        (0.0625 * 0.001_690_873_998_142_639_6, 0.999_828_822_072_874_9),
        (0.0625 * 0.000_513_393_824_067_903_3, 0.999_953_871_005_627_9),
        (0.0625 * 0.000_132_052_341_256_099_76, 0.999_989_482_014_818_5),
        (0.0625 * 0.000_028_110_164_327_940_134, 0.999_998_017_140_595_4),
        (0.0625 * 4.823_718_203_261_55e-6, 0.999_999_698_894_152_6),
        (0.0625 * 6.477_756_603_592_972e-7, 0.999_999_964_239_080_9),
        (0.0625 * 6.583_518_512_718_34e-8, 0.999_999_996_787_199_1),
        (0.0625 * 4.876_006_097_424_062e-9, 0.999_999_999_789_732_9),
        (0.0625 * 2.521_634_791_853_014_7e-10, 0.999_999_999_990_393_9),
        (0.0625 * 8.675_931_414_979_604e-12, 0.999_999_999_999_708_1),
    ],
    &[
        // 6th layer weights and abcissas: transformed 1/32, 3/32, ...
        (0.03125 * 1.567_781_431_307_221_8, 0.049_055_967_305_077_885),
        (0.03125 * 1.543_881_116_176_959_2, 0.146_417_984_290_587_94),
        (0.03125 * 1.497_226_222_541_036_2, 0.241_566_319_538_883_66),
        (0.03125 * 1.430_008_354_872_299_7, 0.333_142_264_577_638_07),
        (0.03125 * 1.345_278_884_766_251_6, 0.419_952_111_278_447_17),
        (0.03125 * 1.246_701_207_451_857_8, 0.501_013_389_379_309_1),
        (0.03125 * 1.138_272_243_376_305_3, 0.575_584_490_635_151_7),
        (0.03125 * 1.024_044_933_111_811_6, 0.643_176_758_985_204_7),
        (0.03125 * 0.907_879_379_154_895_4, 0.703_550_005_147_142),
        (0.03125 * 0.793_242_700_820_516_7, 0.756_693_908_633_73),
        (0.03125 * 0.683_068_516_344_263_8, 0.802_798_741_343_241_3),
        (0.03125 * 0.579_678_103_087_787_7, 0.842_219_246_350_756_8),
        (0.03125 * 0.484_758_091_214_755_4, 0.875_435_397_630_408_7),
        (0.03125 * 0.399_384_741_525_717_1, 0.903_013_281_513_573_9),
        (0.03125 * 0.324_082_539_611_528_9, 0.925_568_634_068_612_7),
        (0.03125 * 0.258_904_639_514_053_5, 0.943_734_786_052_757_2),
        (0.03125 * 0.203_523_998_858_601_76, 0.958_136_022_710_213_7),
        (0.03125 * 0.157_326_203_484_366_16, 0.969_366_732_896_917_3),
        (0.03125 * 0.119_497_411_288_695_93, 0.977_976_235_186_665),
        (0.03125 * 0.089_103_139_240_941_46, 0.984_458_831_167_430_8),
        (0.03125 * 0.065_155_533_432_536_2, 0.989_248_431_090_133_9),
        (0.03125 * 0.046_668_208_054_846_616, 0.992_716_997_196_827_3),
        (0.03125 * 0.032_698_732_726_609_03, 0.995_176_026_155_327_4),
        (0.03125 * 0.022_379_471_063_648_477, 0.996_880_318_128_191_9),
        (0.03125 * 0.014_937_835_096_050_13, 0.998_033_336_315_433_8),
        (0.03125 * 0.009_707_223_739_391_69, 0.998_793_534_298_805_9),
        (0.03125 * 0.006_130_037_632_083_030_5, 0.999_281_111_921_791_9),
        (0.03125 * 0.003_754_250_977_431_834_5, 0.999_584_750_351_517_6),
        (0.03125 * 0.002_225_082_706_478_642_7, 0.999_767_971_599_560_9),
        (0.03125 * 0.001_273_327_944_708_238_2, 0.999_874_865_048_780_3),
        (0.03125 * 0.000_701_859_515_684_242_3, 0.999_935_019_925_082_4),
        (0.03125 * 0.000_371_666_936_216_777_6, 0.999_967_593_067_943_5),
        (0.03125 * 0.000_188_564_429_767_003_2, 0.999_984_519_902_270_8),
        (0.03125 * 0.000_091_390_817_490_710_13, 0.999_992_937_876_662_9),
        (0.03125 * 0.000_042_183_183_841_757_6, 0.999_996_932_449_190_4),
        (0.03125 * 0.000_018_481_813_599_879_218, 0.999_998_735_471_865_9),
        (0.03125 * 7.659_575_852_520_317e-6, 0.999_999_507_005_719_5),
        (0.03125 * 2.991_661_587_813_878_6e-6, 0.999_999_818_893_712_8),
        (0.03125 * 1.096_883_512_590_126_5e-6, 0.999_999_937_554_078_3),
        (0.03125 * 3.759_541_186_236_063e-7, 0.999_999_979_874_503_2),
        (0.03125 * 1.199_244_278_290_277e-7, 0.999_999_993_964_134_2),
        (0.03125 * 3.543_477_717_142_195e-8, 0.999_999_998_323_362),
        (0.03125 * 9.649_888_896_108_964e-9, 0.999_999_999_570_787_8),
        (0.03125 * 2.409_177_325_647_594e-9, 0.999_999_999_899_277_7),
        (0.03125 * 5.482_835_779_709_498e-10, 0.999_999_999_978_455_3),
        (0.03125 * 1.130_605_534_749_468e-10, 0.999_999_999_995_824_6),
        (0.03125 * 2.098_933_540_451_147e-11, 0.999_999_999_999_271_5),
        (0.03125 * 3.484_193_767_026_105_8e-12, 0.999_999_999_999_886_3),
    ],
    &[
        // 7th layer weights and abcissas: transformed 1/64, 3/64, ...
        (0.015625 * 1.570_042_029_279_593_1, 0.024_539_763_574_649_16),
        (0.015625 * 1.564_021_403_773_232, 0.073_525_122_985_671_29),
        (0.015625 * 1.552_053_169_845_412, 0.122_229_122_201_557_64),
        (0.015625 * 1.534_281_738_154_303_5, 0.170_467_972_382_010_53),
        (0.015625 * 1.510_919_723_074_169_8, 0.218_063_473_469_712),
        (0.015625 * 1.482_243_297_885_538, 0.264_845_076_583_447_97),
        (0.015625 * 1.448_586_254_961_322_7, 0.310_651_780_552_846),
        (0.015625 * 1.410_332_971_446_259, 0.355_333_825_165_074_56),
        (0.015625 * 1.367_910_511_680_896_5, 0.398_754_150_467_237_8),
        (0.015625 * 1.321_780_117_443_772_9, 0.440_789_599_033_900_86),
        (0.015625 * 1.272_428_345_537_862_7, 0.481_331_846_116_905_05),
        (0.015625 * 1.220_358_109_579_358_3, 0.520_288_050_691_230_2),
        (0.015625 * 1.166_079_869_932_434_6, 0.557_581_228_260_778_3),
        (0.015625 * 1.110_103_193_965_340_3, 0.593_150_353_591_953_1),
        (0.015625 * 1.052_928_879_955_266_7, 0.626_950_208_051_042_8),
        (0.015625 * 0.995_041_804_046_132_7, 0.658_950_991_743_350_1),
        (0.015625 * 0.936_904_612_745_667_9, 0.689_137_725_061_667_7),
        (0.015625 * 0.878_952_345_552_782_1, 0.717_509_467_487_324_1),
        (0.015625 * 0.821_588_035_266_964_7, 0.744_078_383_547_347_3),
        (0.015625 * 0.765_179_298_908_956_1, 0.768_868_686_768_246_6),
        (0.015625 * 0.710_055_901_205_469, 0.791_915_492_376_142_1),
        (0.015625 * 0.656_508_246_131_627_5, 0.813_263_608_502_973_9),
        (0.015625 * 0.604_786_730_578_403_6, 0.832_966_293_919_410_9),
        (0.015625 * 0.555_101_878_003_633_5, 0.851_084_007_987_848_8),
        (0.015625 * 0.507_625_158_831_908_1, 0.867_683_175_775_646),
        (0.015625 * 0.462_490_398_055_367_74, 0.882_834_988_244_669),
        (0.015625 * 0.419_795_668_445_015_5, 0.896_614_254_280_076),
        (0.015625 * 0.379_605_569_386_651_63, 0.909_098_318_163_020_4),
        (0.015625 * 0.341_953_795_923_016_83, 0.920_366_053_031_952_8),
        (0.015625 * 0.306_845_909_417_916_95, 0.930_496_937_997_153_4),
        (0.015625 * 0.274_262_229_689_068_1, 0.939_570_223_933_274_7),
        (0.015625 * 0.244_160_777_869_839_92, 0.947_664_190_615_153_1),
        (0.015625 * 0.216_480_209_117_296_18, 0.954_855_495_805_022_7),
        (0.015625 * 0.191_142_684_133_427_5, 0.961_218_615_151_116_4),
        (0.015625 * 0.168_056_637_948_269_17, 0.966_825_370_312_355_8),
        (0.015625 * 0.147_119_413_257_856_93, 0.971_744_541_565_487_3),
        (0.015625 * 0.128_219_733_631_200_98, 0.976_041_560_256_576_7),
        (0.015625 * 0.111_239_998_988_744_52, 0.979_778_275_800_615_7),
        (0.015625 * 0.096_058_391_865_189_47, 0.983_012_791_481_101_1),
        (0.015625 * 0.082_550_788_110_701_74, 0.985_799_363_025_283_5),
        (0.015625 * 0.070_592_469_906_867, 0.988_188_353_800_742_7),
        (0.015625 * 0.060_059_642_358_636_3, 0.990_226_240_467_527_7),
        (0.015625 * 0.050_830_757_572_570_474, 0.991_955_663_002_677_6),
        (0.015625 * 0.042_787_652_157_725_675, 0.993_415_513_169_264),
        (0.015625 * 0.035_816_505_604_196_434, 0.994_641_055_712_511_2),
        (0.015625 * 0.029_808_628_117_310_127, 0.995_664_076_816_953_1),
        (0.015625 * 0.024_661_087_314_753_284, 0.996_513_054_640_253_7),
        (0.015625 * 0.020_277_183_817_500_124, 0.997_213_347_043_468_8),
        (0.015625 * 0.016_566_786_254_247_574, 0.997_787_391_958_906_5),
        (0.015625 * 0.013_446_536_605_285_732, 0.998_254_916_171_996_2),
        (0.015625 * 0.010_839_937_168_255_907, 0.998_633_148_640_677_4),
        (0.015625 * 0.008_677_330_749_539_181, 0.998_937_034_833_512_1),
        (0.015625 * 0.006_895_785_969_066_003, 0.999_179_448_934_886),
        (0.015625 * 0.005_438_899_797_623_999, 0.999_371_401_140_937_7),
        (0.015625 * 0.004_256_529_599_017_858, 0.999_522_237_651_217_2),
        (0.015625 * 0.003_304_466_994_034_830_4, 0.999_639_831_345_600_4),
        (0.015625 * 0.002_544_065_767_529_173, 0.999_730_761_519_808_4),
        (0.015625 * 0.001_941_835_775_984_367_5, 0.999_800_481_431_138_4),
        (0.015625 * 0.001_469_014_359_942_979_1, 0.999_853_472_773_111_4),
        (0.015625 * 0.001_101_126_113_451_938_4, 0.999_893_386_547_592_5),
        (0.015625 * 0.000_817_541_013_324_694_9, 0.999_923_170_129_289_3),
        (0.015625 * 0.000_601_039_879_911_474_2, 0.999_945_180_614_458_7),
        (0.015625 * 0.000_437_394_956_159_116_86, 0.999_961_284_807_856_6),
        (0.015625 * 0.000_314_972_091_860_212, 0.999_972_946_425_232_3),
        (0.015625 * 0.000_224_359_652_050_085_5, 0.999_981_301_270_120_7),
        (0.015625 * 0.000_158_027_884_007_011_92, 0.999_987_221_282_000_7),
        (0.015625 * 0.000_110_021_128_466_666_97, 0.999_991_368_448_344_9),
        (0.015625 * 0.000_075_683_996_586_201_48, 0.999_994_239_627_616_7),
        (0.015625 * 0.000_051_421_497_447_658_804, 0.999_996_203_347_166_2),
        (0.015625 * 0.000_034_492_124_759_343_2, 0.999_997_529_623_805_2),
        (0.015625 * 0.000_022_832_118_109_036_146, 0.999_998_413_810_964_8),
        (0.015625 * 0.000_014_908_514_031_870_609, 0.999_998_995_410_689_9),
        (0.015625 * 9.598_194_128_378_471e-6, 0.999_999_372_707_335_4),
        (0.015625 * 6.089_910_032_094_904e-6, 0.999_999_613_988_550_2),
        (0.015625 * 3.806_198_326_464_489_7e-6, 0.999_999_766_023_332_4),
        (0.015625 * 2.342_166_720_852_809_5e-6, 0.999_999_860_371_214_6),
        (0.015625 * 1.418_306_715_549_391_7e-6, 0.999_999_918_004_794_7),
        (0.015625 * 8.447_375_638_485_986e-7, 0.999_999_952_642_664_5),
        (0.015625 * 4.945_828_870_275_42e-7, 0.999_999_973_113_236),
        (0.015625 * 2.844_992_365_915_980_6e-7, 0.999_999_985_003_076_3),
        (0.015625 * 1.606_939_457_907_622_5e-7, 0.999_999_991_786_456),
        (0.015625 * 8.907_139_514_024_239e-8, 0.999_999_995_585_633_6),
        (0.015625 * 4.842_095_019_807_237e-8, 0.999_999_997_673_236_8),
        (0.015625 * 2.579_956_822_953_589_4e-8, 0.999_999_998_797_983_5),
        (0.015625 * 1.346_464_552_230_203_8e-8, 0.999_999_999_391_776_9),
        (0.015625 * 6.878_461_095_589_9e-9, 0.999_999_999_698_754_4),
        (0.015625 * 3.437_185_674_465_009e-9, 0.999_999_999_854_056_1),
        (0.015625 * 1.678_889_768_216_190_6e-9, 0.999_999_999_930_888_4),
        (0.015625 * 8.009_978_447_972_966e-10, 0.999_999_999_968_033_2),
        (0.015625 * 3.729_950_184_305_279e-10, 0.999_999_999_985_568_8),
        (0.015625 * 1.693_945_778_941_164_8e-10, 0.999_999_999_993_646_3),
        (0.015625 * 7.496_739_757_381_822e-11, 0.999_999_999_997_274_1),
        (0.015625 * 3.230_446_433_325_236_6e-11, 0.999_999_999_998_861_2),
        (0.015625 * 1.354_251_291_233_627_5e-11, 0.999_999_999_999_537_3),
        (0.015625 * 5.518_236_946_817_489e-12, 0.999_999_999_999_817_1),
        (0.015625 * 2.183_592_209_923_360_7e-12, 0.999_999_999_999_929_8),
    ],
]; // end weights