Function evaluation for segments with zero extent

[szabo137]

Hi!

I would like to understand the following behavior:

julia> using QuadGK

julia> quadgk_count(x->x,0,0)
(0.0, 0.0, 15)

A 1D integral with the same upper and lower bound should always return zero. For quadgk, I assumed a shortcut that returns just zero without evaluating the integrand, but it is evaluated 15 times instead.

This is also the case for segmented integration:

julia> using QuadGK

julia> quadgk_print(x->x,0,0,1)
f(0.0) = 0.0
f(0.0) = 0.0
f(0.0) = 0.0
f(0.0) = 0.0
f(0.0) = 0.0
f(0.0) = 0.0
f(0.0) = 0.0
f(0.0) = 0.0
f(0.0) = 0.0
f(0.0) = 0.0
f(0.0) = 0.0
f(0.0) = 0.0
f(0.0) = 0.0
f(0.0) = 0.0
f(0.0) = 0.0
f(0.5) = 0.5
f(0.39610752249605075) = 0.39610752249605075
f(0.6038924775039493) = 0.6038924775039493
f(0.025446043828620757) = 0.025446043828620757
f(0.9745539561713792) = 0.9745539561713792
f(0.004272314439593694) = 0.004272314439593694
f(0.9957276855604063) = 0.9957276855604063
f(0.12923440720030277) = 0.12923440720030277
f(0.8707655927996972) = 0.8707655927996972
f(0.06756778832011545) = 0.06756778832011545
f(0.9324322116798845) = 0.9324322116798845
f(0.2970774243113014) = 0.2970774243113014
f(0.7029225756886985) = 0.7029225756886985
f(0.20695638226615443) = 0.20695638226615443
f(0.7930436177338456) = 0.7930436177338456
(0.5, 0.0, 30)

The integrand is evaluated 15 times in the first segment before going to the segment with a finite extent.

In my opinion, a shortcut for segments with no extent could help to save some evaluations.

[stevengj]

Right now this is not special-cased — it always puts N points in between the two endpoints, and doesn't check whether the endpoints are equal. Of course, it gives the correct answer because in this case the weight is zero.

It could be special cased to evaluate the function only once — you would still need one function evaluation in order to determine the return type (e.g. you could be integrating a function that returns a vector, not a scalar).

What is the use-case for this optimization?

[szabo137]

Hi @stevengj

Thank you for the reply. Since I embedded the call of quadgk in a larger system, I can not provide a full example. But the idea is, that I somehow calculate the positions of possible singularities within the boundaries of an integral. Then I integrate over the segments before and after the singularity. However, depening on parameters of the integral, the position of the singularity might be outside of the interval I want to integrate over. In this case, one of the segmented integrals is zero.

Here is some dummy code to illustrate this:

function expensive_integrand(x,s)
    # this has a singularity at s, 
    # which is not necessarily in (a,b)
    # but we could have s<=a. 
end

function something_to_calculate(s,a,b)
    crit_point = max(a,s)
    res,_ = quadgk(t->expensive_integrand(t,s),a,crit_point,b)
    return res
end

If I now call something_to_calculate(-1,0,1), the expensive_integrand is called 15 times on (0,-1). Depending on how expensive one call is, this adds some unnecessary overhead.

Sidenote: I know that I could handle this case by myself, but since this behaviour is not documented in QuadGK (which is obviously fine) it took me a while to realize, where the overhead in my code came from in the first place. Therefore, I thought I could be a good idea to have this special case in QuadGK.

[stevengj]

Should be a straighforward fix to add this special case. After this line

QuadGK.jl/src/evalrule.jl

Line 19 in 89a3b8f

s = convert(eltype(x), 0.5) * (b-a)

just add something like:

if iszero(s)
    I_s = f(a+s) * (w[end] * s) # zero of correct type
    return Segment(oftype(s, a), oftype(s, b), I_s, nrm(I_s))
end

and similarly for the other evalrule methods. Do you want to take a crack at a PR?

[stevengj]

PS. Note that it is better not to call quadgk separately for each segment, but instead to call it on all your segments at once, so that it can do the adaptivity globally.