From 8026cde08089010d4c78651885890d14c2f1b329 Mon Sep 17 00:00:00 2001 From: Bart de Koning Date: Mon, 22 Jul 2024 20:32:56 +0200 Subject: [PATCH 1/3] Set up PR --- GaussTuranExampleTemp/example.jl | 63 ++++ .../high_order_derivatives.jl | 57 ++++ GaussTuranExampleTemp/linear_algebra_fix.jl | 292 ++++++++++++++++++ Project.toml | 3 + ext/IntegralsOptimExt.jl | 286 +++++++++++++++++ src/Integrals.jl | 2 + 6 files changed, 703 insertions(+) create mode 100644 GaussTuranExampleTemp/example.jl create mode 100644 GaussTuranExampleTemp/high_order_derivatives.jl create mode 100644 GaussTuranExampleTemp/linear_algebra_fix.jl create mode 100644 ext/IntegralsOptimExt.jl diff --git a/GaussTuranExampleTemp/example.jl b/GaussTuranExampleTemp/example.jl new file mode 100644 index 0000000..3a8a629 --- /dev/null +++ b/GaussTuranExampleTemp/example.jl @@ -0,0 +1,63 @@ +using Optim +using PreallocationTools +using Integrals +using DoubleFloats + +include("high_order_derivatives.jl") + +# This fix is needed to avoid crashing: https://github.com/JuliaLang/julia/pull/54201 +include("linear_algebra_fix.jl") + +FT = Double64 +n = 5 +s = 1 +a = FT(0.0) +b = FT(1.0) + +""" + Example functions ϕⱼ from [1] +""" +function f(x::T, j::Int)::T where {T<:Number} + pow = if j % 2 == 0 + j / 2 - 1 + else + (j - 1) / 2 - 1 / 3 + end + x^pow +end + +deriv! = @generate_derivs(2) +I = Integrals.GaussTuran( + f, + deriv!, + a, + b, + n, + s; + optimization_options = Optim.Options(; + x_abstol = FT(1e-250), + g_tol = FT(1e-250), + show_trace = true, + show_every = 100, + iterations = 10_000, + ), +) + +# Integration error |I(fⱼ) - ₐ∫ᵇfⱼ(x)dx| first N functions fⱼ +max_int_error_upper = + maximum(abs(I(x -> f(x, j)) - I.cache.rhs_upper[j]) for j = 1:I.cache.N) +@show max_int_error_upper +# max_int_error_upper = 3.0814879110195774e-32 + +# Integration error |I(fⱼ) - ₐ∫ᵇfⱼ(x)dx| last n functions fⱼ +max_int_error_lower = maximum( + abs(I(x -> f(x, j)) - I.cache.rhs_lower[j-I.cache.N]) for + j = (I.cache.N+1):(I.cache.N+I.cache.n) +) +@show max_int_error_lower +# max_int_error_lower = 2.8858134286698342e-30 + +# Example with eˣ +exp_error = abs(I(Base.exp) - (Base.exp(1) - 1)) +@show exp_error; +# exp_error = -6.436621416953051e-18 \ No newline at end of file diff --git a/GaussTuranExampleTemp/high_order_derivatives.jl b/GaussTuranExampleTemp/high_order_derivatives.jl new file mode 100644 index 0000000..8de081a --- /dev/null +++ b/GaussTuranExampleTemp/high_order_derivatives.jl @@ -0,0 +1,57 @@ +using ForwardDiff + +""" + Generate a function `derivs!(out, f, x)` which computes the 0th up to the max_orderth derivative + of the scalar-to-scalar function f at x and stores them in ascending derivative order in `out`. + Hence `out` must be at least of length `max_order + 1`. +""" +macro generate_derivs(max_order::Int) + # Create nested dual number of required depth (arg_0, …, arg_{max_order}) + arg_assignments = [:(arg_0 = x)] + for i = 1:max_order + arg_name = Symbol("arg_$i") + prev_arg_name = Symbol("arg_$(i-1)") + push!( + arg_assignments, + :($arg_name = ForwardDiff.Dual{Val{$i}}($prev_arg_name, one($prev_arg_name))), + ) + end + + # Unpack the results + arg_max = Symbol("arg_$max_order") + res_unpacks = [:(res_0 = f($arg_max))] + for i = 1:max_order + res_name = Symbol("res_$i") + res_prev_name = Symbol("res_$(i-1)") + push!(res_unpacks, :($res_name = only($res_prev_name.partials))) + end + + # Assign the results + out_assignments = Expr[] + for i = 0:max_order + res = Symbol("res_$i") + res_temp = Symbol("$(res)_temp_0") + push!(out_assignments, :($res_temp = $res)) + # Create temporary variables to get to + # res_{i}.value.value.value… + for j = 1:(max_order-i) + res_temp = Symbol("$(res)_temp_$j") + res_temp_prev = Symbol("$(res)_temp_$(j-1)") + push!(out_assignments, :($res_temp = $res_temp_prev.value)) + end + res_temp = Symbol("$(res)_temp_$(max_order - i)") + push!(out_assignments, :(out[$(i + 1)] = $res_temp)) + end + + # Construct the complete function definition + func_def = quote + function derivs!(out, f, x::T)::Nothing where {T<:Number} + $(arg_assignments...) + $(res_unpacks...) + $(out_assignments...) + return nothing + end + end + + return func_def +end \ No newline at end of file diff --git a/GaussTuranExampleTemp/linear_algebra_fix.jl b/GaussTuranExampleTemp/linear_algebra_fix.jl new file mode 100644 index 0000000..1055b8e --- /dev/null +++ b/GaussTuranExampleTemp/linear_algebra_fix.jl @@ -0,0 +1,292 @@ +using LinearAlgebra +using LinearAlgebra: AdjOrTrans, require_one_based_indexing, _ustrip + +# manually hoisting b[j] significantly improves performance as of Dec 2015 +# manually eliding bounds checking significantly improves performance as of Dec 2015 +# replacing repeated references to A.data[j,j] with [Ajj = A.data[j,j] and references to Ajj] +# does not significantly impact performance as of Dec 2015 +# in the transpose and conjugate transpose naive substitution variants, +# accumulating in z rather than b[j,k] significantly improves performance as of Dec 2015 +function LinearAlgebra.generic_trimatdiv!(C::AbstractVecOrMat, uploc, isunitc, tfun::Function, A::AbstractMatrix, B::AbstractVecOrMat) + require_one_based_indexing(C, A, B) + mA, nA = size(A) + m, n = size(B, 1), size(B,2) + if nA != m + throw(DimensionMismatch(lazy"second dimension of left hand side A, $nA, and first dimension of right hand side B, $m, must be equal")) + end + if size(C) != size(B) + throw(DimensionMismatch(lazy"size of output, $(size(C)), does not match size of right hand side, $(size(B))")) + end + iszero(mA) && return C + oA = oneunit(eltype(A)) + @inbounds if uploc == 'U' + if isunitc == 'N' + if tfun === identity + for k in 1:n + amm = A[m,m] + iszero(amm) && throw(SingularException(m)) + Cm = C[m,k] = amm \ B[m,k] + # fill C-column + for i in m-1:-1:1 + C[i,k] = oA \ B[i,k] - _ustrip(A[i,m]) * Cm + end + for j in m-1:-1:1 + ajj = A[j,j] + iszero(ajj) && throw(SingularException(j)) + Cj = C[j,k] = _ustrip(ajj) \ C[j,k] + for i in j-1:-1:1 + C[i,k] -= _ustrip(A[i,j]) * Cj + end + end + end + else # tfun in (adjoint, transpose) + for k in 1:n + for j in 1:m + ajj = A[j,j] + iszero(ajj) && throw(SingularException(j)) + Bj = B[j,k] + for i in 1:j-1 + Bj -= tfun(A[i,j]) * C[i,k] + end + C[j,k] = tfun(ajj) \ Bj + end + end + end + else # isunitc == 'U' + if tfun === identity + for k in 1:n + Cm = C[m,k] = oA \ B[m,k] + # fill C-column + for i in m-1:-1:1 + C[i,k] = oA \ B[i,k] - _ustrip(A[i,m]) * Cm + end + for j in m-1:-1:1 + Cj = C[j,k] + for i in 1:j-1 + C[i,k] -= _ustrip(A[i,j]) * Cj + end + end + end + else # tfun in (adjoint, transpose) + for k in 1:n + for j in 1:m + Bj = B[j,k] + for i in 1:j-1 + Bj -= tfun(A[i,j]) * C[i,k] + end + C[j,k] = oA \ Bj + end + end + end + end + else # uploc == 'L' + if isunitc == 'N' + if tfun === identity + for k in 1:n + a11 = A[1,1] + iszero(a11) && throw(SingularException(1)) + C1 = C[1,k] = a11 \ B[1,k] + # fill C-column + for i in 2:m + C[i,k] = oA \ B[i,k] - _ustrip(A[i,1]) * C1 + end + for j in 2:m + ajj = A[j,j] + iszero(ajj) && throw(SingularException(j)) + Cj = C[j,k] = _ustrip(ajj) \ C[j,k] + for i in j+1:m + C[i,k] -= _ustrip(A[i,j]) * Cj + end + end + end + else # tfun in (adjoint, transpose) + for k in 1:n + for j in m:-1:1 + ajj = A[j,j] + iszero(ajj) && throw(SingularException(j)) + Bj = B[j,k] + for i in j+1:m + Bj -= tfun(A[i,j]) * C[i,k] + end + C[j,k] = tfun(ajj) \ Bj + end + end + end + else # isunitc == 'U' + if tfun === identity + for k in 1:n + C1 = C[1,k] = oA \ B[1,k] + # fill C-column + for i in 2:m + C[i,k] = oA \ B[i,k] - _ustrip(A[i,1]) * C1 + end + for j in 2:m + Cj = C[j,k] + for i in j+1:m + C[i,k] -= _ustrip(A[i,j]) * Cj + end + end + end + else # tfun in (adjoint, transpose) + for k in 1:n + for j in m:-1:1 + Bj = B[j,k] + for i in j+1:m + Bj -= tfun(A[i,j]) * C[i,k] + end + C[j,k] = oA \ Bj + end + end + end + end + end + return C +end +# conjugate cases +function LinearAlgebra.generic_trimatdiv!(C::AbstractVecOrMat, uploc, isunitc, ::Function, xA::AdjOrTrans, B::AbstractVecOrMat) + A = parent(xA) + require_one_based_indexing(C, A, B) + mA, nA = size(A) + m, n = size(B, 1), size(B,2) + if nA != m + throw(DimensionMismatch(lazy"second dimension of left hand side A, $nA, and first dimension of right hand side B, $m, must be equal")) + end + if size(C) != size(B) + throw(DimensionMismatch(lazy"size of output, $(size(C)), does not match size of right hand side, $(size(B))")) + end + iszero(mA) && return C + oA = oneunit(eltype(A)) + @inbounds if uploc == 'U' + if isunitc == 'N' + for k in 1:n + amm = conj(A[m,m]) + iszero(amm) && throw(SingularException(m)) + Cm = C[m,k] = amm \ B[m,k] + # fill C-column + for i in m-1:-1:1 + C[i,k] = oA \ B[i,k] - _ustrip(conj(A[i,m])) * Cm + end + for j in m-1:-1:1 + ajj = conj(A[j,j]) + iszero(ajj) && throw(SingularException(j)) + Cj = C[j,k] = _ustrip(ajj) \ C[j,k] + for i in j-1:-1:1 + C[i,k] -= _ustrip(conj(A[i,j])) * Cj + end + end + end + else # isunitc == 'U' + for k in 1:n + Cm = C[m,k] = oA \ B[m,k] + # fill C-column + for i in m-1:-1:1 + C[i,k] = oA \ B[i,k] - _ustrip(conj(A[i,m])) * Cm + end + for j in m-1:-1:1 + Cj = C[j,k] + for i in 1:j-1 + C[i,k] -= _ustrip(conj(A[i,j])) * Cj + end + end + end + end + else # uploc == 'L' + if isunitc == 'N' + for k in 1:n + a11 = conj(A[1,1]) + iszero(a11) && throw(SingularException(1)) + C1 = C[1,k] = a11 \ B[1,k] + # fill C-column + for i in 2:m + C[i,k] = oA \ B[i,k] - _ustrip(conj(A[i,1])) * C1 + end + for j in 2:m + ajj = conj(A[j,j]) + iszero(ajj) && throw(SingularException(j)) + Cj = C[j,k] = _ustrip(ajj) \ C[j,k] + for i in j+1:m + C[i,k] -= _ustrip(conj(A[i,j])) * Cj + end + end + end + else # isunitc == 'U' + for k in 1:n + C1 = C[1,k] = oA \ B[1,k] + # fill C-column + for i in 2:m + C[i,k] = oA \ B[i,k] - _ustrip(conj(A[i,1])) * C1 + end + for j in 1:m + Cj = C[j,k] + for i in j+1:m + C[i,k] -= _ustrip(conj(A[i,j])) * Cj + end + end + end + end + end + return C +end + +function generic_mattridiv!(C::AbstractMatrix, uploc, isunitc, tfun::Function, A::AbstractMatrix, B::AbstractMatrix) + require_one_based_indexing(C, A, B) + m, n = size(A) + if size(B, 1) != n + throw(DimensionMismatch(lazy"right hand side B needs first dimension of size $n, has size $(size(B,1))")) + end + if size(C) != size(A) + throw(DimensionMismatch(lazy"size of output, $(size(C)), does not match size of left hand side, $(size(A))")) + end + oB = oneunit(eltype(B)) + unit = isunitc == 'U' + @inbounds if uploc == 'U' + if tfun === identity + for i in 1:m + for j in 1:n + Aij = A[i,j] + for k in 1:j - 1 + Aij -= C[i,k]*B[k,j] + end + unit || (iszero(B[j,j]) && throw(SingularException(j))) + C[i,j] = Aij / (unit ? oB : B[j,j]) + end + end + else # tfun in (adjoint, transpose) + for i in 1:m + for j in n:-1:1 + Aij = A[i,j] + for k in j + 1:n + Aij -= C[i,k]*tfun(B[j,k]) + end + unit || (iszero(B[j,j]) && throw(SingularException(j))) + C[i,j] = Aij / (unit ? oB : tfun(B[j,j])) + end + end + end + else # uploc == 'L' + if tfun === identity + for i in 1:m + for j in n:-1:1 + Aij = A[i,j] + for k in j + 1:n + Aij -= C[i,k]*B[k,j] + end + unit || (iszero(B[j,j]) && throw(SingularException(j))) + C[i,j] = Aij / (unit ? oB : B[j,j]) + end + end + else # tfun in (adjoint, transpose) + for i in 1:m + for j in 1:n + Aij = A[i,j] + for k in 1:j - 1 + Aij -= C[i,k]*tfun(B[j,k]) + end + unit || (iszero(B[j,j]) && throw(SingularException(j))) + C[i,j] = Aij / (unit ? oB : tfun(B[j,j])) + end + end + end + end + return C +end \ No newline at end of file diff --git a/Project.toml b/Project.toml index d281ac3..1de41f8 100644 --- a/Project.toml +++ b/Project.toml @@ -21,6 +21,8 @@ Cubature = "667455a9-e2ce-5579-9412-b964f529a492" FastGaussQuadrature = "442a2c76-b920-505d-bb47-c5924d526838" ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" MCIntegration = "ea1e2de9-7db7-4b42-91ee-0cd1bf6df167" +Optim = "429524aa-4258-5aef-a3af-852621145aeb" +PreallocationTools = "d236fae5-4411-538c-8e31-a6e3d9e00b46" Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f" [extensions] @@ -29,6 +31,7 @@ IntegralsCubaExt = "Cuba" IntegralsCubatureExt = "Cubature" IntegralsFastGaussQuadratureExt = "FastGaussQuadrature" IntegralsForwardDiffExt = "ForwardDiff" +IntegralsOptimExt = ["Optim", "PreallocationTools", "ForwardDiff"] IntegralsMCIntegrationExt = "MCIntegration" IntegralsZygoteExt = ["Zygote", "ChainRulesCore"] diff --git a/ext/IntegralsOptimExt.jl b/ext/IntegralsOptimExt.jl new file mode 100644 index 0000000..520ad24 --- /dev/null +++ b/ext/IntegralsOptimExt.jl @@ -0,0 +1,286 @@ + +module IntegralsOptimExt +using Base.Threads +using Integrals +if isdefined(Base, :get_extension) + using Optim + using ForwardDiff + using PreallocationTools +else + using ..Optim + using ..ForwardDiff + using ..PreallocationTools +end + +######################################################## +## Computing Gauss-Turán quadrature rules + +function DEFAULT_w(x::T)::T where {T} + one(T) +end + +""" +Cached data for the `GaussTuranLoss!` call. +""" +struct GaussTuranCache{T} + n::Int + s::Int + N::Int + a::T + b::T + ε::T + rhs_upper::Vector{T} + rhs_lower::Vector{T} + M_upper_buffer::LazyBufferCache{typeof(identity)} + M_lower_buffer::LazyBufferCache{typeof(identity)} + X_buffer::LazyBufferCache{typeof(identity)} + A_buffer::LazyBufferCache{typeof(identity)} + function GaussTuranCache( + n, + s, + N, + a::T, + b::T, + ε::T, + rhs_upper::Vector{T}, + rhs_lower::Vector{T}, + ) where {T} + new{T}( + n, + s, + N, + a, + b, + ε, + rhs_upper, + rhs_lower, + LazyBufferCache(), + LazyBufferCache(), + LazyBufferCache(), + LazyBufferCache(), + ) + end +end + +""" +Function whose root defines the quadrature rule. +""" +function GaussTuranLoss!(f!, ΔX, cache) + (; + n, + s, + N, + a, + rhs_upper, + rhs_lower, + M_upper_buffer, + M_lower_buffer, + A_buffer, + X_buffer, + ) = cache + M_upper = M_upper_buffer[ΔX, (N, N)] + M_lower = M_lower_buffer[ΔX, (n, N)] + A = A_buffer[ΔX, N] + X = X_buffer[ΔX, n] + + # Compute X from ΔX + cumsum!(X, ΔX) + X .+= a + + # Evaluating f! + for (i, x) in enumerate(X) + Threads.@threads for j = 1:N + f!(view(M_upper, j, i:n:N), x, j) + end + Threads.@threads for j = (N+1):(N+n) + f!(view(M_lower, j - N, i:n:N), x, j) + end + end + + # Solving for A + A .= M_upper \ rhs_upper + + # Computing output + out = zero(eltype(ΔX)) + for i in eachindex(X) + out_term = -rhs_lower[i] + for j in eachindex(A) + out_term += A[j] * M_lower[i, j] + end + out += out_term^2 + end + sqrt(out) +end + +""" + Callable result object of the Gauss-Turán quadrature rule + computation algorithm. +""" +struct GaussTuranResult{T,RType,dfType,deriv!Type} + X::Vector{T} + A::Matrix{T} + res::RType + cache::GaussTuranCache + df::dfType + deriv!::deriv!Type + f_tmp::Vector{T} + + function GaussTuranResult(res, cache::GaussTuranCache{T}, df, deriv!) where {T} + (; A_buffer, s, n, N, a) = cache + X = cumsum(res.minimizer) .+ a + df.f(res.minimizer) + A = reshape(A_buffer[T[], N], (n, 2s + 1)) + f_tmp = zeros(T, 2s + 1) + new{T,typeof(res),typeof(df),typeof(deriv!)}(X, A, res, cache, df, deriv!, f_tmp) + end +end + +""" + Input: function f(x, d) which gives the dth derivative of f +""" +function (I::GaussTuranResult{T} where {T})(integrand) + (; X, A, cache, deriv!, f_tmp) = I + out = zero(eltype(X)) + for (i, x) in enumerate(X) + deriv!(f_tmp, integrand, x) + for m = 1:(2*cache.s+1) + out += A[i, m] * f_tmp[m] + end + end + out +end + +""" + GaussTuran(f!, a, b, n, s; w = DEFAULT_w, ε = nothing, X₀ = nothing) + +Determine a quadrature rule + +I(g) = ∑ₘ∑ᵢ Aₘᵢ * ∂ᵐ⁻¹g(xᵢ) (m = 1, … 2s + 1, i = 1, …, n) + +that gives the precise integral ₐ∫ᵇf(x)dx for given linearly independent functions f₁, f₂, …, f₂₍ₛ₊₁₎ₙ. + +Method: + +The equations + +∑ₘ∑ᵢ Aₘᵢ * ∂ᵐ⁻¹fⱼ(xᵢ) = ₐ∫ᵇw(x)fⱼ(x)dx j = 1, …, 2(s+1)n + +define an overdetermined linear system M(X)A = b in the weights Aₘᵢ for a given X = (x₁, x₂, …, xₙ). +We split the matrix M into a square upper part M_upper of size (2s+1)n x (2s+1)n and a lower part M_lower of size n x (2s+1)n, +and the right hand size b analogously. From this we obtain A = M_upper⁻¹ * b_upper. Then we can asses the correctness of X by comparing +M_lower * A to b_lower, i.e. how well the last n equations holds. This yields the loss function + +loss(X) = ||M_lower * A - b_lower||₂ = ||M_lower * M_upper⁻¹ * b_upper - b_lower||₂. + +We have the constraints that we want X to be ordered and in the interval (a, b). To achieve this, we formulate the loss +in terms of ΔX = (Δx₁, Δx₂, …, Δxₙ) = (x₁ - a, x₂ - x₁, …, xₙ - xₙ₋₁) on which we set the constraints + +ε ≤ Δxᵢ ≤ b - a - 2 * ε i = 1, …, n +nε ≤ a + ∑ΔX ≤ b - a - ε + +where ε is an enforced minimum distance between the nodes. This prevents that consecutive nodes converge towards eachother making +M_upper singular. + +## Inputs + + - `f`: Function with signature `f(x::T, j)::T` that returns the d-th derivative of fⱼ at x + - `deriv!`: Function generated with @generate_derivs(2s) for automatically computing the derivatives of the fⱼ + - `a`: Integration lower bound + - `b`: Integration upper bound + - `n`: The number of nodes in the quadrature rule + - `s`: Determines the highest order derivative required from the functions fⱼ, currently 2(s + 1) + +## Keyword Arguments + + - `w`: the integrand weighting function, must have signature w(x::Number)::Number. Defaults to `w(x) = 1`. + - `ε`: the minimum distance between nodes. Defaults to 1e-3 * (b - a) / (n + 1). + - `X₀`: The initial guess for the nodes. Defaults to uniformly distributed over (a, b). + - `integration_kwargs`: The key word arguments passed to `solve` for integrating w * fⱼ + - `optimization_kwargs`: The key word arguments passed to `Optim.Options` for the minization problem + for finding X. +""" +function Integrals.GaussTuran( + f, + deriv!, + a::T, + b::T, + n, + s; + w = DEFAULT_w, + ε = nothing, + X₀ = nothing, + integration_kwargs::NamedTuple = (; reltol = 1e-120), + optimization_options::Optim.Options = Optim.Options(), +) where {T<:AbstractFloat} + # Initial guess + if isnothing(X₀) + X₀ = collect(range(a, b, length = n + 2)[2:(end-1)]) + else + @assert length(X₀) == n + end + ΔX₀ = diff(X₀) + pushfirst!(ΔX₀, X₀[1] - a) + + # Minimum distance between nodes + if isnothing(ε) + ε = 1e-3 * (b - a) / (n + 1) + else + @assert 0 < ε ≤ (b - a) / (n + 1) + end + ε = T(ε) + + # Integrate w * f + integrand = (out, x, j) -> out[] = w(x) * f(x, j) + function integrate(j) + prob = IntegralProblem{true}(integrand, (a, b), j) + res = solve(prob, QuadGKJL(); integration_kwargs...) + res.u[] + end + N = (2s + 1) * n + rhs_upper = [integrate(j) for j = 1:N] + rhs_lower = [integrate(j) for j = (N+1):(N+n)] + + # Solving constrained non linear problem for ΔX, see + # https://julianlsolvers.github.io/Optim.jl/stable/examples/generated/ipnewton_basics/ + + # The cache for evaluating GaussTuranLoss + cache = GaussTuranCache(n, s, N, a, b, ε, rhs_upper, rhs_lower) + + # In-place derivative computation of f + function f!(out, x, j::Int)::Nothing + deriv!(out, x -> f(x, j), x) + end + + # The function whose root defines the quadrature rule + # Note: the optimization method requires a Hessian, + # which brings the highest order derivative required to 2s + 2 + func(ΔX) = GaussTuranLoss!(f!, ΔX, cache) + df = TwiceDifferentiable(func, ΔX₀; autodiff = :forward) + + # The constraints on ΔX + ΔX_lb = fill(ε, length(ΔX₀)) + ΔX_ub = fill(b - a - 2 * ε, length(ΔX₀)) + + # Defining the variable and constraints nε ≤ a + ∑ΔX ≤ b - a - ε + sum_variable!(c, ΔX) = (c[1] = a + sum(ΔX); c) + sum_jacobian!(J, ΔX) = (J[1, :] .= one(eltype(ΔX)); J) + sum_hessian!(H, ΔX, λ) = nothing + sum_lb = [n * ε] + sum_ub = [b - a - ε] + constraints = TwiceDifferentiableConstraints( + sum_variable!, + sum_jacobian!, + sum_hessian!, + ΔX_lb, + ΔX_ub, + sum_lb, + sum_ub, + ) + + # Solve for the quadrature rule by minimizing the loss function + res = Optim.optimize(df, constraints, T.(ΔX₀), IPNewton(), optimization_options) + + GaussTuranResult(res, cache, df, deriv!) +end + +end # module IntegralsOptimExt \ No newline at end of file diff --git a/src/Integrals.jl b/src/Integrals.jl index d11fab9..715b0f1 100644 --- a/src/Integrals.jl +++ b/src/Integrals.jl @@ -20,6 +20,8 @@ include("sampled.jl") include("trapezoidal.jl") include("simpsons.jl") +function GaussTuran end + abstract type QuadSensitivityAlg end struct ReCallVJP{V} vjp::V From b60c6a467ccb8fce4235c5a4422f89651e73c1d2 Mon Sep 17 00:00:00 2001 From: Bart de Koning Date: Mon, 22 Jul 2024 22:36:49 +0200 Subject: [PATCH 2/3] Use TaylorDiff --- GaussTuranExampleTemp/example.jl | 28 ++- .../high_order_derivatives.jl | 57 ----- GaussTuranExampleTemp/linear_algebra_fix.jl | 204 +++++++++--------- Project.toml | 3 +- ...sOptimExt.jl => IntegralsGaussTuranExt.jl} | 126 +++++------ src/algorithms_sampled.jl | 4 +- 6 files changed, 180 insertions(+), 242 deletions(-) delete mode 100644 GaussTuranExampleTemp/high_order_derivatives.jl rename ext/{IntegralsOptimExt.jl => IntegralsGaussTuranExt.jl} (76%) diff --git a/GaussTuranExampleTemp/example.jl b/GaussTuranExampleTemp/example.jl index 3a8a629..c57dc44 100644 --- a/GaussTuranExampleTemp/example.jl +++ b/GaussTuranExampleTemp/example.jl @@ -1,14 +1,13 @@ using Optim +using TaylorDiff using PreallocationTools using Integrals using DoubleFloats -include("high_order_derivatives.jl") - # This fix is needed to avoid crashing: https://github.com/JuliaLang/julia/pull/54201 include("linear_algebra_fix.jl") -FT = Double64 +FT = Float64 # Double64 n = 5 s = 1 a = FT(0.0) @@ -17,7 +16,7 @@ b = FT(1.0) """ Example functions ϕⱼ from [1] """ -function f(x::T, j::Int)::T where {T<:Number} +function f(x::T, j::Int)::T where {T <: Number} pow = if j % 2 == 0 j / 2 - 1 else @@ -26,10 +25,8 @@ function f(x::T, j::Int)::T where {T<:Number} x^pow end -deriv! = @generate_derivs(2) I = Integrals.GaussTuran( f, - deriv!, a, b, n, @@ -39,25 +36,26 @@ I = Integrals.GaussTuran( g_tol = FT(1e-250), show_trace = true, show_every = 100, - iterations = 10_000, - ), + iterations = 10_000 + ) ) # Integration error |I(fⱼ) - ₐ∫ᵇfⱼ(x)dx| first N functions fⱼ -max_int_error_upper = - maximum(abs(I(x -> f(x, j)) - I.cache.rhs_upper[j]) for j = 1:I.cache.N) +max_int_error_upper = maximum(abs(I(x -> f(x, j)) - I.cache.rhs_upper[j]) +for j in 1:(I.cache.N)) @show max_int_error_upper -# max_int_error_upper = 3.0814879110195774e-32 +# max_int_error_upper = 2.220446049250313e-16 # Integration error |I(fⱼ) - ₐ∫ᵇfⱼ(x)dx| last n functions fⱼ max_int_error_lower = maximum( - abs(I(x -> f(x, j)) - I.cache.rhs_lower[j-I.cache.N]) for - j = (I.cache.N+1):(I.cache.N+I.cache.n) + abs(I(x -> f(x, j)) - I.cache.rhs_lower[j - I.cache.N]) +for +j in (I.cache.N + 1):(I.cache.N + I.cache.n) ) @show max_int_error_lower -# max_int_error_lower = 2.8858134286698342e-30 +# max_int_error_lower = 2.0211055051788662e-12 # Example with eˣ exp_error = abs(I(Base.exp) - (Base.exp(1) - 1)) @show exp_error; -# exp_error = -6.436621416953051e-18 \ No newline at end of file +# exp_error = 1.5543122344752192e-15 diff --git a/GaussTuranExampleTemp/high_order_derivatives.jl b/GaussTuranExampleTemp/high_order_derivatives.jl deleted file mode 100644 index 8de081a..0000000 --- a/GaussTuranExampleTemp/high_order_derivatives.jl +++ /dev/null @@ -1,57 +0,0 @@ -using ForwardDiff - -""" - Generate a function `derivs!(out, f, x)` which computes the 0th up to the max_orderth derivative - of the scalar-to-scalar function f at x and stores them in ascending derivative order in `out`. - Hence `out` must be at least of length `max_order + 1`. -""" -macro generate_derivs(max_order::Int) - # Create nested dual number of required depth (arg_0, …, arg_{max_order}) - arg_assignments = [:(arg_0 = x)] - for i = 1:max_order - arg_name = Symbol("arg_$i") - prev_arg_name = Symbol("arg_$(i-1)") - push!( - arg_assignments, - :($arg_name = ForwardDiff.Dual{Val{$i}}($prev_arg_name, one($prev_arg_name))), - ) - end - - # Unpack the results - arg_max = Symbol("arg_$max_order") - res_unpacks = [:(res_0 = f($arg_max))] - for i = 1:max_order - res_name = Symbol("res_$i") - res_prev_name = Symbol("res_$(i-1)") - push!(res_unpacks, :($res_name = only($res_prev_name.partials))) - end - - # Assign the results - out_assignments = Expr[] - for i = 0:max_order - res = Symbol("res_$i") - res_temp = Symbol("$(res)_temp_0") - push!(out_assignments, :($res_temp = $res)) - # Create temporary variables to get to - # res_{i}.value.value.value… - for j = 1:(max_order-i) - res_temp = Symbol("$(res)_temp_$j") - res_temp_prev = Symbol("$(res)_temp_$(j-1)") - push!(out_assignments, :($res_temp = $res_temp_prev.value)) - end - res_temp = Symbol("$(res)_temp_$(max_order - i)") - push!(out_assignments, :(out[$(i + 1)] = $res_temp)) - end - - # Construct the complete function definition - func_def = quote - function derivs!(out, f, x::T)::Nothing where {T<:Number} - $(arg_assignments...) - $(res_unpacks...) - $(out_assignments...) - return nothing - end - end - - return func_def -end \ No newline at end of file diff --git a/GaussTuranExampleTemp/linear_algebra_fix.jl b/GaussTuranExampleTemp/linear_algebra_fix.jl index 1055b8e..1361bb7 100644 --- a/GaussTuranExampleTemp/linear_algebra_fix.jl +++ b/GaussTuranExampleTemp/linear_algebra_fix.jl @@ -7,10 +7,12 @@ using LinearAlgebra: AdjOrTrans, require_one_based_indexing, _ustrip # does not significantly impact performance as of Dec 2015 # in the transpose and conjugate transpose naive substitution variants, # accumulating in z rather than b[j,k] significantly improves performance as of Dec 2015 -function LinearAlgebra.generic_trimatdiv!(C::AbstractVecOrMat, uploc, isunitc, tfun::Function, A::AbstractMatrix, B::AbstractVecOrMat) +function LinearAlgebra.generic_trimatdiv!( + C::AbstractVecOrMat, uploc, isunitc, tfun::Function, + A::AbstractMatrix, B::AbstractVecOrMat) require_one_based_indexing(C, A, B) mA, nA = size(A) - m, n = size(B, 1), size(B,2) + m, n = size(B, 1), size(B, 2) if nA != m throw(DimensionMismatch(lazy"second dimension of left hand side A, $nA, and first dimension of right hand side B, $m, must be equal")) end @@ -23,58 +25,58 @@ function LinearAlgebra.generic_trimatdiv!(C::AbstractVecOrMat, uploc, isunitc, t if isunitc == 'N' if tfun === identity for k in 1:n - amm = A[m,m] + amm = A[m, m] iszero(amm) && throw(SingularException(m)) - Cm = C[m,k] = amm \ B[m,k] + Cm = C[m, k] = amm \ B[m, k] # fill C-column - for i in m-1:-1:1 - C[i,k] = oA \ B[i,k] - _ustrip(A[i,m]) * Cm + for i in (m - 1):-1:1 + C[i, k] = oA \ B[i, k] - _ustrip(A[i, m]) * Cm end - for j in m-1:-1:1 - ajj = A[j,j] + for j in (m - 1):-1:1 + ajj = A[j, j] iszero(ajj) && throw(SingularException(j)) - Cj = C[j,k] = _ustrip(ajj) \ C[j,k] - for i in j-1:-1:1 - C[i,k] -= _ustrip(A[i,j]) * Cj + Cj = C[j, k] = _ustrip(ajj) \ C[j, k] + for i in (j - 1):-1:1 + C[i, k] -= _ustrip(A[i, j]) * Cj end end end else # tfun in (adjoint, transpose) for k in 1:n for j in 1:m - ajj = A[j,j] + ajj = A[j, j] iszero(ajj) && throw(SingularException(j)) - Bj = B[j,k] - for i in 1:j-1 - Bj -= tfun(A[i,j]) * C[i,k] + Bj = B[j, k] + for i in 1:(j - 1) + Bj -= tfun(A[i, j]) * C[i, k] end - C[j,k] = tfun(ajj) \ Bj + C[j, k] = tfun(ajj) \ Bj end end end else # isunitc == 'U' if tfun === identity for k in 1:n - Cm = C[m,k] = oA \ B[m,k] + Cm = C[m, k] = oA \ B[m, k] # fill C-column - for i in m-1:-1:1 - C[i,k] = oA \ B[i,k] - _ustrip(A[i,m]) * Cm + for i in (m - 1):-1:1 + C[i, k] = oA \ B[i, k] - _ustrip(A[i, m]) * Cm end - for j in m-1:-1:1 - Cj = C[j,k] - for i in 1:j-1 - C[i,k] -= _ustrip(A[i,j]) * Cj + for j in (m - 1):-1:1 + Cj = C[j, k] + for i in 1:(j - 1) + C[i, k] -= _ustrip(A[i, j]) * Cj end end end else # tfun in (adjoint, transpose) for k in 1:n for j in 1:m - Bj = B[j,k] - for i in 1:j-1 - Bj -= tfun(A[i,j]) * C[i,k] + Bj = B[j, k] + for i in 1:(j - 1) + Bj -= tfun(A[i, j]) * C[i, k] end - C[j,k] = oA \ Bj + C[j, k] = oA \ Bj end end end @@ -83,58 +85,58 @@ function LinearAlgebra.generic_trimatdiv!(C::AbstractVecOrMat, uploc, isunitc, t if isunitc == 'N' if tfun === identity for k in 1:n - a11 = A[1,1] + a11 = A[1, 1] iszero(a11) && throw(SingularException(1)) - C1 = C[1,k] = a11 \ B[1,k] + C1 = C[1, k] = a11 \ B[1, k] # fill C-column for i in 2:m - C[i,k] = oA \ B[i,k] - _ustrip(A[i,1]) * C1 + C[i, k] = oA \ B[i, k] - _ustrip(A[i, 1]) * C1 end for j in 2:m - ajj = A[j,j] + ajj = A[j, j] iszero(ajj) && throw(SingularException(j)) - Cj = C[j,k] = _ustrip(ajj) \ C[j,k] - for i in j+1:m - C[i,k] -= _ustrip(A[i,j]) * Cj + Cj = C[j, k] = _ustrip(ajj) \ C[j, k] + for i in (j + 1):m + C[i, k] -= _ustrip(A[i, j]) * Cj end end end else # tfun in (adjoint, transpose) for k in 1:n for j in m:-1:1 - ajj = A[j,j] + ajj = A[j, j] iszero(ajj) && throw(SingularException(j)) - Bj = B[j,k] - for i in j+1:m - Bj -= tfun(A[i,j]) * C[i,k] + Bj = B[j, k] + for i in (j + 1):m + Bj -= tfun(A[i, j]) * C[i, k] end - C[j,k] = tfun(ajj) \ Bj + C[j, k] = tfun(ajj) \ Bj end end end else # isunitc == 'U' if tfun === identity for k in 1:n - C1 = C[1,k] = oA \ B[1,k] + C1 = C[1, k] = oA \ B[1, k] # fill C-column for i in 2:m - C[i,k] = oA \ B[i,k] - _ustrip(A[i,1]) * C1 + C[i, k] = oA \ B[i, k] - _ustrip(A[i, 1]) * C1 end for j in 2:m - Cj = C[j,k] - for i in j+1:m - C[i,k] -= _ustrip(A[i,j]) * Cj + Cj = C[j, k] + for i in (j + 1):m + C[i, k] -= _ustrip(A[i, j]) * Cj end end end else # tfun in (adjoint, transpose) for k in 1:n for j in m:-1:1 - Bj = B[j,k] - for i in j+1:m - Bj -= tfun(A[i,j]) * C[i,k] + Bj = B[j, k] + for i in (j + 1):m + Bj -= tfun(A[i, j]) * C[i, k] end - C[j,k] = oA \ Bj + C[j, k] = oA \ Bj end end end @@ -143,11 +145,12 @@ function LinearAlgebra.generic_trimatdiv!(C::AbstractVecOrMat, uploc, isunitc, t return C end # conjugate cases -function LinearAlgebra.generic_trimatdiv!(C::AbstractVecOrMat, uploc, isunitc, ::Function, xA::AdjOrTrans, B::AbstractVecOrMat) +function LinearAlgebra.generic_trimatdiv!(C::AbstractVecOrMat, uploc, isunitc, ::Function, + xA::AdjOrTrans, B::AbstractVecOrMat) A = parent(xA) require_one_based_indexing(C, A, B) mA, nA = size(A) - m, n = size(B, 1), size(B,2) + m, n = size(B, 1), size(B, 2) if nA != m throw(DimensionMismatch(lazy"second dimension of left hand side A, $nA, and first dimension of right hand side B, $m, must be equal")) end @@ -159,33 +162,33 @@ function LinearAlgebra.generic_trimatdiv!(C::AbstractVecOrMat, uploc, isunitc, : @inbounds if uploc == 'U' if isunitc == 'N' for k in 1:n - amm = conj(A[m,m]) + amm = conj(A[m, m]) iszero(amm) && throw(SingularException(m)) - Cm = C[m,k] = amm \ B[m,k] + Cm = C[m, k] = amm \ B[m, k] # fill C-column - for i in m-1:-1:1 - C[i,k] = oA \ B[i,k] - _ustrip(conj(A[i,m])) * Cm + for i in (m - 1):-1:1 + C[i, k] = oA \ B[i, k] - _ustrip(conj(A[i, m])) * Cm end - for j in m-1:-1:1 - ajj = conj(A[j,j]) + for j in (m - 1):-1:1 + ajj = conj(A[j, j]) iszero(ajj) && throw(SingularException(j)) - Cj = C[j,k] = _ustrip(ajj) \ C[j,k] - for i in j-1:-1:1 - C[i,k] -= _ustrip(conj(A[i,j])) * Cj + Cj = C[j, k] = _ustrip(ajj) \ C[j, k] + for i in (j - 1):-1:1 + C[i, k] -= _ustrip(conj(A[i, j])) * Cj end end end else # isunitc == 'U' for k in 1:n - Cm = C[m,k] = oA \ B[m,k] + Cm = C[m, k] = oA \ B[m, k] # fill C-column - for i in m-1:-1:1 - C[i,k] = oA \ B[i,k] - _ustrip(conj(A[i,m])) * Cm + for i in (m - 1):-1:1 + C[i, k] = oA \ B[i, k] - _ustrip(conj(A[i, m])) * Cm end - for j in m-1:-1:1 - Cj = C[j,k] - for i in 1:j-1 - C[i,k] -= _ustrip(conj(A[i,j])) * Cj + for j in (m - 1):-1:1 + Cj = C[j, k] + for i in 1:(j - 1) + C[i, k] -= _ustrip(conj(A[i, j])) * Cj end end end @@ -193,33 +196,33 @@ function LinearAlgebra.generic_trimatdiv!(C::AbstractVecOrMat, uploc, isunitc, : else # uploc == 'L' if isunitc == 'N' for k in 1:n - a11 = conj(A[1,1]) + a11 = conj(A[1, 1]) iszero(a11) && throw(SingularException(1)) - C1 = C[1,k] = a11 \ B[1,k] + C1 = C[1, k] = a11 \ B[1, k] # fill C-column for i in 2:m - C[i,k] = oA \ B[i,k] - _ustrip(conj(A[i,1])) * C1 + C[i, k] = oA \ B[i, k] - _ustrip(conj(A[i, 1])) * C1 end for j in 2:m - ajj = conj(A[j,j]) + ajj = conj(A[j, j]) iszero(ajj) && throw(SingularException(j)) - Cj = C[j,k] = _ustrip(ajj) \ C[j,k] - for i in j+1:m - C[i,k] -= _ustrip(conj(A[i,j])) * Cj + Cj = C[j, k] = _ustrip(ajj) \ C[j, k] + for i in (j + 1):m + C[i, k] -= _ustrip(conj(A[i, j])) * Cj end end end else # isunitc == 'U' for k in 1:n - C1 = C[1,k] = oA \ B[1,k] + C1 = C[1, k] = oA \ B[1, k] # fill C-column for i in 2:m - C[i,k] = oA \ B[i,k] - _ustrip(conj(A[i,1])) * C1 + C[i, k] = oA \ B[i, k] - _ustrip(conj(A[i, 1])) * C1 end for j in 1:m - Cj = C[j,k] - for i in j+1:m - C[i,k] -= _ustrip(conj(A[i,j])) * Cj + Cj = C[j, k] + for i in (j + 1):m + C[i, k] -= _ustrip(conj(A[i, j])) * Cj end end end @@ -228,7 +231,8 @@ function LinearAlgebra.generic_trimatdiv!(C::AbstractVecOrMat, uploc, isunitc, : return C end -function generic_mattridiv!(C::AbstractMatrix, uploc, isunitc, tfun::Function, A::AbstractMatrix, B::AbstractMatrix) +function generic_mattridiv!(C::AbstractMatrix, uploc, isunitc, tfun::Function, + A::AbstractMatrix, B::AbstractMatrix) require_one_based_indexing(C, A, B) m, n = size(A) if size(B, 1) != n @@ -243,23 +247,23 @@ function generic_mattridiv!(C::AbstractMatrix, uploc, isunitc, tfun::Function, A if tfun === identity for i in 1:m for j in 1:n - Aij = A[i,j] - for k in 1:j - 1 - Aij -= C[i,k]*B[k,j] + Aij = A[i, j] + for k in 1:(j - 1) + Aij -= C[i, k] * B[k, j] end - unit || (iszero(B[j,j]) && throw(SingularException(j))) - C[i,j] = Aij / (unit ? oB : B[j,j]) + unit || (iszero(B[j, j]) && throw(SingularException(j))) + C[i, j] = Aij / (unit ? oB : B[j, j]) end end else # tfun in (adjoint, transpose) for i in 1:m for j in n:-1:1 - Aij = A[i,j] - for k in j + 1:n - Aij -= C[i,k]*tfun(B[j,k]) + Aij = A[i, j] + for k in (j + 1):n + Aij -= C[i, k] * tfun(B[j, k]) end - unit || (iszero(B[j,j]) && throw(SingularException(j))) - C[i,j] = Aij / (unit ? oB : tfun(B[j,j])) + unit || (iszero(B[j, j]) && throw(SingularException(j))) + C[i, j] = Aij / (unit ? oB : tfun(B[j, j])) end end end @@ -267,26 +271,26 @@ function generic_mattridiv!(C::AbstractMatrix, uploc, isunitc, tfun::Function, A if tfun === identity for i in 1:m for j in n:-1:1 - Aij = A[i,j] - for k in j + 1:n - Aij -= C[i,k]*B[k,j] + Aij = A[i, j] + for k in (j + 1):n + Aij -= C[i, k] * B[k, j] end - unit || (iszero(B[j,j]) && throw(SingularException(j))) - C[i,j] = Aij / (unit ? oB : B[j,j]) + unit || (iszero(B[j, j]) && throw(SingularException(j))) + C[i, j] = Aij / (unit ? oB : B[j, j]) end end else # tfun in (adjoint, transpose) for i in 1:m for j in 1:n - Aij = A[i,j] - for k in 1:j - 1 - Aij -= C[i,k]*tfun(B[j,k]) + Aij = A[i, j] + for k in 1:(j - 1) + Aij -= C[i, k] * tfun(B[j, k]) end - unit || (iszero(B[j,j]) && throw(SingularException(j))) - C[i,j] = Aij / (unit ? oB : tfun(B[j,j])) + unit || (iszero(B[j, j]) && throw(SingularException(j))) + C[i, j] = Aij / (unit ? oB : tfun(B[j, j])) end end end end return C -end \ No newline at end of file +end diff --git a/Project.toml b/Project.toml index 1de41f8..25a276b 100644 --- a/Project.toml +++ b/Project.toml @@ -23,6 +23,7 @@ ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" MCIntegration = "ea1e2de9-7db7-4b42-91ee-0cd1bf6df167" Optim = "429524aa-4258-5aef-a3af-852621145aeb" PreallocationTools = "d236fae5-4411-538c-8e31-a6e3d9e00b46" +TaylorDiff = "b36ab563-344f-407b-a36a-4f200bebf99c" Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f" [extensions] @@ -31,8 +32,8 @@ IntegralsCubaExt = "Cuba" IntegralsCubatureExt = "Cubature" IntegralsFastGaussQuadratureExt = "FastGaussQuadrature" IntegralsForwardDiffExt = "ForwardDiff" -IntegralsOptimExt = ["Optim", "PreallocationTools", "ForwardDiff"] IntegralsMCIntegrationExt = "MCIntegration" +IntegralsGaussTuranExt = ["Optim", "TaylorDiff", "PreallocationTools"] IntegralsZygoteExt = ["Zygote", "ChainRulesCore"] [compat] diff --git a/ext/IntegralsOptimExt.jl b/ext/IntegralsGaussTuranExt.jl similarity index 76% rename from ext/IntegralsOptimExt.jl rename to ext/IntegralsGaussTuranExt.jl index 520ad24..b06e693 100644 --- a/ext/IntegralsOptimExt.jl +++ b/ext/IntegralsGaussTuranExt.jl @@ -1,14 +1,14 @@ -module IntegralsOptimExt +module IntegralsGaussTuranExt using Base.Threads using Integrals if isdefined(Base, :get_extension) using Optim - using ForwardDiff + using TaylorDiff using PreallocationTools else using ..Optim - using ..ForwardDiff + using ..TaylorDiff using ..PreallocationTools end @@ -36,14 +36,14 @@ struct GaussTuranCache{T} X_buffer::LazyBufferCache{typeof(identity)} A_buffer::LazyBufferCache{typeof(identity)} function GaussTuranCache( - n, - s, - N, - a::T, - b::T, - ε::T, - rhs_upper::Vector{T}, - rhs_lower::Vector{T}, + n, + s, + N, + a::T, + b::T, + ε::T, + rhs_upper::Vector{T}, + rhs_lower::Vector{T} ) where {T} new{T}( n, @@ -57,7 +57,7 @@ struct GaussTuranCache{T} LazyBufferCache(), LazyBufferCache(), LazyBufferCache(), - LazyBufferCache(), + LazyBufferCache() ) end end @@ -65,19 +65,19 @@ end """ Function whose root defines the quadrature rule. """ -function GaussTuranLoss!(f!, ΔX, cache) +function GaussTuranLoss!(f, ΔX::AbstractVector{T}, cache) where {T} (; - n, - s, - N, - a, - rhs_upper, - rhs_lower, - M_upper_buffer, - M_lower_buffer, - A_buffer, - X_buffer, - ) = cache + n, + s, + N, + a, + rhs_upper, + rhs_lower, + M_upper_buffer, + M_lower_buffer, + A_buffer, + X_buffer +) = cache M_upper = M_upper_buffer[ΔX, (N, N)] M_lower = M_lower_buffer[ΔX, (n, N)] A = A_buffer[ΔX, N] @@ -89,11 +89,12 @@ function GaussTuranLoss!(f!, ΔX, cache) # Evaluating f! for (i, x) in enumerate(X) - Threads.@threads for j = 1:N - f!(view(M_upper, j, i:n:N), x, j) + # Threads.@threads for j = 1:N + for j in 1:N + M_upper[j, i:n:N] .= derivatives(x -> f(x, j), x, one(T), Val(2s + 1)).value end - Threads.@threads for j = (N+1):(N+n) - f!(view(M_lower, j - N, i:n:N), x, j) + Threads.@threads for j in (N + 1):(N + n) + M_lower[j - N, i:n:N] .= derivatives(x -> f(x, j), x, one(T), Val(2s + 1)).value end end @@ -116,22 +117,19 @@ end Callable result object of the Gauss-Turán quadrature rule computation algorithm. """ -struct GaussTuranResult{T,RType,dfType,deriv!Type} +struct GaussTuranResult{T, RType, dfType} X::Vector{T} A::Matrix{T} res::RType cache::GaussTuranCache df::dfType - deriv!::deriv!Type - f_tmp::Vector{T} - function GaussTuranResult(res, cache::GaussTuranCache{T}, df, deriv!) where {T} + function GaussTuranResult(res, cache::GaussTuranCache{T}, df) where {T} (; A_buffer, s, n, N, a) = cache X = cumsum(res.minimizer) .+ a df.f(res.minimizer) A = reshape(A_buffer[T[], N], (n, 2s + 1)) - f_tmp = zeros(T, 2s + 1) - new{T,typeof(res),typeof(df),typeof(deriv!)}(X, A, res, cache, df, deriv!, f_tmp) + new{T, typeof(res), typeof(df)}(X, A, res, cache, df) end end @@ -139,19 +137,20 @@ end Input: function f(x, d) which gives the dth derivative of f """ function (I::GaussTuranResult{T} where {T})(integrand) - (; X, A, cache, deriv!, f_tmp) = I + (; X, A, cache) = I + (; s) = cache out = zero(eltype(X)) for (i, x) in enumerate(X) - deriv!(f_tmp, integrand, x) - for m = 1:(2*cache.s+1) - out += A[i, m] * f_tmp[m] + derivs = derivatives(integrand, x, 1.0, Val(2s + 1)).value + for (m, deriv) in enumerate(derivs) + out += A[i, m] * deriv end end out end """ - GaussTuran(f!, a, b, n, s; w = DEFAULT_w, ε = nothing, X₀ = nothing) + GaussTuran(f, a, b, n, s; w = DEFAULT_w, ε = nothing, X₀ = nothing) Determine a quadrature rule @@ -161,7 +160,7 @@ that gives the precise integral ₐ∫ᵇf(x)dx for given linearly independent f Method: -The equations +The equations ∑ₘ∑ᵢ Aₘᵢ * ∂ᵐ⁻¹fⱼ(xᵢ) = ₐ∫ᵇw(x)fⱼ(x)dx j = 1, …, 2(s+1)n @@ -183,8 +182,7 @@ M_upper singular. ## Inputs - - `f`: Function with signature `f(x::T, j)::T` that returns the d-th derivative of fⱼ at x - - `deriv!`: Function generated with @generate_derivs(2s) for automatically computing the derivatives of the fⱼ + - `f`: Function with signature `f(x::T, j)::T` that returns fⱼ at x - `a`: Integration lower bound - `b`: Integration upper bound - `n`: The number of nodes in the quadrature rule @@ -197,24 +195,23 @@ M_upper singular. - `X₀`: The initial guess for the nodes. Defaults to uniformly distributed over (a, b). - `integration_kwargs`: The key word arguments passed to `solve` for integrating w * fⱼ - `optimization_kwargs`: The key word arguments passed to `Optim.Options` for the minization problem - for finding X. + for finding X. """ function Integrals.GaussTuran( - f, - deriv!, - a::T, - b::T, - n, - s; - w = DEFAULT_w, - ε = nothing, - X₀ = nothing, - integration_kwargs::NamedTuple = (; reltol = 1e-120), - optimization_options::Optim.Options = Optim.Options(), -) where {T<:AbstractFloat} + f, + a::T, + b::T, + n, + s; + w = DEFAULT_w, + ε = nothing, + X₀ = nothing, + integration_kwargs::NamedTuple = (; reltol = 1e-120), + optimization_options::Optim.Options = Optim.Options() +) where {T <: AbstractFloat} # Initial guess if isnothing(X₀) - X₀ = collect(range(a, b, length = n + 2)[2:(end-1)]) + X₀ = collect(range(a, b, length = n + 2)[2:(end - 1)]) else @assert length(X₀) == n end @@ -237,8 +234,8 @@ function Integrals.GaussTuran( res.u[] end N = (2s + 1) * n - rhs_upper = [integrate(j) for j = 1:N] - rhs_lower = [integrate(j) for j = (N+1):(N+n)] + rhs_upper = [integrate(j) for j in 1:N] + rhs_lower = [integrate(j) for j in (N + 1):(N + n)] # Solving constrained non linear problem for ΔX, see # https://julianlsolvers.github.io/Optim.jl/stable/examples/generated/ipnewton_basics/ @@ -246,15 +243,10 @@ function Integrals.GaussTuran( # The cache for evaluating GaussTuranLoss cache = GaussTuranCache(n, s, N, a, b, ε, rhs_upper, rhs_lower) - # In-place derivative computation of f - function f!(out, x, j::Int)::Nothing - deriv!(out, x -> f(x, j), x) - end - # The function whose root defines the quadrature rule # Note: the optimization method requires a Hessian, # which brings the highest order derivative required to 2s + 2 - func(ΔX) = GaussTuranLoss!(f!, ΔX, cache) + func(ΔX) = GaussTuranLoss!(f, ΔX, cache) df = TwiceDifferentiable(func, ΔX₀; autodiff = :forward) # The constraints on ΔX @@ -274,13 +266,13 @@ function Integrals.GaussTuran( ΔX_lb, ΔX_ub, sum_lb, - sum_ub, + sum_ub ) # Solve for the quadrature rule by minimizing the loss function res = Optim.optimize(df, constraints, T.(ΔX₀), IPNewton(), optimization_options) - GaussTuranResult(res, cache, df, deriv!) + GaussTuranResult(res, cache, df) end -end # module IntegralsOptimExt \ No newline at end of file +end # module IntegralsGaussTuranExt diff --git a/src/algorithms_sampled.jl b/src/algorithms_sampled.jl index bc1085f..38ad0c3 100644 --- a/src/algorithms_sampled.jl +++ b/src/algorithms_sampled.jl @@ -10,7 +10,7 @@ Example with sampled data: ```julia using Integrals f = x -> x^2 -x = range(0, 1, length=20) +x = range(0, 1, length = 20) y = f.(x) problem = SampledIntegralProblem(y, x) method = TrapezoidalRule() @@ -31,7 +31,7 @@ Example with equidistant data: ```julia using Integrals f = x -> x^2 -x = range(0, 1, length=20) +x = range(0, 1, length = 20) y = f.(x) problem = SampledIntegralProblem(y, x) method = SimpsonsRule() From 7e0b9ab38629c49a3887ace7fe8c3460cf2a58d6 Mon Sep 17 00:00:00 2001 From: Bart de Koning Date: Tue, 23 Jul 2024 00:03:23 +0200 Subject: [PATCH 3/3] Fix using TaylorDiff with DoubleFloats --- GaussTuranExampleTemp/example.jl | 10 +++++----- Project.toml | 2 +- ext/IntegralsGaussTuranExt.jl | 6 +++--- 3 files changed, 9 insertions(+), 9 deletions(-) diff --git a/GaussTuranExampleTemp/example.jl b/GaussTuranExampleTemp/example.jl index c57dc44..2098d52 100644 --- a/GaussTuranExampleTemp/example.jl +++ b/GaussTuranExampleTemp/example.jl @@ -1,13 +1,13 @@ using Optim +using DoubleFloats using TaylorDiff using PreallocationTools using Integrals -using DoubleFloats # This fix is needed to avoid crashing: https://github.com/JuliaLang/julia/pull/54201 include("linear_algebra_fix.jl") -FT = Float64 # Double64 +FT = Double64 n = 5 s = 1 a = FT(0.0) @@ -44,7 +44,7 @@ I = Integrals.GaussTuran( max_int_error_upper = maximum(abs(I(x -> f(x, j)) - I.cache.rhs_upper[j]) for j in 1:(I.cache.N)) @show max_int_error_upper -# max_int_error_upper = 2.220446049250313e-16 +# max_int_error_upper = 2.465190328815662e-32 # Integration error |I(fⱼ) - ₐ∫ᵇfⱼ(x)dx| last n functions fⱼ max_int_error_lower = maximum( @@ -53,9 +53,9 @@ for j in (I.cache.N + 1):(I.cache.N + I.cache.n) ) @show max_int_error_lower -# max_int_error_lower = 2.0211055051788662e-12 +# max_int_error_lower = 9.079315240327527e-30 # Example with eˣ exp_error = abs(I(Base.exp) - (Base.exp(1) - 1)) @show exp_error; -# exp_error = 1.5543122344752192e-15 +# exp_error = 6.43662141695295e-18 diff --git a/Project.toml b/Project.toml index 25a276b..03e0ca0 100644 --- a/Project.toml +++ b/Project.toml @@ -32,8 +32,8 @@ IntegralsCubaExt = "Cuba" IntegralsCubatureExt = "Cubature" IntegralsFastGaussQuadratureExt = "FastGaussQuadrature" IntegralsForwardDiffExt = "ForwardDiff" -IntegralsMCIntegrationExt = "MCIntegration" IntegralsGaussTuranExt = ["Optim", "TaylorDiff", "PreallocationTools"] +IntegralsMCIntegrationExt = "MCIntegration" IntegralsZygoteExt = ["Zygote", "ChainRulesCore"] [compat] diff --git a/ext/IntegralsGaussTuranExt.jl b/ext/IntegralsGaussTuranExt.jl index b06e693..8b54e84 100644 --- a/ext/IntegralsGaussTuranExt.jl +++ b/ext/IntegralsGaussTuranExt.jl @@ -87,14 +87,14 @@ function GaussTuranLoss!(f, ΔX::AbstractVector{T}, cache) where {T} cumsum!(X, ΔX) X .+= a - # Evaluating f! + # Evaluating f for (i, x) in enumerate(X) # Threads.@threads for j = 1:N for j in 1:N - M_upper[j, i:n:N] .= derivatives(x -> f(x, j), x, one(T), Val(2s + 1)).value + M_upper[j, i:n:N] .= derivatives(x -> f(x, j), x, 1, Val(2s + 1)).value end Threads.@threads for j in (N + 1):(N + n) - M_lower[j - N, i:n:N] .= derivatives(x -> f(x, j), x, one(T), Val(2s + 1)).value + M_lower[j - N, i:n:N] .= derivatives(x -> f(x, j), x, 1, Val(2s + 1)).value end end