Skip to content

Commit eda5c03

Browse files
lxvmlxvm
committed
strip units in handle_infinities
1 parent d672c1e commit eda5c03

File tree

4 files changed

+72
-50
lines changed

4 files changed

+72
-50
lines changed

‎src/QuadGK.jl‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -40,14 +40,14 @@ struct InplaceIntegrand{F,R,RI}
4040
Ig::R
4141
Ik::R
4242
fx::R
43-
Idiff::RI
43+
Idiff::R
4444

4545
# final result array
4646
I::RI
4747
end
4848

4949
InplaceIntegrand(f!::F, I::RI, fx::R) where {F,RI,R} =
50-
InplaceIntegrand{F,R,RI}(f!, similar(fx), similar(fx), similar(fx), similar(fx), fx, similar(I), I)
50+
InplaceIntegrand{F,R,RI}(f!, similar(fx), similar(fx), similar(fx), similar(fx), fx, similar(fx), I)
5151

5252
include("gausskronrod.jl")
5353
include("evalrule.jl")

‎src/adapt.jl‎

Lines changed: 46 additions & 29 deletions
Original file line numberDiff line numberDiff line change
@@ -95,7 +95,7 @@ end
9595
# re-sum (paranoia about accumulated roundoff)
9696
function resum(f, segs)
9797
if f isa InplaceIntegrand
98-
I = f.I .= segs[1].I
98+
I = f.Ik .= segs[1].I
9999
E = segs[1].E
100100
for i in 2:length(segs)
101101
I .+= segs[i].I
@@ -119,61 +119,77 @@ realone(x::Number) = one(x) isa Real
119119
# and pass transformed data to workfunc(f, s, tfunc)
120120
function handle_infinities(workfunc, f, s)
121121
s1, s2 = s[1], s[end]
122+
u = float(real(oneunit(s1))) # the units of the segment
122123
if realone(s1) && realone(s2) # check for infinite or semi-infinite intervals
123124
inf1, inf2 = isinf(s1), isinf(s2)
124125
if inf1 || inf2
125126
if inf1 && inf2 # x = t/(1-t^2) coordinate transformation
126-
return workfunc(u -> begin t = u/oneunit(u); t2 = t*t; den = 1 / (1 - t2);
127-
f(u*den) * (1+t2)*den*den; end,
128-
map(x -> isinf(x) ? (signbit(x) ? -oneunit(x) : oneunit(x)) : 2x / (one(x)+hypot(one(x),2x/oneunit(x))), s),
129-
t -> oneunit(s1) * t / (1 - t^2))
127+
I, E = workfunc(t -> begin t2 = t*t; den = 1 / (1 - t2);
128+
f(u*t*den) * (1+t2)*den*den; end,
129+
map(x -> isinf(x) ? (signbit(x) ? -one(x) : one(x)) : 2x / (oneunit(x)+hypot(oneunit(x),2x)), s),
130+
t -> u * t / (1 - t^2))
131+
return u * I, u * E
130132
end
131133
let (s0,si) = inf1 ? (s2,s1) : (s1,s2) # let is needed for JuliaLang/julia#15276
132134
if si < zero(si) # x = s0 - t/(1-t)
133-
return workfunc(u -> begin t = u/oneunit(u); den = 1 / (1 - t);
134-
f(s0 - u*den) * den*den; end,
135-
reverse(map(x -> oneunit(x) / (1 + oneunit(x) / (s0 - x)), s)),
136-
t -> s0 - oneunit(s1)*t/(1-t))
135+
I, E = workfunc(t -> begin den = 1 / (1 - t);
136+
f(s0 - u*t*den) * den*den; end,
137+
reverse(map(x -> 1 / (1 + oneunit(x) / (s0 - x)), s)),
138+
t -> s0 - u*t/(1-t))
139+
return u * I, u * E
137140
else # x = s0 + t/(1-t)
138-
return workfunc(u -> begin t = u/oneunit(u); den = 1 / (1 - t);
139-
f(s0 + u*den) * den*den; end,
140-
map(x -> oneunit(x) / (1 + oneunit(x) / (x - s0)), s),
141-
t -> s0 + oneunit(s1)*t/(1-t))
141+
I, E = workfunc(t -> begin den = 1 / (1 - t);
142+
f(s0 + u*t*den) * den*den; end,
143+
map(x -> 1 / (1 + oneunit(x) / (x - s0)), s),
144+
t -> s0 + u*t/(1-t))
145+
return u * I, u * E
142146
end
143147
end
144148
end
145149
end
146-
return workfunc(f, s, identity)
150+
I, E = workfunc(f, map(x -> x/oneunit(x), s), identity)
151+
return u * I, u * E
147152
end
148153

149154
function handle_infinities(workfunc, f::InplaceIntegrand, s)
150155
s1, s2 = s[1], s[end]
156+
result = f.I
157+
u = float(real(oneunit(s1))) # the units of the segment
151158
if realone(s1) && realone(s2) # check for infinite or semi-infinite intervals
152159
inf1, inf2 = isinf(s1), isinf(s2)
153160
if inf1 || inf2
154-
ftmp = f.fx # original integrand may have different type
155161
if inf1 && inf2 # x = t/(1-t^2) coordinate transformation
156-
return workfunc(InplaceIntegrand((v, u) -> begin t = u/oneunit(u); t2 = t*t; den = 1 / (1 - t2);
157-
f.f!(ftmp, u*den); v .= ftmp .* ((1+t2)*den*den); end, f.I, f.fx * float(one(s1))),
158-
map(x -> isinf(x) ? (signbit(x) ? -oneunit(x) : oneunit(x)) : 2x / (one(x)+hypot(one(x),2x/oneunit(x))), s),
159-
t -> oneunit(s1) * t / (1 - t^2))
162+
I, E = workfunc(InplaceIntegrand((v, t) -> begin t2 = t*t; den = 1 / (1 - t2);
163+
f.f!(v, u*t*den); v .*= ((1+t2)*den*den); end, f.fg, f.fk, f.Ig, f.Ik, f.fx, f.Idiff, similar(f.fx)),
164+
map(x -> isinf(x) ? (signbit(x) ? -one(x) : one(x)) : 2x / (oneunit(x)+hypot(oneunit(x),2x)), s),
165+
t -> u * t / (1 - t^2))
166+
result .= u .* I
167+
return result, u * E
160168
end
161169
let (s0,si) = inf1 ? (s2,s1) : (s1,s2) # let is needed for JuliaLang/julia#15276
162170
if si < zero(si) # x = s0 - t/(1-t)
163-
return workfunc(InplaceIntegrand((v, u) -> begin t = u/oneunit(u); den = 1 / (1 - t);
164-
f.f!(ftmp, s0 - u*den); v .= ftmp .* (den * den); end, f.I, f.fx * float(one(s1))),
165-
reverse(map(x -> oneunit(x) / (1 + oneunit(x) / (s0 - x)), s)),
166-
t -> s0 - oneunit(s1)*t/(1-t))
171+
I, E = workfunc(InplaceIntegrand((v, t) -> begin den = 1 / (1 - t);
172+
f.f!(v, s0 - u*t*den); v .*= (den * den); end, f.fg, f.fk, f.Ig, f.Ik, f.fx, f.Idiff, similar(f.fx)),
173+
reverse(map(x -> 1 / (1 + oneunit(x) / (s0 - x)), s)),
174+
t -> s0 - u*t/(1-t))
175+
result .= u .* I
176+
return result, u * E
167177
else # x = s0 + t/(1-t)
168-
return workfunc(InplaceIntegrand((v, u) -> begin t = u/oneunit(u); den = 1 / (1 - t);
169-
f.f!(ftmp, s0 + u*den); v .= ftmp .* (den * den); end, f.I, f.fx * float(one(s1))),
170-
map(x -> oneunit(x) / (1 + oneunit(x) / (x - s0)), s),
171-
t -> s0 + oneunit(s1)*t/(1-t))
178+
I, E = workfunc(InplaceIntegrand((v, t) -> begin den = 1 / (1 - t);
179+
f.f!(v, s0 + u*t*den); v .*= (den * den); end, f.fg, f.fk, f.Ig, f.Ik, f.fx, f.Idiff, similar(f.fx)),
180+
map(x -> 1 / (1 + oneunit(x) / (x - s0)), s),
181+
t -> s0 + u*t/(1-t))
182+
result .= u .* I
183+
return result, u * E
172184
end
173185
end
174186
end
175187
end
176-
return workfunc(f, s, identity)
188+
I, E = workfunc(InplaceIntegrand(f.f!, f.fg, f.fk, f.Ig, f.Ik, f.fx, f.Idiff, similar(f.fx)),
189+
map(x -> x/oneunit(x), s),
190+
identity)
191+
result .= u .* I
192+
return result, u * E
177193
end
178194

179195
function check_endpoint_roundoff(a, b, x; throw_error::Bool=false)
@@ -253,8 +269,9 @@ quadgk(f, segs...; kws...) =
253269

254270
function quadgk(f, segs::T...;
255271
atol=nothing, rtol=nothing, maxevals=10^7, order=7, norm=norm, segbuf=nothing) where {T}
272+
utol = isnothing(atol) ? atol : atol/float(real(oneunit(T))) # remove units of domain
256273
handle_infinities(f, segs) do f, s, _
257-
do_quadgk(f, s, order, atol, rtol, maxevals, norm, segbuf)
274+
do_quadgk(f, s, order, utol, rtol, maxevals, norm, segbuf)
258275
end
259276
end
260277

‎src/batch.jl‎

Lines changed: 22 additions & 17 deletions
Original file line numberDiff line numberDiff line change
@@ -144,35 +144,40 @@ end
144144

145145
function handle_infinities(workfunc, f::BatchIntegrand, s)
146146
s1, s2 = s[1], s[end]
147+
u = float(real(oneunit(s1))) # the units of the segment
148+
tbuf = similar(f.x, typeof(s1/oneunit(s1)))
147149
if realone(s1) && realone(s2) # check for infinite or semi-infinite intervals
148150
inf1, inf2 = isinf(s1), isinf(s2)
149151
if inf1 || inf2
150-
xbuf = f.x # buffer to store evaluation points
151-
ytmp = f.y # original integrand may have different type
152-
xtmp = similar(xbuf, typeof(one(eltype(f.x))))
153-
ybuf = similar(ytmp, typeof(oneunit(eltype(f.y))*float(one(s1))))
154152
if inf1 && inf2 # x = t/(1-t^2) coordinate transformation
155-
return workfunc(BatchIntegrand((v, u) -> begin resize!(xtmp, length(u)); resize!(ytmp, length(v)); t = xtmp .= u ./ oneunit(s1);
156-
f.f!(ytmp, u .= oneunit(s1) .* t ./ (1 .- t .* t)); v .= ytmp .* (1 .+ t .* t) ./ (1 .- t .* t) .^ 2; end, ybuf, xbuf, f.max_batch),
157-
map(x -> isinf(x) ? (signbit(x) ? -oneunit(x) : oneunit(x)) : 2x / (one(x)+hypot(one(x),2x/oneunit(x))), s),
158-
t -> oneunit(s1) * t / (1 - t^2))
153+
I, E = workfunc(BatchIntegrand((v, t) -> begin resize!(f.x, length(t));
154+
f.f!(v, f.x .= u .* t ./ (1 .- t .* t)); v .*= (1 .+ t .* t) ./ (1 .- t .* t) .^ 2; end, f.y, tbuf, f.max_batch),
155+
map(x -> isinf(x) ? (signbit(x) ? -one(x) : one(x)) : 2x / (oneunit(x)+hypot(oneunit(x),2x)), s),
156+
t -> u * t / (1 - t^2))
157+
return u * I, u * E
159158
end
160159
let (s0,si) = inf1 ? (s2,s1) : (s1,s2) # let is needed for JuliaLang/julia#15276
161160
if si < zero(si) # x = s0 - t/(1-t)
162-
return workfunc(BatchIntegrand((v, u) -> begin resize!(xtmp, length(u)); resize!(ytmp, length(v)); t = xtmp .= u ./ oneunit(s1);
163-
f.f!(ytmp, u .= s0 .- oneunit(s1) .* t ./ (1 .- t)); v .= ytmp ./ (1 .- t) .^ 2; end, ybuf, xbuf, f.max_batch),
164-
reverse(map(x -> oneunit(x) / (1 + oneunit(x) / (s0 - x)), s)),
165-
t -> s0 - oneunit(s1)*t/(1-t))
161+
I, E = workfunc(BatchIntegrand((v, t) -> begin resize!(f.x, length(t));
162+
f.f!(v, f.x .= s0 .- u .* t ./ (1 .- t)); v ./= (1 .- t) .^ 2; end, f.y, tbuf, f.max_batch),
163+
reverse(map(x -> 1 / (1 + oneunit(x) / (s0 - x)), s)),
164+
t -> s0 - u*t/(1-t))
165+
return u * I, u * E
166166
else # x = s0 + t/(1-t)
167-
return workfunc(BatchIntegrand((v, u) -> begin resize!(xtmp, length(u)); resize!(ytmp, length(v)); t = xtmp .= u ./ oneunit(s1);
168-
f.f!(ytmp, u .= s0 .+ oneunit(s1) .* t ./ (1 .- t)); v .= ytmp ./ (1 .- t) .^ 2; end, ybuf, xbuf, f.max_batch),
169-
map(x -> oneunit(x) / (1 + oneunit(x) / (x - s0)), s),
170-
t -> s0 + oneunit(s1)*t/(1-t))
167+
I, E = workfunc(BatchIntegrand((v, t) -> begin resize!(f.x, length(t));
168+
f.f!(v, f.x .= s0 .+ u .* t ./ (1 .- t)); v ./= (1 .- t) .^ 2; end, f.y, tbuf, f.max_batch),
169+
map(x -> 1 / (1 + oneunit(x) / (x - s0)), s),
170+
t -> s0 + u*t/(1-t))
171+
return u * I, u * E
171172
end
172173
end
173174
end
174175
end
175-
return workfunc(f, s, identity)
176+
I, E = workfunc(BatchIntegrand((y, t) -> begin resize!(f.x, length(t));
177+
f.f!(y, f.x .= u .* t); end, f.y, tbuf, f.max_batch),
178+
map(x -> x/oneunit(x), s),
179+
identity)
180+
return u * I, u * E
176181
end
177182

178183
"""

‎test/runtests.jl‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -70,12 +70,12 @@ module Test19626
7070
for (lb, ub) in ulims
7171
## function
7272
f = x -> 1/(1+(x/oneunit(x))^2)
73-
buf = QuadGK.alloc_segbuf(MockQuantity, MockQuantity, MockQuantity)
73+
buf = QuadGK.alloc_segbuf(MockQuantity, Float64, MockQuantity)
7474
@test QuadGK.quadgk(f, lb, ub, atol=MockQuantity(0.0))[1] ≈
7575
QuadGK.quadgk(f, lb, ub, atol=MockQuantity(0.0), segbuf=buf)[1]
7676
## inplace
7777
fiip = (y, x) -> y[1] = 1/(1+(x/oneunit(x))^2)
78-
ibuf = QuadGK.alloc_segbuf(MockQuantity, Array{MockQuantity,1}, MockQuantity)
78+
ibuf = QuadGK.alloc_segbuf(MockQuantity, Array{Float64,1}, MockQuantity)
7979
@test QuadGK.quadgk!(fiip, [MockQuantity(0.0)], lb, ub, atol=MockQuantity(0.0), norm=abs∘first)[1][] ≈
8080
QuadGK.quadgk!(fiip, [MockQuantity(0.0)], lb, ub, atol=MockQuantity(0.0), norm=abs∘getindex, segbuf=ibuf)[1][]
8181
## batch

0 commit comments

Comments
 (0)