From de3c0ae39ccbd70f6ddf26b1e8d3948c17ec1b12 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tam=C3=A1s=20K=2E=20Papp?= Date: Tue, 21 May 2024 17:11:26 +0200 Subject: [PATCH 01/14] WIP: rewrite univariate and partials spec --- Project.toml | 1 + src/SpectralKit.jl | 13 +- src/derivatives.jl | 772 +++++++++++++++++++++++++-------------- src/utilities.jl | 15 + test/Project.toml | 1 + test/runtests.jl | 12 +- test/test_derivatives.jl | 261 ++++++++----- test/utilities.jl | 10 +- 8 files changed, 700 insertions(+), 385 deletions(-) diff --git a/Project.toml b/Project.toml index dc3fde2..4840493 100644 --- a/Project.toml +++ b/Project.toml @@ -6,6 +6,7 @@ version = "0.15.0" [deps] ArgCheck = "dce04be8-c92d-5529-be00-80e4d2c0e197" DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" +OrderedCollections = "bac558e1-5e72-5ebc-8fee-abe8a469f55d" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" [compat] diff --git a/src/SpectralKit.jl b/src/SpectralKit.jl index 0ac0159..f2d6000 100644 --- a/src/SpectralKit.jl +++ b/src/SpectralKit.jl @@ -2,15 +2,16 @@ module SpectralKit using ArgCheck: @argcheck using DocStringExtensions: FUNCTIONNAME, SIGNATURES, TYPEDEF +using OrderedCollections: OrderedSet using StaticArrays: MVector, SVector, sacollect include("utilities.jl") include("derivatives.jl") -include("domains.jl") -include("transformations.jl") -include("generic_api.jl") -include("chebyshev.jl") -include("smolyak_traversal.jl") -include("smolyak_api.jl") +# include("domains.jl") +# include("transformations.jl") +# include("generic_api.jl") +# include("chebyshev.jl") +# include("smolyak_traversal.jl") +# include("smolyak_api.jl") end # module diff --git a/src/derivatives.jl b/src/derivatives.jl index 3c07113..a790902 100644 --- a/src/derivatives.jl +++ b/src/derivatives.jl @@ -2,30 +2,46 @@ ##### Internal implementation for derivatives. See docs of [`Derivatives`](@ref). ##### -export derivatives, ∂ +export 𝑑, ∂ #### -#### derivatives API +#### Operations we support, with fallbacks defined below. This is deliberately kept +#### minimal, now all versions are defined for commutative operators. #### -""" -A small AD framework used *internally*, for calculating derivatives. +_one(::Type{T}) where {T<:Real} = one(T) -Supports only the operations required by this module. The restriction is deliberate, -should not be used for arithmetic operators outside this package. +_add(x::Real, y::Real) = x + y -See [`derivatives`](@ref) for the exposed API. -""" -struct Derivatives{N,T} - "The function value and derivatives." - derivatives::NTuple{N,T} - function Derivatives(derivatives::NTuple{N,T}) where {N,T} - new{N,T}(derivatives) +_sub(x::Real, y::Real) = x - y + +_mul(x::Real, y::Real) = x * y + +_mul(x, y, z) = _mul(_mul(x, y), z) + +#### +#### Univariate expansions and derivatives +#### + + +# """ +# A small AD framework used *internally*, for calculating derivatives. + +# Supports only the operations required by this module. The restriction is deliberate, +# should not be used for arithmetic operators outside this package. + +# See [`∂`](@ref) for the exposed API. +# """ +struct 𝑑Expansion{Dp1,T} + "The function value and derivatives. `Dp1` is the degree of the last derivative + 1." + coefficients::SVector{Dp1,T} + function 𝑑Expansion(coefficients::SVector{Dp1,T}) where {Dp1,T} + new{Dp1,T}(coefficients) end end -function Base.show(io::IO, x::Derivatives) - for (i, d) in enumerate(x.derivatives) +function Base.show(io::IO, x::𝑑Expansion) + for (i, d) in enumerate(x.coefficients) i ≠ 1 && print(io, " + ") print(io, d) i ≥ 2 && print(io, "⋅Δ") @@ -33,338 +49,530 @@ function Base.show(io::IO, x::Derivatives) end end -Base.eltype(::Type{Derivatives{N,T}}) where {N,T} = T +Base.eltype(::Type{𝑑Expansion{Dp1,T}}) where {Dp1,T} = T -@inline Base.getindex(x::Derivatives, i::Int) = x.derivatives[i + 1] - -"Types accepted as scalars in this package." -const Scalar = Union{Real,Derivatives} +@inline Base.getindex(x::𝑑Expansion, i::Int) = x.coefficients[i + 1] """ - derivatives(x, ::Val(N) = Val(1)) +$(SIGNATURES) -Obtain `N` derivatives (and the function value) at a scalar `x`. The `i`th derivative -can be accessed with `[i]` from results, with `[0]` for the function value. +Requests the calculation of `D ≥ 0` derivatives within this package. -# Important note about transformations +The preferred syntax is `𝑑^D`, which is type-stable. +""" +struct 𝑑Derivatives{D} + function 𝑑Derivatives{D}() where D + @argcheck D isa Int + @argcheck D ≥ 0 + new{D}() + end +end -Always use `derivatives` *before* a transformation for correct results. For example, for -some transformation `t` and value `x` in the transformed domain, +Base.show(io::IO, ::𝑑Derivatives{D}) where D = print(io, "𝑑^$(D)") -```julia -# right -linear_combination(basis, θ, transform_to(domain(basis), t, derivatives(x))) -# right (convenience form) -(linear_combination(basis, θ) ∘ t)(derivatives(x)) -``` +""" +A callable that calculates the value and the derivative of the argument. Higher-order +derivatives can be obtained by using an exponent, or multiplication. -instead of +```jldoctest +julia> 𝑑(2.0) +2.0 + 1.0⋅Δ -```julia -# WRONG -linear_combination(basis, θ, derivatives(transform_to(domain(basis), t, x))) +julia> (𝑑^3)(2.0) +2.0 + 1.0⋅Δ + 0.0⋅Δ² + 0.0⋅Δ³ ``` -For multivariate calculations, use the [`∂`](@ref) interface. +Note that non-literal exponentiation requires `^Val(y)`, for type stability. +""" +const 𝑑 = 𝑑Derivatives{1}() -# Example +Base.:*(::𝑑Derivatives{D1}, ::𝑑Derivatives{D2}) where {D1,D2} = 𝑑Derivatives{D1+D2}() -```jldoctest -julia> basis = Chebyshev(InteriorGrid(), 3) -Chebyshev polynomials (1st kind), InteriorGrid(), dimension: 3 +function Base.literal_pow(::typeof(^), ::𝑑Derivatives{D}, ::Val{Y}) where {D,Y} + @argcheck Y isa Int && Y ≥ 0 + 𝑑Derivatives{D*Y}() +end -julia> C = collect(basis_at(basis, derivatives(0.1))) -3-element Vector{SpectralKit.Derivatives{2, Float64}}: - 1.0 + 0.0⋅Δ - 0.1 + 1.0⋅Δ - -0.98 + 0.4⋅Δ +# ^::Int deliberately unsupported +Base.:^(d::𝑑Derivatives, y::Val{Y}) where Y = Base.literal_pow(^, d, y) -julia> C[1][1] # 1st derivative of the linear term is 1 -0.0 -``` -""" -function derivatives(x::T, ::Val{N} = Val(1)) where {N, T <: Real} - Derivatives((x, ntuple(i -> i == 1 ? one(T) : zero(T), Val(N))...)) +function (::𝑑Derivatives{D})(x::T) where {D, T <: Real} + 𝑑Expansion(SVector(x, ntuple(i -> i == 1 ? one(T) : zero(T), Val(D))...)) end -#### -#### partial derivatives API -#### +function _one(::Type{𝑑Expansion{Dp1,T}}) where {Dp1,T} + 𝑑Expansion(SVector(ntuple(i -> i == 1 ? _one(T) : zero(T), Val(Dp1)))) +end -""" -A specification for partial derivatives. Each element of `lookup` is an N-dimensional -tuple, the elements of which determine the order of the derivative along that -coordinate, eg `(1, 2, 0)` means 1st derivative along coordinate 1, 2nd derivative along -coordinate 2. +function _add(x::𝑑Expansion{Dp1}, y::𝑑Expansion{Dp1}) where Dp1 + 𝑑Expansion(map(_add, x.coefficients, y.coefficients)) +end -The elementwise maximum is stored in `M`. It is a type parameter because it needs to be -available as such for `derivatives`. The implementation than calculates derivatives -along coordinates, and combines them according to `lookups`. +function _sub(x::𝑑Expansion, y::𝑑Expansion) + 𝑑Expansion(map(_sub, x.coefficients, y.coefficients)) +end -Internal. -""" -struct ∂Specification{M,L} - lookups::L - function ∂Specification{M}(lookups::L) where {M,L} - @argcheck M isa Tuple{Vararg{Int}} - @argcheck !isempty(lookups) - @argcheck lookups isa Tuple{Vararg{typeof(M)}} - new{M,L}(lookups) - end +function _mul(x::Real, y::𝑑Expansion) + 𝑑Expansion(map(y -> _mul(x, y), y.coefficients)) end -function Base.show(io::IO, ∂specification::∂Specification) - (; lookups) = ∂specification - print(io, "partial derivatives") - for (j, lookup) in enumerate(lookups) - print(io, "\n[$j] ") - if all(iszero, lookup) - print(io, "f") - else - s = sum(lookup) - print(io, "∂") - s ≠ 1 && print(io, SuperScript(s)) - print(io, "f/") - for (i, l) in enumerate(lookup) - l == 0 && continue # don't print ∂⁰ - print(io, "∂") - if l ≠ 1 - print(io, SuperScript(l)) - end - print(io, "x", SubScript(i)) - end - end +@generated function _mul(x::𝑑Expansion{Dp1}, y::𝑑Expansion{Dp1}) where {Dp1} + _sum_terms(k) = mapreduce(i -> :(_mul($(binomial(k, i)), xd[$(i + 1)], yd[$(k - i + 1)])), + (a, b) -> :(_add($(a), $(b))), 0:k) + _derivatives(k) = mapfoldl(_sum_terms, (a, b) -> :($(a)..., $(b)), 0:(Dp1-1); init = ()) + quote + xd = x.coefficients + yd = y.coefficients + 𝑑Expansion(SVector($(_derivatives(Dp1)))) end end -""" -$(SIGNATURES) +"Types accepted as scalars in this package." +const Scalar = Union{Real,𝑑Expansion} -Convert a partial specification to a lookup, for `N`-dimensional arguments. Eg +# """ +# derivatives(x, ::Val(N) = Val(1)) -```jldoctest -julia> SpectralKit._partial_to_lookup(Val(3), (1, 3, 1)) -(2, 0, 1) -``` -""" -function _partial_to_lookup(::Val{N}, partial::Tuple{Vararg{Int}}) where N - @argcheck all(p -> 0 < p ≤ N, partial) - ntuple(d -> sum(p -> p == d, partial; init = 0), Val(N)) -end +# Obtain `N` derivatives (and the function value) at a scalar `x`. The `i`th derivative +# can be accessed with `[i]` from results, with `[0]` for the function value. -""" -$(SIGNATURES) +# # Important note about transformations -Partial derivative specification. The first argument is `Val(::Int)` or simply an `Int` -(for convenience, using constant folding), determining the dimension of the argument. +# Always use `derivatives` *before* a transformation for correct results. For example, for +# some transformation `t` and value `x` in the transformed domain, -Subsequent arguments are indices of the input variable. +# ```julia +# # right +# linear_combination(basis, θ, transform_to(domain(basis), t, derivatives(x))) +# # right (convenience form) +# (linear_combination(basis, θ) ∘ t)(derivatives(x)) +# ``` -```jldoctest -julia> ∂(3, (), (1, 1), (2, 3)) -partial derivatives -[1] f -[2] ∂²f/∂²x₁ -[3] ∂²f/∂x₂∂x₃ -``` -""" -@inline function ∂(::Val{N}, partials...) where N - @argcheck N ≥ 1 "Needs at least one dimension." - @argcheck !isempty(partials) "Empty partial derivative specification." - lookups = map(p -> _partial_to_lookup(Val(N), p), partials) - M = ntuple(d -> maximum(l -> l[d], lookups), Val(N)) - ∂Specification{M}(lookups) -end +# instead of -@inline ∂(N::Integer, partials...) = ∂(Val(Int(N)), partials...) +# ```julia +# # WRONG +# linear_combination(basis, θ, derivatives(transform_to(domain(basis), t, x))) +# ``` -""" -Partial derivatives to be evaluated at some `x`. These need to be [`_lift`](@ref)ed, -then combined with [`_product`](@ref) from bases. Internal, use `∂(specification, x)` to -construct. -""" -struct ∂Input{TS<:∂Specification,TX<:SVector} - ∂specification::TS - x::TX - function ∂Input(∂specification::TS, x::TX) where {M,N,TS<:∂Specification{M},TX<:SVector{N}} - @argcheck length(M) == N - new{TS,TX}(∂specification, x) - end -end +# For multivariate calculations, use the [`∂`](@ref) interface. -function Base.show(io::IO, ∂x::∂Input) - show(io, ∂x.∂specification) - print(io, "\nat ", ∂x.x) -end +# # Example -""" -$(SIGNATURES) +# ```jldoctest +# julia> basis = Chebyshev(InteriorGrid(), 3) +# Chebyshev polynomials (1st kind), InteriorGrid(), dimension: 3 -Input wrappert type for evaluating partial derivatives `∂specification` at `x`. +# julia> C = collect(basis_at(basis, derivatives(0.1))) +# 3-element Vector{SpectralKit.Derivatives{2, Float64}}: +# 1.0 + 0.0⋅Δ +# 0.1 + 1.0⋅Δ +# -0.98 + 0.4⋅Δ -```jldoctest -julia> using StaticArrays - -julia> s = ∂(Val(2), (), (1,), (2,), (1, 2)) -partial derivatives -[1] f -[2] ∂f/∂x₁ -[3] ∂f/∂x₂ -[4] ∂²f/∂x₁∂x₂ - -julia> ∂(s, SVector(1, 2)) -partial derivatives -[1] f -[2] ∂f/∂x₁ -[3] ∂f/∂x₂ -[4] ∂²f/∂x₁∂x₂ -at [1, 2] -``` -""" -function ∂(∂specification::∂Specification{M}, x::Union{AbstractVector,Tuple}) where M - N = length(M) - ∂Input(∂specification, SVector{N}(x)) -end +# julia> C[1][1] # 1st derivative of the linear term is 1 +# 0.0 +# ``` +# """ -""" -$(SIGNATURES) +#### +#### partial derivatives API +#### -Shorthand for `∂(x, ∂(Val(length(x)), partials...))`. Ideally needs an `SVector` or a -`Tuple` so that size information can be obtained statically. -""" -@inline function ∂(x::SVector{N}, partials...) where N - ∂specification = ∂(Val(N), partials...) - ∂Input(∂specification, x) +struct Partials{N} + I::NTuple{N,Int} + function Partials(I::NTuple{N,Int}) where N + @argcheck all(i -> i ≥ 0, I) + @argcheck I ≡ () || last(I) ≠ 0 + new{N}(I) + end end -@inline ∂(x::Tuple, partials...) = ∂(SVector(x), partials...) - -""" -See [`_lift`](@ref). Internal. -""" -struct ∂InputLifted{TS<:∂Specification,TL<:Tuple} - ∂specification::TS - lifted_x::TL +function Partials(I::Integer...) + N = length(I) + while N > 0 && I[N] == 0 + @show N + N -= 1 + end + Partials(ntuple(i -> Int(I[i]), N)) end -""" -$(SIGNATURES) - -Lift a partial derivative calculation into a tuple of `Derivatives`. Internal. -""" -@generated function _lift(∂x::∂Input{<:∂Specification{M}}) where M - _lifted_x = [:(derivatives(x[$(i)], Val($(m)))) for (i, m) in enumerate(M)] - quote - x = ∂x.x - lifted_x = ($(_lifted_x...),) - ∂InputLifted(∂x.∂specification, lifted_x) +function _is_strict_subset(p1::Partials{N1}, p2::Partials{N2}) where {N1,N2} + N1 > N2 && return false # valid because derivatives are always positive + I1 = p1.I + I2 = p2.I + strict = false + for i in 1:min(N1, N2) + I1[i] > I2[i] && return false + strict |= I1[i] < I2[i] end + strict || N1 < N2 end -#### -#### products (used by tensor / Smolyak bases) -#### - -""" -$(SIGNATURES) - -Conceptually equivalent to `prod(getindex.(sources, indices))`, which it returns when -`kind` is `nothing`, a placeholder calculating any derivatives. Internal. -""" -_product(kind::Nothing, sources, indices) = mapreduce(getindex, *, sources, indices) +function Base.isless(p1::Partials{N1}, p2::Partials{N2}) where {N1, N2} + _is_strict_subset(p1, p2) && return true + _is_strict_subset(p2, p1) && return false + p1.I < p2.I +end -""" -$(SIGNATURES) +Base.isequal(p1::Partials, p2::Partials) = p1.I == p2.I -Type that is returnedby [`_product`](@ref). -""" -function _product_type(::Type{Nothing}, source_eltypes) - mapfoldl(eltype, promote_type, source_eltypes) +function _partials_minimal_representation(partials) + descending_partials = sort!(collect(Partials, partials); rev = true) + minimal_partials = Vector{Partials}() + for p in descending_partials + if isempty(minimal_partials) || !_is_strict_subset(p, minimal_partials[end]) + push!(minimal_partials, p) + end + end + minimal_partials end -""" -Container for output of evaluating partial derivatives. Each corresponds to an -specification in a [`∂Specification`](@ref). Can be indexed with integers, iterated, or -converted to a `Tuple`. -""" -struct ∂Output{N,T} - values::NTuple{N,T} +function _partials_canonical_expansion(::Val{N}, partials) where N + result = OrderedSet{NTuple{N,Int}}() + function _plus1_cartesian_indices(p::Partials{M}) where M + (; I) = p + @argcheck M ≤ N + CartesianIndices(ntuple(i -> i ≤ M ? I[i] + 1 : 1, Val(N))) + end + for p in partials + for ι in _plus1_cartesian_indices(p) + i = map(i -> i - 1, Tuple(ι)) + if !(i ∈ result) + push!(result, i) + end + end + end + result end -function Base.show(io::IO, ∂output::∂Output) - print(io, "SpectralKit.∂Output(") - join(io, ∂output.values, ", ") - print(io, ")") -end -@inline Base.Tuple(∂output::∂Output) = ∂output.values -@inline Base.length(∂output::∂Output) = length(∂output.values) +# const _PARTIALS = Tuple{Vararg{Int}} + +# function _partials_max(D1::NTuple{N1,Int}, D2::NTuple{N2,Int}) where {N1,N2} +# ntuple(Val(max(N1,N2))) do i +# if i ≤ N1 +# if i ≤ N2 +# max(D1[i], D2[i]) +# else +# D1[i] +# end +# else +# D2[i] +# end +# end +# end + +# function _partials_isless_or_eq(D1::NTuple{N1,Int}, D2::NTuple{N2,Int}) where {N1,N2} +# N1 ≤ N2 || return false +# for i in 1:N1 +# D1[i] ≤ D2[i] || return false +# end +# true +# end + +# function _normalize_partials(D::Tuple) +# i = something(findlast(d -> d > 0, D), 0) +# D[1:i] +# end + +# function _expand_M_Ds_union(k_d_pairs) +# UDs = Vector{_PARTIALS}() +# M = () +# for D in Ds +# M = _partials_max(M, D) +# for i in CartesianIndices(map(d -> 0:d, D)) +# _sorted_insert_new!(UDs, _normalize_partials(Tuple(i))) +# end +# end +# M, UDs +# end + +# @generated function _check_M_Ds_consistency(::Val{N}, ::Val{Ds}) where {N,Ds} +# @argcheck Ds isa Tuple{Vararg{_PARTIALS}} && all(D -> all(d -> d ≥ 0, D), Ds) +# @argcheck M isa NTuple{K,Int} && all(m -> m ≥ 0, M) +# M′, Ds′ = _calculate_M_Ds_union(Ds) +# C = M′ ≡ M && Ds′ ≡ Ds +# quote +# Val($(C)) +# end +# end + +# struct ∂Derivatives{K,M,Ds} +# function ∂Derivatives{K,M,Ds}() where {K,M,Ds} +# @argcheck K isa Int && K ≥ 0 +# @argcheck _check_M_Ds_consistency(Val(N), Val(Ds)) ≡ Val(true) +# new{K,M,Ds}() +# end +# function ∂Derivatives(dimension::Val{K}, degree::Val{m}) where {K,m} +# @argcheck K isa Int && K ≥ 1 +# @argcheck m isa Int && m ≥ 0 +# m == 0 && return new{0,(),((),)}() +# z = ntuple(_ -> 0, Val(K - 1)) +# M = (z..., m) +# Ds = ntuple(i -> (z..., i - 1), Val(m + 1)) +# new{K,M,Ds}() +# end +# end + +# ∂(k::Val{K}, m::Val{M}) where {K,M} = ∂Derivatives(k, m) + +# @inline ∂(dimension::Int, degree::Int = 1) = ∂(Val(dimension), Val(degree)) + +# @generated function Base.:∪(∂ds::∂Derivatives...) +# _get_Ds(::∂Derivatives{K,N,Ds}) where {K,N,Ds} = Ds +# M, Ds = _calculate_M_Ds_union(Iterators.flatten((_get_Ds(∂d) for ∂d in ∂ds))) +# K = length(M) +# quote +# ∂Derivatives{$(K), $(M), $(Ds)}() +# end +# end + +# function Base.:^(::∂Derivatives{K,M,D}, ::Val{Y}) where {K,M,D,Y} +# @argcheck Y isa Int && Y ≥ 0 +# (dhead..., dlast) = D +# @argcheck(all(iszero, dhead...), +# "The ^ operator can only be applied for ∂ derivatives along one dimension.") +# ∂Derivatives(Val(K), Val(M * Y)) +# end + +# function _collapse_D_union!(D::Vector{_PARTIALS}) +# D = collect(D) +# U = Vector{Pair{Int,Int}}() +# while !isempty(D) +# M = reduce(_partials_max, D) +# k = findfirst(d -> d > 0, M) +# k ≡ nothing && break # we are done +# m = M[k] +# push!(U, k => m) +# _d = (ntuple(_ -> 0, k - 1)..., m) +# filter!(d -> !_partials_isless_or_eq(d, _d), D) +# end +# U +# end + +# function Base.show(io::IO, ::∂Derivatives{K,M,D}) where {K,M,D} +# print(io, "partial derivatives (at least $(K) dimensions), up to ") +# join(io, ("∂($(k), $(m))" for (k, m) in _collapse_D_union!(collect(_PARTIALS, D))), " ∪ ") +# end + +# struct ∂DerivativesInput{P<:∂Derivatives,N,T} +# partial_derivatives::P +# x::SVector{N,T} +# function ∂DerivativesInput(∂derivatives::P, +# x::SVector{N,T}) where {K,P<:∂Derivatives{K},N,T} +# @argcheck N ≤ K "Can't differentiate the $(K)th element of a $(N)-element input." +# new{P,N,T}(∂derivatives, x) +# end +# end + +# (∂derivatives::∂Derivatives)(x::SVector) = ∂DerivativesInput(partial_derivatives, x) + +# """ +# See [`_lift`](@ref). Internal. +# """ +# struct ∂DerivativesInputLifted{P<:∂Derivatives,L<:Tuple} +# partial_derivatives::P +# lifted_x::L +# end + +# """ +# $(SIGNATURES) + +# Lift a partial derivative calculation into a tuple of `Derivatives`. Internal. +# """ +# @generated function _lift(∂x::∂DerivativesInput{<:∂Derivatives{K,M,N}}) where {K,M,N} +# _lifted_x = [:(∂(Val($(i ≤ K ? M[i] : 0)))(x[$(i)])) for i in 1:N] +# quote +# x = ∂x.x +# lifted_x = ($(_lifted_x...),) +# ∂PartialDerivativesInputLifted(∂x.∂derivatives, lifted_x) +# end +# end + +# struct ∂Expansion{P<:∂Derivatives,N,T} +# ∂derivatives::P +# coefficients::SVector{N,T} +# function ∂Expansion(∂derivatives::∂Derivatives{K,M,D}, +# coefficients::SVector{N,T}) where {K,M,D,N,T} +# @argcheck length(D) == N +# new{typeof(P),N,T}(∂derivatives, coefficients) +# end +# end + +# function Base.show(io::IO, expansion::∂Expansion{<:∂Derivatives{K,M,D}}) where {K,M,D} +# (; coefficients) = expansion +# print(io, "multivariate expansion:") +# for (c, d) in enumerate(zip(coefficiends, D)) +# print(io, "\n") +# _print_partial_notation(io, d) +# print(io, " ", c) +# end +# end + +# """ +# $(SIGNATURES) + +# Conceptually equivalent to `prod(getindex.(sources, indices))`, which it returns when +# `kind` is `nothing`, a placeholder calculating any derivatives. Internal. +# """ +# _product(kind::Nothing, sources, indices) = mapreduce(getindex, *, sources, indices) + +# """ +# $(SIGNATURES) + +# Type that is returnedby [`_product`](@ref). +# """ +# function _product_type(::Type{Nothing}, source_eltypes) +# mapfoldl(eltype, promote_type, source_eltypes) +# end + +# @generated function _product(∂derivatives::∂Derivatives{K,M,D}, +# sources::NTuple{N}, +# indices) where {K,M,D,N} +# xs = [gensym(:x) for _ in 1:N] +# assignments = [:($(xs[i]) = sources[indices[$(i)]]) for i in 1:N] +# products = [begin +# ι = NTuple(i -> i ≤ K ? d[i] + 1 : 1, Val(N)) +# mapreduce(i -> :($(xs[i])[$(ι[i])]), +# (a, b) -> :($(a) * $(b)), +# 1:N) +# end for d in D] +# quote +# $(assignments...) +# ∂Expansion(∂derivatives, SVector($(products)...)) +# end +# end + +# function _product_type(::Type{∂Derivatives{M,L}}, source_eltypes) where {M,L} +# T = _product_type(Nothing, map(eltype, source_eltypes)) +# N = length(fieldtypes(L)) +# ∂Output{N,T} +# end + +# function _add(x::∂Expansion{P}, y::∂Expansion{P}) where P +# ∂Expansion(x.∂derivatives, map(+, x.values, y.values)) +# end + +# _mul(x::Real, y::∂Expansion) = ∂Expansion(y.∂derivatives, map(y -> _mul(x, y), y.values)) -@inline Base.getindex(∂output::∂Output, i) = ∂output.values[i] +##### +##### FIXME revise and move documentation below +##### -@inline Base.iterate(∂output::∂Output, i...) = Base.iterate(∂output.values, i...) -function _product(∂specification::∂Specification, sources, indices) - (; lookups) = ∂specification - p = map(lookups) do lookup - mapreduce((l, s, i) -> s[i][l], *, lookup, sources, indices) - end - ∂Output(p) -end +# """ +# $(SIGNATURES) + +# Partial derivative specification. The first argument is `Val(::Int)` or simply an `Int` +# (for convenience, using constant folding), determining the dimension of the argument. + +# Subsequent arguments are indices of the input variable. + +# ```jldoctest +# julia> ∂(3, (), (1, 1), (2, 3)) +# partial derivatives +# [1] f +# [2] ∂²f/∂²x₁ +# [3] ∂²f/∂x₂∂x₃ +# ``` +# """ +# @inline function ∂(::Val{N}, partials...) where N +# @argcheck N ≥ 1 "Needs at least one dimension." +# @argcheck !isempty(partials) "Empty partial derivative specification." +# lookups = map(p -> _partial_to_lookup(Val(N), p), partials) +# M = ntuple(d -> maximum(l -> l[d], lookups), Val(N)) +# ∂Specification{M}(lookups) +# end + +# @inline ∂(N::Integer, partials...) = ∂(Val(Int(N)), partials...) + +# """ +# Partial derivatives to be evaluated at some `x`. These need to be [`_lift`](@ref)ed, +# then combined with [`_product`](@ref) from bases. Internal, use `∂(specification, x)` to +# construct. +# """ +# struct ∂Input{TS<:∂Specification,TX<:SVector} +# ∂specification::TS +# x::TX +# function ∂Input(∂specification::TS, x::TX) where {M,N,TS<:∂Specification{M},TX<:SVector{N}} +# @argcheck length(M) == N +# new{TS,TX}(∂specification, x) +# end +# end + +# function Base.show(io::IO, ∂x::∂Input) +# show(io, ∂x.∂specification) +# print(io, "\nat ", ∂x.x) +# end + +# """ +# $(SIGNATURES) + +# Input wrappert type for evaluating partial derivatives `∂specification` at `x`. + +# ```jldoctest +# julia> using StaticArrays + +# julia> s = ∂(Val(2), (), (1,), (2,), (1, 2)) +# partial derivatives +# [1] f +# [2] ∂f/∂x₁ +# [3] ∂f/∂x₂ +# [4] ∂²f/∂x₁∂x₂ + +# julia> ∂(s, SVector(1, 2)) +# partial derivatives +# [1] f +# [2] ∂f/∂x₁ +# [3] ∂f/∂x₂ +# [4] ∂²f/∂x₁∂x₂ +# at [1, 2] +# ``` +# """ +# function ∂(∂specification::∂Specification{M}, x::Union{AbstractVector,Tuple}) where M +# N = length(M) +# ∂Input(∂specification, SVector{N}(x)) +# end + +# """ +# $(SIGNATURES) + +# Shorthand for `∂(x, ∂(Val(length(x)), partials...))`. Ideally needs an `SVector` or a +# `Tuple` so that size information can be obtained statically. +# """ +# @inline function ∂(x::SVector{N}, partials...) where N +# ∂specification = ∂(Val(N), partials...) +# ∂Input(∂specification, x) +# end + +# @inline ∂(x::Tuple, partials...) = ∂(SVector(x), partials...) -function _product_type(::Type{∂Specification{M,L}}, source_eltypes) where {M,L} - T = _product_type(Nothing, map(eltype, source_eltypes)) - N = length(fieldtypes(L)) - ∂Output{N,T} -end #### -#### operations we support -#### -#### This is deliberately kept minimal, now all versions are defined for commutative -#### operators. +#### products (used by tensor / Smolyak bases) #### -_one(::Type{T}) where {T<:Real} = one(T) -function _one(::Type{Derivatives{N,T}}) where {N,T} - Derivatives(ntuple(i -> i == 1 ? _one(T) : zero(T), Val(N))) -end +# """ +# Container for output of evaluating partial derivatives. Each corresponds to an +# specification in a [`∂Specification`](@ref). Can be indexed with integers, iterated, or +# converted to a `Tuple`. +# """ +# struct ∂Output{N,T} +# values::NTuple{N,T} +# end -_add(x::Real, y::Real) = x + y +# function Base.show(io::IO, ∂output::∂Output) +# print(io, "SpectralKit.∂Output(") +# join(io, ∂output.values, ", ") +# print(io, ")") +# end -function _add(x::Derivatives, y::Derivatives) - Derivatives(map(_add, x.derivatives, y.derivatives)) -end +# @inline Base.Tuple(∂output::∂Output) = ∂output.values -function _add(x::∂Output, y::∂Output) - ∂Output(map(+, x.values, y.values)) -end - -_sub(x::Real, y::Real) = x - y - -function _sub(x::Derivatives, y::Derivatives) - Derivatives(map(_sub, x.derivatives, y.derivatives)) -end - -_mul(x::Real, y::Real) = x * y - -_mul(x, y, z) = _mul(_mul(x, y), z) - -function _mul(x::Real, y::Derivatives) - Derivatives(map(y -> _mul(x, y), y.derivatives)) -end +# @inline Base.length(∂output::∂Output) = length(∂output.values) -_mul(x::Real, y::∂Output) = ∂Output(map(y -> _mul(x, y), y.values)) +# @inline Base.getindex(∂output::∂Output, i) = ∂output.values[i] -@generated function _mul(x::Derivatives{N}, y::Derivatives{N}) where {N} - _sum_terms(k) = mapreduce(i -> :(_mul($(binomial(k, i)), xd[$(i + 1)], yd[$(k - i + 1)])), - (a, b) -> :(_add($(a), $(b))), 0:k) - _derivatives(k) = mapfoldl(_sum_terms, (a, b) -> :($(a)..., $(b)), 0:(N-1); init = ()) - quote - xd = x.derivatives - yd = y.derivatives - Derivatives($(_derivatives(N))) - end -end +# @inline Base.iterate(∂output::∂Output, i...) = Base.iterate(∂output.values, i...) diff --git a/src/utilities.jl b/src/utilities.jl index 84e2f49..ef76c98 100644 --- a/src/utilities.jl +++ b/src/utilities.jl @@ -39,6 +39,21 @@ Base.print(io::IO, s::SubScript) = print_number(io, _SUBSCRIPT_DIGITS, s.i) """ $(SIGNATURES) +Print notation for partial derivatives, where `d[i]` stands for ``∂ⁱ/∂x[d]ⁱ``. +""" +function _print_partial_notation(io::IO, d) + if d ≡ () + print(io, "value") + else + for (i, p) in enumerate(d) + print(io, "∂", SubScript(i), SuperScript(p)) + end + end +end + +""" +$(SIGNATURES) + If `T <: NTuple{N}`, convert `v` into an `NTuple{N}`. Used for ingesting `::AbstractVector` arguments in contexts where an `NTuple` or diff --git a/test/Project.toml b/test/Project.toml index 9e95a9e..2bd97b0 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -2,6 +2,7 @@ BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf" DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" FiniteDifferences = "26cc04aa-876d-5657-8c51-4c34ba976000" +Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" Sobol = "ed01d8cd-4d21-5b2a-85b4-cc3bdc58bad4" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" diff --git a/test/runtests.jl b/test/runtests.jl index b6c78fe..f1791ac 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,13 +1,13 @@ using SpectralKit -using SpectralKit: PM1 +# using SpectralKit: PM1 using Test, DocStringExtensions, StaticArrays, BenchmarkTools, FiniteDifferences include("utilities.jl") include("test_utilities.jl") -include("test_domains.jl") include("test_derivatives.jl") -include("test_transformations.jl") -include("test_generic_api.jl") -include("test_chebyshev.jl") -include("test_smolyak.jl") +# include("test_domains.jl") +# include("test_transformations.jl") +# include("test_generic_api.jl") +# include("test_chebyshev.jl") +# include("test_smolyak.jl") diff --git a/test/test_derivatives.jl b/test/test_derivatives.jl index 8681a37..84852aa 100644 --- a/test/test_derivatives.jl +++ b/test/test_derivatives.jl @@ -1,103 +1,190 @@ -using SpectralKit: ∂Specification, ∂Output -using SpectralKit: _add # used to form a scalar for testing -using SpectralKit: Derivatives # test internals +using SpectralKit, Test +# test internals +using SpectralKit: 𝑑Derivatives, _one, _add, _sub, _mul, + Partials, _is_strict_subset, _partials_minimal_representation, _partials_canonical_expansion + # _partials_max, +# _normalize_partials, _calculate_M_Ds_union, _collapse_D_union!, ∂Derivatives +using Random: randperm +using StaticArrays: SVector -@testset "partial derivatives interface" begin - @test ∂(2, (1, 2), (2, 2)) == ∂Specification{(1, 2)}(((1, 1), (0, 2))) - @test ∂(3, (), (1, 1, 2), (3, 2)) == ∂Specification{(2, 1, 1)}(((0, 0, 0), (2, 1, 0), (0, 1, 1))) - @test_throws ArgumentError ∂(3, (4, )) - @test_throws ArgumentError ∂(3, (-1, )) - ∂spec = ∂(2, (1, 2), (2, 2)) - @test repr(∂spec) == "partial derivatives\n[1] ∂²f/∂x₁∂x₂\n[2] ∂²f/∂²x₂" - @test repr(∂(∂spec, [1.0, 2.0])) == - "partial derivatives\n[1] ∂²f/∂x₁∂x₂\n[2] ∂²f/∂²x₂\nat [1.0, 2.0]" - @test ∂((1.0, 2.0), (1, 2)) ≡ ∂(SVector(1.0, 2.0), (1, 2)) ≡ ∂(∂(Val(2), (1, 2)), (1.0, 2.0)) - ∂o = ∂Output((1.0, 2.0)) - @test repr(∂o) == "SpectralKit.∂Output(1.0, 2.0)" - @test ∂o[1] == 1.0 - @test length(∂o) == 2 - @test Tuple(∂o) == (1.0, 2.0) +@testset "𝑑Derivatives" begin + @test_throws ArgumentError 𝑑^-1 + @test 𝑑^3 == 𝑑Derivatives{3}() + @test 𝑑^3 * 𝑑 == 𝑑Derivatives{4}() end -@testset "univariate derivatives check" begin - z = 0.6 - x = derivatives(z, Val(2)) - @test x[0] == z - @test x[1] == 1 - @test x[2] == 0 - b = Chebyshev(EndpointGrid(), 4) - f(z) = sum(basis_at(b, z)) - bz = basis_at(b, derivatives(z, Val(2))) - iterator_sanity_checks(bz) - ∑b = reduce(_add, bz) - @test ∑b[0] ≈ f(z) - @test ∑b[1] ≈ DD(f, z) atol = 1e-8 - @test ∑b[2] ≈ DD(f, z, 2) atol = 1e-8 -end - -@testset "univariate derivatives sanity checks" begin - d = derivatives(0.5, Val(2)) - @test eltype(d) ≡ Float64 - @test repr(d) == "0.5 + 1.0⋅Δ + 0.0⋅Δ²" +@testset "𝑑Expansion" begin + N = 3 + d = 𝑑^3 + x = 0.5 + 𝑑x = @inferred d(x) + @test 𝑑x.coefficients == SVector(x, 1.0, ntuple(_ -> 0.0, N - 1)...) + f(x) = _sub(_add(_one(typeof(x)), _mul(2.0, x)), _mul(x, x)) + f𝑑x = @inferred f(𝑑x) + for i in 0:N + @test f𝑑x[i] ≈ DD(f, x, i) atol = 1e-8 + end end -@testset "univariate transformed derivatives" begin - D = 1 # derivatives up to this one - b = Chebyshev(EndpointGrid(), 4) - θ = randn(dimension(b)) - for t in (BoundedLinear(-2.0, 7.0), ) - for i in 1:100 - z = rand_pm1(i) - x = transform_from(PM1(), t, z) - ℓ = linear_combination(b ∘ t, θ) - Dx = @inferred ℓ(derivatives(x, Val(D))) - @test Dx[0] == ℓ(x) - for d in 1:D - @test DD(ℓ, x, d) ≈ Dx[d] atol = 1e-8 +function rand_Partials(max_length = 4, max_d = 5, zero_prob = 0.2) + l = rand(0:max_length) + if l == 0 + Partials(()) + else + d = rand(0:max_d, l) + for i in 1:l + if rand() < zero_prob + d[i] = 0 end end + if d[end] == 0 + d[end] += 1 + end + Partials(Tuple(d)) end end -@testset "Smolyak derivatives check" begin - N = 3 - b = smolyak_basis(Chebyshev, InteriorGrid(), SmolyakParameters(3), N) - t = coordinate_transformations((BoundedLinear(1.0, 2.7), - SemiInfRational(3.0, 0.7), - InfRational(1.7, 0.3))) - D = ∂(Val(3), (), (1,), (2, ), (3, ), (1, 2)) - D̃ = [(f, x) -> f(x), - (f, x) -> DD(x1 -> f((x1, x[2], x[3])), x[1]), - (f, x) -> DD(x2 -> f((x[1], x2, x[3])), x[2]), - (f, x) -> DD(x3 -> f((x[1], x[2], x3)), x[3]), - (f, x) -> DD(x1 -> DD(x2 -> f((x1, x2, x[3])), x[2]), x[1])] - bt = b ∘ t - θ = randn(dimension(bt)) - ℓ = linear_combination(bt, θ) - d = domain(b) - for i in 1:100 - z = [rand_pm1() for _ in 1:N] - x = transform_from(d, t, z) - if i ≤ 5 # just a few allocation tests - @test @ballocated($ℓ(∂($D, $x))) == 0 - end - ℓDx = @inferred ℓ(∂(D, x)) - for (i, a) in enumerate(Tuple(ℓDx)) - @test a ≈ D̃[i](ℓ, x) atol = 1e-3 # cross-derivatives: lower tolerance - end +@testset "Partials sanity checks and subset tests" begin + p123 = Partials((1, 2, 3)) + @test !_is_strict_subset(p123, p123) + @test _is_strict_subset(Partials((1, 2)), p123) + @test _is_strict_subset(Partials((1, 2, 2)), p123) + @test _is_strict_subset(Partials((1, 0, 2)), p123) + @test !_is_strict_subset(Partials((1, 0, 2, 1)), p123) + @test _is_strict_subset(Partials(()), p123) + @test !SpectralKit._is_strict_subset(Partials((1,3)), Partials((3,0,1))) + @test_throws ArgumentError Partials((1, 2, -1)) + + for _ in 1:100 + p1 = rand_Partials() + p2 = rand_Partials() + S1 = Set(_partials_canonical_expansion(Val(5), [p1])) + S2 = Set(_partials_canonical_expansion(Val(5), [p2])) + @test _is_strict_subset(p1, p2) == (S1 ⊆ S2 && S1 ≠ S2) end end -@testset "iteration for ∂Output" begin - o_src = (1.0, 2.0, 3.0) - ∂o = SpectralKit.∂Output(o_src) - for (o1, o2) in zip(o_src, ∂o) - @test o1 ≡ o2 +@testset "Partials total order and strict subset" begin + for _ in 1:1000 + p1 = rand_Partials() + p2 = rand_Partials() + @test isless(p1, p2) + isless(p2, p1) + isequal(p1, p2) == 1 + if _is_strict_subset(p1, p2) + @test isless(p1, p2) + end end - function f(o) - o1, o2, o3 = o - o1 + o2 + o3 +end + +@testset "random partial derivatives minimal canonical roundtrip" begin + for _ in 1:1000 + p_rand = [rand_Partials(3, 3, 0.9) for _ in 1:rand(1:20)] + p_min = _partials_minimal_representation(p_rand) + for _ in 1:50 + @test p_min == _partials_minimal_representation(p_min[randperm(length(p_min))]) + @test p_min == _partials_minimal_representation(p_rand[randperm(length(p_rand))]) + end + p_exp = Set(_partials_canonical_expansion(Val(5), p_min)) + p_exp2 = mapreduce(p -> Set(_partials_canonical_expansion(Val(5), [p])), ∪, p_rand) + @test p_exp == p_exp2 end - @test @inferred(f(∂o)) == sum(o_src) - @test @ballocated($f($∂o)) == 0 end + +# @testset "partial derivatives interface" begin +# @test ∂(2, (1, 2), (2, 2)) == ∂Specification{(1, 2)}(((1, 1), (0, 2))) +# @test ∂(3, (), (1, 1, 2), (3, 2)) == ∂Specification{(2, 1, 1)}(((0, 0, 0), (2, 1, 0), (0, 1, 1))) +# @test_throws ArgumentError ∂(3, (4, )) +# @test_throws ArgumentError ∂(3, (-1, )) +# ∂spec = ∂(2, (1, 2), (2, 2)) +# @test repr(∂spec) == "partial derivatives\n[1] ∂²f/∂x₁∂x₂\n[2] ∂²f/∂²x₂" +# @test repr(∂(∂spec, [1.0, 2.0])) == +# "partial derivatives\n[1] ∂²f/∂x₁∂x₂\n[2] ∂²f/∂²x₂\nat [1.0, 2.0]" +# @test ∂((1.0, 2.0), (1, 2)) ≡ ∂(SVector(1.0, 2.0), (1, 2)) ≡ ∂(∂(Val(2), (1, 2)), (1.0, 2.0)) +# ∂o = ∂Output((1.0, 2.0)) +# @test repr(∂o) == "SpectralKit.∂Output(1.0, 2.0)" +# @test ∂o[1] == 1.0 +# @test length(∂o) == 2 +# @test Tuple(∂o) == (1.0, 2.0) +# end + +# @testset "univariate derivatives check" begin +# z = 0.6 +# x = derivatives(z, Val(2)) +# @test x[0] == z +# @test x[1] == 1 +# @test x[2] == 0 +# b = Chebyshev(EndpointGrid(), 4) +# f(z) = sum(basis_at(b, z)) +# bz = basis_at(b, derivatives(z, Val(2))) +# iterator_sanity_checks(bz) +# ∑b = reduce(_add, bz) +# @test ∑b[0] ≈ f(z) +# @test ∑b[1] ≈ DD(f, z) atol = 1e-8 +# @test ∑b[2] ≈ DD(f, z, 2) atol = 1e-8 +# end + +# @testset "univariate derivatives sanity checks" begin +# d = derivatives(0.5, Val(2)) +# @test eltype(d) ≡ Float64 +# @test repr(d) == "0.5 + 1.0⋅Δ + 0.0⋅Δ²" +# end + +# @testset "univariate transformed derivatives" begin +# D = 1 # derivatives up to this one +# b = Chebyshev(EndpointGrid(), 4) +# θ = randn(dimension(b)) +# for t in (BoundedLinear(-2.0, 7.0), ) +# for i in 1:100 +# z = rand_pm1(i) +# x = transform_from(PM1(), t, z) +# ℓ = linear_combination(b ∘ t, θ) +# Dx = @inferred ℓ(derivatives(x, Val(D))) +# @test Dx[0] == ℓ(x) +# for d in 1:D +# @test DD(ℓ, x, d) ≈ Dx[d] atol = 1e-8 +# end +# end +# end +# end + +# @testset "Smolyak derivatives check" begin +# N = 3 +# b = smolyak_basis(Chebyshev, InteriorGrid(), SmolyakParameters(3), N) +# t = coordinate_transformations((BoundedLinear(1.0, 2.7), +# SemiInfRational(3.0, 0.7), +# InfRational(1.7, 0.3))) +# D = ∂(Val(3), (), (1,), (2, ), (3, ), (1, 2)) +# D̃ = [(f, x) -> f(x), +# (f, x) -> DD(x1 -> f((x1, x[2], x[3])), x[1]), +# (f, x) -> DD(x2 -> f((x[1], x2, x[3])), x[2]), +# (f, x) -> DD(x3 -> f((x[1], x[2], x3)), x[3]), +# (f, x) -> DD(x1 -> DD(x2 -> f((x1, x2, x[3])), x[2]), x[1])] +# bt = b ∘ t +# θ = randn(dimension(bt)) +# ℓ = linear_combination(bt, θ) +# d = domain(b) +# for i in 1:100 +# z = [rand_pm1() for _ in 1:N] +# x = transform_from(d, t, z) +# if i ≤ 5 # just a few allocation tests +# @test @ballocated($ℓ(∂($D, $x))) == 0 +# end +# ℓDx = @inferred ℓ(∂(D, x)) +# for (i, a) in enumerate(Tuple(ℓDx)) +# @test a ≈ D̃[i](ℓ, x) atol = 1e-3 # cross-derivatives: lower tolerance +# end +# end +# end + +# @testset "iteration for ∂Output" begin +# o_src = (1.0, 2.0, 3.0) +# ∂o = SpectralKit.∂Output(o_src) +# for (o1, o2) in zip(o_src, ∂o) +# @test o1 ≡ o2 +# end +# function f(o) +# o1, o2, o3 = o +# o1 + o2 + o3 +# end +# @test @inferred(f(∂o)) == sum(o_src) +# @test @ballocated($f($∂o)) == 0 +# end diff --git a/test/utilities.jl b/test/utilities.jl index 3bcc5ec..eec3601 100644 --- a/test/utilities.jl +++ b/test/utilities.jl @@ -76,9 +76,11 @@ function iterator_sanity_checks(itr) @test count(_ -> true, itr) == length(itr) end -"Derivative of f at x." +"nth derivative of f at x." function DD(f, x, n = 1; p = 10) - # workaround for https://github.com/JuliaDiff/FiniteDifferences.jl/issues/224 - # remove anonymous call when that is fixed - central_fdm(p, n)(x -> f(x), x) + if n == 0 + f(x) + else + central_fdm(p, n)(f, x) + end end From b797dba36c585aa4790ad81835966470dd65612f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tam=C3=A1s=20K=2E=20Papp?= Date: Tue, 21 May 2024 17:12:02 +0200 Subject: [PATCH 02/14] remove unused code --- src/derivatives.jl | 53 ---------------------------------------------- 1 file changed, 53 deletions(-) diff --git a/src/derivatives.jl b/src/derivatives.jl index a790902..29a21d1 100644 --- a/src/derivatives.jl +++ b/src/derivatives.jl @@ -245,59 +245,6 @@ function _partials_canonical_expansion(::Val{N}, partials) where N result end - - -# const _PARTIALS = Tuple{Vararg{Int}} - -# function _partials_max(D1::NTuple{N1,Int}, D2::NTuple{N2,Int}) where {N1,N2} -# ntuple(Val(max(N1,N2))) do i -# if i ≤ N1 -# if i ≤ N2 -# max(D1[i], D2[i]) -# else -# D1[i] -# end -# else -# D2[i] -# end -# end -# end - -# function _partials_isless_or_eq(D1::NTuple{N1,Int}, D2::NTuple{N2,Int}) where {N1,N2} -# N1 ≤ N2 || return false -# for i in 1:N1 -# D1[i] ≤ D2[i] || return false -# end -# true -# end - -# function _normalize_partials(D::Tuple) -# i = something(findlast(d -> d > 0, D), 0) -# D[1:i] -# end - -# function _expand_M_Ds_union(k_d_pairs) -# UDs = Vector{_PARTIALS}() -# M = () -# for D in Ds -# M = _partials_max(M, D) -# for i in CartesianIndices(map(d -> 0:d, D)) -# _sorted_insert_new!(UDs, _normalize_partials(Tuple(i))) -# end -# end -# M, UDs -# end - -# @generated function _check_M_Ds_consistency(::Val{N}, ::Val{Ds}) where {N,Ds} -# @argcheck Ds isa Tuple{Vararg{_PARTIALS}} && all(D -> all(d -> d ≥ 0, D), Ds) -# @argcheck M isa NTuple{K,Int} && all(m -> m ≥ 0, M) -# M′, Ds′ = _calculate_M_Ds_union(Ds) -# C = M′ ≡ M && Ds′ ≡ Ds -# quote -# Val($(C)) -# end -# end - # struct ∂Derivatives{K,M,Ds} # function ∂Derivatives{K,M,Ds}() where {K,M,Ds} # @argcheck K isa Int && K ≥ 0 From 49fa3b81f8d1a486fc74dff2aecfbc75a09db096 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tam=C3=A1s=20K=2E=20Papp?= Date: Wed, 22 May 2024 12:25:28 +0200 Subject: [PATCH 03/14] re-enable rest of package, minor rearrangements --- src/SpectralKit.jl | 12 +- src/derivatives.jl | 249 ++++++++++++++++----------------- src/generic_api.jl | 2 +- src/smolyak_api.jl | 165 +--------------------- src/smolyak_traversal.jl | 157 +++++++++++++++++++++ src/transformations.jl | 30 ++-- test/runtests.jl | 13 +- test/test_chebyshev.jl | 2 +- test/test_derivatives.jl | 24 +++- test/test_domains.jl | 7 +- test/test_smolyak.jl | 123 +--------------- test/test_smolyak_traversal.jl | 127 +++++++++++++++++ 12 files changed, 463 insertions(+), 448 deletions(-) create mode 100644 test/test_smolyak_traversal.jl diff --git a/src/SpectralKit.jl b/src/SpectralKit.jl index f2d6000..e794d09 100644 --- a/src/SpectralKit.jl +++ b/src/SpectralKit.jl @@ -7,11 +7,11 @@ using StaticArrays: MVector, SVector, sacollect include("utilities.jl") include("derivatives.jl") -# include("domains.jl") -# include("transformations.jl") -# include("generic_api.jl") -# include("chebyshev.jl") -# include("smolyak_traversal.jl") -# include("smolyak_api.jl") +include("domains.jl") +include("transformations.jl") +include("generic_api.jl") +include("chebyshev.jl") +include("smolyak_traversal.jl") +include("smolyak_api.jl") end # module diff --git a/src/derivatives.jl b/src/derivatives.jl index 29a21d1..f557515 100644 --- a/src/derivatives.jl +++ b/src/derivatives.jl @@ -227,14 +227,14 @@ function _partials_minimal_representation(partials) minimal_partials end -function _partials_canonical_expansion(::Val{N}, partials) where N +function _partials_canonical_expansion(::Val{N}, Ps) where N result = OrderedSet{NTuple{N,Int}}() function _plus1_cartesian_indices(p::Partials{M}) where M (; I) = p @argcheck M ≤ N CartesianIndices(ntuple(i -> i ≤ M ? I[i] + 1 : 1, Val(N))) end - for p in partials + for p in Ps for ι in _plus1_cartesian_indices(p) i = map(i -> i - 1, Tuple(ι)) if !(i ∈ result) @@ -245,107 +245,106 @@ function _partials_canonical_expansion(::Val{N}, partials) where N result end -# struct ∂Derivatives{K,M,Ds} -# function ∂Derivatives{K,M,Ds}() where {K,M,Ds} -# @argcheck K isa Int && K ≥ 0 -# @argcheck _check_M_Ds_consistency(Val(N), Val(Ds)) ≡ Val(true) -# new{K,M,Ds}() -# end -# function ∂Derivatives(dimension::Val{K}, degree::Val{m}) where {K,m} -# @argcheck K isa Int && K ≥ 1 -# @argcheck m isa Int && m ≥ 0 -# m == 0 && return new{0,(),((),)}() -# z = ntuple(_ -> 0, Val(K - 1)) -# M = (z..., m) -# Ds = ntuple(i -> (z..., i - 1), Val(m + 1)) -# new{K,M,Ds}() -# end -# end +function _partials_expansion_degrees(::Val{N}, partials) where N + degrees = zero(MVector{N,Int}) + for P in partials + for (j, i) in enumerate(P.I) + degrees[j] = max(degrees[j], i) + end + end + Tuple(degrees) +end -# ∂(k::Val{K}, m::Val{M}) where {K,M} = ∂Derivatives(k, m) +_partials_minimum_input_dimension(partials) = maximum(P -> length(P.I), partials) -# @inline ∂(dimension::Int, degree::Int = 1) = ∂(Val(dimension), Val(degree)) +function _is_minimal_representation(::Val{Ps}) where Ps + _Ps = fieldtypes(Ps) + _Ps isa Tuple{Vararg{Partials}} || return false + _partials_minimal_representation(_Ps) == collect(_Ps) +end -# @generated function Base.:∪(∂ds::∂Derivatives...) -# _get_Ds(::∂Derivatives{K,N,Ds}) where {K,N,Ds} = Ds -# M, Ds = _calculate_M_Ds_union(Iterators.flatten((_get_Ds(∂d) for ∂d in ∂ds))) -# K = length(M) -# quote -# ∂Derivatives{$(K), $(M), $(Ds)}() -# end -# end +struct ∂Derivatives{Ps} + function ∂Derivatives{Ps}() where {Ps} + if Ps isa Partials + new{Tuple{Ps}}() + else + @argcheck _is_minimal_representation(Val(Ps)) + new{Ps}() + end + end +end -# function Base.:^(::∂Derivatives{K,M,D}, ::Val{Y}) where {K,M,D,Y} -# @argcheck Y isa Int && Y ≥ 0 -# (dhead..., dlast) = D -# @argcheck(all(iszero, dhead...), -# "The ^ operator can only be applied for ∂ derivatives along one dimension.") -# ∂Derivatives(Val(K), Val(M * Y)) -# end +∂(I::Tuple{Vararg{Int}}) = ∂Derivatives{Partials(I)}() -# function _collapse_D_union!(D::Vector{_PARTIALS}) -# D = collect(D) -# U = Vector{Pair{Int,Int}}() -# while !isempty(D) -# M = reduce(_partials_max, D) -# k = findfirst(d -> d > 0, M) -# k ≡ nothing && break # we are done -# m = M[k] -# push!(U, k => m) -# _d = (ntuple(_ -> 0, k - 1)..., m) -# filter!(d -> !_partials_isless_or_eq(d, _d), D) -# end -# U -# end +∂(I::Integer...) = ∂Derivatives{Partials(I...)}() -# function Base.show(io::IO, ::∂Derivatives{K,M,D}) where {K,M,D} -# print(io, "partial derivatives (at least $(K) dimensions), up to ") -# join(io, ("∂($(k), $(m))" for (k, m) in _collapse_D_union!(collect(_PARTIALS, D))), " ∪ ") -# end +∂() = ∂Derivatives{Tuple{}}() -# struct ∂DerivativesInput{P<:∂Derivatives,N,T} -# partial_derivatives::P -# x::SVector{N,T} -# function ∂DerivativesInput(∂derivatives::P, -# x::SVector{N,T}) where {K,P<:∂Derivatives{K},N,T} -# @argcheck N ≤ K "Can't differentiate the $(K)th element of a $(N)-element input." -# new{P,N,T}(∂derivatives, x) -# end -# end +function Base.:<<(d::𝑑Derivatives{D}, ::Val{N}) where {D,N} + @argcheck N isa Int && N ≥ 1 + @argcheck D ≥ 1 + ∂Derivatives{Partials(ntuple(i -> 0, Val(N - 1))..., D)}() +end -# (∂derivatives::∂Derivatives)(x::SVector) = ∂DerivativesInput(partial_derivatives, x) +@generated function Base.:∪(∂ds::∂Derivatives...) + _get_Ps(::Type{<:∂Derivatives{Ps}}) where {Ps} = [fieldtypes(Ps)...] + Ps = mapreduce(_get_Ps, vcat, ∂ds) + quote + ∂Derivatives{Tuple{$(_partials_minimal_representation(Ps))...}}() + end +end -# """ -# See [`_lift`](@ref). Internal. -# """ -# struct ∂DerivativesInputLifted{P<:∂Derivatives,L<:Tuple} -# partial_derivatives::P -# lifted_x::L -# end +function Base.show(io::IO, ::∂Derivatives{Ps}) where {Ps} + _Ps = fieldtypes(Ps) + _repr(P::Partials) = "∂($(join(P.I, ", ")))" + if isempty(_Ps) + print(io, "∂()") + elseif length(_Ps) == 1 + print(io, _repr(_Ps[1])) + else + print(io, "union(") + join(io, (_repr(P) for P in fieldtypes(Ps)), ", ") + print(io, ")") + end +end -# """ -# $(SIGNATURES) +@generated function _expand_coordinates(::∂Derivatives{Ps}, x::NTuple{N}) where {Ps,N} + _Ps = fieldtypes(Ps) + @argcheck N ≥ _partials_minimum_input_dimension(_Ps) + M = _partials_expansion_degrees(Val(N), _Ps) + x = [:(𝑑Derivatives{$(m)}()(x[$(j)])) for (j, m) in enumerate(M)] + quote + tuple($(x...)) + end +end -# Lift a partial derivative calculation into a tuple of `Derivatives`. Internal. -# """ -# @generated function _lift(∂x::∂DerivativesInput{<:∂Derivatives{K,M,N}}) where {K,M,N} -# _lifted_x = [:(∂(Val($(i ≤ K ? M[i] : 0)))(x[$(i)])) for i in 1:N] -# quote -# x = ∂x.x -# lifted_x = ($(_lifted_x...),) -# ∂PartialDerivativesInputLifted(∂x.∂derivatives, lifted_x) -# end -# end +struct ∂CoordinateExpansion{D<:∂Derivatives,S<:Tuple} + ∂D::D + x::S + function ∂CoordinateExpansion(∂D::D, x::S) where {D<:∂Derivatives,S<:Tuple} + new{D,S}(∂D, x) + end +end -# struct ∂Expansion{P<:∂Derivatives,N,T} -# ∂derivatives::P -# coefficients::SVector{N,T} -# function ∂Expansion(∂derivatives::∂Derivatives{K,M,D}, -# coefficients::SVector{N,T}) where {K,M,D,N,T} -# @argcheck length(D) == N -# new{typeof(P),N,T}(∂derivatives, coefficients) -# end -# end +(∂D::∂Derivatives)(x::Tuple) = ∂CoordinateExpansion(∂D, _expand_coordinates(∂D, x)) + +(∂D::∂Derivatives)(x::AbstractVector) = ∂CoordinateExpansion(∂D, Tuple(x)) + +struct ∂Expansion{D,N,T} + ∂D::D + coefficients::SVector{N,T} + function ∂Expansion(∂D::D, coefficients::SVector{N,T}) where {D<:∂Derivatives,N,T} + new{D,N,T}(coefficients) + end +end + +function _add(x::∂Expansion{D,N}, y::∂Expansion{D,N}) where {D,N} + ∂Expansion(x.D, map(+, x.values, y.values)) +end + +function _mul(x::Real, y::∂Expansion) + ∂Expansion(y.D, map(y -> _mul(x, y), y.coefficients)) +end # function Base.show(io::IO, expansion::∂Expansion{<:∂Derivatives{K,M,D}}) where {K,M,D} # (; coefficients) = expansion @@ -357,51 +356,39 @@ end # end # end -# """ -# $(SIGNATURES) - -# Conceptually equivalent to `prod(getindex.(sources, indices))`, which it returns when -# `kind` is `nothing`, a placeholder calculating any derivatives. Internal. -# """ -# _product(kind::Nothing, sources, indices) = mapreduce(getindex, *, sources, indices) - -# """ -# $(SIGNATURES) +""" +$(SIGNATURES) -# Type that is returnedby [`_product`](@ref). -# """ -# function _product_type(::Type{Nothing}, source_eltypes) -# mapfoldl(eltype, promote_type, source_eltypes) -# end +Conceptually equivalent to `prod(x))`, which it returns when `kind` is `nothing`, a +placeholder calculating any derivatives. Internal. +""" +_product(kind::Nothing, x::Tuple) = prod(x) -# @generated function _product(∂derivatives::∂Derivatives{K,M,D}, -# sources::NTuple{N}, -# indices) where {K,M,D,N} -# xs = [gensym(:x) for _ in 1:N] -# assignments = [:($(xs[i]) = sources[indices[$(i)]]) for i in 1:N] -# products = [begin -# ι = NTuple(i -> i ≤ K ? d[i] + 1 : 1, Val(N)) -# mapreduce(i -> :($(xs[i])[$(ι[i])]), -# (a, b) -> :($(a) * $(b)), -# 1:N) -# end for d in D] -# quote -# $(assignments...) -# ∂Expansion(∂derivatives, SVector($(products)...)) -# end -# end +""" +$(SIGNATURES) -# function _product_type(::Type{∂Derivatives{M,L}}, source_eltypes) where {M,L} -# T = _product_type(Nothing, map(eltype, source_eltypes)) -# N = length(fieldtypes(L)) -# ∂Output{N,T} -# end +Type that is returnedby [`_product`](@ref). +""" +function _product_type(::Type{Nothing}, source_eltypes) + mapfoldl(eltype, promote_type, source_eltypes) +end -# function _add(x::∂Expansion{P}, y::∂Expansion{P}) where P -# ∂Expansion(x.∂derivatives, map(+, x.values, y.values)) -# end +@generated function _product(∂D::∂Derivatives{Ps}, x::NTuple{N,𝑑Expansion}) where {Ps,N} + function _product(d) + # FIXME could skip bounds checking if verified at the beginning + mapreduce(i -> :(x[$i].coefficients[$(d[i]) + 1]), (a, b) -> :($(a) * $(b)), 1:N) + end + products = [_product(d) for d in _partials_canonical_expansion(Val(N), fieldtypes(Ps))] + quote + ∂Expansion(∂derivatives, SVector($(products)...)) + end +end -# _mul(x::Real, y::∂Expansion) = ∂Expansion(y.∂derivatives, map(y -> _mul(x, y), y.values)) +function _product_type(::Type{D}, source_eltypes) where {Ps,D<:∂Derivatives{Ps}} + T = _product_type(Nothing, map(eltype, source_eltypes)) + N = length(_partials_canonical_expansion(Val(N), fieldtypes(Ps))) + ∂Expansion{D,N,T} +end ##### ##### FIXME revise and move documentation below diff --git a/src/generic_api.jl b/src/generic_api.jl index ca41226..0eecf2d 100644 --- a/src/generic_api.jl +++ b/src/generic_api.jl @@ -2,7 +2,7 @@ ##### Generic API ##### -export is_function_basis, domain, dimension, basis_at, linear_combination, InteriorGrid, +export is_function_basis, dimension, basis_at, linear_combination, InteriorGrid, InteriorGrid2, EndpointGrid, grid, collocation_matrix, augment_coefficients, is_subset_basis diff --git a/src/smolyak_api.jl b/src/smolyak_api.jl index f231647..bb4c110 100644 --- a/src/smolyak_api.jl +++ b/src/smolyak_api.jl @@ -2,158 +2,7 @@ ##### Smolyak bases ##### -export SmolyakParameters, smolyak_basis - -struct SmolyakParameters{B,M} - function SmolyakParameters{B,M}() where {B,M} - @argcheck B isa Int && B ≥ 0 - @argcheck M isa Int && M ≥ 0 - M > B && @warn "M > B replaced with M = B" M B - new{B,min(B,M)}() # maintain M ≤ B - end -end - -function Base.show(io::IO, ::SmolyakParameters{B,M}) where {B,M} - print(io, "Smolyak parameters, ∑bᵢ ≤ $(B), all bᵢ ≤ $(M)") -end - -""" -$(SIGNATURES) - -Parameters for Smolyak grids that are independent of the dimension of the domain. - -Polynomials are organized into blocks (of eg `1, 2, 2, 4, 8, 16, …`) polynomials (and -corresponding gridpoints), indexed with a *block index* `b` that starts at `0`. `B ≥ ∑ -bᵢ` and `0 ≤ bᵢ ≤ M` constrain the number of blocks along each dimension `i`. - -`M > B` is not an error, but will be normalized to `M = B` with a warning. -""" -@inline function SmolyakParameters(B::Integer, M::Integer = B) - SmolyakParameters{Int(B),Int(M)}() -end - -""" -$(TYPEDEF) - -Indexing specification in a Smolyak basis/interpolation. - -# Type parameters - -- `N`: the dimension of indices - -- `H`: highest index visited for all dimensions - -- `B ≥ 0`: sum of block indices, starting from `0` (ie `B = 0` has just one element), - -- `M`: upper bound on each block index - -# Constructor - -Takes the dimension `N` as a parameter, `grid_kind`, and a `SmolyakParameters` object, -calculating everything else. - -# Details - -Consider positive integer indices `(i1, …, iN)`, each starting at one. - -Let `ℓ(b) = nesting_total_length(Chebyshev, grid_knid, kind, b)`, and `b1` denote the -smallest integer such that `i1 ≤ ℓ(b1)`, and similarly for `i2, …, iN`. Extend this with -`ℓ(-1) = 0` for the purposes of notation. - -An index `(i1, …, iN)` is visited iff all of the following hold: - -1. `1 ≤ i1 ≤ ℓ(M)`, …, `1 ≤ iN ≤ ℓ(M)`, -2. `0 ≤ b1 ≤ M`, …, `1 ≤ bN ≤ M`, -3. `b1 + … + bN ≤ B` - -Visited indexes are in *column-major* order. -""" -struct SmolyakIndices{N,H,B,M,Mp1} - "number of coefficients (cached)" - len::Int - "nesting total lengths (cached)" - nesting_total_lengths::NTuple{Mp1,Int} - function SmolyakIndices{N}(grid_kind::AbstractGrid, - smolyak_parameters::SmolyakParameters{B,M}) where {N,B,M} - @argcheck N ≥ 1 - Mp1 = M + 1 - len = __smolyak_length(grid_kind, Val(N), Val(B), M) - first_block_length = nesting_total_length(Chebyshev, grid_kind, 0) - nesting_total_lengths = ntuple(bp1 -> nesting_total_length(Chebyshev, grid_kind, bp1 - 1), - Val(Mp1)) - H = last(nesting_total_lengths) - new{N,H,B,M,Mp1}(len, nesting_total_lengths) - end -end - -function Base.show(io::IO, smolyak_indices::SmolyakIndices{N,H,B,M}) where {N,H,B,M} - (; len) = smolyak_indices - print(io, "Smolyak indexing, ∑bᵢ ≤ $(B), all bᵢ ≤ $(M), dimension $(len)") -end - -@inline highest_visited_index(::SmolyakIndices{N,H}) where {N,H} = H - -Base.eltype(::Type{<:SmolyakIndices{N}}) where N = NTuple{N,Int} - -@inline Base.length(ι::SmolyakIndices) = ι.len - -@inline function Base.iterate(ι::SmolyakIndices{N,H,B}) where {N,H,B} - slack, indices, blocks, limits = __inc_init(ι.nesting_total_lengths, Val(N), Val(B)) - indices, (slack, indices, blocks, limits) -end - -@inline function Base.iterate(ι::SmolyakIndices, (slack, indices, blocks, limits)) - valid, Δ, indices′, blocks′, limits′ = __inc(ι.nesting_total_lengths, slack, indices, - blocks, limits) - valid || return nothing - slack′ = slack + Δ - indices′, (slack′, indices′, blocks′, limits′) -end - -#### -#### product traversal -#### - -struct SmolyakProduct{I<:SmolyakIndices,S<:Tuple,P} - smolyak_indices::I - sources::S - product_kind::P - @doc """ - $(SIGNATURES) - - An iterator conceptually equivalent to - - ``` - [prod(getindex.(sources, indices)) for indices in smolyak_indices] - ``` - - using [`_product`](@ref) instead to account for derivatives. Detailed docs of the - arguments are in [`SmolyakIndices`](@ref). - - Caller should arrange the elements of `sources` in the correct order, see - [`nested_extrema_indices`](@ref). Each element in `sources` should have at least - `H` elements (cf type parameters of [`SmolyakIndices`](@ref)), this is not checked. - """ - function SmolyakProduct(smolyak_indices::I, sources::S, - product_kind::P) where {N,I<:SmolyakIndices{N},S,P} - @argcheck length(sources) == N - new{I,S,P}(smolyak_indices, sources, product_kind) - end -end - -Base.length(smolyak_product::SmolyakProduct) = length(smolyak_product.smolyak_indices) - -function Base.eltype(::Type{SmolyakProduct{I,S,P}}) where {I,S,P} - _product_type(P, fieldtypes(S)) -end - -@inline function Base.iterate(smolyak_product::SmolyakProduct, state...) - (; smolyak_indices, sources, product_kind) = smolyak_product - itr = iterate(smolyak_indices, state...) - itr ≡ nothing && return nothing - indices, state′ = itr - _product(product_kind, sources, indices), state′ -end +export smolyak_basis struct SmolyakBasis{I<:SmolyakIndices,U<:UnivariateBasis} <: MultivariateBasis smolyak_indices::I @@ -252,16 +101,14 @@ function basis_at(smolyak_basis::SmolyakBasis{<:SmolyakIndices{N}}, end function basis_at(smolyak_basis::SmolyakBasis{<:SmolyakIndices{N}}, - Lx::∂InputLifted) where {N} - (; ∂specification, lifted_x) = Lx - @argcheck length(lifted_x) == N + Dx::∂CoordinateExpansion) where {N} + (; ∂D, x) = Dx + @argcheck length(x) == N SmolyakProduct(smolyak_basis.smolyak_indices, - _univariate_bases_at(smolyak_basis, lifted_x), - ∂specification) + _univariate_bases_at(smolyak_basis, x), + ∂D) end -basis_at(smolyak_basis::SmolyakBasis, ∂x::∂Input) = basis_at(smolyak_basis, _lift(∂x)) - struct SmolyakGridIterator{T,I,S} smolyak_indices::I sources::S diff --git a/src/smolyak_traversal.jl b/src/smolyak_traversal.jl index 5118bf6..13b665a 100644 --- a/src/smolyak_traversal.jl +++ b/src/smolyak_traversal.jl @@ -2,6 +2,40 @@ ##### Smolyak implementation details ##### +### +### Blocking parameters +### + +export SmolyakParameters + +struct SmolyakParameters{B,M} + function SmolyakParameters{B,M}() where {B,M} + @argcheck B isa Int && B ≥ 0 + @argcheck M isa Int && M ≥ 0 + M > B && @warn "M > B replaced with M = B" M B + new{B,min(B,M)}() # maintain M ≤ B + end +end + +function Base.show(io::IO, ::SmolyakParameters{B,M}) where {B,M} + print(io, "Smolyak parameters, ∑bᵢ ≤ $(B), all bᵢ ≤ $(M)") +end + +""" +$(SIGNATURES) + +Parameters for Smolyak grids that are *independent of the dimension of the domain*. + +Polynomials are organized into blocks (of eg `1, 2, 2, 4, 8, 16, …`) polynomials (and +corresponding gridpoints), indexed with a *block index* `b` that starts at `0`. `B ≥ ∑ +bᵢ` and `0 ≤ bᵢ ≤ M` constrain the number of blocks along each dimension `i`. + +`M > B` is not an error, but will be normalized to `M = B` with a warning. +""" +@inline function SmolyakParameters(B::Integer, M::Integer = B) + SmolyakParameters{Int(B),Int(M)}() +end + #### #### Nesting sizes and shuffling #### @@ -220,3 +254,126 @@ function __smolyak_length(grid_kind::AbstractGrid, ::Val{N}, ::Val{B}, M::Int) w end sum(c) end + +""" +$(TYPEDEF) + +Indexing specification in a Smolyak basis/interpolation. + +# Type parameters + +- `N`: the dimension of indices + +- `H`: highest index visited for all dimensions + +- `B ≥ 0`: sum of block indices, starting from `0` (ie `B = 0` has just one element), + +- `M`: upper bound on each block index + +# Constructor + +Takes the dimension `N` as a parameter, `grid_kind`, and a `SmolyakParameters` object, +calculating everything else. + +# Details + +Consider positive integer indices `(i1, …, iN)`, each starting at one. + +Let `ℓ(b) = nesting_total_length(Chebyshev, grid_knid, kind, b)`, and `b1` denote the +smallest integer such that `i1 ≤ ℓ(b1)`, and similarly for `i2, …, iN`. Extend this with +`ℓ(-1) = 0` for the purposes of notation. + +An index `(i1, …, iN)` is visited iff all of the following hold: + +1. `1 ≤ i1 ≤ ℓ(M)`, …, `1 ≤ iN ≤ ℓ(M)`, +2. `0 ≤ b1 ≤ M`, …, `1 ≤ bN ≤ M`, +3. `b1 + … + bN ≤ B` + +Visited indexes are in *column-major* order. +""" +struct SmolyakIndices{N,H,B,M,Mp1} + "number of coefficients (cached)" + len::Int + "nesting total lengths (cached)" + nesting_total_lengths::NTuple{Mp1,Int} + function SmolyakIndices{N}(grid_kind::AbstractGrid, + smolyak_parameters::SmolyakParameters{B,M}) where {N,B,M} + @argcheck N ≥ 1 + Mp1 = M + 1 + len = __smolyak_length(grid_kind, Val(N), Val(B), M) + first_block_length = nesting_total_length(Chebyshev, grid_kind, 0) + nesting_total_lengths = ntuple(bp1 -> nesting_total_length(Chebyshev, grid_kind, bp1 - 1), + Val(Mp1)) + H = last(nesting_total_lengths) + new{N,H,B,M,Mp1}(len, nesting_total_lengths) + end +end + +function Base.show(io::IO, smolyak_indices::SmolyakIndices{N,H,B,M}) where {N,H,B,M} + (; len) = smolyak_indices + print(io, "Smolyak indexing, ∑bᵢ ≤ $(B), all bᵢ ≤ $(M), dimension $(len)") +end + +@inline highest_visited_index(::SmolyakIndices{N,H}) where {N,H} = H + +Base.eltype(::Type{<:SmolyakIndices{N}}) where N = NTuple{N,Int} + +@inline Base.length(ι::SmolyakIndices) = ι.len + +@inline function Base.iterate(ι::SmolyakIndices{N,H,B}) where {N,H,B} + slack, indices, blocks, limits = __inc_init(ι.nesting_total_lengths, Val(N), Val(B)) + indices, (slack, indices, blocks, limits) +end + +@inline function Base.iterate(ι::SmolyakIndices, (slack, indices, blocks, limits)) + valid, Δ, indices′, blocks′, limits′ = __inc(ι.nesting_total_lengths, slack, indices, + blocks, limits) + valid || return nothing + slack′ = slack + Δ + indices′, (slack′, indices′, blocks′, limits′) +end + +#### +#### product traversal +#### + +struct SmolyakProduct{I<:SmolyakIndices,S<:Tuple,P} + smolyak_indices::I + sources::S + product_kind::P + @doc """ + $(SIGNATURES) + + An iterator conceptually equivalent to + + ``` + [prod(getindex.(sources, indices)) for indices in smolyak_indices] + ``` + + using [`_product`](@ref) instead to account for derivatives. Detailed docs of the + arguments are in [`SmolyakIndices`](@ref). + + Caller should arrange the elements of `sources` in the correct order, see + [`nested_extrema_indices`](@ref). Each element in `sources` should have at least + `H` elements (cf type parameters of [`SmolyakIndices`](@ref)), this is not checked. + """ + function SmolyakProduct(smolyak_indices::I, sources::S, + product_kind::P) where {N,I<:SmolyakIndices{N},S,P} + @argcheck length(sources) == N + new{I,S,P}(smolyak_indices, sources, product_kind) + end +end + +Base.length(smolyak_product::SmolyakProduct) = length(smolyak_product.smolyak_indices) + +function Base.eltype(::Type{SmolyakProduct{I,S,P}}) where {I,S,P} + _product_type(P, fieldtypes(S)) +end + +@inline function Base.iterate(smolyak_product::SmolyakProduct, state...) + (; smolyak_indices, sources, product_kind) = smolyak_product + itr = iterate(smolyak_indices, state...) + itr ≡ nothing && return nothing + indices, state′ = itr + _product(product_kind, map(getindex, sources, indices)), state′ +end diff --git a/src/transformations.jl b/src/transformations.jl index a52b617..961e73d 100644 --- a/src/transformations.jl +++ b/src/transformations.jl @@ -2,10 +2,9 @@ ##### transformations ##### -export transform_to, transform_from, coordinate_transformations, +export domain, domain_kind, transform_to, transform_from, coordinate_transformations, BoundedLinear, InfRational, SemiInfRational - #### #### generic api #### @@ -118,14 +117,9 @@ function transform_to(domain::CoordinateDomains{T}, ct::CoordinateTransformation SVector(transform_to(domain, ct, _ntuple_like(T, x))) end -function transform_to(domain::CoordinateDomains, ct::CoordinateTransformations, ∂x::∂Input) - transform_to(domain, ct, _lift(∂x)) -end - function transform_to(domain::CoordinateDomains, ct::CoordinateTransformations, - ∂x::∂InputLifted) - (; ∂specification, lifted_x) = ∂x - ∂InputLifted(∂specification, transform_to(domain, ct, lifted_x)) + Dx::∂CoordinateExpansion) + ∂CoordinateExpansion(x.∂D, transform_to(domain, ct, x.x)) end function transform_from(domain::CoordinateDomains, ct::CoordinateTransformations, x::Tuple) @@ -187,13 +181,13 @@ function transform_to(::PM1, t::BoundedLinear, y::Real) (y - m) / s end -function transform_to(domain::PM1, t::BoundedLinear, y::Derivatives{N}) where N +function transform_to(domain::PM1, t::BoundedLinear, y::𝑑Expansion{Dp1}) where Dp1 (; m, s) = t (; derivatives) = y y0, yD... = derivatives x0 = transform_to(domain, t, y0) xD = map(y -> y / s, yD) - Derivatives((x0, xD...)) + 𝑑Expansion((x0, xD...)) end function domain(t::BoundedLinear) @@ -256,16 +250,16 @@ function transform_to(::PM1, t::SemiInfRational, y::Real) end end -function transform_to(domain::PM1, t::SemiInfRational, y::Derivatives{N}) where N +function transform_to(domain::PM1, t::SemiInfRational, y::𝑑Expansion{Dp1}) where Dp1 (; A, L) = t (; derivatives) = y x0 = transform_to(domain, t, derivatives[1]) - N == 1 && return Derivatives((x0, )) + Dp1 == 1 && return 𝑑Expansion(SVector(x0)) # based on Boyd (2001), Table E.7 Q = abs2(x0 - 1) x1 = (derivatives[2] * Q) / (2*L) - N == 2 && return Derivatives((x0, x1)) - error("$(N-1)th derivative not implemented yet, open an issue.") + Dp1 == 2 && return 𝑑Expansion(SVector(x0, x1)) + error("$(Dp1-1)th derivative not implemented yet, open an issue.") end function domain(t::SemiInfRational) @@ -321,16 +315,16 @@ function transform_to(::PM1, t::InfRational, y::Real) end end -function transform_to(domain::PM1, t::InfRational, y::Derivatives{N}) where N +function transform_to(domain::PM1, t::InfRational, y::𝑑Expansion{Dp1}) where Dp1 (; A, L) = t (; derivatives) = y x0 = transform_to(domain, t, derivatives[1]) - N == 1 && return Derivatives((x0, )) + N == 1 && return Derivatives(SVector(x0)) # based on Boyd (2001), Table E.5 Q = 1 - abs2(x0) sQ = √Q x1 = (derivatives[2] * Q * sQ) / L - N == 2 && return Derivatives((x0, x1)) + N == 2 && return 𝑑Expansion(SVector(x0, x1)) error("$(N-1)th derivative not implemented yet, open an issue.") end diff --git a/test/runtests.jl b/test/runtests.jl index f1791ac..fec177f 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,13 +1,14 @@ using SpectralKit -# using SpectralKit: PM1 using Test, DocStringExtensions, StaticArrays, BenchmarkTools, FiniteDifferences include("utilities.jl") include("test_utilities.jl") include("test_derivatives.jl") -# include("test_domains.jl") -# include("test_transformations.jl") -# include("test_generic_api.jl") -# include("test_chebyshev.jl") -# include("test_smolyak.jl") +include("test_domains.jl") +include("test_transformations.jl") +include("test_chebyshev.jl") +include("test_smolyak_traversal.jl") +include("test_smolyak.jl") + +include("test_generic_api.jl") diff --git a/test/test_chebyshev.jl b/test/test_chebyshev.jl index c3e44d4..2da50bb 100644 --- a/test/test_chebyshev.jl +++ b/test/test_chebyshev.jl @@ -27,7 +27,7 @@ θ = rand(N) @test linear_combination(basis, θ, x) ≈ sum(chebyshev_cos(x, i) * θ for (i,θ) in enumerate(θ)) - @test linear_combination(basis, θ, derivatives(x))[1] ≈ + @test linear_combination(basis, θ, 𝑑(x))[1] ≈ sum(chebyshev_cos_deriv(x, i) * θ for (i,θ) in enumerate(θ)) end diff --git a/test/test_derivatives.jl b/test/test_derivatives.jl index 84852aa..394438c 100644 --- a/test/test_derivatives.jl +++ b/test/test_derivatives.jl @@ -1,9 +1,8 @@ using SpectralKit, Test # test internals using SpectralKit: 𝑑Derivatives, _one, _add, _sub, _mul, - Partials, _is_strict_subset, _partials_minimal_representation, _partials_canonical_expansion - # _partials_max, -# _normalize_partials, _calculate_M_Ds_union, _collapse_D_union!, ∂Derivatives + Partials, _is_strict_subset, _partials_minimal_representation, + _partials_canonical_expansion, ∂Derivatives, ∂CoordinateExpansion using Random: randperm using StaticArrays: SVector @@ -89,6 +88,25 @@ end end end +@testset "partial derivatives API entry points" begin + @test ∂() ≡ ∂Derivatives{Tuple{}}() + @test ∂(2, 2) ≡ ∂Derivatives{Tuple{Partials(2, 2)}}() + @test 𝑑^2 << Val(2) ≡ ∂Derivatives{Tuple{Partials(0, 2)}}() + @test ∂(1, 2) ∪ ∂(2, 1) ≡ + ∂Derivatives{Tuple{_partials_minimal_representation([Partials(1, 2), + Partials(2, 1)])...}}() + @test repr(∂()) == "∂()" + @test repr(∂(1, 2)) == "∂(1, 2)" +end + +@testset "partials derivatives expansions" begin + D = ∂(2, 2) + x = SVector(1.0, 2.0) + Dx = @inferred(D(x)) + @test Dx isa ∂CoordinateExpansion{<:typeof(D)} +end + + # @testset "partial derivatives interface" begin # @test ∂(2, (1, 2), (2, 2)) == ∂Specification{(1, 2)}(((1, 1), (0, 2))) # @test ∂(3, (), (1, 1, 2), (3, 2)) == ∂Specification{(2, 1, 1)}(((0, 0, 0), (2, 1, 0), (0, 1, 1))) diff --git a/test/test_domains.jl b/test/test_domains.jl index bde3ca5..4bac7c2 100644 --- a/test/test_domains.jl +++ b/test/test_domains.jl @@ -1,5 +1,8 @@ +using SpectralKit, Test +using SpectralKit: PM1, UnivariateDomain + @testset "domain API" begin - d1 = SpectralKit.PM1() + d1 = PM1() @test repr(d1) == "[-1,1]" @test extrema(d1) == (-1, 1) @test (minimum(d1), maximum(d1)) == extrema(d1) @@ -11,6 +14,6 @@ @test length(dn) == 2 @test dn[1] == d1 @test Tuple(dn) ≡ (d1, d1) - d2 = coordinate_domains(SpectralKit.PM1(), SpectralKit.UnivariateDomain(-3.0, Inf)) + d2 = coordinate_domains(PM1(), SpectralKit.UnivariateDomain(-3.0, Inf)) @test repr(d2) == "[-1,1]×[-3.0,∞]" end diff --git a/test/test_smolyak.jl b/test/test_smolyak.jl index c006ad5..8ae0296 100644 --- a/test/test_smolyak.jl +++ b/test/test_smolyak.jl @@ -1,124 +1,5 @@ -using SpectralKit: nesting_total_length, nesting_block_length, SmolyakIndices, - __smolyak_length, SmolyakGridShuffle, SmolyakProduct - -"grids we test on" -GRIDS = (EndpointGrid(), InteriorGrid(), InteriorGrid2()) - -@testset "printing SmolyakParameters" begin - @test repr(SmolyakParameters(3, 2)) == "Smolyak parameters, ∑bᵢ ≤ 3, all bᵢ ≤ 2" -end - -#### -#### blocks -#### - -@testset "block length" begin - for grid_kind in GRIDS - for b in 0:5 - nA = nesting_total_length(Chebyshev, grid_kind, b) - gA = grid(Chebyshev(grid_kind, nA)) - gB = grid(Chebyshev(grid_kind, nesting_total_length(Chebyshev, grid_kind, b + 1))) - - @test is_approximate_subset(gA, gB) - @test sum(b -> nesting_block_length(Chebyshev, grid_kind, b), 0:b) == nA - end - end -end - -""" -Collect shuffled indices for the given `grid_kind` for block indices `0, …, b`, returned as -a `Vector{Vector{Int}}`. For testing. -""" -function shuffled_indices_upto_b(grid_kind, b) - _grid(b) = grid(Chebyshev(grid_kind, nesting_total_length(Chebyshev, grid_kind, b))) - g0 = _grid(b) - indices = Vector{Vector{Int}}() - for b in b:(-1):1 - in_b = is_approximately_in(g0, _grid(b)) - notin_bm1 = .!is_approximately_in(g0, _grid(b - 1)) - mask = in_b .& notin_bm1 - push!(indices, findall(mask)) - end - push!(indices, [(length(g0) + 1) ÷ 2]) - reverse(indices) -end - -@testset "block shuffle" begin - @testset "endpoint" begin - for grid_kind in GRIDS - for b in 0:6 - len = nesting_total_length(Chebyshev, grid_kind, b) - ι = SmolyakGridShuffle(grid_kind, len) - @test length(ι) == len - @test @inferred eltype(ι) == Int - @test collect(ι) == reduce(vcat, shuffled_indices_upto_b(grid_kind, b)) - end - end - end -end - -#### -#### traversal -#### - -""" -Naive implementation of Smolyan index iteration, traversing a `CartesianIndices` and keeping -valid indexes. For testing/comparison. Returns a vector of `indexes => blocks` pairs. -""" -function smolyak_indices_check(grid_kind, N, B, M) - m = nesting_total_length(Chebyshev, grid_kind, M) - b_table = fill(M, m) - for b in (M-1):(-1):0 - b_table[1:nesting_total_length(Chebyshev, grid_kind, b)] .= b - end - T = NTuple{N,Int} - result = Vector{Pair{T,T}}() - for ι in CartesianIndices(ntuple(_ -> 1:m, N)) - ix = Tuple(ι) - blocks = map(i -> b_table[i], ix) - if sum(blocks) ≤ B - push!(result, ix => blocks) - end - end - result -end - -@testset "Smolyak indices" begin - for grid_kind in GRIDS - for B in 0:3 - for M in 0:B - for N in 1:4 - ι = SmolyakIndices{N}(grid_kind, SmolyakParameters(B, M)) - x1 = @inferred collect(ι) - x2 = first.(smolyak_indices_check(grid_kind, N, B, M)) - len = @inferred __smolyak_length(grid_kind, Val(N), Val(B), M) - @test x1 == x2 - @test len == length(x1) == length(ι) == length(x2) - end - end - end - end -end - -@testset "Smolyak product primitives" begin - for grid_kind in GRIDS - for B in 0:3 - for M in 0:B - for N in 1:4 - ι = SmolyakIndices{N}(grid_kind, SmolyakParameters(B, M)) - ℓ = nesting_total_length(Chebyshev, grid_kind, min(B,M)) - sources = ntuple(_ -> rand(SVector{ℓ, Float64}), Val(N)) - P = SmolyakProduct(ι, sources, nothing) - @test length(ι) == length(P) - @test eltype(P) == Float64 - for (i, p) in zip(ι, P) - @test prod(getindex.(sources, i)) ≈ p - end - end - end - end - end -end +using SpectralKit, Test +using SpectralKit: PM1 #### #### api diff --git a/test/test_smolyak_traversal.jl b/test/test_smolyak_traversal.jl new file mode 100644 index 0000000..ef28c0e --- /dev/null +++ b/test/test_smolyak_traversal.jl @@ -0,0 +1,127 @@ +using SpectralKit: nesting_total_length, nesting_block_length, SmolyakIndices, + __smolyak_length, SmolyakGridShuffle, SmolyakProduct + +"grids we test on" +GRIDS = (EndpointGrid(), InteriorGrid(), InteriorGrid2()) + +### +### block parameters +### + + +@testset "printing SmolyakParameters" begin + @test repr(SmolyakParameters(3, 2)) == "Smolyak parameters, ∑bᵢ ≤ 3, all bᵢ ≤ 2" +end + + +#### +#### blocks +#### + +@testset "block length" begin + for grid_kind in GRIDS + for b in 0:5 + nA = nesting_total_length(Chebyshev, grid_kind, b) + gA = grid(Chebyshev(grid_kind, nA)) + gB = grid(Chebyshev(grid_kind, nesting_total_length(Chebyshev, grid_kind, b + 1))) + + @test is_approximate_subset(gA, gB) + @test sum(b -> nesting_block_length(Chebyshev, grid_kind, b), 0:b) == nA + end + end +end + +""" +Collect shuffled indices for the given `grid_kind` for block indices `0, …, b`, returned as +a `Vector{Vector{Int}}`. For testing. +""" +function shuffled_indices_upto_b(grid_kind, b) + _grid(b) = grid(Chebyshev(grid_kind, nesting_total_length(Chebyshev, grid_kind, b))) + g0 = _grid(b) + indices = Vector{Vector{Int}}() + for b in b:(-1):1 + in_b = is_approximately_in(g0, _grid(b)) + notin_bm1 = .!is_approximately_in(g0, _grid(b - 1)) + mask = in_b .& notin_bm1 + push!(indices, findall(mask)) + end + push!(indices, [(length(g0) + 1) ÷ 2]) + reverse(indices) +end + +@testset "block shuffle" begin + @testset "endpoint" begin + for grid_kind in GRIDS + for b in 0:6 + len = nesting_total_length(Chebyshev, grid_kind, b) + ι = SmolyakGridShuffle(grid_kind, len) + @test length(ι) == len + @test @inferred eltype(ι) == Int + @test collect(ι) == reduce(vcat, shuffled_indices_upto_b(grid_kind, b)) + end + end + end +end + +#### +#### traversal +#### + +""" +Naive implementation of Smolyan index iteration, traversing a `CartesianIndices` and keeping +valid indexes. For testing/comparison. Returns a vector of `indexes => blocks` pairs. +""" +function smolyak_indices_check(grid_kind, N, B, M) + m = nesting_total_length(Chebyshev, grid_kind, M) + b_table = fill(M, m) + for b in (M-1):(-1):0 + b_table[1:nesting_total_length(Chebyshev, grid_kind, b)] .= b + end + T = NTuple{N,Int} + result = Vector{Pair{T,T}}() + for ι in CartesianIndices(ntuple(_ -> 1:m, N)) + ix = Tuple(ι) + blocks = map(i -> b_table[i], ix) + if sum(blocks) ≤ B + push!(result, ix => blocks) + end + end + result +end + +@testset "Smolyak indices" begin + for grid_kind in GRIDS + for B in 0:3 + for M in 0:B + for N in 1:4 + ι = SmolyakIndices{N}(grid_kind, SmolyakParameters(B, M)) + x1 = @inferred collect(ι) + x2 = first.(smolyak_indices_check(grid_kind, N, B, M)) + len = @inferred __smolyak_length(grid_kind, Val(N), Val(B), M) + @test x1 == x2 + @test len == length(x1) == length(ι) == length(x2) + end + end + end + end +end + +@testset "Smolyak product primitives" begin + for grid_kind in GRIDS + for B in 0:3 + for M in 0:B + for N in 1:4 + ι = SmolyakIndices{N}(grid_kind, SmolyakParameters(B, M)) + ℓ = nesting_total_length(Chebyshev, grid_kind, min(B,M)) + sources = ntuple(_ -> rand(SVector{ℓ, Float64}), Val(N)) + P = SmolyakProduct(ι, sources, nothing) + @test length(ι) == length(P) + @test eltype(P) == Float64 + for (i, p) in zip(ι, P) + @test prod(getindex.(sources, i)) ≈ p + end + end + end + end + end +end From c755cfb0d794cec41a894d626134a1eae1a4fa98 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tam=C3=A1s=20K=2E=20Papp?= Date: Wed, 22 May 2024 15:40:11 +0200 Subject: [PATCH 04/14] fix partial derivatives calculations --- src/derivatives.jl | 10 +++++----- test/test_smolyak.jl | 14 +++++++++++++- 2 files changed, 18 insertions(+), 6 deletions(-) diff --git a/src/derivatives.jl b/src/derivatives.jl index f557515..6616932 100644 --- a/src/derivatives.jl +++ b/src/derivatives.jl @@ -328,22 +328,22 @@ end (∂D::∂Derivatives)(x::Tuple) = ∂CoordinateExpansion(∂D, _expand_coordinates(∂D, x)) -(∂D::∂Derivatives)(x::AbstractVector) = ∂CoordinateExpansion(∂D, Tuple(x)) +(∂D::∂Derivatives)(x::AbstractVector) = ∂D(Tuple(x)) struct ∂Expansion{D,N,T} ∂D::D coefficients::SVector{N,T} function ∂Expansion(∂D::D, coefficients::SVector{N,T}) where {D<:∂Derivatives,N,T} - new{D,N,T}(coefficients) + new{D,N,T}(∂D, coefficients) end end function _add(x::∂Expansion{D,N}, y::∂Expansion{D,N}) where {D,N} - ∂Expansion(x.D, map(+, x.values, y.values)) + ∂Expansion(x.∂D, map(+, x.coefficients, y.coefficients)) end function _mul(x::Real, y::∂Expansion) - ∂Expansion(y.D, map(y -> _mul(x, y), y.coefficients)) + ∂Expansion(y.∂D, map(y -> _mul(x, y), y.coefficients)) end # function Base.show(io::IO, expansion::∂Expansion{<:∂Derivatives{K,M,D}}) where {K,M,D} @@ -380,7 +380,7 @@ end end products = [_product(d) for d in _partials_canonical_expansion(Val(N), fieldtypes(Ps))] quote - ∂Expansion(∂derivatives, SVector($(products)...)) + ∂Expansion(∂D, SVector($(products...))) end end diff --git a/test/test_smolyak.jl b/test/test_smolyak.jl index 8ae0296..c1863a1 100644 --- a/test/test_smolyak.jl +++ b/test/test_smolyak.jl @@ -1,5 +1,5 @@ using SpectralKit, Test -using SpectralKit: PM1 +using SpectralKit: PM1, ∂Expansion #### #### api @@ -28,6 +28,18 @@ end @test linear_combination(basis, θ, y) ≈ f(y) end end + + @testset "sanity check for derivatives" begin + # NOTE this just checks that it runs and is inferred, but does not check + # correctness, derivatives derived below should be compared + # x[1] * x[2] + 5 * x[1] - 3 * x[2] + 5 + # f1(x) = x[2] + 5 + # f2(x) = x[1] - 3 + # f12(x) = 1 + D = ∂(1, 1) + y = SVector(1.0, 2.0) + @test @inferred(linear_combination(basis, θ, D(y))) isa ∂Expansion + end end @testset "Smolyak API allocations" begin From 208630d417621c9375622a68eb7e48bfabdb2693 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tam=C3=A1s=20K=2E=20Papp?= Date: Wed, 22 May 2024 16:21:06 +0200 Subject: [PATCH 05/14] documentation, first pass --- src/derivatives.jl | 105 +++++++++++++++++++++++++++++++++++++++------ src/utilities.jl | 28 ++++++------ 2 files changed, 108 insertions(+), 25 deletions(-) diff --git a/src/derivatives.jl b/src/derivatives.jl index 6616932..3260aa3 100644 --- a/src/derivatives.jl +++ b/src/derivatives.jl @@ -1,12 +1,13 @@ ##### -##### Internal implementation for derivatives. See docs of [`Derivatives`](@ref). +##### Internal implementation for derivatives. ##### export 𝑑, ∂ #### -#### Operations we support, with fallbacks defined below. This is deliberately kept -#### minimal, now all versions are defined for commutative operators. +#### Operations we support for calculating (untransformed) basis functions and their +#### products, with fallbacks defined below. This is deliberately kept minimal, not all +#### versions are defined for commutative operators. #### _one(::Type{T}) where {T<:Real} = one(T) @@ -23,18 +24,20 @@ _mul(x, y, z) = _mul(_mul(x, y), z) #### Univariate expansions and derivatives #### - -# """ -# A small AD framework used *internally*, for calculating derivatives. - -# Supports only the operations required by this module. The restriction is deliberate, -# should not be used for arithmetic operators outside this package. - -# See [`∂`](@ref) for the exposed API. -# """ struct 𝑑Expansion{Dp1,T} "The function value and derivatives. `Dp1` is the degree of the last derivative + 1." coefficients::SVector{Dp1,T} + @doc """ + $(SIGNATURES) + + Taylor expansion around a given value. + + Implements small AD framework used *internally*, for calculating derivatives, that + supports only the operations required by this module. The restriction is deliberate, + should not be used for arithmetic operators outside this package. + + See [`∂`](@ref) for the exposed API. This type supports `eltype` and `getindex`. + """ function 𝑑Expansion(coefficients::SVector{Dp1,T}) where {Dp1,T} new{Dp1,T}(coefficients) end @@ -130,6 +133,7 @@ end "Types accepted as scalars in this package." const Scalar = Union{Real,𝑑Expansion} +# FIXME incorporate into docs # """ # derivatives(x, ::Val(N) = Val(1)) @@ -180,6 +184,23 @@ const Scalar = Union{Real,𝑑Expansion} struct Partials{N} I::NTuple{N,Int} + """ + $(SIGNATURES) + + Partial derivatives up to given indices. Eg `Partials((1, 2))` would contain + - the value, + - first derivatives ``∂_1`, ``∂_2``, + - cross derivatives ``∂_1 ∂_2``, + - second derivative ``∂^2_2`` + + This is just a building block used internally by `∂Derivatives`. The actual order in + containers is determined by [`_partials_canonical_expansion`](@ref). See also + [`_partials_minimal_representation`](@ref). + + This constructor enforces that the last index is non-zero. Use the other + `Partials(I...)` constructor to strip trailing zeros. Note that the number indices + just determines the *minimum* length for the multivariate arguments to be differentiated. + """ function Partials(I::NTuple{N,Int}) where N @argcheck all(i -> i ≥ 0, I) @argcheck I ≡ () || last(I) ≠ 0 @@ -196,6 +217,12 @@ function Partials(I::Integer...) Partials(ntuple(i -> Int(I[i]), N)) end +""" +$(SIGNATURES) + +True iff the partial derivatives contained in the first argument is a *strict* subset of +those in the second. +""" function _is_strict_subset(p1::Partials{N1}, p2::Partials{N2}) where {N1,N2} N1 > N2 && return false # valid because derivatives are always positive I1 = p1.I @@ -208,6 +235,13 @@ function _is_strict_subset(p1::Partials{N1}, p2::Partials{N2}) where {N1,N2} strict || N1 < N2 end +""" +$(SIGNATURES) + +Imposes a total ordering on `Partials` with the property that a strict subset implies an +order, but not necessarily vice versa. This determines the traversal and help with +eliminating nested specifications. +""" function Base.isless(p1::Partials{N1}, p2::Partials{N2}) where {N1, N2} _is_strict_subset(p1, p2) && return true _is_strict_subset(p2, p1) && return false @@ -216,6 +250,17 @@ end Base.isequal(p1::Partials, p2::Partials) = p1.I == p2.I +""" +$(SIGNATURES) + +Collapse specification of partial derivatives (an iterable of `Partials`) so a canonical +“minimal” representation with respect to the total ordering. Eliminates nested +specifications and duplicates. + +!!! NOTE + Takes any iterable of `Partials`, returns a `Vector`. Allocates, for use in + generated functions. +""" function _partials_minimal_representation(partials) descending_partials = sort!(collect(Partials, partials); rev = true) minimal_partials = Vector{Partials}() @@ -227,6 +272,13 @@ function _partials_minimal_representation(partials) minimal_partials end +""" +$(SIGNATURES) + +The ordering of partial derivatives contains for containers, for `N`-dimensional +arguments. Returns an iterable of `N`-tuples that contain integer indices of partial +derivatives, eg `(0, 1, 2)` for ``∂_2 ∂_3^2``. +""" function _partials_canonical_expansion(::Val{N}, Ps) where N result = OrderedSet{NTuple{N,Int}}() function _plus1_cartesian_indices(p::Partials{M}) where M @@ -245,6 +297,12 @@ function _partials_canonical_expansion(::Val{N}, Ps) where N result end +""" +$(SIGNATURES) + +Elementwise maximum of the iterable from `_partials_canonical_expansion`, an `N`-tuple +of nonnegative integers. +""" function _partials_expansion_degrees(::Val{N}, partials) where N degrees = zero(MVector{N,Int}) for P in partials @@ -255,8 +313,18 @@ function _partials_expansion_degrees(::Val{N}, partials) where N Tuple(degrees) end +""" +$(SIGNATURES) + +The smallest input dimension (length) a partial derivative specification can support. +""" _partials_minimum_input_dimension(partials) = maximum(P -> length(P.I), partials) +""" +$(SIGNATURES) + +Check if `Ps` is a minimal representation, ie a valid type parameter for `∂Derivatives`. +""" function _is_minimal_representation(::Val{Ps}) where Ps _Ps = fieldtypes(Ps) _Ps isa Tuple{Vararg{Partials}} || return false @@ -264,10 +332,23 @@ function _is_minimal_representation(::Val{Ps}) where Ps end struct ∂Derivatives{Ps} + @doc """ + $(SIGNATURES) + + A callable that requests that the given partial derivatives of its argument are + evaluated. + + The partial derivatives are encoded in the `Ps` as a `Tuple{...}` of `Partials`. + They are checked to be “minimal”, see [`_is_minimal_representation`](@ref), except + when `Ps` is a `Partials`, then it is wrapped and used as is. + + The API entry point is `∂`s, combined with `<<` and `∪`/`union`. + """ function ∂Derivatives{Ps}() where {Ps} if Ps isa Partials new{Tuple{Ps}}() else + # FIXME check that this does not allocate @argcheck _is_minimal_representation(Val(Ps)) new{Ps}() end diff --git a/src/utilities.jl b/src/utilities.jl index ef76c98..80d3dbc 100644 --- a/src/utilities.jl +++ b/src/utilities.jl @@ -36,20 +36,22 @@ end Base.print(io::IO, s::SubScript) = print_number(io, _SUBSCRIPT_DIGITS, s.i) -""" -$(SIGNATURES) +# FIXME function below is currently unused, decide if we need to keep it after +# refactoring Base.show(::IO, ::∂Expansion) +# """ +# $(SIGNATURES) -Print notation for partial derivatives, where `d[i]` stands for ``∂ⁱ/∂x[d]ⁱ``. -""" -function _print_partial_notation(io::IO, d) - if d ≡ () - print(io, "value") - else - for (i, p) in enumerate(d) - print(io, "∂", SubScript(i), SuperScript(p)) - end - end -end +# Print notation for partial derivatives, where `d[i]` stands for ``∂ⁱ/∂x[d]ⁱ``. +# """ +# function _print_partial_notation(io::IO, d) +# if d ≡ () +# print(io, "value") +# else +# for (i, p) in enumerate(d) +# print(io, "∂", SubScript(i), SuperScript(p)) +# end +# end +# end """ $(SIGNATURES) From cee43e402cc9e458fa400b1f5f3f8e949495bf65 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tam=C3=A1s=20K=2E=20Papp?= Date: Wed, 22 May 2024 16:44:53 +0200 Subject: [PATCH 06/14] add a shortcut --- src/generic_api.jl | 8 ++++++++ src/transformations.jl | 4 ++++ test/runtests.jl | 3 +-- test/test_generic_api.jl | 18 ++++++++++++++++++ 4 files changed, 31 insertions(+), 2 deletions(-) diff --git a/src/generic_api.jl b/src/generic_api.jl index 0eecf2d..e238256 100644 --- a/src/generic_api.jl +++ b/src/generic_api.jl @@ -297,3 +297,11 @@ function Base.getindex(basis::TransformedBasis{<:MultivariateBasis}, i::Int) end # FIXME add augmentation for transformed bases + +function transform_to(basis::FunctionBasis, transformation, x) + transform_to(domain(basis), transformation, x) +end + +function transform_from(basis::FunctionBasis, transformation, x) + transform_from(domain(basis), transformation, x) +end diff --git a/src/transformations.jl b/src/transformations.jl index 961e73d..4451c17 100644 --- a/src/transformations.jl +++ b/src/transformations.jl @@ -31,6 +31,8 @@ domain_kind(::Type{<:AbstractUnivariateTransformation}) = :univariate Transform `x` to `domain` using `transformation`. +`domain` can be replaced by `basis` for a shortcut which uses `domain(basis)`. + !!! FIXME document, especially differentiability requirements at infinite endpoints """ @@ -41,6 +43,8 @@ function transform_to end Transform `x` from `domain` using `transformation`. +`domain` can be replaced by `basis` for a shortcut which uses `domain(basis)`. + !!! FIXME document, especially differentiability requirements at infinite endpoints """ diff --git a/test/runtests.jl b/test/runtests.jl index fec177f..b414af3 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -10,5 +10,4 @@ include("test_transformations.jl") include("test_chebyshev.jl") include("test_smolyak_traversal.jl") include("test_smolyak.jl") - -include("test_generic_api.jl") +include("test_generic_api.jl") # NOTE moved last as it used constructs from above diff --git a/test/test_generic_api.jl b/test/test_generic_api.jl index d3a7ccb..378668b 100644 --- a/test/test_generic_api.jl +++ b/test/test_generic_api.jl @@ -48,6 +48,15 @@ end end end +@testset "transform_to and transform_from univariate shortcuts" begin + basis = Chebyshev(EndpointGrid(), 5) + t = BoundedLinear(1.0, 2.0) + y = rand_pm1() + x = transform_to(domain(basis), t, y) + @test transform_to(basis, t, y) == x + @test transform_from(basis, t, x) == transform_from(domain(basis), t, x) +end + @testset "transformed bases and linear combinations (bivariate)" begin basis0 = smolyak_basis(Chebyshev, InteriorGrid(), SmolyakParameters(2, 2), Val(2)) t = coordinate_transformations(BoundedLinear(1.0, 2.0), SemiInfRational(0, 1)) @@ -70,6 +79,15 @@ end end end +@testset "transform_to and transform_from bivariate shortcuts" begin + basis = smolyak_basis(Chebyshev, InteriorGrid(), SmolyakParameters(2, 2), Val(2)) + t = coordinate_transformations(BoundedLinear(1.0, 2.0), SemiInfRational(0, 1)) + y = SVector(rand_pm1(), rand_pm1()) + x = transform_to(domain(basis), t, y) + @test transform_to(basis, t, y) == x + @test transform_from(basis, t, x) == transform_from(domain(basis), t, x) +end + @testset "subset fallback" begin @test !is_subset_basis(Chebyshev(InteriorGrid(), 4), # just test the fallback method smolyak_basis(Chebyshev, InteriorGrid(), SmolyakParameters(2, 2), 2)) From 507b06ae0eb314b226d631a4eb15f3e951335880 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tam=C3=A1s=20K=2E=20Papp?= Date: Wed, 22 May 2024 16:44:59 +0200 Subject: [PATCH 07/14] add compat for new dependency --- Project.toml | 1 + 1 file changed, 1 insertion(+) diff --git a/Project.toml b/Project.toml index 4840493..966168f 100644 --- a/Project.toml +++ b/Project.toml @@ -12,6 +12,7 @@ StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" [compat] ArgCheck = "1, 2" DocStringExtensions = "0.8, 0.9" +OrderedCollections = "1" StaticArrays = "1" julia = "1.9" From a88c2a6f8cb474bb417db37b1216b2b346b8ce13 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tam=C3=A1s=20K=2E=20Papp?= Date: Wed, 22 May 2024 16:52:40 +0200 Subject: [PATCH 08/14] add JET and Aqua testing, remove cruft from Project.toml --- Project.toml | 9 --------- README.md | 1 + test/Project.toml | 2 ++ test/runtests.jl | 12 ++++++++++++ 4 files changed, 15 insertions(+), 9 deletions(-) diff --git a/Project.toml b/Project.toml index 966168f..d4c33a2 100644 --- a/Project.toml +++ b/Project.toml @@ -15,12 +15,3 @@ DocStringExtensions = "0.8, 0.9" OrderedCollections = "1" StaticArrays = "1" julia = "1.9" - -[extras] -BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf" -ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" -Sobol = "ed01d8cd-4d21-5b2a-85b4-cc3bdc58bad4" -Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" - -[targets] -test = ["BenchmarkTools", "Test", "ForwardDiff", "Sobol"] diff --git a/README.md b/README.md index d01c241..e7656a3 100644 --- a/README.md +++ b/README.md @@ -6,6 +6,7 @@ [![Documentation](https://img.shields.io/badge/docs-stable-blue.svg)](https://tpapp.github.io/SpectralKit.jl/stable) [![Documentation](https://img.shields.io/badge/docs-master-blue.svg)](https://tpapp.github.io/SpectralKit.jl/dev) [![DOI](https://zenodo.org/badge/220448027.svg)](https://zenodo.org/badge/latestdoi/220448027) +[![Aqua QA](https://raw.githubusercontent.com/JuliaTesting/Aqua.jl/master/badge.svg)](https://github.com/JuliaTesting/Aqua.jl) Building blocks of spectral methods for Julia. Currently includes Chebyshev polynomials on univariate and Smolyak (multivariate) grids, with domain transformations to semi-infinite and infinite domains. diff --git a/test/Project.toml b/test/Project.toml index 2bd97b0..ddfa4a9 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -1,7 +1,9 @@ [deps] +Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595" BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf" DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" FiniteDifferences = "26cc04aa-876d-5657-8c51-4c34ba976000" +JET = "c3a54625-cd67-489e-a8e7-0a5a0ff4e31b" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" Sobol = "ed01d8cd-4d21-5b2a-85b4-cc3bdc58bad4" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" diff --git a/test/runtests.jl b/test/runtests.jl index b414af3..315e867 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -11,3 +11,15 @@ include("test_chebyshev.jl") include("test_smolyak_traversal.jl") include("test_smolyak.jl") include("test_generic_api.jl") # NOTE moved last as it used constructs from above + +using JET +@testset "static analysis with JET.jl" begin + @test isempty(JET.get_reports(report_package(SpectralKit, target_modules=(SpectralKit,)))) +end + +@testset "QA with Aqua" begin + import Aqua + Aqua.test_all(SpectralKit; ambiguities = false) + # testing separately, cf https://github.com/JuliaTesting/Aqua.jl/issues/77 + Aqua.test_ambiguities(SpectralKit) +end From 7d25c4866c426a6f036ba5ca9bdb3d607e811452 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tam=C3=A1s=20K=2E=20Papp?= Date: Thu, 23 May 2024 07:29:30 +0200 Subject: [PATCH 09/14] make linear_combination work with SVector of coefficients --- src/derivatives.jl | 4 ++-- src/generic_api.jl | 4 ++++ test/test_generic_api.jl | 9 +++++++++ 3 files changed, 15 insertions(+), 2 deletions(-) diff --git a/src/derivatives.jl b/src/derivatives.jl index 3260aa3..fd43851 100644 --- a/src/derivatives.jl +++ b/src/derivatives.jl @@ -12,11 +12,11 @@ export 𝑑, ∂ _one(::Type{T}) where {T<:Real} = one(T) -_add(x::Real, y::Real) = x + y +_add(x::T, y::T) where {T <: Union{Real,SVector}} = x + y _sub(x::Real, y::Real) = x - y -_mul(x::Real, y::Real) = x * y +_mul(x::Union{Real,SVector}, y::Real) = x * y _mul(x, y, z) = _mul(_mul(x, y), z) diff --git a/src/generic_api.jl b/src/generic_api.jl index e238256..396ad8e 100644 --- a/src/generic_api.jl +++ b/src/generic_api.jl @@ -106,6 +106,10 @@ $(SIGNATURES) Helper function for linear combinations of basis elements at `x`. Always checks that `θ` and `basis` have compatible dimensions. + +!!! NOTE + `x` and `θ` can be anything that supports `_mul(θ[i], b[i])` and `_add` on the + result of this. """ @inline function _linear_combination(basis, θ, x) # an implementation of mapreduce, to work around diff --git a/test/test_generic_api.jl b/test/test_generic_api.jl index 378668b..7f8909c 100644 --- a/test/test_generic_api.jl +++ b/test/test_generic_api.jl @@ -88,6 +88,15 @@ end @test transform_from(basis, t, x) == transform_from(domain(basis), t, x) end +@testset "linear combination SVector passthrough" begin + N = 10 + M = 3 + basis = Chebyshev(InteriorGrid(), N) ∘ SemiInfRational(0.0, 1.0) + θ = rand(SVector{M,Float64}, N) + x = 2.0 + @test linear_combination(basis, θ, x) ≡ SVector{M}(linear_combination(basis, map(x -> x[i], θ), x) for i in 1:M) +end + @testset "subset fallback" begin @test !is_subset_basis(Chebyshev(InteriorGrid(), 4), # just test the fallback method smolyak_basis(Chebyshev, InteriorGrid(), SmolyakParameters(2, 2), 2)) From eb6f2f52b6e30ca5a5807ed9be2ef6bddcf4f047 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tam=C3=A1s=20K=2E=20Papp?= Date: Mon, 17 Jun 2024 17:31:39 +0200 Subject: [PATCH 10/14] add simple test for Chebyshev derivatives --- src/transformations.jl | 20 ++++++++++---------- test/test_chebyshev.jl | 17 +++++++++++++++++ test/test_derivatives.jl | 5 +++++ 3 files changed, 32 insertions(+), 10 deletions(-) diff --git a/src/transformations.jl b/src/transformations.jl index 4451c17..715a41e 100644 --- a/src/transformations.jl +++ b/src/transformations.jl @@ -187,11 +187,11 @@ end function transform_to(domain::PM1, t::BoundedLinear, y::𝑑Expansion{Dp1}) where Dp1 (; m, s) = t - (; derivatives) = y - y0, yD... = derivatives + (; coefficients) = y + y0, yD... = coefficients x0 = transform_to(domain, t, y0) xD = map(y -> y / s, yD) - 𝑑Expansion((x0, xD...)) + 𝑑Expansion(SVector(x0, xD...)) end function domain(t::BoundedLinear) @@ -256,12 +256,12 @@ end function transform_to(domain::PM1, t::SemiInfRational, y::𝑑Expansion{Dp1}) where Dp1 (; A, L) = t - (; derivatives) = y - x0 = transform_to(domain, t, derivatives[1]) + (; coefficients) = y + x0 = transform_to(domain, t, coefficients[1]) Dp1 == 1 && return 𝑑Expansion(SVector(x0)) # based on Boyd (2001), Table E.7 Q = abs2(x0 - 1) - x1 = (derivatives[2] * Q) / (2*L) + x1 = (coefficients[2] * Q) / (2*L) Dp1 == 2 && return 𝑑Expansion(SVector(x0, x1)) error("$(Dp1-1)th derivative not implemented yet, open an issue.") end @@ -321,13 +321,13 @@ end function transform_to(domain::PM1, t::InfRational, y::𝑑Expansion{Dp1}) where Dp1 (; A, L) = t - (; derivatives) = y - x0 = transform_to(domain, t, derivatives[1]) - N == 1 && return Derivatives(SVector(x0)) + (; coefficients) = y + x0 = transform_to(domain, t, coefficients[1]) + N == 1 && return Coefficients(SVector(x0)) # based on Boyd (2001), Table E.5 Q = 1 - abs2(x0) sQ = √Q - x1 = (derivatives[2] * Q * sQ) / L + x1 = (coefficients[2] * Q * sQ) / L N == 2 && return 𝑑Expansion(SVector(x0, x1)) error("$(N-1)th derivative not implemented yet, open an issue.") end diff --git a/test/test_chebyshev.jl b/test/test_chebyshev.jl index 2da50bb..d2756b4 100644 --- a/test/test_chebyshev.jl +++ b/test/test_chebyshev.jl @@ -72,3 +72,20 @@ @test_throws ArgumentError augment_coefficients(basis, basis, randn(4)) end end + +@testset "univariate derivatives" begin + basis = Chebyshev(InteriorGrid(), 5) + N = 5 + D = 𝑑^Val(N) + for transformation in (BoundedLinear(-2, 3), ) + transformed_basis = basis ∘ transformation + f = linear_combination(transformed_basis, randn(dimension(transformed_basis))) + for _ in 1:100 + x = transform_from(basis, transformation, rand_pm1()) + y = f(D(x)) + for i in 0:N + @test y[i] ≈ DD(f, x, i) atol = 1e-6 + end + end + end +end diff --git a/test/test_derivatives.jl b/test/test_derivatives.jl index 394438c..f626fb4 100644 --- a/test/test_derivatives.jl +++ b/test/test_derivatives.jl @@ -1,3 +1,8 @@ +##### +##### NOTE: we only test building blocks here, derivative calculations exposed in the +##### API are tested in test_chebyshev.jl and test_smolyak.jl +##### + using SpectralKit, Test # test internals using SpectralKit: 𝑑Derivatives, _one, _add, _sub, _mul, From 2d3c0ac44e9b253d8ad29e478a0bdc587bd2ae8a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tam=C3=A1s=20K=2E=20Papp?= Date: Mon, 17 Jun 2024 17:35:28 +0200 Subject: [PATCH 11/14] add tests for other univariate transformations --- src/transformations.jl | 6 +++--- test/test_chebyshev.jl | 9 +++++---- 2 files changed, 8 insertions(+), 7 deletions(-) diff --git a/src/transformations.jl b/src/transformations.jl index 715a41e..f7a148f 100644 --- a/src/transformations.jl +++ b/src/transformations.jl @@ -323,13 +323,13 @@ function transform_to(domain::PM1, t::InfRational, y::𝑑Expansion{Dp1}) where (; A, L) = t (; coefficients) = y x0 = transform_to(domain, t, coefficients[1]) - N == 1 && return Coefficients(SVector(x0)) + Dp1 == 1 && return Coefficients(SVector(x0)) # based on Boyd (2001), Table E.5 Q = 1 - abs2(x0) sQ = √Q x1 = (coefficients[2] * Q * sQ) / L - N == 2 && return 𝑑Expansion(SVector(x0, x1)) - error("$(N-1)th derivative not implemented yet, open an issue.") + Dp1 == 2 && return 𝑑Expansion(SVector(x0, x1)) + error("$(Dp1-1)th derivative not implemented yet, open an issue.") end domain(::InfRational) = UnivariateDomain(-Inf, Inf) diff --git a/test/test_chebyshev.jl b/test/test_chebyshev.jl index d2756b4..0730e88 100644 --- a/test/test_chebyshev.jl +++ b/test/test_chebyshev.jl @@ -75,12 +75,13 @@ end @testset "univariate derivatives" begin basis = Chebyshev(InteriorGrid(), 5) - N = 5 - D = 𝑑^Val(N) - for transformation in (BoundedLinear(-2, 3), ) + for (transformation, N) in ((BoundedLinear(-2, 3), 5), + (SemiInfRational(0.7, 0.3), 1), + (InfRational(0.4, 0.9), 1)) + D = 𝑑^Val(N) transformed_basis = basis ∘ transformation f = linear_combination(transformed_basis, randn(dimension(transformed_basis))) - for _ in 1:100 + for _ in 1:50 x = transform_from(basis, transformation, rand_pm1()) y = f(D(x)) for i in 0:N From 4cf67a8471cdc6744b39e49722c528033796cfd8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tam=C3=A1s=20K=2E=20Papp?= Date: Mon, 17 Jun 2024 17:59:11 +0200 Subject: [PATCH 12/14] clean up tests a bit --- test/test_chebyshev.jl | 10 ++++---- test/test_generic_api.jl | 4 +-- test/test_transformations.jl | 48 ++++++++++++++++++++++-------------- test/utilities.jl | 28 ++++++++++++++++----- 4 files changed, 59 insertions(+), 31 deletions(-) diff --git a/test/test_chebyshev.jl b/test/test_chebyshev.jl index 0730e88..7d99aeb 100644 --- a/test/test_chebyshev.jl +++ b/test/test_chebyshev.jl @@ -16,8 +16,8 @@ @test dimension(basis) == N # check linear combinations - for i in 1:100 - x = rand_pm1(i) + for _ in 1:100 + x = rand_in_domain(basis) bx = @inferred basis_at(basis, x) @test length(bx) == N @@ -44,8 +44,8 @@ end # augmented coefficients - for i in 1:100 - x = rand_pm1(i) + for _ in 1:100 + x = rand_in_domain(basis) θ = rand(N) destination_basis = Chebyshev(grid_kind, N + 5) destination_θ = augment_coefficients(basis, destination_basis, θ) @@ -82,7 +82,7 @@ end transformed_basis = basis ∘ transformation f = linear_combination(transformed_basis, randn(dimension(transformed_basis))) for _ in 1:50 - x = transform_from(basis, transformation, rand_pm1()) + x = transform_from(basis, transformation, rand_in_domain(basis)) y = f(D(x)) for i in 0:N @test y[i] ≈ DD(f, x, i) atol = 1e-6 diff --git a/test/test_generic_api.jl b/test/test_generic_api.jl index 7f8909c..b7f6ed3 100644 --- a/test/test_generic_api.jl +++ b/test/test_generic_api.jl @@ -51,7 +51,7 @@ end @testset "transform_to and transform_from univariate shortcuts" begin basis = Chebyshev(EndpointGrid(), 5) t = BoundedLinear(1.0, 2.0) - y = rand_pm1() + y = rand_in_domain(basis) x = transform_to(domain(basis), t, y) @test transform_to(basis, t, y) == x @test transform_from(basis, t, x) == transform_from(domain(basis), t, x) @@ -82,7 +82,7 @@ end @testset "transform_to and transform_from bivariate shortcuts" begin basis = smolyak_basis(Chebyshev, InteriorGrid(), SmolyakParameters(2, 2), Val(2)) t = coordinate_transformations(BoundedLinear(1.0, 2.0), SemiInfRational(0, 1)) - y = SVector(rand_pm1(), rand_pm1()) + y = rand_in_domain(basis) x = transform_to(domain(basis), t, y) @test transform_to(basis, t, y) == x @test transform_from(basis, t, x) == transform_from(domain(basis), t, x) diff --git a/test/test_transformations.jl b/test/test_transformations.jl index 9e44e64..898318d 100644 --- a/test/test_transformations.jl +++ b/test/test_transformations.jl @@ -1,4 +1,4 @@ -@testset "BoundedLinear" begin +@testset "bounded linear domain transformations" begin @test_throws DomainError BoundedLinear(-1.0, Inf) @test_throws DomainError BoundedLinear(-1.0, -2.0) @@ -7,17 +7,21 @@ @test extrema(domain(trans)) == (A, B) - for i in 1:100 - x = rand_pm1(i) + for _ in 1:100 + x = rand_pm1() y = transform_from(PM1(), trans, x) - i == 1 && @test y ≈ A - i == 2 && @test y ≈ B - i > 2 && @test A < y < B + if x == -1 + @test y ≈ A + elseif x == 1 + @test y ≈ B + else + @test A < y < B + end @test transform_to(PM1(), trans, y) ≈ x end end -@testset "Chebyshev semi-infinite" begin +@testset "semi-infinite domain transformations" begin @test_throws DomainError SemiInfRational(-1.0, Inf) @test_throws DomainError SemiInfRational(-1.0, 0.0) @test_throws DomainError SemiInfRational(NaN, 2.0) @@ -28,17 +32,21 @@ end @test extrema(domain(trans)) == (A, Inf) - for i in 1:100 - x = rand_pm1(i) + for _ in 1:100 + x = rand_pm1() y = transform_from(PM1(), trans, x) - i == 1 && @test y ≈ A - i == 2 && @test y ≈ Inf - i > 2 && @test A < y < Inf + if x == -1 + @test y ≈ A + elseif x == 1 + @test y ≈ Inf + else + @test A < y < Inf + end @test transform_to(PM1(), trans, y) ≈ x end end -@testset "Chebyshev infinite" begin +@testset "infinite domain transformations" begin @test_throws DomainError InfRational(1.0, Inf) @test_throws DomainError InfRational(1.0, 0.0) @test_throws DomainError InfRational(1.0, -2.0) @@ -50,12 +58,16 @@ end @test extrema(domain(trans)) == (-Inf, Inf) - for i in 1:100 - x = rand_pm1(i) + for _ in 1:100 + x = rand_pm1() y = transform_from(PM1(), trans, x) - i == 1 && @test y == -Inf - i == 2 && @test y == Inf - i > 2 && @test isfinite(y) + if x == -1 + @test y == -Inf + elseif x == 1 + @test y == Inf + else + @test isfinite(y) + end @test transform_to(PM1(), trans, y) ≈ x end end diff --git a/test/utilities.jl b/test/utilities.jl index eec3601..c75778c 100644 --- a/test/utilities.jl +++ b/test/utilities.jl @@ -2,6 +2,8 @@ ##### utility functions for tests ##### +using SpectralKit: TransformedBasis, SmolyakBasis, SmolyakIndices # dispatch for rand_in_domain + struct KroneckerVector{T} <: AbstractVector{T} i::Int len::Int @@ -31,15 +33,29 @@ end """ $(SIGNATURES) -Return a random scalar in `(-1.0, 1.0)` (the default), except when `i == 1`, return `-1.0`, -when `i == 2`, return `1.0`. +Return a random value in [-1,1], putting an atomic mass on endpoints. The intention is to provide comprehensive testing for endpoints. """ -function rand_pm1(i = 3) - i == 1 && return -1.0 - i == 2 && return 1.0 - rand() * 2 - 1.0 +rand_pm1() = clamp((rand() - 0.5) * 2.5, -1, 1) + +""" +$(SIGNATURES) + +Return a random value in the domain of the given basis, putting an atomic mass on endpoints. + +The intention is to provide comprehensive testing for endpoints. +""" +rand_in_domain(::Chebyshev) = rand_pm1() + +function rand_in_domain(basis::SmolyakBasis{<:SmolyakIndices{N}}) where N + (; univariate_parent) = basis + SVector(ntuple(_ -> rand_in_domain(univariate_parent), Val(N))) +end + +function rand_in_domain(basis::TransformedBasis) + (; parent, transformation) + transform_from(parent, transformation, rand_in_domain(parent)) end "Flags (`true`) for elements in `a` that are within `atol` of some element in `b`." From 0a6d2884731addce6987bfc5f2edd7efdfbcf18a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tam=C3=A1s=20K=2E=20Papp?= Date: Tue, 18 Jun 2024 10:52:38 +0200 Subject: [PATCH 13/14] re-enable endpoint continuity for derivatives --- src/transformations.jl | 2 +- test/test_chebyshev.jl | 39 +++++++++++++++++++++++++++++++++++ test/test_composite.jl | 40 ------------------------------------ test/test_transformations.jl | 20 ++++++++++++++++++ test/utilities.jl | 12 ----------- 5 files changed, 60 insertions(+), 53 deletions(-) diff --git a/src/transformations.jl b/src/transformations.jl index f7a148f..b11fceb 100644 --- a/src/transformations.jl +++ b/src/transformations.jl @@ -247,7 +247,7 @@ function transform_to(::PM1, t::SemiInfRational, y::Real) (; A, L) = t z = y - A x = (z - L) / (z + L) - if (y == Inf && L > 0) || (y == -Inf && L < 0) + if y == Inf || y == -Inf one(x) else x diff --git a/test/test_chebyshev.jl b/test/test_chebyshev.jl index 7d99aeb..d9d5593 100644 --- a/test/test_chebyshev.jl +++ b/test/test_chebyshev.jl @@ -90,3 +90,42 @@ end end end end + +@testset "endpoint continuity for derivatives" begin + N = 10 + basis = Chebyshev(InteriorGrid(), N) + + # NOTE here we are checking that in some sense, derivatives give the right limit at + # endpoints for transformations to ∞. We use the analytical derivatives for + # comparison, based on the chain rule. + x_pinf = 𝑑(Inf) + x_minf = 𝑑(-Inf) + + @testset "SemiInfRational endpoints continuity" begin + trans = SemiInfRational(2.3, 0.7) + + for i in 1:N + θ = e_i(basis ∘ trans, i) + y_pinf = @inferred linear_combination(basis ∘ trans, θ, x_pinf) + @test y_pinf[0] == 1 + @test y_pinf[1] == 0 + y_minf = @inferred linear_combination(basis ∘ trans, θ, x_minf) + @test y_minf[0] == 1 + @test y_minf[1] == 0 + end + end + + @testset "InfRational endpoints continuity" begin + trans = InfRational(0.3, 3.0) + + for i in 1:N + θ = e_i(basis ∘ trans, i) + y_pinf = @inferred linear_combination(basis ∘ trans, θ, x_pinf) + @test y_pinf[0] == 1 + @test y_pinf[1] == 0 + y_minf = @inferred linear_combination(basis ∘ trans, θ, x_minf) + @test y_minf[0] == (-1)^(i+1) + @test y_minf[1] == 0 + end + end +end diff --git a/test/test_composite.jl b/test/test_composite.jl index 8680032..37768b7 100644 --- a/test/test_composite.jl +++ b/test/test_composite.jl @@ -44,43 +44,3 @@ end @test ForwardDiff.gradient(F, x) ≈ ForwardDiff.gradient(x -> f(from_pm1(ct, x)), x) end end - -# @testset "endpoint continuity for derivatives" begin -# N = 10 -# basis = Chebyshev(InteriorGrid(), N) - -# # NOTE here we are checking that in some sense, derivatives give the right limit at -# # endpoints for transformations to ∞. We use the derivative of the inverse rule (h′ -# # below) for comparison. - -# @testset "SemiInfRational endpoints continuity" begin -# trans = SemiInfRational(1.0, 1.0) -# x = 1.0 -# y = from_pm1(trans, x) -# h′ = ForwardDiff.derivative(x -> from_pm1(trans, x), x) - -# for i in 1:N -# θ = e_i(basis, i) -# f = y -> linear_combination(basis, θ, to_pm1(trans, y)) -# fy, f′y = f_f′(f, y) -# @test fy ≈ linear_combination(basis, θ, x) -# @test f′y ≈ ForwardDiff.derivative(x -> linear_combination(basis, θ, x), x) / h′ -# end -# end - -# @testset "InfRational endpoints continuity" begin -# trans = InfRational(0.0, 1.0) -# for x ∈ (-1.0, 1.0) -# y = from_pm1(trans, x) -# h′ = ForwardDiff.derivative(x -> from_pm1(trans, x), x) - -# for i in 1:N -# θ = e_i(basis, i) -# f = y -> linear_combination(basis, θ, to_pm1(trans, y)) -# fy, f′y = f_f′(f, y) -# @test fy ≈ linear_combination(basis, θ, x) -# @test f′y ≈ ForwardDiff.derivative(x -> linear_combination(basis, θ, x), x) / h′ -# end -# end -# end -# end diff --git a/test/test_transformations.jl b/test/test_transformations.jl index 898318d..b37d22d 100644 --- a/test/test_transformations.jl +++ b/test/test_transformations.jl @@ -1,3 +1,5 @@ +using SpectralKit: PM1 + @testset "bounded linear domain transformations" begin @test_throws DomainError BoundedLinear(-1.0, Inf) @test_throws DomainError BoundedLinear(-1.0, -2.0) @@ -44,6 +46,15 @@ end end @test transform_to(PM1(), trans, y) ≈ x end + + # compare to analytical limits NOTE extend when we add more derivatives + y_pinf = @inferred transform_to(PM1(), trans, 𝑑(Inf)) + @test y_pinf[0] == 1 == @inferred transform_to(PM1(), trans, Inf) + @test y_pinf[1] == 0 + + y_minf = @inferred transform_to(PM1(), trans, 𝑑(-Inf)) + @test y_minf[0] == 1 == @inferred transform_to(PM1(), trans, -Inf) + @test y_minf[1] == 0 end @testset "infinite domain transformations" begin @@ -70,6 +81,15 @@ end end @test transform_to(PM1(), trans, y) ≈ x end + + # compare to analytical limits NOTE extend when we add more derivatives + y_pinf = @inferred transform_to(PM1(), trans, 𝑑(Inf)) + @test y_pinf[0] == 1 == transform_to(PM1(), trans, Inf) + @test y_pinf[1] == 0 + + y_minf = @inferred transform_to(PM1(), trans, 𝑑(-Inf)) + @test y_minf[0] == -1 == transform_to(PM1(), trans, -Inf) + @test y_minf[1] == 0 end @testset "coordinate transformations" begin diff --git a/test/utilities.jl b/test/utilities.jl index c75778c..c24ed19 100644 --- a/test/utilities.jl +++ b/test/utilities.jl @@ -4,18 +4,6 @@ using SpectralKit: TransformedBasis, SmolyakBasis, SmolyakIndices # dispatch for rand_in_domain -struct KroneckerVector{T} <: AbstractVector{T} - i::Int - len::Int -end - -Base.size(v::KroneckerVector) = (v.len, ) - -@inline function Base.getindex(v::KroneckerVector{T}, i::Int) where {T} - @boundscheck 1 ≤ i ≤ v.len - i == v.i ? one(T) : zero(T) -end - chebyshev_cos(x, n) = cos((n - 1) * acos(x)) function chebyshev_cos_deriv(x, n) From c3cb53fe5a055ef5946a1584791406af8050a495 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tam=C3=A1s=20K=2E=20Papp?= Date: Tue, 18 Jun 2024 12:06:46 +0200 Subject: [PATCH 14/14] try to fix documenter workflow --- .github/workflows/CI.yml | 23 --------------- .github/workflows/documentation.yml | 27 ++++++++++++++++++ docs/src/index.md | 6 ++-- src/derivatives.jl | 43 ++++++++++++++++++++++++++++- 4 files changed, 72 insertions(+), 27 deletions(-) create mode 100644 .github/workflows/documentation.yml diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index 53decd2..af924b8 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -45,26 +45,3 @@ jobs: - uses: codecov/codecov-action@v1 with: file: lcov.info - docs: - name: Documentation - runs-on: ubuntu-latest - steps: - - uses: actions/checkout@v2 - - uses: julia-actions/setup-julia@v1 - with: - version: '1' - - run: | - julia --project=docs -e ' - using Pkg - Pkg.develop(PackageSpec(path=pwd())) - Pkg.instantiate()' - - run: | - julia --project=docs -e ' - using Documenter: doctest, DocMeta - using SpectralKit - DocMeta.setdocmeta!(SpectralKit, :DocTestSetup, :(using SpectralKit); recursive=true) - doctest(SpectralKit)' - - run: julia --project=docs docs/make.jl - env: - GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} - DOCUMENTER_KEY: ${{ secrets.DOCUMENTER_KEY }} diff --git a/.github/workflows/documentation.yml b/.github/workflows/documentation.yml new file mode 100644 index 0000000..d1f5fb2 --- /dev/null +++ b/.github/workflows/documentation.yml @@ -0,0 +1,27 @@ +# see https://juliadocs.github.io/Documenter.jl/dev/man/hosting/#GitHub-Actions +# uncomment file below once you set up authentication as described in the Documenter manual + +name: Documentation + +# on: +# push: +# branches: +# - master +# tags: '*' +# pull_request: + +# jobs: +# build: +# runs-on: ubuntu-latest +# steps: +# - uses: actions/checkout@v2 +# - uses: julia-actions/setup-julia@latest +# with: +# version: '1.6' +# - name: Install dependencies +# run: julia --project=docs/ -e 'using Pkg; Pkg.develop(PackageSpec(path=pwd())); Pkg.instantiate()' +# - name: Build and deploy +# env: +# GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} # If authenticating with GitHub Actions token +# DOCUMENTER_KEY: ${{ secrets.DOCUMENTER_KEY }} # If authenticating with SSH deploy key +# run: julia --project=docs/ docs/make.jl diff --git a/docs/src/index.md b/docs/src/index.md index a7a6a25..99323d9 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -10,7 +10,7 @@ This package was designed primarily for solving functional equations, as usually 4. allocation-free, thread safe linear combinations for the above with a given set of coefficients, 5. using [static arrays](https://github.com/JuliaArrays/StaticArrays.jl) extensively to avoid allocation and unroll *some* loops. -While there is some functionality in this package to *fit* approximations to existing functions, it is not ideal for that, as it was optimized for mapping a set of coefficients to residuals of functional equations at gridpoints. +While there is some functionality in this package to *fit* approximations to existing functions, it does not use optimized algorithms (DCT) for that, as it was optimized for mapping a set of coefficients to residuals of functional equations at gridpoints. Also, while the package should interoperate seamlessly with most AD frameworks, only the derivative API (explained below) is guaranteed to have correct derivatives of limits near infinity. @@ -169,10 +169,10 @@ is_subset_basis !!! note API for derivatives is still experimental and subject to change. -For univariate functions, use [`derivatives`](@ref). For multivariate functions, use partial derivatives with `∂`. +For univariate functions, use [`𝑑`](@ref). For multivariate functions, use partial derivatives with [`∂`](@ref). ```@docs -derivatives +𝑑 ∂ ``` diff --git a/src/derivatives.jl b/src/derivatives.jl index fd43851..d36eef5 100644 --- a/src/derivatives.jl +++ b/src/derivatives.jl @@ -86,6 +86,9 @@ julia> (𝑑^3)(2.0) ``` Note that non-literal exponentiation requires `^Val(y)`, for type stability. + +See [`linear_combination`](@ref) for examples of evaluating derivatives of basis +functions and linear combinations. """ const 𝑑 = 𝑑Derivatives{1}() @@ -211,7 +214,6 @@ end function Partials(I::Integer...) N = length(I) while N > 0 && I[N] == 0 - @show N N -= 1 end Partials(ntuple(i -> Int(I[i]), N)) @@ -355,6 +357,45 @@ struct ∂Derivatives{Ps} end end +""" +$(SIGNATURES) + +Partial derivatives along the given coordinates. + +The following are equivalent, and represent ``\\partial_1 \\partial^2_2``, ie the first +derivative along the first axis, and the second partial derivative along the second +axis. + +```@jldoctest +julia> ∂(1, 2) +∂(1, 2) + +julia> ∂((1, 2)) +∂(1, 2) +``` + +Only the vararg form allows trailing zeros, which are stripped: +```@jldoctest +julia> ∂(1, 0) +∂(1) + +julia> ∂((1, 0)) +ERROR: ArgumentError: I ≡ () || last(I) ≠ 0 must hold. +``` + +Use the empty form for no derivatives: +```@jldoctest +julia> ∂() +∂() +``` + +Combine derivatives using `union` or `∪`: + +```jldoctest +julia> ∂(1, 2) ∪ ∂(2, 1) +union(∂(2, 1), ∂(1, 2)) +``` +""" ∂(I::Tuple{Vararg{Int}}) = ∂Derivatives{Partials(I)}() ∂(I::Integer...) = ∂Derivatives{Partials(I...)}()