avoid integrand evaluations at endpoints

lxvm commented 11 months ago

This pr concerns issues #86, #82, and #39, and addresses errors arising from evaluation of the integrand at the endpoints of the interval. This pr handles the endpoints by:

throwing a DomainError if the roundoff happens on an initial segment
silently returning the best integral and error estimate if the roundoff happens during refinement

I think it would be better to warn the user of the following situations in which quadgk silently returns a possibly unconverged result

endpoint roundoff error
maxevals reached

Should we do this and add return codes?

codecov[bot] commented 11 months ago

Codecov Report

Patch coverage: 100.00% and project coverage change: +0.04% :tada:

Comparison is base (10bc1a1) 98.16% compared to head (173eff3) 98.20%.

Additional details and impacted files

```diff @@ Coverage Diff @@ ## master #91 +/- ## ========================================== + Coverage 98.16% 98.20% +0.04% ========================================== Files 6 6 Lines 599 614 +15 ========================================== + Hits 588 603 +15 Misses 11 11 ``` | [Files Changed](https://app.codecov.io/gh/JuliaMath/QuadGK.jl/pull/91?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=JuliaMath) | Coverage Δ | | |---|---|---| | [src/adapt.jl](https://app.codecov.io/gh/JuliaMath/QuadGK.jl/pull/91?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=JuliaMath#diff-c3JjL2FkYXB0Lmps) | `96.99% <100.00%> (+0.32%)` | :arrow_up: | | [src/batch.jl](https://app.codecov.io/gh/JuliaMath/QuadGK.jl/pull/91?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=JuliaMath#diff-c3JjL2JhdGNoLmps) | `99.17% <100.00%> (+0.01%)` | :arrow_up: |

:umbrella: View full report in Codecov by Sentry.
:loudspeaker: Have feedback on the report? Share it here.

stevengj commented 11 months ago

Should we do this and add return codes?

The problem is that adding return codes to quadgk is breaking.

Of course, you can always tell from the error estimate — if the final error estimate is greater than the requested tolerance, you know that it must have returned early.

stevengj commented 11 months ago

This introduces a type-instability into evalrule. Does that affect performance?

lxvm commented 11 months ago

Of course, you can always tell from the error estimate — if the final error estimate is greater than the requested tolerance, you know that it must have returned early.

Right, that's a good point

This introduces a type-instability into evalrule. Does that affect performance?

Thanks for pointing this out, so I refactored a bit and put the type instability into the refine routine.

As for performance, the endpoint check has an overhead of about 3 ns, which should be negligible. Here is a benchmark:

master

julia> using BenchmarkTools

julia> using QuadGK

julia> @benchmark quadgk(x -> x, 0, 1) # zero-cost integrand, no subdivisions
BenchmarkTools.Trial: 10000 samples with 993 evaluations.
 Range (min … max):  34.597 ns … 133.580 ns  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     35.161 ns               ┊ GC (median):    0.00%
 Time  (mean ± σ):   36.762 ns ±   6.744 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

  █▅▆                                                          ▁
  ███▆▅▅▃▅▅▆▇▇█▆▆▆▄▆▆▆▆▆▆▆▅▆█▇▅▄▇█▅▅▄▆▅▅▄▆▅▄▄▄▅▅▄▁▅▅▅▅▆▆▅▅█▅▄▃ █
  34.6 ns       Histogram: log(frequency) by time      72.5 ns <

 Memory estimate: 0 bytes, allocs estimate: 0.

julia> @benchmark quadgk(x -> sin(100x), 0, 1) # low-cost integrand with subdivisions
BenchmarkTools.Trial: 10000 samples with 1 evaluation.
 Range (min … max):  11.684 μs … 63.263 μs  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     11.858 μs              ┊ GC (median):    0.00%
 Time  (mean ± σ):   12.611 μs ±  1.834 μs  ┊ GC (mean ± σ):  0.00% ± 0.00%

  ▆█▅▂             ▃▆▆▄▃▁             ▁▂▁                     ▂
  █████▆▅▄▁▁▁▁▁▁▄▁▁███████▇▆▆▇▇▇▇▅▆▆▄▅███▇▇▆▄▁▁▃▁▃▃▄▅▄▁▃▅▁▄▄▅ █
  11.7 μs      Histogram: log(frequency) by time      16.9 μs <

 Memory estimate: 1.69 KiB, allocs estimate: 3.

this pr

julia> using BenchmarkTools

julia> using QuadGK

julia> @benchmark quadgk(x -> x, 0, 1) # zero-cost integrand, no subdivisions
BenchmarkTools.Trial: 10000 samples with 992 evaluations.
 Range (min … max):  37.661 ns … 132.844 ns  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     37.885 ns               ┊ GC (median):    0.00%
 Time  (mean ± σ):   39.574 ns ±   7.290 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

  █  ▁                                                         ▁
  █▇▃█▃▃▃▄▅▄▄▅▅▅▄▂▃▄▂▂▄▆▆▅▅▄▅▅▅▄▅▄▅▅▄▄▄▅▅▅▅▄▄▃▅▅▅▆▆▇▆▆▆▅▅▄▄▄▄▄ █
  37.7 ns       Histogram: log(frequency) by time      75.2 ns <

 Memory estimate: 0 bytes, allocs estimate: 0.

julia> @benchmark quadgk(x -> sin(100x), 0, 1) # low-cost integrand with subdivisions
BenchmarkTools.Trial: 10000 samples with 1 evaluation.
 Range (min … max):  11.616 μs … 79.306 μs  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     11.986 μs              ┊ GC (median):    0.00%
 Time  (mean ± σ):   12.831 μs ±  2.304 μs  ┊ GC (mean ± σ):  0.00% ± 0.00%

  ▄██▄▃     ▂▆▇▄▃▁      ▂▂                                    ▂
  █████▇▁▃▃▄██████▇███▅▇███▇▆▆▅▆▆▄▄▅▄▃▄▅▆▇▇▅▆▃▄▄▄▁▄▁▄▄▄▃▁▄▃▁▄ █
  11.6 μs      Histogram: log(frequency) by time      20.5 μs <

 Memory estimate: 1.69 KiB, allocs estimate: 3.

JuliaMath / QuadGK.jl