From 82fbf9a5a9ef25cd9aa25cd5323e1e6c66da96fa Mon Sep 17 00:00:00 2001 From: lxvm Date: Mon, 13 Nov 2023 12:48:49 -0500 Subject: [PATCH 1/9] make handle_infinities consistently use unitful limits --- src/adapt.jl | 38 +++++++++++++++++++------------------- src/batch.jl | 26 +++++++++++++------------- 2 files changed, 32 insertions(+), 32 deletions(-) diff --git a/src/adapt.jl b/src/adapt.jl index 0a79d15..ec8cf1a 100644 --- a/src/adapt.jl +++ b/src/adapt.jl @@ -123,21 +123,21 @@ function handle_infinities(workfunc, f, s) inf1, inf2 = isinf(s1), isinf(s2) if inf1 || inf2 if inf1 && inf2 # x = t/(1-t^2) coordinate transformation - return workfunc(t -> begin t2 = t*t; den = 1 / (1 - t2); - f(oneunit(s1)*t*den) * (1+t2)*den*den*oneunit(s1); end, - map(x -> isinf(x) ? (signbit(x) ? -one(x) : one(x)) : 2x / (oneunit(x)+hypot(oneunit(x),2x)), s), + return workfunc(u -> begin t = u/oneunit(u); t2 = t*t; den = 1 / (1 - t2); + f(u*den) * (1+t2)*den*den; end, + map(x -> isinf(x) ? (signbit(x) ? -oneunit(x) : oneunit(x)) : 2x / (one(x)+hypot(one(x),2x/oneunit(x))), s), t -> oneunit(s1) * t / (1 - t^2)) end let (s0,si) = inf1 ? (s2,s1) : (s1,s2) # let is needed for JuliaLang/julia#15276 if si < zero(si) # x = s0 - t/(1-t) - return workfunc(t -> begin den = 1 / (1 - t); - f(s0 - oneunit(s1)*t*den) * den*den*oneunit(s1); end, - reverse(map(x -> 1 / (1 + oneunit(x) / (s0 - x)), s)), + return workfunc(u -> begin t = u/oneunit(u); den = 1 / (1 - t); + f(s0 - u*den) * den*den; end, + reverse(map(x -> oneunit(x) / (1 + oneunit(x) / (s0 - x)), s)), t -> s0 - oneunit(s1)*t/(1-t)) else # x = s0 + t/(1-t) - return workfunc(t -> begin den = 1 / (1 - t); - f(s0 + oneunit(s1)*t*den) * den*den*oneunit(s1); end, - map(x -> 1 / (1 + oneunit(x) / (x - s0)), s), + return workfunc(u -> begin t = u/oneunit(u); den = 1 / (1 - t); + f(s0 + u*den) * den*den; end, + map(x -> oneunit(x) / (1 + oneunit(x) / (x - s0)), s), t -> s0 + oneunit(s1)*t/(1-t)) end end @@ -151,23 +151,23 @@ function handle_infinities(workfunc, f::InplaceIntegrand, s) if realone(s1) && realone(s2) # check for infinite or semi-infinite intervals inf1, inf2 = isinf(s1), isinf(s2) if inf1 || inf2 - ftmp = f.fx # original integrand may have different units + ftmp = f.fx # original integrand may have different type if inf1 && inf2 # x = t/(1-t^2) coordinate transformation - return workfunc(InplaceIntegrand((v, t) -> begin t2 = t*t; den = 1 / (1 - t2); - f.f!(ftmp, oneunit(s1)*t*den); v .= ftmp .* ((1+t2)*den*den*oneunit(s1)); end, f.I, f.fx * oneunit(s1)), - map(x -> isinf(x) ? (signbit(x) ? -one(x) : one(x)) : 2x / (oneunit(x)+hypot(oneunit(x),2x)), s), + return workfunc(InplaceIntegrand((v, u) -> begin t = u/oneunit(u); t2 = t*t; den = 1 / (1 - t2); + f.f!(ftmp, u*den); v .= ftmp .* ((1+t2)*den*den); end, f.I, f.fx * float(one(s1))), + map(x -> isinf(x) ? (signbit(x) ? -oneunit(x) : oneunit(x)) : 2x / (one(x)+hypot(one(x),2x/oneunit(x))), s), t -> oneunit(s1) * t / (1 - t^2)) end let (s0,si) = inf1 ? (s2,s1) : (s1,s2) # let is needed for JuliaLang/julia#15276 if si < zero(si) # x = s0 - t/(1-t) - return workfunc(InplaceIntegrand((v, t) -> begin den = 1 / (1 - t); - f.f!(ftmp, s0 - oneunit(s1)*t*den); v .= ftmp .* (den * den * oneunit(s1)); end, f.I, f.fx * oneunit(s1)), - reverse(map(x -> 1 / (1 + oneunit(x) / (s0 - x)), s)), + return workfunc(InplaceIntegrand((v, u) -> begin t = u/oneunit(u); den = 1 / (1 - t); + f.f!(ftmp, s0 - u*den); v .= ftmp .* (den * den); end, f.I, f.fx * float(one(s1))), + reverse(map(x -> oneunit(x) / (1 + oneunit(x) / (s0 - x)), s)), t -> s0 - oneunit(s1)*t/(1-t)) else # x = s0 + t/(1-t) - return workfunc(InplaceIntegrand((v, t) -> begin den = 1 / (1 - t); - f.f!(ftmp, s0 + oneunit(s1)*t*den); v .= ftmp .* (den * den * oneunit(s1)); end, f.I, f.fx * oneunit(s1)), - map(x -> 1 / (1 + oneunit(x) / (x - s0)), s), + return workfunc(InplaceIntegrand((v, u) -> begin t = u/oneunit(u); den = 1 / (1 - t); + f.f!(ftmp, s0 + u*den); v .= ftmp .* (den * den); end, f.I, f.fx * float(one(s1))), + map(x -> oneunit(x) / (1 + oneunit(x) / (x - s0)), s), t -> s0 + oneunit(s1)*t/(1-t)) end end diff --git a/src/batch.jl b/src/batch.jl index d4221b7..fc53756 100644 --- a/src/batch.jl +++ b/src/batch.jl @@ -147,26 +147,26 @@ function handle_infinities(workfunc, f::BatchIntegrand, s) if realone(s1) && realone(s2) # check for infinite or semi-infinite intervals inf1, inf2 = isinf(s1), isinf(s2) if inf1 || inf2 - xtmp = f.x # buffer to store evaluation points - ytmp = f.y # original integrand may have different units - xbuf = similar(xtmp, typeof(one(eltype(f.x)))) - ybuf = similar(ytmp, typeof(oneunit(eltype(f.y))*oneunit(s1))) + xbuf = f.x # buffer to store evaluation points + ytmp = f.y # original integrand may have different type + xtmp = similar(xbuf, typeof(one(eltype(f.x)))) + ybuf = similar(ytmp, typeof(oneunit(eltype(f.y))*float(one(s1)))) if inf1 && inf2 # x = t/(1-t^2) coordinate transformation - return workfunc(BatchIntegrand((v, t) -> begin resize!(xtmp, length(t)); resize!(ytmp, length(v)); - f.f!(ytmp, xtmp .= oneunit(s1) .* t ./ (1 .- t .* t)); v .= ytmp .* (1 .+ t .* t) .* oneunit(s1) ./ (1 .- t .* t) .^ 2; end, ybuf, xbuf, f.max_batch), - map(x -> isinf(x) ? (signbit(x) ? -one(x) : one(x)) : 2x / (oneunit(x)+hypot(oneunit(x),2x)), s), + return workfunc(BatchIntegrand((v, u) -> begin resize!(xtmp, length(u)); resize!(ytmp, length(v)); t = xtmp .= u ./ oneunit(s1); + 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), + map(x -> isinf(x) ? (signbit(x) ? -oneunit(x) : oneunit(x)) : 2x / (one(x)+hypot(one(x),2x/oneunit(x))), s), t -> oneunit(s1) * t / (1 - t^2)) end let (s0,si) = inf1 ? (s2,s1) : (s1,s2) # let is needed for JuliaLang/julia#15276 if si < zero(si) # x = s0 - t/(1-t) - return workfunc(BatchIntegrand((v, t) -> begin resize!(xtmp, length(t)); resize!(ytmp, length(v)); - f.f!(ytmp, xtmp .= s0 .- oneunit(s1) .* t ./ (1 .- t)); v .= ytmp .* oneunit(s1) ./ (1 .- t) .^ 2; end, ybuf, xbuf, f.max_batch), - reverse(map(x -> 1 / (1 + oneunit(x) / (s0 - x)), s)), + return workfunc(BatchIntegrand((v, u) -> begin resize!(xtmp, length(u)); resize!(ytmp, length(v)); t = xtmp .= u ./ oneunit(s1); + f.f!(ytmp, u .= s0 .- oneunit(s1) .* t ./ (1 .- t)); v .= ytmp ./ (1 .- t) .^ 2; end, ybuf, xbuf, f.max_batch), + reverse(map(x -> oneunit(x) / (1 + oneunit(x) / (s0 - x)), s)), t -> s0 - oneunit(s1)*t/(1-t)) else # x = s0 + t/(1-t) - return workfunc(BatchIntegrand((v, t) -> begin resize!(xtmp, length(t)); resize!(ytmp, length(v)); - f.f!(ytmp, xtmp .= s0 .+ oneunit(s1) .* t ./ (1 .- t)); v .= ytmp .* oneunit(s1) ./ (1 .- t) .^ 2; end, ybuf, xbuf, f.max_batch), - map(x -> 1 / (1 + oneunit(x) / (x - s0)), s), + return workfunc(BatchIntegrand((v, u) -> begin resize!(xtmp, length(u)); resize!(ytmp, length(v)); t = xtmp .= u ./ oneunit(s1); + f.f!(ytmp, u .= s0 .+ oneunit(s1) .* t ./ (1 .- t)); v .= ytmp ./ (1 .- t) .^ 2; end, ybuf, xbuf, f.max_batch), + map(x -> oneunit(x) / (1 + oneunit(x) / (x - s0)), s), t -> s0 + oneunit(s1)*t/(1-t)) end end From f205c64cb5b283bb8c77a3efd28c97e099394db2 Mon Sep 17 00:00:00 2001 From: lxvm Date: Mon, 13 Nov 2023 12:48:57 -0500 Subject: [PATCH 2/9] add tests --- test/runtests.jl | 31 ++++++++++++++++++++++++++++++- 1 file changed, 30 insertions(+), 1 deletion(-) diff --git a/test/runtests.jl b/test/runtests.jl index 8c56014..d9aed5a 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -41,15 +41,23 @@ module Test19626 end # Following definitions needed for quadgk to work with MockQuantity - import Base: +, -, *, abs, isnan, isinf, isless, float + import Base: +, -, *, /, abs, isnan, isinf, isless, float, one, real, signbit +(a::MockQuantity, b::MockQuantity) = MockQuantity(a.val+b.val) -(a::MockQuantity, b::MockQuantity) = MockQuantity(a.val-b.val) + -(a::MockQuantity) = MockQuantity(-a.val) *(a::MockQuantity, b::Number) = MockQuantity(a.val*b) + *(a::Number, b::MockQuantity) = MockQuantity(a*b.val) + /(a::MockQuantity, b::Number) = MockQuantity(a.val/b) + /(a::MockQuantity, b::MockQuantity) = a.val/b.val abs(a::MockQuantity) = MockQuantity(abs(a.val)) float(a::MockQuantity) = a isnan(a::MockQuantity) = isnan(a.val) isinf(a::MockQuantity) = isinf(a.val) isless(a::MockQuantity, b::MockQuantity) = isless(a.val, b.val) + one(::Type{MockQuantity}) = one(fieldtype(MockQuantity, :val)) + one(a::MockQuantity) = one(a.val) + real(a::MockQuantity) = MockQuantity(real(a.val)) + signbit(a::MockQuantity) = signbit(a.val) # isapprox only needed for test purposes Base.isapprox(a::MockQuantity, b::MockQuantity) = isapprox(a.val, b.val) @@ -57,6 +65,27 @@ module Test19626 # Test physical quantity-valued functions @test QuadGK.quadgk(x->MockQuantity(x), 0.0, 1.0, atol=MockQuantity(0.0))[1] ≈ MockQuantity(0.5) + + # Test that quantities work with infinity transformations and segbufs + # for all quadgk interfaces + lims = [(-1.0, 1.0), (-Inf, 0.0), (0.0, Inf), (-Inf, Inf)] + ulims = map(x -> map(MockQuantity, x), lims) + for (lb, ub) in ulims + ## function + f = x -> 1/(1+(x/oneunit(x))^2) + buf = QuadGK.alloc_segbuf(MockQuantity, MockQuantity, MockQuantity) + @test QuadGK.quadgk(f, lb, ub, atol=MockQuantity(0.0))[1] ≈ + QuadGK.quadgk(f, lb, ub, atol=MockQuantity(0.0), segbuf=buf)[1] + ## inplace + fiip = (y, x) -> y[] = 1/(1+(x/oneunit(x))^2) + ibuf = QuadGK.alloc_segbuf(MockQuantity, Array{MockQuantity,0}, MockQuantity) + @test QuadGK.quadgk!(fiip, fill(MockQuantity(0.0)), lb, ub, atol=MockQuantity(0.0), norm=abs∘only)[1][] ≈ + QuadGK.quadgk!(fiip, fill(MockQuantity(0.0)), lb, ub, atol=MockQuantity(0.0), norm=abs∘only, segbuf=ibuf)[1][] + ## batch + fbatch = BatchIntegrand{Float64}((y, x) -> y .= 1 ./ (1 .+ (x ./ oneunit.(x)) .^ 2)) + @test QuadGK.quadgk(fbatch, lb, ub, atol=MockQuantity(0.0))[1] ≈ + QuadGK.quadgk(fbatch, lb, ub, atol=MockQuantity(0.0), segbuf=buf)[1] + end end @testset "inference" begin From 06f4f25f591fe98ca86787ad5d4d0fd933f72586 Mon Sep 17 00:00:00 2001 From: lxvm Date: Mon, 13 Nov 2023 13:02:00 -0500 Subject: [PATCH 3/9] fix CI tests --- test/runtests.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/runtests.jl b/test/runtests.jl index d9aed5a..fcea846 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -79,8 +79,8 @@ module Test19626 ## inplace fiip = (y, x) -> y[] = 1/(1+(x/oneunit(x))^2) ibuf = QuadGK.alloc_segbuf(MockQuantity, Array{MockQuantity,0}, MockQuantity) - @test QuadGK.quadgk!(fiip, fill(MockQuantity(0.0)), lb, ub, atol=MockQuantity(0.0), norm=abs∘only)[1][] ≈ - QuadGK.quadgk!(fiip, fill(MockQuantity(0.0)), lb, ub, atol=MockQuantity(0.0), norm=abs∘only, segbuf=ibuf)[1][] + @test QuadGK.quadgk!(fiip, fill(MockQuantity(0.0)), lb, ub, atol=MockQuantity(0.0), norm=abs∘getindex)[1][] ≈ + QuadGK.quadgk!(fiip, fill(MockQuantity(0.0)), lb, ub, atol=MockQuantity(0.0), norm=abs∘getindex, segbuf=ibuf)[1][] ## batch fbatch = BatchIntegrand{Float64}((y, x) -> y .= 1 ./ (1 .+ (x ./ oneunit.(x)) .^ 2)) @test QuadGK.quadgk(fbatch, lb, ub, atol=MockQuantity(0.0))[1] ≈ From 1764eee26d67e46ce2acec8f9c25f1a6fbfa371e Mon Sep 17 00:00:00 2001 From: lxvm Date: Mon, 13 Nov 2023 13:19:15 -0500 Subject: [PATCH 4/9] fix another CI problem --- test/runtests.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/test/runtests.jl b/test/runtests.jl index fcea846..66dc398 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -77,10 +77,10 @@ module Test19626 @test QuadGK.quadgk(f, lb, ub, atol=MockQuantity(0.0))[1] ≈ QuadGK.quadgk(f, lb, ub, atol=MockQuantity(0.0), segbuf=buf)[1] ## inplace - fiip = (y, x) -> y[] = 1/(1+(x/oneunit(x))^2) - ibuf = QuadGK.alloc_segbuf(MockQuantity, Array{MockQuantity,0}, MockQuantity) - @test QuadGK.quadgk!(fiip, fill(MockQuantity(0.0)), lb, ub, atol=MockQuantity(0.0), norm=abs∘getindex)[1][] ≈ - QuadGK.quadgk!(fiip, fill(MockQuantity(0.0)), lb, ub, atol=MockQuantity(0.0), norm=abs∘getindex, segbuf=ibuf)[1][] + fiip = (y, x) -> y[1] = 1/(1+(x/oneunit(x))^2) + ibuf = QuadGK.alloc_segbuf(MockQuantity, Array{MockQuantity,1}, MockQuantity) + @test QuadGK.quadgk!(fiip, [MockQuantity(0.0)], lb, ub, atol=MockQuantity(0.0), norm=abs∘first)[1][] ≈ + QuadGK.quadgk!(fiip, [MockQuantity(0.0)], lb, ub, atol=MockQuantity(0.0), norm=abs∘getindex, segbuf=ibuf)[1][] ## batch fbatch = BatchIntegrand{Float64}((y, x) -> y .= 1 ./ (1 .+ (x ./ oneunit.(x)) .^ 2)) @test QuadGK.quadgk(fbatch, lb, ub, atol=MockQuantity(0.0))[1] ≈ From b709b8b9120c06526921c504cb916c14cd1f9119 Mon Sep 17 00:00:00 2001 From: lxvm Date: Mon, 20 Nov 2023 09:51:19 -0500 Subject: [PATCH 5/9] make alloc_segbuf domain unitless --- src/adapt.jl | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/adapt.jl b/src/adapt.jl index ec8cf1a..75884c7 100644 --- a/src/adapt.jl +++ b/src/adapt.jl @@ -270,7 +270,9 @@ starting with the given `size`. The buffer can then be reused across multiple compatible calls to `quadgk(...)` to avoid repeated allocation. """ function alloc_segbuf(domain_type=Float64, range_type=Float64, error_type=Float64; size=1) - Vector{Segment{domain_type, range_type, error_type}}(undef, size) + x = float(real(one(domain_type))) # quadrature point type + err = oneunit(error_type)/oneunit(domain_type) # error has units of the integral + Vector{Segment{typeof(x), range_type, typeof(err)}}(undef, size) end """ From c7cbc9cd5c90d46a4cc5272358b06b2173d16a43 Mon Sep 17 00:00:00 2001 From: lxvm Date: Mon, 20 Nov 2023 11:44:59 -0500 Subject: [PATCH 6/9] strip units in handle_infinities --- src/QuadGK.jl | 4 +-- src/adapt.jl | 75 +++++++++++++++++++++++++++++------------------- src/batch.jl | 39 ++++++++++++++----------- test/runtests.jl | 4 +-- 4 files changed, 72 insertions(+), 50 deletions(-) diff --git a/src/QuadGK.jl b/src/QuadGK.jl index 90c8681..8878092 100644 --- a/src/QuadGK.jl +++ b/src/QuadGK.jl @@ -40,14 +40,14 @@ struct InplaceIntegrand{F,R,RI} Ig::R Ik::R fx::R - Idiff::RI + Idiff::R # final result array I::RI end InplaceIntegrand(f!::F, I::RI, fx::R) where {F,RI,R} = - InplaceIntegrand{F,R,RI}(f!, similar(fx), similar(fx), similar(fx), similar(fx), fx, similar(I), I) + InplaceIntegrand{F,R,RI}(f!, similar(fx), similar(fx), similar(fx), similar(fx), fx, similar(fx), I) include("gausskronrod.jl") include("evalrule.jl") diff --git a/src/adapt.jl b/src/adapt.jl index 75884c7..c8009e4 100644 --- a/src/adapt.jl +++ b/src/adapt.jl @@ -95,7 +95,7 @@ end # re-sum (paranoia about accumulated roundoff) function resum(f, segs) if f isa InplaceIntegrand - I = f.I .= segs[1].I + I = f.Ik .= segs[1].I E = segs[1].E for i in 2:length(segs) I .+= segs[i].I @@ -119,61 +119,77 @@ realone(x::Number) = one(x) isa Real # and pass transformed data to workfunc(f, s, tfunc) function handle_infinities(workfunc, f, s) s1, s2 = s[1], s[end] + u = float(real(oneunit(s1))) # the units of the segment if realone(s1) && realone(s2) # check for infinite or semi-infinite intervals inf1, inf2 = isinf(s1), isinf(s2) if inf1 || inf2 if inf1 && inf2 # x = t/(1-t^2) coordinate transformation - return workfunc(u -> begin t = u/oneunit(u); t2 = t*t; den = 1 / (1 - t2); - f(u*den) * (1+t2)*den*den; end, - map(x -> isinf(x) ? (signbit(x) ? -oneunit(x) : oneunit(x)) : 2x / (one(x)+hypot(one(x),2x/oneunit(x))), s), - t -> oneunit(s1) * t / (1 - t^2)) + I, E = workfunc(t -> begin t2 = t*t; den = 1 / (1 - t2); + f(u*t*den) * (1+t2)*den*den; end, + map(x -> isinf(x) ? (signbit(x) ? -one(x) : one(x)) : 2x / (oneunit(x)+hypot(oneunit(x),2x)), s), + t -> u * t / (1 - t^2)) + return u * I, u * E end let (s0,si) = inf1 ? (s2,s1) : (s1,s2) # let is needed for JuliaLang/julia#15276 if si < zero(si) # x = s0 - t/(1-t) - return workfunc(u -> begin t = u/oneunit(u); den = 1 / (1 - t); - f(s0 - u*den) * den*den; end, - reverse(map(x -> oneunit(x) / (1 + oneunit(x) / (s0 - x)), s)), - t -> s0 - oneunit(s1)*t/(1-t)) + I, E = workfunc(t -> begin den = 1 / (1 - t); + f(s0 - u*t*den) * den*den; end, + reverse(map(x -> 1 / (1 + oneunit(x) / (s0 - x)), s)), + t -> s0 - u*t/(1-t)) + return u * I, u * E else # x = s0 + t/(1-t) - return workfunc(u -> begin t = u/oneunit(u); den = 1 / (1 - t); - f(s0 + u*den) * den*den; end, - map(x -> oneunit(x) / (1 + oneunit(x) / (x - s0)), s), - t -> s0 + oneunit(s1)*t/(1-t)) + I, E = workfunc(t -> begin den = 1 / (1 - t); + f(s0 + u*t*den) * den*den; end, + map(x -> 1 / (1 + oneunit(x) / (x - s0)), s), + t -> s0 + u*t/(1-t)) + return u * I, u * E end end end end - return workfunc(f, s, identity) + I, E = workfunc(f, map(x -> x/oneunit(x), s), identity) + return u * I, u * E end function handle_infinities(workfunc, f::InplaceIntegrand, s) s1, s2 = s[1], s[end] + result = f.I + u = float(real(oneunit(s1))) # the units of the segment if realone(s1) && realone(s2) # check for infinite or semi-infinite intervals inf1, inf2 = isinf(s1), isinf(s2) if inf1 || inf2 - ftmp = f.fx # original integrand may have different type if inf1 && inf2 # x = t/(1-t^2) coordinate transformation - return workfunc(InplaceIntegrand((v, u) -> begin t = u/oneunit(u); t2 = t*t; den = 1 / (1 - t2); - f.f!(ftmp, u*den); v .= ftmp .* ((1+t2)*den*den); end, f.I, f.fx * float(one(s1))), - map(x -> isinf(x) ? (signbit(x) ? -oneunit(x) : oneunit(x)) : 2x / (one(x)+hypot(one(x),2x/oneunit(x))), s), - t -> oneunit(s1) * t / (1 - t^2)) + I, E = workfunc(InplaceIntegrand((v, t) -> begin t2 = t*t; den = 1 / (1 - t2); + 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)), + map(x -> isinf(x) ? (signbit(x) ? -one(x) : one(x)) : 2x / (oneunit(x)+hypot(oneunit(x),2x)), s), + t -> u * t / (1 - t^2)) + result .= u .* I + return result, u * E end let (s0,si) = inf1 ? (s2,s1) : (s1,s2) # let is needed for JuliaLang/julia#15276 if si < zero(si) # x = s0 - t/(1-t) - return workfunc(InplaceIntegrand((v, u) -> begin t = u/oneunit(u); den = 1 / (1 - t); - f.f!(ftmp, s0 - u*den); v .= ftmp .* (den * den); end, f.I, f.fx * float(one(s1))), - reverse(map(x -> oneunit(x) / (1 + oneunit(x) / (s0 - x)), s)), - t -> s0 - oneunit(s1)*t/(1-t)) + I, E = workfunc(InplaceIntegrand((v, t) -> begin den = 1 / (1 - t); + 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)), + reverse(map(x -> 1 / (1 + oneunit(x) / (s0 - x)), s)), + t -> s0 - u*t/(1-t)) + result .= u .* I + return result, u * E else # x = s0 + t/(1-t) - return workfunc(InplaceIntegrand((v, u) -> begin t = u/oneunit(u); den = 1 / (1 - t); - f.f!(ftmp, s0 + u*den); v .= ftmp .* (den * den); end, f.I, f.fx * float(one(s1))), - map(x -> oneunit(x) / (1 + oneunit(x) / (x - s0)), s), - t -> s0 + oneunit(s1)*t/(1-t)) + I, E = workfunc(InplaceIntegrand((v, t) -> begin den = 1 / (1 - t); + 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)), + map(x -> 1 / (1 + oneunit(x) / (x - s0)), s), + t -> s0 + u*t/(1-t)) + result .= u .* I + return result, u * E end end end end - return workfunc(f, s, identity) + I, E = workfunc(InplaceIntegrand(f.f!, f.fg, f.fk, f.Ig, f.Ik, f.fx, f.Idiff, similar(f.fx)), + map(x -> x/oneunit(x), s), + identity) + result .= u .* I + return result, u * E end function check_endpoint_roundoff(a, b, x; throw_error::Bool=false) @@ -254,8 +270,9 @@ quadgk(f, segs...; kws...) = function quadgk(f, segs::T...; atol=nothing, rtol=nothing, maxevals=10^7, order=7, norm=norm, segbuf=nothing) where {T} + utol = isnothing(atol) ? atol : atol/float(real(oneunit(T))) # remove units of domain handle_infinities(f, segs) do f, s, _ - do_quadgk(f, s, order, atol, rtol, maxevals, norm, segbuf) + do_quadgk(f, s, order, utol, rtol, maxevals, norm, segbuf) end end diff --git a/src/batch.jl b/src/batch.jl index fc53756..f366620 100644 --- a/src/batch.jl +++ b/src/batch.jl @@ -144,35 +144,40 @@ end function handle_infinities(workfunc, f::BatchIntegrand, s) s1, s2 = s[1], s[end] + u = float(real(oneunit(s1))) # the units of the segment + tbuf = similar(f.x, typeof(s1/oneunit(s1))) if realone(s1) && realone(s2) # check for infinite or semi-infinite intervals inf1, inf2 = isinf(s1), isinf(s2) if inf1 || inf2 - xbuf = f.x # buffer to store evaluation points - ytmp = f.y # original integrand may have different type - xtmp = similar(xbuf, typeof(one(eltype(f.x)))) - ybuf = similar(ytmp, typeof(oneunit(eltype(f.y))*float(one(s1)))) if inf1 && inf2 # x = t/(1-t^2) coordinate transformation - return workfunc(BatchIntegrand((v, u) -> begin resize!(xtmp, length(u)); resize!(ytmp, length(v)); t = xtmp .= u ./ oneunit(s1); - 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), - map(x -> isinf(x) ? (signbit(x) ? -oneunit(x) : oneunit(x)) : 2x / (one(x)+hypot(one(x),2x/oneunit(x))), s), - t -> oneunit(s1) * t / (1 - t^2)) + I, E = workfunc(BatchIntegrand((v, t) -> begin resize!(f.x, length(t)); + 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), + map(x -> isinf(x) ? (signbit(x) ? -one(x) : one(x)) : 2x / (oneunit(x)+hypot(oneunit(x),2x)), s), + t -> u * t / (1 - t^2)) + return u * I, u * E end let (s0,si) = inf1 ? (s2,s1) : (s1,s2) # let is needed for JuliaLang/julia#15276 if si < zero(si) # x = s0 - t/(1-t) - return workfunc(BatchIntegrand((v, u) -> begin resize!(xtmp, length(u)); resize!(ytmp, length(v)); t = xtmp .= u ./ oneunit(s1); - f.f!(ytmp, u .= s0 .- oneunit(s1) .* t ./ (1 .- t)); v .= ytmp ./ (1 .- t) .^ 2; end, ybuf, xbuf, f.max_batch), - reverse(map(x -> oneunit(x) / (1 + oneunit(x) / (s0 - x)), s)), - t -> s0 - oneunit(s1)*t/(1-t)) + I, E = workfunc(BatchIntegrand((v, t) -> begin resize!(f.x, length(t)); + f.f!(v, f.x .= s0 .- u .* t ./ (1 .- t)); v ./= (1 .- t) .^ 2; end, f.y, tbuf, f.max_batch), + reverse(map(x -> 1 / (1 + oneunit(x) / (s0 - x)), s)), + t -> s0 - u*t/(1-t)) + return u * I, u * E else # x = s0 + t/(1-t) - return workfunc(BatchIntegrand((v, u) -> begin resize!(xtmp, length(u)); resize!(ytmp, length(v)); t = xtmp .= u ./ oneunit(s1); - f.f!(ytmp, u .= s0 .+ oneunit(s1) .* t ./ (1 .- t)); v .= ytmp ./ (1 .- t) .^ 2; end, ybuf, xbuf, f.max_batch), - map(x -> oneunit(x) / (1 + oneunit(x) / (x - s0)), s), - t -> s0 + oneunit(s1)*t/(1-t)) + I, E = workfunc(BatchIntegrand((v, t) -> begin resize!(f.x, length(t)); + f.f!(v, f.x .= s0 .+ u .* t ./ (1 .- t)); v ./= (1 .- t) .^ 2; end, f.y, tbuf, f.max_batch), + map(x -> 1 / (1 + oneunit(x) / (x - s0)), s), + t -> s0 + u*t/(1-t)) + return u * I, u * E end end end end - return workfunc(f, s, identity) + I, E = workfunc(BatchIntegrand((y, t) -> begin resize!(f.x, length(t)); + f.f!(y, f.x .= u .* t); end, f.y, tbuf, f.max_batch), + map(x -> x/oneunit(x), s), + identity) + return u * I, u * E end """ diff --git a/test/runtests.jl b/test/runtests.jl index 66dc398..5ff0c5a 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -73,12 +73,12 @@ module Test19626 for (lb, ub) in ulims ## function f = x -> 1/(1+(x/oneunit(x))^2) - buf = QuadGK.alloc_segbuf(MockQuantity, MockQuantity, MockQuantity) + buf = QuadGK.alloc_segbuf(MockQuantity, Float64, MockQuantity) @test QuadGK.quadgk(f, lb, ub, atol=MockQuantity(0.0))[1] ≈ QuadGK.quadgk(f, lb, ub, atol=MockQuantity(0.0), segbuf=buf)[1] ## inplace fiip = (y, x) -> y[1] = 1/(1+(x/oneunit(x))^2) - ibuf = QuadGK.alloc_segbuf(MockQuantity, Array{MockQuantity,1}, MockQuantity) + ibuf = QuadGK.alloc_segbuf(MockQuantity, Array{Float64,1}, MockQuantity) @test QuadGK.quadgk!(fiip, [MockQuantity(0.0)], lb, ub, atol=MockQuantity(0.0), norm=abs∘first)[1][] ≈ QuadGK.quadgk!(fiip, [MockQuantity(0.0)], lb, ub, atol=MockQuantity(0.0), norm=abs∘getindex, segbuf=ibuf)[1][] ## batch From 119956da52136c1cb47c577b93992e9253868816 Mon Sep 17 00:00:00 2001 From: lxvm Date: Mon, 20 Nov 2023 11:57:18 -0500 Subject: [PATCH 7/9] add unitless buffer to BatchIntegrand --- src/batch.jl | 37 +++++++++++++++++++------------------ test/runtests.jl | 2 +- 2 files changed, 20 insertions(+), 19 deletions(-) diff --git a/src/batch.jl b/src/batch.jl index f366620..23b8e60 100644 --- a/src/batch.jl +++ b/src/batch.jl @@ -19,17 +19,19 @@ If `x` or `X` are not specified, `quadgk` internally creates a new `BatchIntegra user-supplied `y` buffer and a freshly-allocated `x` buffer based on the domain types. So, if you want to reuse the `x` buffer between calls, supply `{Y,X}` or pass `y,x` explicitly. """ -struct BatchIntegrand{Y,X,Ty<:AbstractVector{Y},Tx<:AbstractVector{X},F} +struct BatchIntegrand{Y,X,Ty<:AbstractVector{Y},Tx<:AbstractVector{X},F,T} # in-place function f!(y, x) that takes an array of x values and outputs an array of results in-place f!::F y::Ty x::Tx + t::T max_batch::Int # maximum number of x to supply in parallel end function BatchIntegrand(f!, y::AbstractVector, x::AbstractVector=similar(y, Nothing); max_batch::Integer=typemax(Int)) max_batch > 0 || throw(ArgumentError("max_batch must be positive")) - return BatchIntegrand(f!, y, x, max_batch) + X = eltype(x) + return BatchIntegrand(f!, y, x, X <: Nothing ? nothing : similar(x, typeof(one(X))), max_batch) end BatchIntegrand{Y,X}(f!; kws...) where {Y,X} = BatchIntegrand(f!, Y[], X[]; kws...) BatchIntegrand{Y}(f!; kws...) where {Y} = BatchIntegrand(f!, Y[]; kws...) @@ -65,20 +67,20 @@ function evalrules(f::BatchIntegrand, s::NTuple{N}, x,w,gw, nrm) where {N} l = length(x) m = 2l-1 # evaluations per segment n = (N-1)*m # total evaluations - resize!(f.x, n) + resize!(f.t, n) resize!(f.y, n) for i in 1:(N-1) # fill buffer with evaluation points a = s[i]; b = s[i+1] check_endpoint_roundoff(a, b, x, throw_error=true) c = convert(eltype(x), 0.5) * (b-a) o = (i-1)*m - f.x[l+o] = a + c + f.t[l+o] = a + c for j in 1:l-1 - f.x[j+o] = a + (1 + x[j]) * c - f.x[m+1-j+o] = a + (1 - x[j]) * c + f.t[j+o] = a + (1 + x[j]) * c + f.t[m+1-j+o] = a + (1 - x[j]) * c end end - f.f!(f.y, f.x) # evaluate integrand + f.f!(f.y, f.t) # evaluate integrand return ntuple(Val(N-1)) do i return batchevalrule(view(f.y, (1+(i-1)*m):(i*m)), s[i], s[i+1], x,w,gw, nrm) end @@ -105,7 +107,7 @@ function refine(f::BatchIntegrand, segs::Vector{T}, I, E, numevals, x,w,gw,n, at len > nsegs && DataStructures.percolate_down!(segs, 1, y, Reverse, len-nsegs) end - resize!(f.x, 2m*nsegs) + resize!(f.t, 2m*nsegs) resize!(f.y, 2m*nsegs) for i in 1:nsegs # fill buffer with evaluation points s = segs[len-i+1] @@ -114,15 +116,15 @@ function refine(f::BatchIntegrand, segs::Vector{T}, I, E, numevals, x,w,gw,n, at check_endpoint_roundoff(a, b, x) && return segs c = convert(eltype(x), 0.5) * (b-a) o = (2i-j)*m - f.x[l+o] = a + c + f.t[l+o] = a + c for k in 1:l-1 # early return if integrand evaluated at endpoints - f.x[k+o] = a + (1 + x[k]) * c - f.x[m+1-k+o] = a + (1 - x[k]) * c + f.t[k+o] = a + (1 + x[k]) * c + f.t[m+1-k+o] = a + (1 - x[k]) * c end end end - f.f!(f.y, f.x) + f.f!(f.y, f.t) resize!(segs, len+nsegs) for i in 1:nsegs # evaluate segments and update estimates & heap @@ -145,13 +147,12 @@ end function handle_infinities(workfunc, f::BatchIntegrand, s) s1, s2 = s[1], s[end] u = float(real(oneunit(s1))) # the units of the segment - tbuf = similar(f.x, typeof(s1/oneunit(s1))) if realone(s1) && realone(s2) # check for infinite or semi-infinite intervals inf1, inf2 = isinf(s1), isinf(s2) if inf1 || inf2 if inf1 && inf2 # x = t/(1-t^2) coordinate transformation I, E = workfunc(BatchIntegrand((v, t) -> begin resize!(f.x, length(t)); - 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), + f.f!(v, f.x .= u .* t ./ (1 .- t .* t)); v .*= (1 .+ t .* t) ./ (1 .- t .* t) .^ 2; end, f.y, f.x, f.t, f.max_batch), map(x -> isinf(x) ? (signbit(x) ? -one(x) : one(x)) : 2x / (oneunit(x)+hypot(oneunit(x),2x)), s), t -> u * t / (1 - t^2)) return u * I, u * E @@ -159,13 +160,13 @@ function handle_infinities(workfunc, f::BatchIntegrand, s) let (s0,si) = inf1 ? (s2,s1) : (s1,s2) # let is needed for JuliaLang/julia#15276 if si < zero(si) # x = s0 - t/(1-t) I, E = workfunc(BatchIntegrand((v, t) -> begin resize!(f.x, length(t)); - f.f!(v, f.x .= s0 .- u .* t ./ (1 .- t)); v ./= (1 .- t) .^ 2; end, f.y, tbuf, f.max_batch), + f.f!(v, f.x .= s0 .- u .* t ./ (1 .- t)); v ./= (1 .- t) .^ 2; end, f.y, f.x, f.t, f.max_batch), reverse(map(x -> 1 / (1 + oneunit(x) / (s0 - x)), s)), t -> s0 - u*t/(1-t)) return u * I, u * E else # x = s0 + t/(1-t) I, E = workfunc(BatchIntegrand((v, t) -> begin resize!(f.x, length(t)); - f.f!(v, f.x .= s0 .+ u .* t ./ (1 .- t)); v ./= (1 .- t) .^ 2; end, f.y, tbuf, f.max_batch), + f.f!(v, f.x .= s0 .+ u .* t ./ (1 .- t)); v ./= (1 .- t) .^ 2; end, f.y, f.x, f.t, f.max_batch), map(x -> 1 / (1 + oneunit(x) / (x - s0)), s), t -> s0 + u*t/(1-t)) return u * I, u * E @@ -174,7 +175,7 @@ function handle_infinities(workfunc, f::BatchIntegrand, s) end end I, E = workfunc(BatchIntegrand((y, t) -> begin resize!(f.x, length(t)); - f.f!(y, f.x .= u .* t); end, f.y, tbuf, f.max_batch), + f.f!(y, f.x .= u .* t); end, f.y, f.x, f.t, f.max_batch), map(x -> x/oneunit(x), s), identity) return u * I, u * E @@ -199,6 +200,6 @@ simultaneously. In particular, there are two differences from `quadgk` """ function quadgk(f::BatchIntegrand{Y,Nothing}, segs::T...; kws...) where {Y,T} FT = float(T) # the gk points are floating-point - g = BatchIntegrand(f.f!, f.y, similar(f.x, FT), f.max_batch) + g = BatchIntegrand(f.f!, f.y, similar(f.x, FT), similar(f.x, typeof(float(one(FT)))), f.max_batch) return quadgk(g, segs...; kws...) end diff --git a/test/runtests.jl b/test/runtests.jl index 5ff0c5a..26d537f 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -341,7 +341,7 @@ end end # test constructors - ref = BatchIntegrand(f!, Float64[], Nothing[], typemax(Int)) + ref = BatchIntegrand(f!, Float64[], Nothing[], nothing, typemax(Int)) for b in ( BatchIntegrand(f!, Float64[]), BatchIntegrand(f!, Float64[], Nothing[]), From bfd48b5443ef67edefe9526964e3d2f97caf8337 Mon Sep 17 00:00:00 2001 From: lxvm Date: Sat, 30 Dec 2023 00:32:33 -0800 Subject: [PATCH 8/9] ensure unit gets applied in ordinary case --- src/adapt.jl | 4 ++-- test/runtests.jl | 6 +++--- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/src/adapt.jl b/src/adapt.jl index c8009e4..7da5191 100644 --- a/src/adapt.jl +++ b/src/adapt.jl @@ -147,7 +147,7 @@ function handle_infinities(workfunc, f, s) end end end - I, E = workfunc(f, map(x -> x/oneunit(x), s), identity) + I, E = workfunc(t -> f(t*u), map(x -> x/oneunit(x), s), identity) return u * I, u * E end @@ -185,7 +185,7 @@ function handle_infinities(workfunc, f::InplaceIntegrand, s) end end end - I, E = workfunc(InplaceIntegrand(f.f!, f.fg, f.fk, f.Ig, f.Ik, f.fx, f.Idiff, similar(f.fx)), + I, E = workfunc(InplaceIntegrand((y,t) -> f.f!(y, t*u), f.fg, f.fk, f.Ig, f.Ik, f.fx, f.Idiff, similar(f.fx)), map(x -> x/oneunit(x), s), identity) result .= u .* I diff --git a/test/runtests.jl b/test/runtests.jl index 26d537f..43f915d 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -72,17 +72,17 @@ module Test19626 ulims = map(x -> map(MockQuantity, x), lims) for (lb, ub) in ulims ## function - f = x -> 1/(1+(x/oneunit(x))^2) + f = x -> 1/(1+(x/MockQuantity(1.0))^2) buf = QuadGK.alloc_segbuf(MockQuantity, Float64, MockQuantity) @test QuadGK.quadgk(f, lb, ub, atol=MockQuantity(0.0))[1] ≈ QuadGK.quadgk(f, lb, ub, atol=MockQuantity(0.0), segbuf=buf)[1] ## inplace - fiip = (y, x) -> y[1] = 1/(1+(x/oneunit(x))^2) + fiip = (y, x) -> y[1] = 1/(1+(x/MockQuantity(1.0))^2) ibuf = QuadGK.alloc_segbuf(MockQuantity, Array{Float64,1}, MockQuantity) @test QuadGK.quadgk!(fiip, [MockQuantity(0.0)], lb, ub, atol=MockQuantity(0.0), norm=abs∘first)[1][] ≈ QuadGK.quadgk!(fiip, [MockQuantity(0.0)], lb, ub, atol=MockQuantity(0.0), norm=abs∘getindex, segbuf=ibuf)[1][] ## batch - fbatch = BatchIntegrand{Float64}((y, x) -> y .= 1 ./ (1 .+ (x ./ oneunit.(x)) .^ 2)) + fbatch = BatchIntegrand{Float64}((y, x) -> y .= 1 ./ (1 .+ (x ./ MockQuantity(1.0)) .^ 2)) @test QuadGK.quadgk(fbatch, lb, ub, atol=MockQuantity(0.0))[1] ≈ QuadGK.quadgk(fbatch, lb, ub, atol=MockQuantity(0.0), segbuf=buf)[1] end From f328bd7b2cb5c21e795a815f61b09f6ba9ddb51d Mon Sep 17 00:00:00 2001 From: lxvm Date: Sat, 30 Dec 2023 02:06:57 -0800 Subject: [PATCH 9/9] restore resum --- src/adapt.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/adapt.jl b/src/adapt.jl index 7da5191..0459bb5 100644 --- a/src/adapt.jl +++ b/src/adapt.jl @@ -95,7 +95,7 @@ end # re-sum (paranoia about accumulated roundoff) function resum(f, segs) if f isa InplaceIntegrand - I = f.Ik .= segs[1].I + I = f.I .= segs[1].I E = segs[1].E for i in 2:length(segs) I .+= segs[i].I