From b60c6a467ccb8fce4235c5a4422f89651e73c1d2 Mon Sep 17 00:00:00 2001 From: Bart de Koning Date: Mon, 22 Jul 2024 22:36:49 +0200 Subject: [PATCH] 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()