diff --git a/Project.toml b/Project.toml index 6f84304a..02426731 100644 --- a/Project.toml +++ b/Project.toml @@ -8,10 +8,10 @@ Bessels = "0e736298-9ec6-45e8-9647-e4fc86a2fe38" ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" CommonSolve = "38540f10-b2f7-11e9-35d8-d573e4eb0ff2" DataInterpolations = "82cc6244-b520-54b8-b5a6-8a565e85f1d0" +FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341" ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" ForwardDiffChainRules = "c9556dd2-1aed-4cfe-8560-1557cf593001" GraphRecipes = "bd48cda9-67a9-57be-86fa-5b3c104eda73" -Integrals = "de52edbc-65ea-441a-8357-d3a637375a31" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae" MatterPower = "743831d1-0e6a-4e86-a0c4-1a098a68aca3" @@ -51,10 +51,10 @@ Bessels = "0.2.8" ChainRulesCore = "1.25.1" CommonSolve = "0.2.4" DataInterpolations = "8.0.0" +FFTW = "1.10.0" ForwardDiff = "0.10, 1.0.0" ForwardDiffChainRules = "0.3.0" GraphRecipes = "0.5.13" -Integrals = "4.5.0, 5.0.0" LinearAlgebra = "1.10.0" LinearSolve = "3.19.0" MatterPower = "0.5.0" diff --git a/docs/src/comparison.md b/docs/src/comparison.md index 7621b323..7539bfc8 100644 --- a/docs/src/comparison.md +++ b/docs/src/comparison.md @@ -412,8 +412,9 @@ end l = 20:20:2000 # CLASS default is lmax = 2500 Dl1 = Dl_class([:TT, :TE, :EE, :phiphi, :TPhi, :Ephi], l, pars) -jl = SphericalBesselCache(l) -Dl2 = Dl_symboltz([:TT, :TE, :EE, :ψψ, :ψT, :ψE], jl, pars) +jl = SphericalBesselCache(l) # precompute +τquad = ClenshawCurtisQuadrature(500) # precompute +Dl2 = Dl_symboltz([:TT, :TE, :EE, :ψψ, :ψT, :ψE], jl, pars; τquad) plot_compare(l, l, Dl1[:, 1], Dl2[:, 1], "l", "Dₗ(TT)"; tol = 2e-12) ``` ```@example class diff --git a/src/SymBoltz.jl b/src/SymBoltz.jl index 27e5b754..1b5a986a 100644 --- a/src/SymBoltz.jl +++ b/src/SymBoltz.jl @@ -40,6 +40,7 @@ include("models/curvature.jl") include("models/inflation.jl") include("models/cosmologies.jl") include("solve.jl") +include("quadrature.jl") include("observables/distances.jl") include("observables/fourier.jl") include("observables/angular.jl") @@ -53,6 +54,7 @@ export solve, solvebg, solvept, remake, issuccess, parameter_updater export parameters_Planck18 export spectrum_primordial, spectrum_matter, spectrum_matter_nonlinear, spectrum_cmb, correlation_function, variance_matter, stddev_matter, los_integrate, source_grid, source_grid_adaptive, sound_horizon, distance_luminosity, SphericalBesselCache export express_derivatives +export Quadrature, TrapezoidalQuadrature, ClenshawCurtisQuadrature, GaussQuadrature, GaussKronrodQuadrature, nodes, weights, integrate, transform using PrecompileTools: @compile_workload @compile_workload begin diff --git a/src/observables/angular.jl b/src/observables/angular.jl index d700b372..270cd75b 100644 --- a/src/observables/angular.jl +++ b/src/observables/angular.jl @@ -1,4 +1,3 @@ -using Integrals using Bessels: besselj!, sphericalbesselj using DataInterpolations using MatterPower @@ -16,7 +15,7 @@ struct SphericalBesselCache{Tdy <: Union{Matrix{Float64}, Nothing}} x::Vector{Float64} end -function SphericalBesselCache(ls::AbstractVector; xmax = 10*ls[end], dx = 2π/15, hermite = true) +function SphericalBesselCache(ls::AbstractVector; dx = 2π/15, xmax = 10*ls[end]+dx, hermite = true) # add one dx to xmax to prevent floating-point bounds checking issues xmin = 0.0 xs = range(xmin, xmax, length = trunc(Int, (xmax - xmin) / dx)) # fixed length (so endpoints are exact) that gives step as close to dx as possible invdx = 1.0 / step(xs) # using the resulting step, which need not be exactly dx @@ -102,7 +101,7 @@ ChainRulesCore.frule((_, _, Δx), ::typeof(jl), l, x) = jl(l, x), jl′(l, x) * # TODO: use u = k*χ as integration variable, so oscillations of Bessel functions are the same for every k? # TODO: define and document symbolic dispatch! """ - los_integrate(Ss::AbstractMatrix{T}, ls::AbstractVector, τs::AbstractVector, ks::AbstractVector, jl::SphericalBesselCache; l_limber = typemax(Int), integrator = TrapezoidalRule(), thread = true, verbose = false) where {T <: Real} + los_integrate(Ss::AbstractMatrix{T}, ls::AbstractVector, τs::AbstractVector, ks::AbstractVector, jl::SphericalBesselCache[, quad::Quadrature]; l_limber = typemax(Int), thread = true, verbose = false) where {T <: Real} For the given `ls` and `ks`, compute the line-of-sight-integrals ```math @@ -116,7 +115,7 @@ Iₗ ≈ √(π/(2l+1)) S(τ₀-(l+1/2)/k, k) ``` is used for `l ≥ l_limber`. """ -function los_integrate(Ss::AbstractMatrix{T}, ls::AbstractVector, τs::AbstractVector, ks::AbstractVector, jl::SphericalBesselCache; l_limber = typemax(Int), integrator = TrapezoidalRule(), thread = true, verbose = false) where {T <: Real} +function los_integrate(Ss::AbstractMatrix{T}, ls::AbstractVector, τs::AbstractVector, ks::AbstractVector, jl::SphericalBesselCache, quad::Quadrature; l_limber = typemax(Int), thread = true, verbose = false) where {T <: Real} # Julia is column-major; make sure innermost loop indices appear first in slice expressions (https://docs.julialang.org/en/v1/manual/performance-tips/#man-performance-column-major) @assert size(Ss, 1) == length(τs) "size(Ss, 1) = $(size(Ss, 1)) and length(τs) = $(length(τs)) differ" @assert size(Ss, 2) == length(ks) "size(Ss, 2) = $(size(Ss, 2)) and length(ks) = $(length(ks)) differ" @@ -126,7 +125,6 @@ function los_integrate(Ss::AbstractMatrix{T}, ls::AbstractVector, τs::AbstractV τs = collect(τs) # force array to avoid floating point errors with ranges in following χs due to (e.g. tiny negative χ) τ0 = τs[end] χs = τ0 .- τs - halfdτs = 0.5 .* (τs[begin+1:end] .- τs[begin:end-1]) # precompute before loops Is = zeros(T, length(ks), length(ls)) verbose && l_limber < typemax(Int) && println("Using Limber approximation for l ≥ $l_limber") @@ -134,6 +132,7 @@ function los_integrate(Ss::AbstractMatrix{T}, ls::AbstractVector, τs::AbstractV # TODO: skip and set jl to zero if l ≳ kτ0 or another cutoff? @fastmath @inbounds @tasks for il in eachindex(ls) # parallellize independent loop iterations @set scheduler = thread ? :dynamic : :serial + @local integrand = zeros(T, length(τs)) # local workspace for LOS integrand per task l = ls[il] verbose && print("\rLOS integrating with l = $l") for ik in reverse(eachindex(ks)) @@ -151,17 +150,8 @@ function los_integrate(Ss::AbstractMatrix{T}, ls::AbstractVector, τs::AbstractV I = √(π/(2l+1)) * S / k end else - prev = Ss[1, ik] * jl(l, k*χs[1]) # set up first point - for iτ in 2:length(τs) - kχ = k * χs[iτ] - halfdτ = halfdτs[iτ-1] - _jl = jl(l, kχ) - curr = Ss[iτ, ik] * _jl - dI = halfdτ * (curr + prev) - I += dI - kχ < l && abs(_jl) < 1e-20 && break # time cut approximation - prev = curr - end + integrand .= (@view Ss[:, ik]) .* jl.(l, k .* χs) + I = quad(integrand) end Is[ik, il] = I k*τ0 < l && abs(I) < 1e-20 && break # multipole cut approximation @@ -182,7 +172,7 @@ end # TODO: integrate splines instead of trapz! https://discourse.julialang.org/t/how-to-speed-up-the-numerical-integration-with-interpolation/96223/5 @doc raw""" - spectrum_cmb(ΘlAs::AbstractMatrix, ΘlBs::AbstractMatrix, P0s::AbstractVector, ls::AbstractVector, ks::AbstractVector; integrator = TrapezoidalRule(), normalization = :Cl, thread = true) + spectrum_cmb(ΘlAs::AbstractMatrix, ΘlBs::AbstractMatrix, P0s::AbstractVector, ls::AbstractVector, ks::AbstractVector[, quad::Quadrature]; normalization = :Cl, thread = true) Compute the angular power spectrum ```math @@ -191,22 +181,20 @@ Cₗᴬᴮ = (2/π) ∫\mathrm{d}k \, k² P₀(k) Θₗᴬ(k,τ₀) Θₗᴮ(k, for the given `ls`. If `normalization == :Dl`, compute ``Dₗ = Cₗ l (l+1) / 2π`` instead. """ -function spectrum_cmb(ΘlAs::AbstractMatrix, ΘlBs::AbstractMatrix, P0s::AbstractVector, ls::AbstractVector, ks::AbstractVector; integrator = TrapezoidalRule(), normalization = :Cl, thread = true) +function spectrum_cmb(ΘlAs::AbstractMatrix, ΘlBs::AbstractMatrix, P0s::AbstractVector, ls::AbstractVector, ks::AbstractVector, quad::Quadrature; normalization = :Cl, thread = true) size(ΘlAs) == size(ΘlBs) || error("ΘlAs and ΘlBs have different sizes") eltype(ΘlAs) == eltype(ΘlBs) || error("ΘlAs and ΘlBs have different types") Cls = similar(ΘlAs, length(ls)) - ks_with0 = [0.0; ks] # add dummy value with k=0 for integration @tasks for il in eachindex(ls) # TODO: skip kτ0 ≲ l? @set scheduler = thread ? :dynamic : :static - @local dCl_dks_with0 = zeros(eltype(ΘlAs), length(ks_with0)) # local task workspace (must zero first element) + @local dCl_dks = zeros(eltype(ΘlAs), length(ks)) # local task workspace (must zero first element) ΘlA = @view ΘlAs[:, il] ΘlB = @view ΘlBs[:, il] - @. dCl_dks_with0[2:end] = 2/π * ks^2 * P0s * ΘlA * ΘlB - spline = CubicSpline(dCl_dks_with0, ks_with0) - Cls[il] = DataInterpolations.integral(spline, ks_with0[begin], ks_with0[end]) # integrate over k (_with0 adds one additional point at (0,0)) + dCl_dks .= 2/π .* ks .^ 2 .* P0s .* ΘlA .* ΘlB + Cls[il] = quad(dCl_dks) # integrate over k (_with0 adds one additional point at (0,0)) end if normalization == :Cl @@ -219,7 +207,7 @@ function spectrum_cmb(ΘlAs::AbstractMatrix, ΘlBs::AbstractMatrix, P0s::Abstrac end """ - spectrum_cmb(modes::AbstractVector{<:Symbol}, prob::CosmologyProblem, jl::SphericalBesselCache; normalization = :Cl, unit = nothing, kτ0s = 0.1*jl.l[begin]:2π/2:10*jl.l[end], xs = 0.0:0.0008:1.0, l_limber = 10, integrator = TrapezoidalRule(), bgopts = (alg = bgalg(prob), reltol = 1e-7, abstol = 1e-7), ptopts = (alg = ptalg(prob),, reltol = 1e-5, abstol = 1e-5), sourceopts = (rtol = 1e-3, atol = 0.9), downsampleopts = (Ttol = 4e-3, Etol = 2e-4, ψtol = 1e-3), coarse_length = 9, thread = true, verbose = false, kwargs...) + spectrum_cmb(modes::AbstractVector{<:Symbol}, prob::CosmologyProblem, jl::SphericalBesselCache; normalization = :Cl, unit = nothing, kτ0_l_min = 0.1, kτ0_l_max = 10.0, τquad = ClenshawCurtisQuadrature(500), kquad = TrapezoidalQuadrature(8500), l_limber = 10, bgopts = (alg = bgalg(prob), reltol = 1e-7, abstol = 1e-7), ptopts = (alg = ptalg(prob), reltol = 1e-5, abstol = 1e-5), sourceopts = (rtol = 1e-3, atol = 0.9), coarse_length = 9, thread = true, verbose = false, kwargs...) Compute angular CMB power spectra ``Cₗᴬᴮ`` at angular wavenumbers `ls` from the cosmological problem `prob`. The requested `modes` are specified as a vector of symbols in the form `:AB`, where `A` and `B` are `T` (temperature), `E` (E-mode polarization) or `ψ` (lensing). @@ -245,19 +233,18 @@ modes = [:TT, :TE, :ψψ, :ψT] Dls = spectrum_cmb(modes, prob, jl; normalization = :Dl, unit = u"μK") ``` """ -function spectrum_cmb(modes::AbstractVector{<:Symbol}, prob::CosmologyProblem, jl::SphericalBesselCache; normalization = :Cl, unit = nothing, kτ0s = 0.1*jl.l[begin]:2π/2:10*jl.l[end], xs = 0.0:0.0008:1.0, l_limber = 10, integrator = TrapezoidalRule(), bgopts = (alg = bgalg(prob), reltol = 1e-7, abstol = 1e-7), ptopts = (alg = ptalg(prob), reltol = 1e-5, abstol = 1e-5), sourceopts = (rtol = 1e-3, atol = 0.9), downsampleopts = (Ttol = 4e-3, Etol = 2e-4, ψtol = 1e-3), coarse_length = 9, thread = true, verbose = false, kwargs...) +function spectrum_cmb(modes::AbstractVector{<:Symbol}, prob::CosmologyProblem, jl::SphericalBesselCache; normalization = :Cl, unit = nothing, kτ0_l_min = 0.1, kτ0_l_max = 10.0, τquad = ClenshawCurtisQuadrature(500), kquad = TrapezoidalQuadrature(8500), l_limber = 10, bgopts = (alg = bgalg(prob), reltol = 1e-7, abstol = 1e-7), ptopts = (alg = ptalg(prob), reltol = 1e-5, abstol = 1e-5), sourceopts = (rtol = 1e-3, atol = 0.9), coarse_length = 9, thread = true, verbose = false, kwargs...) ls = jl.l + lmin = ls[begin] + lmax = ls[end] sol = solve(prob; bgopts, verbose) τ0 = getsym(sol, prob.M.τ0)(sol) - ks_fine = collect(kτ0s ./ τ0) + kquad = transform(kquad, (kτ0_l_min*lmin/τ0, kτ0_l_max*lmax/τ0)) # map from [-1, +1] to [kmin, kmax] + ks_fine = nodes(kquad) τs = sol.bg.t # by default, use background (thermodynamics) time points for line of sight integration - if !isnothing(xs) - # use user's array of x = (τ-τi)/(τ0-τi) - xs[begin] == 0 || error("xs begins with $(xs[begin]), but should begin with 0") - xs[end] == 1 || error("xs ends with $(xs[end]), but should end with 1") - τs = τs[begin] .+ (τs[end] .- τs[begin]) .* xs - end + τquad = transform(τquad, (τs[begin], τs[end])) # map from [-1, 1] to [τi, τ0] + τs = nodes(τquad) # Integrate perturbations to calculate source function on coarse k-grid iT = 'T' in join(modes) ? 1 : 0 @@ -269,25 +256,22 @@ function spectrum_cmb(modes::AbstractVector{<:Symbol}, prob::CosmologyProblem, j Θls = zeros(eltype(Ss), max(iT, iE, iψ), length(ks_fine), length(ls)) if iT > 0 - STs, τTs = source_grid_downsample(Ss[1, :, :], τs; tol = downsampleopts.Ttol) # downsample in τ - verbose && println("Downsampled T source function from ", length(τs), " to ", length(τTs), " time points") + STs = Ss[1, :, :] STs = source_grid(STs, ks_coarse, ks_fine; thread) # upsample in k - Θls[iT, :, :] .= los_integrate(STs, ls, τTs, ks_fine, jl; integrator, verbose, thread, kwargs...) + Θls[iT, :, :] .= los_integrate(STs, ls, τs, ks_fine, jl, τquad; verbose, thread, kwargs...) end if iE > 0 - SEs, τEs = source_grid_downsample(Ss[2, :, :], τs; tol = downsampleopts.Etol) # downsample in τ - verbose && println("Downsampled E source function from ", length(τs), " to ", length(τEs), " time points") - @. SEs ./= (ks_coarse' * (τ0-τEs))^2 + SEs = Ss[2, :, :] + @. SEs ./= (ks_coarse' * (τ0-τs))^2 SEs[end, :] .= 0.0 # can be Inf, but is always weighted by zero-valued spherical Bessel function in LOS integration SEs = source_grid(SEs, ks_coarse, ks_fine; thread) # upsample in k - Θls[iE, :, :] .= transpose(@. √((ls+2)*(ls+1)*(ls+0)*(ls-1))) .* los_integrate(SEs, ls, τEs, ks_fine, jl; integrator, verbose, thread, kwargs...) + Θls[iE, :, :] .= transpose(@. √((ls+2)*(ls+1)*(ls+0)*(ls-1))) .* los_integrate(SEs, ls, τs, ks_fine, jl, τquad; verbose, thread, kwargs...) end if iψ > 0 - Sψs, τψs = source_grid_downsample(Ss[3, :, :], τs; tol = downsampleopts.ψtol) # downsample in τ - verbose && println("Downsampled ψ source function from ", length(τs), " to ", length(τψs), " time points") + Sψs = Ss[3, :, :] Sψs = source_grid(Sψs, ks_coarse, ks_fine; thread) # upsample in k Sψs[end, :] .= 0.0 # contains Inf/NaN, but will be weighted by 0 from jl in LOS integral - Θls[iψ, :, :] .= los_integrate(Sψs, ls, τψs, ks_fine, jl; l_limber, integrator, verbose, thread, kwargs...) + Θls[iψ, :, :] .= los_integrate(Sψs, ls, τs, ks_fine, jl, τquad; l_limber, verbose, thread, kwargs...) end P0s = spectrum_primordial(ks_fine, sol) # more accurate @@ -314,7 +298,7 @@ function spectrum_cmb(modes::AbstractVector{<:Symbol}, prob::CosmologyProblem, j iB = geti(Symbol(mode[lastindex(mode)])) ΘlAs = @view(Θls[iA, :, :]) ΘlBs = @view(Θls[iB, :, :]) - spectrum = spectrum_cmb(ΘlAs, ΘlBs, P0s, ls, ks_fine; integrator, normalization, thread) + spectrum = spectrum_cmb(ΘlAs, ΘlBs, P0s, ls, ks_fine, kquad; normalization, thread) spectrum *= factor^2 # possibly make dimensionful spectra[:, i] .= spectrum end diff --git a/src/quadrature.jl b/src/quadrature.jl new file mode 100644 index 00000000..fb89c2d8 --- /dev/null +++ b/src/quadrature.jl @@ -0,0 +1,89 @@ +using FFTW +using QuadGK + +struct Quadrature{T, U} + x::Vector{T} # integration points on [-1, +1] + w::Vector{T} # integration weights on [-1, +1] + a::U + b::U + name::Symbol + + function Quadrature(x::AbstractVector{T}, w::AbstractVector{T}, interval = (-1.0, 1.0); name = Symbol()) where {T} + all(-1 .≤ x .≤ 1) || throw(ArgumentError("Quadrature points must be on the canonical interval [-1, 1]")) + issorted(x) || throw(ArgumentError("Quadrature points must be sorted in ascending order")) + sum(w) ≈ 2 || throw(ArgumentError("Quadrature weights must sum to 2, but sums to $(sum(w))")) + length(x) == length(w) || throw(ArgumentError("Quadrature nodes and weights must have same length")) + new{T, eltype(interval)}(x, w, interval[begin], interval[end], name) + end +end + +function TrapezoidalQuadrature(N::Integer, args...) + N ≥ 2 || throw(ArgumentError("Trapezoidal quadrature needs at least 2 points")) + x = collect(range(-1, 1, length = N)) + w = fill(2/(N-1), N) + w[begin] = 1/(N-1) + w[end] = 1/(N-1) + return Quadrature(x, w, args...; name = Symbol("Trapezoidal")) +end + +function TrapezoidalQuadrature(x::AbstractArray, args...) + N = length(x) + N ≥ 2 || throw(ArgumentError("Trapezoidal quadrature needs at least 2 points")) + xmin, xmax = extrema(x) + x = collect(x) + x .= -1 .+ 2 .* (x .- xmin) / (xmax - xmin) # normalize to canonical [-1, 1] + w = zeros(N) + # each interval contributes (y1+y2)*(x2-x1)/2 + for i in 1:N-1 + dx = x[i+1] - x[i] + w[i] += dx / 2 + w[i+1] += dx / 2 + end + return Quadrature(x, w, (xmin, xmax), args...; name = Symbol("Trapezoidal")) +end + + +function ClenshawCurtisQuadrature(N::Integer, args...) + N ≥ 2 || throw(ArgumentError("Clenshaw-Curtis quadrature needs at least 2 points")) + n = N - 1 # number of FFT points + x = [-cos(π*m/n) for m in 0:n] # Chebyshev nodes in ascending order + v = [1 / (1 - 4*min(m, n-m)^2) for m in 0:n-1] # 1/(1-4k^2) and it's mirror image + w = similar(x) + w[begin:end-1] .= 2/n .* real(fft(v)) # Discrete Cosine Transform (DCT) for O(N log N) time instead of O(N²) + w[begin] = w[end] = w[begin]/2 # modify endpoint factors + return Quadrature(x, w, args...; name = Symbol("Clenshaw-Curtis")) +end + +function GaussQuadrature(N::Integer, args...) + N ≥ 1 || throw(ArgumentError("The number of Gauss quadrature points must be positive")) + x, w = QuadGK.gauss(N) + return Quadrature(x, w, args...; name = Symbol("Gauss")) +end + +function GaussKronrodQuadrature(N::Integer, args...) + N ≥ 3 && isodd(N) || throw(ArgumentError("The number of Gauss-Kronrod quadrature points must be at least 3 and odd")) + x, w = QuadGK.kronrod(div(N-1, 2)) # only returns left half + x = [x; -reverse(x[begin:end-1])] # add right half (antisymmetric) + w = [w; reverse(w[begin:end-1])] # add right half (symmetric) + return Quadrature(x, w, args...; name = Symbol("Gauss-Kronrod")) +end + +transform(q::Quadrature, interval) = Quadrature(q.x, q.w, interval; name = q.name) # arrays are reuse and shared! +Base.eltype(::Quadrature{T, U}) where {T, U} = Base.promote_type(T, U) +Base.nameof(q::Quadrature) = q.name +Base.show(io::IO, q::Quadrature) = print(io, q.name == Symbol() ? "Q" : "$(q.name) q", "uadrature: on [$(q.a), $(q.b)], $(length(q)) points, eltype = $(eltype(q))") +Base.length(q::Quadrature) = length(q.x) +Base.eachindex(q::Quadrature) = eachindex(q.x) +Base.:(==)(q1::Quadrature, q2::Quadrature) = q1.x == q2.x && q1.w == q2.w && q1.a == q2.a && q1.b == q2.b +Base.:(≈)(q1::Quadrature, q2::Quadrature) = q1.x ≈ q2.x && q1.w ≈ q2.w && q1.a ≈ q2.a && q1.b ≈ q2.b +node(q::Quadrature, i) = (q.b+q.a)/2 + (q.b-q.a)/2 * q.x[i] +nodes(q::Quadrature) = (q.b+q.a)/2 .+ (q.b-q.a)/2 .* q.x +weights(q::Quadrature) = (q.b-q.a)/2 .* q.w # weights for integration on [a, b] +@inbounds @fastmath function integrate(q::Quadrature, f::AbstractVector) + length(f) == length(q) || throw(ArgumentError("$(length(f))-point f is incompatible with $(length(q))-point quadrature")) + return (q.b-q.a)/2 * sum(q.w[i] * f[i] for i in eachindex(q)) # integrate samples of f on [a, b] +end +@inbounds @fastmath function integrate(q::Quadrature, f::Function) + return (q.b-q.a)/2 * sum(q.w[i] * f(node(q, i)) for i in eachindex(q)) # integrate f(x) on [a, b] +end +(q::Quadrature)(args...) = integrate(q, args...) # calling is equivalent to integrating diff --git a/src/solve.jl b/src/solve.jl index 344a22b8..45260520 100644 --- a/src/solve.jl +++ b/src/solve.jl @@ -584,11 +584,6 @@ function issuccess(sol::CosmologySolution) return successful_retcode(sol.bg) && (isnothing(sol.pts) || all(successful_retcode(pt) for pt in sol.pts)) end -function integrate(xs, ys; integrator = TrapezoidalRule()) - prob = SampledIntegralProblem(ys, xs) - sol = solve(prob, integrator) - return sol.u -end integrate_cumulative(sol::CosmologySolution, x, y) = cumul_integrate(sol[x], sol[y]) integrate_cumulative(sol::CosmologySolution, y) = integrate_cumulative(sol, sol.prob.M.τ, y) diff --git a/test/runtests.jl b/test/runtests.jl index 65a6f11e..2ad2714e 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -435,8 +435,12 @@ end jl = SphericalBesselCache(ls) kτ0s_coarse, kτ0s_fine = SymBoltz.cmb_kτ0s(ls[begin], ls[end]) ks_coarse, ks_fine = kτ0s_coarse / τ0, kτ0s_fine / τ0 - ks_fine = clamp.(ks_fine, ks_coarse[begin], ks_coarse[end]) # numerics can change endpoint slightly - τs = range(τi, τ0, length = 768) + τquad = ClenshawCurtisQuadrature(500, (τi, τ0)) + kquad = ClenshawCurtisQuadrature(1000, (ks_coarse[begin], ks_coarse[end])) + τs = nodes(τquad) + ks_fine = nodes(kquad) + ks_fine[begin] = max(ks_fine[begin], ks_coarse[begin]) # possible floating point roundoff error + ks_fine[end] = min(ks_fine[end], ks_coarse[end]) # possible floating point roundoff error #sol = solve(prob, ks_coarse) #Ss = sol(ks_fine, τs, M.ST) ptopts = (alg = SymBoltz.ptalg(prob), reltol = 1e-5, abstol = 1e-5) @@ -447,9 +451,9 @@ end Ss = @inferred source_grid(sol, [M.ST], τs) # TODO: save allocation time with out-of-place version? @test Ss == source_grid(prob, [M.ST], τs, ks_coarse; ptopts) Ss = @inferred source_grid(Ss[1, :, :], ks_coarse, ks_fine) # TODO: save allocation time with out-of-place version? - Θ0s = @inferred los_integrate(Ss, ls, τs, ks_fine, jl) # TODO: sequential along τ? # TODO: cache kτ0 and x=τ/τ0 (only depends on l) + Θ0s = @inferred los_integrate(Ss, ls, τs, ks_fine, jl, τquad) # TODO: sequential along τ? # TODO: cache kτ0 and x=τ/τ0 (only depends on l) P0s = @inferred spectrum_primordial(ks_fine, pars[M.g.h], prob.bg.ps[M.I.As]) - Cls = @inferred spectrum_cmb(Θ0s, Θ0s, P0s, ls, ks_fine) + Cls = @inferred spectrum_cmb(Θ0s, Θ0s, P0s, ls, ks_fine, kquad) end @testset "Background differentiation test" begin @@ -733,6 +737,48 @@ end @test isequal(expandeq(eqs, D(ρ); protect = Set(ℋ)), -4ℋ*ρ) end +@testset "Quadrature rules" begin + f = x -> log(1+x) / (1+x^2) + I = π/8*log(2) # analytical result for ∫f(x)dx from 0 to 1 (https://math.stackexchange.com/a/220754/932179) + + @test_throws ArgumentError TrapezoidalQuadrature(1) + @test TrapezoidalQuadrature(2) == Quadrature([-1.0, 1.0], [1/1, 1/1]) + @test TrapezoidalQuadrature(3) == Quadrature([-1.0, 0.0, 1.0], [1/2, 2/2, 1/2]) + @test TrapezoidalQuadrature(4) == Quadrature([-1.0, -1/3, 1/3, 1.0], [1/3, 2/3, 2/3, 1/3]) + @test TrapezoidalQuadrature(2^13, (0.0, 1.0))(f) ≈ I + + @test_throws ArgumentError TrapezoidalQuadrature([1.0]) + @test_throws ArgumentError TrapezoidalQuadrature([2.0, 1.0]) + @test TrapezoidalQuadrature([1.0, 2.0, 3.0, 4.0]) ≈ TrapezoidalQuadrature(4, (1.0, 4.0)) + @test TrapezoidalQuadrature([-1.0, -0.5, 1.0]) ≈ Quadrature([-1.0, -0.5, 1.0], [0.5/2, 0.5/2+1.5/2, 1.5/2]) + + @test_throws ArgumentError ClenshawCurtisQuadrature(1) + @test ClenshawCurtisQuadrature(2) ≈ Quadrature([-1.0, 1.0], [1.0, 1.0]) + @test ClenshawCurtisQuadrature(3) ≈ Quadrature([-1.0, 0.0, 1.0], [1/3, 4/3, 1/3]) + @test ClenshawCurtisQuadrature(4) ≈ Quadrature([-1, -1/2, 1/2, 1], [1/9, 8/9, 8/9, 1/9]) + @test ClenshawCurtisQuadrature(5) ≈ Quadrature([-1, -√2/2, 0, √2/2, 1], [1/15, 8/15, 4/5, 8/15, 1/15]) + @test ClenshawCurtisQuadrature(2^4, (0.0, 1.0))(f) ≈ I + + # https://numfactory.upc.edu/web/Calculo2/P2_Integracio/html/Gauss1D.html + @test_throws ArgumentError GaussQuadrature(0) + @test GaussQuadrature(1) ≈ Quadrature([0.0], [2.0]) + @test GaussQuadrature(2) ≈ Quadrature([-1/√(3), 1/√(3)], [1.0, 1.0]) + @test GaussQuadrature(3) ≈ Quadrature([-√(3/5), 0, √(3/5)], [5/9, 8/9, 5/9]) + @test GaussQuadrature(2^3, (0.0, 1.0))(f) ≈ I + + @test_throws ArgumentError GaussKronrodQuadrature(0) + @test_throws ArgumentError GaussKronrodQuadrature(1) + @test_throws ArgumentError GaussKronrodQuadrature(2) + @test GaussKronrodQuadrature(3) ≈ GaussQuadrature(3) + @test_throws ArgumentError GaussKronrodQuadrature(4) + @test_nowarn GaussKronrodQuadrature(5) + @test GaussKronrodQuadrature(2^3+1, (0.0, 1.0))(f) ≈ I + + # autodiff + quad = GaussQuadrature(2^3) # TODO: Clenshaw-Curtis? + @test ForwardDiff.derivative(x -> transform(quad, (0.0, x))(f), 1.0) ≈ f(1.0) # derivative of F(x) = ∫₀ˣ dy f(y) = f(x), using fundamental theorem +end + @testset "High lmax" begin M = ΛCDM(lmax = 32) prob = CosmologyProblem(M, pars)