Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Type stability #21

Closed
wants to merge 5 commits into from
Closed

Type stability #21

Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
68f31e3
fix struct member
matbesancon Jun 2, 2018
5e64cf6
start fixing type instability
matbesancon Jun 3, 2018
d13fcb1
heapify segs
matbesancon Jun 3, 2018
d9070e4
fixed heapify
matbesancon Jun 4, 2018
0f54b7d
inplace heapify
matbesancon Jun 6, 2018
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
98 changes: 53 additions & 45 deletions src/QuadGK.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -14,7 +14,7 @@ end

using Compat
using DataStructures, Compat.LinearAlgebra
import Base: isless, Order.Reverse, AnyDict
import Base: isless, Order.Reverse

# Adaptive Gauss-Kronrod quadrature routines (arbitrary precision),
# written and contributed to Julia by Steven G. Johnson, 2013.
Expand All @@ -24,7 +24,7 @@ import Base: isless, Order.Reverse, AnyDict
# cache of (T,n) -> (x,w,gw) Kronrod rules, to avoid recomputing them
# unnecessarily for repeated integration. We initialize it with the
# default n=7 rule for double-precision calculations.
const rulecache = AnyDict( (Float64,7) => # precomputed in 100-bit arith.
const rulecache = Dict( (Float64,precision(Float64),7) => # precomputed in 100-bit arith.
([-9.9145537112081263920685469752598e-01,
-9.4910791234275852452618968404809e-01,
-8.6486442335976907278971278864098e-01,
Expand All @@ -47,11 +47,11 @@ const rulecache = AnyDict( (Float64,7) => # precomputed in 100-bit arith.
4.1795918367346938775510204081658e-01]) )

# integration segment (a,b), estimated integral I, and estimated error E
struct Segment
a::Number
b::Number
I
E
struct Segment{TB<:Number,TI,TE}
a::TB
b::TB
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Note that this is a different constraint: currently we're saying that the types of a and b are both subtypes of Number, but here you're saying they have to be the same subtype of Number. Is this intentional?

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

yes indeed, I considered the boundaries of the domain would be of the same type, so this was intentional. Would there be some use cases where this is not the case?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

No, not necessarily, I just wanted to make sure the change was intentional.

I::TI
E::TE
end
isless(i::Segment, j::Segment) = isless(i.E, j.E)

Expand All @@ -62,7 +62,7 @@ function evalrule(f, a,b, x,w,gw, nrm)
# Ik and Ig are integrals via Kronrod and Gauss rules, respectively
s = convert(eltype(x), 0.5) * (b-a)
n1 = 1 - (length(x) & 1) # 0 if even order, 1 if odd order
# unroll first iterationof loop to get correct type of Ik and Ig
# unroll first iteration of loop to get correct type of Ik and Ig
fg = f(a + (1+x[2])*s) + f(a + (1-x[2])*s)
fk = f(a + (1+x[1])*s) + f(a + (1-x[1])*s)
Ig = fg * gw[1]
Expand Down Expand Up @@ -91,45 +91,28 @@ function evalrule(f, a,b, x,w,gw, nrm)
end

rulekey(::Type{BigFloat}, n) = (BigFloat, precision(BigFloat), n)
rulekey(T,n) = (T,n)
rulekey(::Type{T},n) where {T<:AbstractFloat} = (T,precision(T),n)


function do_quadgk(f, s::TS, n, ::Type{Tw}, abstol, reltol, maxevals, nrm) where {Tw,R<:Real,TS<:AbstractArray{R}}
preprocess_real(f,s,n,Tw,abstol, reltol, maxevals, nrm)
end

function do_quadgk(f, s::TS, n, ::Type{Tw}, abstol, reltol, maxevals, nrm) where {Tw,TS<:AbstractArray}
do_quadgk_common(f, s, n, Tw, abstol, reltol, maxevals, nrm)
end

# Internal routine: integrate f over the union of the open intervals
# (s[1],s[2]), (s[2],s[3]), ..., (s[end-1],s[end]), using h-adaptive
# integration with the order-n Kronrod rule and weights of type Tw,
# with absolute tolerance abstol and relative tolerance reltol,
# with maxevals an approximate maximum number of f evaluations.
function do_quadgk(f, s, n, ::Type{Tw}, abstol, reltol, maxevals, nrm) where Tw
if eltype(s) <: Real # check for infinite or semi-infinite intervals
s1 = s[1]; s2 = s[end]; inf1 = isinf(s1); inf2 = isinf(s2)
if inf1 || inf2
if inf1 && inf2 # x = t/(1-t^2) coordinate transformation
return do_quadgk(t -> begin t2 = t*t; den = 1 / (1 - t2);
f(t*den) * (1+t2)*den*den; end,
map(x -> isinf(x) ? copysign(one(x), x) : 2x / (1+hypot(1,2x)), s),
n, Tw, abstol, reltol, maxevals, nrm)
end
s0,si = inf1 ? (s2,s1) : (s1,s2)
if si < 0 # x = s0 - t/(1-t)
return do_quadgk(t -> begin den = 1 / (1 - t);
f(s0 - t*den) * den*den; end,
reverse!(map(x -> 1 / (1 + 1 / (s0 - x)), s)),
n, Tw, abstol, reltol, maxevals, nrm)
else # x = s0 + t/(1-t)
return do_quadgk(t -> begin den = 1 / (1 - t);
f(s0 + t*den) * den*den; end,
map(x -> 1 / (1 + 1 / (x - s0)), s),
n, Tw, abstol, reltol, maxevals, nrm)
end
end
end
function do_quadgk_common(f, s::TS, n, ::Type{Tw}, abstol, reltol, maxevals, nrm) where {Tw,TS<:AbstractArray}

key = rulekey(Tw,n)
x,w,gw = haskey(rulecache, key) ? rulecache[key] :
(rulecache[key] = kronrod(Tw, n))
segs = Segment[]
for i in 1:length(s) - 1
heappush!(segs, evalrule(f, s[i],s[i+1], x,w,gw, nrm), Reverse)
end
segs = heapify!([evalrule(f, s[i], s[i+1], x, w, gw, nrm) for i in 1:length(s)-1], Reverse)
numevals = (2n+1) * length(segs)
I = segs[1].I
E = segs[1].E
Expand All @@ -140,14 +123,14 @@ function do_quadgk(f, s, n, ::Type{Tw}, abstol, reltol, maxevals, nrm) where Tw
# Pop the biggest-error segment and subdivide (h-adaptation)
# until convergence is achieved or maxevals is exceeded.
while E > abstol && E > reltol * nrm(I) && numevals < maxevals
s = heappop!(segs, Reverse)
mid = (s.a + s.b) * 0.5
s1 = evalrule(f, s.a, mid, x,w,gw, nrm)
s2 = evalrule(f, mid, s.b, x,w,gw, nrm)
heappush!(segs, s1, Reverse)
heappush!(segs, s2, Reverse)
I = (I - s.I) + s1.I + s2.I
E = (E - s.E) + s1.E + s2.E
seg = heappop!(segs, Reverse)
mid = (seg.a + seg.b) * 0.5
seg1 = evalrule(f, seg.a, mid, x,w,gw, nrm)
seg2 = evalrule(f, mid, seg.b, x,w,gw, nrm)
heappush!(segs, seg1, Reverse)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Seems to need some kind of check here for seg1 isa eltype(segs) (similarly for seg2), and reallocate segs if not (will also require a refactored do_quadgk_common routine that can take an existing segs heap as a starting point).

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

this would mean all segments don't have the same type? Wouldn't that be strange, boundaries all have the same type?

heappush!(segs, seg2, Reverse)
I = (I - seg.I) + seg1.I + seg2.I
E = (E - seg.E) + seg1.E + seg2.E
numevals += 4n+2
end
# re-sum (paranoia about accumulated roundoff)
Expand All @@ -160,6 +143,31 @@ function do_quadgk(f, s, n, ::Type{Tw}, abstol, reltol, maxevals, nrm) where Tw
return (I, E)
end

function preprocess_real(f, s::TS, n, ::Type{Tw}, abstol, reltol, maxevals, nrm) where {Tw,ET<:Real,TS<:AbstractVector{ET}}
s1 = s[1]; s2 = s[end]; inf1 = isinf(s1); inf2 = isinf(s2)
if inf1 || inf2
if inf1 && inf2 # x = t/(1-t^2) coordinate transformation
return do_quadgk_common(t -> begin t2 = t*t; den = 1 / (1 - t2);
f(t*den) * (1+t2)*den*den; end,
[isinf(x) ? copysign(one(x), x) : 2x / (1+hypot(1,2x)) for x in s],
n, Tw, abstol, reltol, maxevals, nrm)
end
(s0,si) = inf1 ? (s2,s1) : (s1,s2)
if si < 0 # x = s0 - t/(1-t)
return do_quadgk_common(t -> begin den = 1 / (1 - t);
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

You can use a do block rather than t -> begin.

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

it was written like this, I found it a bit dense, I'll rewrite this part

f(s0 - t*den) * den*den; end,
reverse!(map(x -> 1 / (1 + 1 / (s0 - x)), s)),
n, Tw, abstol, reltol, maxevals, nrm)
else # x = s0 + t/(1-t)
return do_quadgk_common(t -> begin den = 1 / (1 - t);
f(s0 + t*den) * den*den; end,
map(x -> 1 / (1 + 1 / (x - s0)), s),
n, Tw, abstol, reltol, maxevals, nrm)
end
end
do_quadgk_common(f, s, n, Tw, abstol, reltol, maxevals, nrm)
end

# Gauss-Kronrod quadrature of f from a to b to c...

function quadgk(f, a::T,b::T,c::T...;
Expand Down