diff --git a/CHANGELOG.md b/CHANGELOG.md index 3f643e8..6687ccb 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,5 +1,25 @@ + # Unreleased + +## Unreleased + +### Added +### Changed + +- change transformation interface, informatively named keyword arguments +- rename grid kinds +- + +### Deprecated +### Removed +### Fixed +### Security + +# v0.16.2 + +No changelog kept between 0.4.0 and 0.16.2, sorry. + # v0.4.0 - add `linear_combination`, rename `evaluate` to `basis_function` [[#6](https://github.com/tpapp/SpectralKit.jl/pull/6)] diff --git a/Project.toml b/Project.toml index a646fac..12d2a7d 100644 --- a/Project.toml +++ b/Project.toml @@ -9,6 +9,8 @@ projects = ["test", "docs"] [deps] ArgCheck = "dce04be8-c92d-5529-be00-80e4d2c0e197" Compat = "34da2185-b29b-5c13-b0c7-acf172513d20" +ConcreteStructs = "2569d6c7-a4a2-43d3-a901-331e8e4be471" +ConstructionBase = "187b0558-2788-49d3-abe0-74a17ed4e7c9" DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" InverseFunctions = "3587e190-3f89-42d0-90ee-14403ec27112" OrderedCollections = "bac558e1-5e72-5ebc-8fee-abe8a469f55d" @@ -17,6 +19,8 @@ StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" [compat] ArgCheck = "1, 2" Compat = "4.18.0" +ConcreteStructs = "0.2.7" +ConstructionBase = "1.6.0" DocStringExtensions = "0.8, 0.9" InverseFunctions = "0.1" OrderedCollections = "1, 2" diff --git a/src/SpectralKit.jl b/src/SpectralKit.jl index 4757c2c..ea2198d 100644 --- a/src/SpectralKit.jl +++ b/src/SpectralKit.jl @@ -1,6 +1,7 @@ module SpectralKit using ArgCheck: @argcheck +using ConcreteStructs: @concrete using DocStringExtensions: FUNCTIONNAME, SIGNATURES, TYPEDEF using OrderedCollections: OrderedSet using StaticArrays: MVector, SVector, sacollect @@ -10,9 +11,9 @@ include("derivatives.jl") include("domains.jl") include("transformations.jl") include("generic_api.jl") -include("chebyshev.jl") +include("univariate.jl") include("smolyak_traversal.jl") -include("smolyak_api.jl") -include("experimental.jl") # experimental code is not part of the API, see its module docstring +include("smolyak.jl") +# include("experimental.jl") # experimental code is not part of the API, see its module docstring end # module diff --git a/src/chebyshev.jl b/src/chebyshev.jl deleted file mode 100644 index e499faf..0000000 --- a/src/chebyshev.jl +++ /dev/null @@ -1,142 +0,0 @@ -#### -#### Chebyshev polynomials on [-1,1] -#### - -export Chebyshev - -""" -$(TYPEDEF) - -The first `N` Chebyhev polynomials of the first kind, defined on `[-1,1]`. -""" -struct Chebyshev{K<:AbstractGrid} <: UnivariateBasis - "Grid specification." - grid_kind::K - "The number of basis functions." - N::Int - @doc """ - $(SIGNATURES) - `N` Chebyshev polynomials (of the first kind) on ``[-1, 1]``, with the associated - grid of `grid_kind`. - - # Example - - ```julia - basis = Chebyshev(InteriorGrid(), 10) - ``` - """ - function Chebyshev(grid_kind::K, N::Int) where K - @argcheck N ≥ 1 - new{K}(grid_kind, N) - end -end - -@inline dimension(basis::Chebyshev) = basis.N - -@inline domain(::Chebyshev) = PM1() - -function Base.show(io::IO, chebyshev::Chebyshev) - (; grid_kind, N) = chebyshev - print(io, "Chebyshev polynomials (1st kind), ", grid_kind, ", dimension: ", N) -end - -#### -#### basis function iterator -#### - -struct ChebyshevIterator{T} - x::T - N::Int -end - -Base.eltype(::Type{<:ChebyshevIterator{T}}) where {T} = T - -Base.length(itr::ChebyshevIterator) = itr.N - -function basis_at(basis::Chebyshev, x::Scalar) - ChebyshevIterator(x, basis.N) -end - -function Base.iterate(itr::ChebyshevIterator{T}) where T - (; x) = itr - _one(T), (2, _one(T), x) -end - -function Base.iterate(itr::ChebyshevIterator{T}, (i, fp, fpp)) where T - (; x, N) = itr - i > N && return nothing - f = _sub(_mul(2, x, fp), fpp) - f::T, (i + 1, f, fp) -end - -#### -#### grids -#### - -""" -$(SIGNATURES) - -Return a gridpoint for collocation, with `1 ≤ i ≤ dimension(basis)`. - -`T` is used *as a hint* for the element type of grid coordinates. The actual type can be -broadened as required. Methods are type stable. - -!!! note - Not all grids have this method defined, especially if it is impractical. See - [`grid`](@ref), which is part of the API, this function isn't. -""" -function gridpoint(::Type{T}, basis::Chebyshev{InteriorGrid}, i::Integer) where {T <: Real} - (; N) = basis - @argcheck 1 ≤ i ≤ N # FIXME use boundscheck - sinpi((N - 2 * i + 1) / T(2 * N))::T # use formula from Xu (2016) -end - -function gridpoint(::Type{T}, basis::Chebyshev{EndpointGrid}, i::Integer) where {T <: Real} - (; N)= basis - @argcheck 1 ≤ i ≤ N # FIXME use boundscheck - if N == 1 - cospi(1/T(2))::T # 0.0 as a fallback, even though it does not have endpoints - else - cospi((N - i) ./ T(N - 1))::T - end -end - -function gridpoint(::Type{T}, basis::Chebyshev{InteriorGrid2}, i::Integer) where {T <: Real} - (; N)= basis - @argcheck 1 ≤ i ≤ N # FIXME use boundscheck - cospi(((N - i + 1) ./ T(N + 1)))::T -end - -struct ChebyshevGridIterator{T,B} - basis::B -end - -Base.eltype(::Type{<:ChebyshevGridIterator{T}}) where {T} = T - -Base.length(itr::ChebyshevGridIterator) = dimension(itr.basis) - -function Base.iterate(itr::ChebyshevGridIterator{T}, i = 1) where {T} - (; basis) = itr - if i ≤ dimension(basis) - gridpoint(T, basis, i), i + 1 - else - nothing - end -end - -grid(::Type{T}, basis::B) where {T<:Real,B<:Chebyshev} = ChebyshevGridIterator{T,B}(basis) - -#### -#### augmenting -#### - -function augment_coefficients(basis1::Chebyshev{K1}, basis2::Chebyshev{K2}, - θ1::AbstractVector) where {K1,K2} - @argcheck is_subset_basis(basis1, basis2) - @argcheck length(θ1) == dimension(basis1) - vcat(θ1, zeros(eltype(θ1), basis2.N - basis1.N)) -end - -function is_subset_basis(basis1::Chebyshev{K1}, basis2::Chebyshev{K2}) where {K1,K2} - basis2.N ≥ basis1.N -end diff --git a/src/derivatives.jl b/src/derivatives.jl index b4e6b26..6b74b77 100644 --- a/src/derivatives.jl +++ b/src/derivatives.jl @@ -485,7 +485,7 @@ end """ $(SIGNATURES) -Conceptually equivalent to `prod(x))`, which it returns when `kind` is `nothing`, a +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) diff --git a/src/generic_api.jl b/src/generic_api.jl index 1b00d76..64dcfce 100644 --- a/src/generic_api.jl +++ b/src/generic_api.jl @@ -2,9 +2,8 @@ ##### Generic API ##### -export is_function_basis, dimension, basis_at, linear_combination, InteriorGrid, - InteriorGrid2, EndpointGrid, grid, collocation_matrix, augment_coefficients, - is_subset_basis, transformation +export is_function_basis, dimension, basis_at, linear_combination, grid, collocation_matrix, + adjust_coefficients """ $(TYPEDEF) @@ -17,8 +16,6 @@ abstract type FunctionBasis end Broadcast.broadcastable(basis::FunctionBasis) = Ref(basis) -abstract type UnivariateBasis <: FunctionBasis end - abstract type MultivariateBasis <: FunctionBasis end """ @@ -31,7 +28,7 @@ interface: - [`domain`](@ref) for querying the domain, -- [`dimension`](@ref) for the dimension, +- [`parameter_dimension`](@ref) for the dimension, - [`basis_at`](@ref) for function evaluation, @@ -69,7 +66,7 @@ function domain end """ `$(FUNCTIONNAME)(basis)` -Return the dimension of `basis`, a positive `Int`. +Return the dimension parameters of `basis`, a positive `Int`. """ function dimension end @@ -159,38 +156,6 @@ though working with `basis ∘ transformation` may be preferred. """ linear_combination(basis, θ) = LinearCombination(basis, θ) -""" -$(TYPEDEF) - -Abstract type for all grid specifications. -""" -abstract type AbstractGrid end - -""" -$(TYPEDEF) - -Grid with interior points (eg Gauss-Chebyshev). -""" -struct InteriorGrid <: AbstractGrid end - -""" -$(TYPEDEF) - -Grid that includes endpoints (eg Gauss-Lobatto). - -!!! note - For small dimensions may fall back to a grid that does not contain endpoints. -""" -struct EndpointGrid <: AbstractGrid end - -""" -$(TYPEDEF) - -Grid with interior points that results in smaller grids than `InteriorGrid` when nested. -Equivalent to an `EndpointGrid` with endpoints dropped. -""" -struct InteriorGrid2 <: AbstractGrid end - """ `$(FUNCTIONNAME)([T], basis)` @@ -227,104 +192,92 @@ function collocation_matrix(basis, x = grid(basis)) end """ -`$(FUNCTIONNAME)(basis1, basis2, θ1)` +`$(FUNCTIONNAME)(θ1, basis1, basis2)` -Return a set of coefficients `θ2` for `basis2` such that +Return a set of coefficients `θ2` for `basis2` such that is either augmentes ```julia -linear_combination(basis1, θ1, x) == linear_combination(basis2, θ2, x) +linear_combination(basis1, θ1, x) ≈ linear_combination(basis2, θ2, x) ``` -for any `x` in the domain. In practice this means padding with zeros. +for any `x` in the domain. In practice this means either padding with zeros, or +truncating the coefficients. -Throw a `ArgumentError` if the bases are incompatible with each other or `x`, or this is not -possible. Methods may not be defined for incompatible bases, compatibility between bases can -be checked with [`is_subset_basis`](@ref). +Throw an error if this is not possible. """ -function augment_coefficients end +function adjust_coefficients end -""" -$(SIGNATURES) +# FIXME remove this +# #### +# #### transformed basis +# #### -Return a `Bool` indicating whether coefficients in `basis1` can be augmented to `basis2` -with [`augment_coefficients`](@ref). +# """ +# Transform the domain of a basis. +# """ +# struct TransformedBasis{B,T} <: FunctionBasis +# parent::B +# transformation::T +# function TransformedBasis(parent::B, transformation::T) where {B,T} +# @argcheck domain_kind(domain(parent)) ≡ domain_kind(T) +# new{B,T}(parent, transformation) +# end +# end -!!! note - `true` does not mean that coefficients from `basis1` can just be padded with zeros, - since they may be in different positions. Always use [`augment_coefficients`](@ref). -""" -is_subset_basis(basis1::FunctionBasis, basis2::FunctionBasis) = false +# function Base.:(∘)(parent::FunctionBasis, transformation) +# TransformedBasis(parent, transformation) +# end -#### -#### transformed basis -#### +# Base.parent(basis::TransformedBasis) = basis.parent -""" -Transform the domain of a basis. -""" -struct TransformedBasis{B,T} <: FunctionBasis - parent::B - transformation::T - function TransformedBasis(parent::B, transformation::T) where {B,T} - @argcheck domain_kind(domain(parent)) ≡ domain_kind(T) - new{B,T}(parent, transformation) - end -end +# """ +# $(SIGNATURES) -function Base.:(∘)(parent::FunctionBasis, transformation) - TransformedBasis(parent, transformation) -end +# Return the transformation of transformed bases, or `nothing` it not applicable. +# """ +# transformation(basis::TransformedBasis) = basis.transformation -Base.parent(basis::TransformedBasis) = basis.parent +# transformation(::FunctionBasis) = nothing -""" -$(SIGNATURES) +# domain(basis::TransformedBasis) = domain(basis.transformation) -Return the transformation of transformed bases, or `nothing` it not applicable. -""" -transformation(basis::TransformedBasis) = basis.transformation +# dimension(basis::TransformedBasis) = dimension(basis.parent) -transformation(::FunctionBasis) = nothing +# function basis_at(basis::TransformedBasis, x) +# (; parent, transformation) = basis +# basis_at(parent, transform_to(domain(parent), transformation, x)) +# end -domain(basis::TransformedBasis) = domain(basis.transformation) +# function grid(::Type{T}, basis::TransformedBasis) where T +# (; parent, transformation) = basis +# d = domain(parent) +# Iterators.map(x -> transform_from(d, transformation, x), grid(T, parent)) +# end -dimension(basis::TransformedBasis) = dimension(basis.parent) +# function Base.:(∘)(linear_combination::LinearCombination, transformation) +# (; basis, θ) = linear_combination +# LinearCombination(basis ∘ transformation, θ) +# end -function basis_at(basis::TransformedBasis, x) - (; parent, transformation) = basis - basis_at(parent, transform_to(domain(parent), transformation, x)) -end - -function grid(::Type{T}, basis::TransformedBasis) where T - (; parent, transformation) = basis - d = domain(parent) - Iterators.map(x -> transform_from(d, transformation, x), grid(T, parent)) -end +# Base.length(basis::TransformedBasis{<:MultivariateBasis}) = length(basis.parent) -function Base.:(∘)(linear_combination::LinearCombination, transformation) - (; basis, θ) = linear_combination - LinearCombination(basis ∘ transformation, θ) -end +# function Base.getindex(basis::TransformedBasis{<:MultivariateBasis}, i::Int) +# (; parent, transformation) = basis +# TransformedBasis(parent[i], Tuple(transformation)[i]) +# end -Base.length(basis::TransformedBasis{<:MultivariateBasis}) = length(basis.parent) +# function is_subset_basis(basis1::TransformedBasis, basis2::TransformedBasis) +# basis1.transformation ≡ basis2.transformation && +# is_subset_basis(basis1.parent, basis2.parent) +# end -function Base.getindex(basis::TransformedBasis{<:MultivariateBasis}, i::Int) - (; parent, transformation) = basis - TransformedBasis(parent[i], Tuple(transformation)[i]) -end +# function augment_coefficients(basis1::TransformedBasis, basis2::TransformedBasis, θ1) +# @argcheck is_subset_basis(basis1, basis2) +# augment_coefficients(basis1.parent, basis2.parent, θ1) +# end -function is_subset_basis(basis1::TransformedBasis, basis2::TransformedBasis) - basis1.transformation ≡ basis2.transformation && - is_subset_basis(basis1.parent, basis2.parent) -end +# function transform_to(basis::FunctionBasis, transformation, x) +# transform_to(domain(basis), transformation, x) +# end -function augment_coefficients(basis1::TransformedBasis, basis2::TransformedBasis, θ1) - @argcheck is_subset_basis(basis1, basis2) - augment_coefficients(basis1.parent, basis2.parent, θ1) -end - -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 +# function transform_from(basis::FunctionBasis, transformation, x) +# transform_from(domain(basis), transformation, x) +# end diff --git a/src/smolyak.jl b/src/smolyak.jl new file mode 100644 index 0000000..be4e219 --- /dev/null +++ b/src/smolyak.jl @@ -0,0 +1,280 @@ +##### +##### Smolyak bases +##### + +export SmolyakLevel, SmolyakBasis + +import ConstructionBase + +struct SmolyakLevel + total::Int + each::Int + @doc """ + $(SIGNATURES) + + Level specification for a Smolyak basis. + + `total` constrains the *sum* of levels in all axes. + + `each` constrains the level along *each* axis. + + Formally, ``0 ≤ ℓᵢ ≤ each ∀ i; ∑ᵢ ℓᵢ ≤ total`` where each `i` is an axis. + + If `each > total`, the normalization `each = total` is set with a warning. + """ + function SmolyakLevel(; total::Int, each::Int = total) + @argcheck total ≥ 0 + @argcheck each ≥ 0 + if each > total + @warn "‘each’ normalized to ‘total’" each total + each = total + end + new(total, each) + end +end + +function Base.show(io::IO, level::SmolyakLevel) + (; total, each) = level + print(io, "SmolyakLevel(total = $(total), each = $(each))") + if !get(io, :compact, false) + printstyled(io, " #= 0 ≤ ∑ℓᵢ ≤ $(total), all 0 ≤ ℓᵢ ≤ $(each) =#"; + color = INFO_COLOR) + end +end + +function ConstructionBase.setproperties(level::SmolyakLevel, patch::NamedTuple) + ConstructionBase.setproperties(level; patch...) +end + +function ConstructionBase.setproperties(level::SmolyakLevel; + total = level.total, each = level.each) + SmolyakLevel(; total = Int(total), each = Int(each)) +end + +""" +$(SIGNATURES) + +A wrapper for iterating through Smolyak indices. See [`__smolyak__init`](@ref) and +[`__smolyak_step`](@ref). Can be plugged straight into `iterate`. + +`g` transforms the value. +""" +function __smolyak_iterate(family, kind, level::SmolyakLevel, f::F, itrs, g::G, + state = nothing) where {F,G} + if state ≡ nothing + accum, state... = __smolyak_init(family, kind, level.total, f, itrs) + g(accum), (accum, state...) + else + accum, Δ, rest... = __smolyak_step(family, kind, level.each, f, itrs, state...) + if Δ == Δ_DONE + nothing + else + slack = state[2] + g(accum), (accum, slack + Δ, rest...) + end + end +end + +struct SmolyakBasis{F,K,D} <: MultivariateBasis + family::F + kind::K + domain_transformations::D + level::SmolyakLevel + grid_level::SmolyakLevel + @doc """ + $(SIGNATURES) + + Smolyak (sparse) basis from `family`, using the given `kind`. + + `domain_transformations::Tuple` determine the domain and the relevant transformation + to the domain of `family` along each coordinate. + + `level` is a [`SmolyakLevel`](@ref), specifying the number of *blocks* used to build + the grid. + + ### Miscellaneous + + As a convenience, `smolyak_basis[i]` returns an univariate basis along coordinate `i`. + """ + function SmolyakBasis(family::F, kind::K, domain_transformations::D, + level::SmolyakLevel, + grid_level::SmolyakLevel = level) where {F,K,D<:Tuple} + new{F,K,D}(family, kind, domain_transformations, level, grid_level) + end +end + +function Base.show(io::IO, basis::SmolyakBasis) + (; family, kind, domain_transformations, level, grid_level) = basis + _print_dimensions(io, basis) + lead = "SmolyakBasis(" + next = ",\n" * ' '^length(lead) + print(io, lead, family, ", ", kind, next, + domain_transformations, next, + level) + if grid_level ≠ level + print(io, next, grid_level) + end + print(io, ")") +end + +function Base.getindex(basis::SmolyakBasis, i::Integer) + (; family, kind, domain_transformations, level, grid_level) = basis + UnivariateBasis(family, kind, domain_transformations[i], level.each, grid_level.each) +end + +# ## Example + +# FIXME these examples need to be updated +# ```jldoctest +# julia> basis = smolyak_basis(Chebyshev, InteriorGrid(), SmolyakParameters(3), 2) +# Sparse multivariate basis on ℝ² +# Smolyak indexing, ∑bᵢ ≤ 3, all bᵢ ≤ 3, dimension 81 +# using Chebyshev polynomials (1st kind), InteriorGrid(), dimension: 27 + +# julia> dimension(basis) +# 81 + +# julia> domain(basis) +# [-1,1]² +# ``` + +# ## Properties + +# *Grids nest*: increasing arguments of `SmolyakParameters` result in a refined grid that +# contains points of the cruder grid. +# """ + +function domain(smolyak_basis::SmolyakBasis) + map(domain, smolyak_basis.domain_transformations) +end + +function dimension(smolyak_basis::SmolyakBasis) + (; family, kind, domain_transformations, level) = smolyak_basis + N = length(domain_transformations) + __smolyak_length(family, kind, Val(N), level.total, level.each) +end + +struct SmolyakBasisAt{I,P,F,K,L<:SmolyakLevel} + family::F + kind::K + level::L + itrs::I + product_kind::P +end + +function Base.eltype(::Type{<:SmolyakBasisAt{I,P}}) where {I,P} + _product_type(P, map(_eltype, fieldtypes(I))) +end + +function Base.length(itr::SmolyakBasisAt) + (; family, kind, level, itrs, product_kind) = itr + N = length(itrs) + __smolyak_length(family, kind, Val(N), level.total, level.each) +end + +_prod(a::Real, b::Tuple{}) = (a, ) + +_prod(a::Real, b::Tuple) = (a * first(b), b...) + +function Base.iterate(itr::SmolyakBasisAt, state = nothing) + (; family, kind, level, itrs, product_kind) = itr + # FIXME this does not yet work for derivatives + __smolyak_iterate(family, kind, level, _prod, itrs, first, state) +end + +function basis_at(smolyak_basis::SmolyakBasis, x::Tuple) + (; family, kind, domain_transformations, level) = smolyak_basis + @argcheck length(x) == length(domain_transformations) + itrs = map((x, d) -> _univariate_basis_itr(family, d, x), x, domain_transformations) + SmolyakBasisAt(family, kind, level, itrs, nothing) +end + +function basis_at(smolyak_basis::SmolyakBasis, x::SVector) + basis_at(smolyak_basis, Tuple(x)) +end + +# function basis_at(smolyak_basis::SmolyakBasis, Dx::∂CoordinateExpansion) +# (; family, kind, domain_transformations, level) = smolyak_basis +# (; ∂D, x) = Dx +# itrs = map((x, d) -> basis_at(UnivariateBasis(family, kind, d, each), x), +# x, domain_transformations) +# BasisAt(family, kind, total, each, itrs, ∂D) +# end + +struct SmolyakGrid{I,F,K,S<:SmolyakLevel} + family::F + kind::K + level::S + itrs::I +end + +function Base.eltype(::Type{<:SmolyakGrid{I}}) where {I} + Tuple{map(_eltype, fieldtypes(I))...} +end + +function Base.length(itr::SmolyakGrid) + __smolyak_length(itr.family, itr.kind, Val(length(itr.itrs)), itr.level.total, itr.level.each) +end + +function Base.iterate(itr::SmolyakGrid, state = nothing) + (; family, kind, level, itrs) = itr + __smolyak_iterate(family, kind, level, (a, b) -> (a, b...), itrs, identity, state) +end + +function grid(::Type{T}, smolyak_basis::SmolyakBasis) where {T<:AbstractFloat} + (; family, kind, domain_transformations, grid_level) = smolyak_basis + (; each) = grid_level + itrs = map(d -> grid(T, UnivariateBasis(family, kind, d, each)), domain_transformations) + SmolyakGrid(family, kind, grid_level, itrs) +end + +#### +#### index traversal for adjustment +#### + +@concrete struct SmolyakIndices{I} + family + kind + level + itrs::I +end + +Base.eltype(::Type{<:SmolyakIndices{I}}) where I = NTuple{length(fieldtypes(I)),Int} + +Base.IteratorSize(::SmolyakIndices) = Base.SizeUnknown() + +""" +$(SIGNATURES) + +Iterate through indices of a Smolyak basis. For each value, the indices in the tuples +correspond to the index of that univariate basis function. Not part of the API. + +Note: only used as a building block for [`adjust_coefficients`](@ref). +""" +function SmolyakIndices(basis::SmolyakBasis) + (; family, kind, domain_transformations, level) = basis + itrs = ntuple(_ -> Counting(), Val(length(domain_transformations))) + SmolyakIndices(family, kind, level, itrs) +end + +function Base.iterate(itr::SmolyakIndices, state = nothing) + (; family, kind, level, itrs) = itr + __smolyak_iterate(family, kind, level, (a, b) -> (a, b...), itrs, identity, state) +end + +function adjust_coefficients(θ1::AbstractVector{T}, basis1::SmolyakBasis, + basis2::SmolyakBasis) where T + (; family, kind, level, domain_transformations) = basis1 + @argcheck family ≡ basis2.family + @argcheck domain_transformations == basis2.domain_transformations + if kind ≡ basis2.kind && level == basis2.level + # these are the same bases + return copy(θ1) + end + θ = Dict{NTuple{length(domain_transformations),Int},T}() + for (x, ι) in zip(θ1, SmolyakIndices(basis1)) + θ[ι] = x + end + z = zero(T) + [get(θ, ι, z) for ι in SmolyakIndices(basis2)] +end diff --git a/src/smolyak_api.jl b/src/smolyak_api.jl deleted file mode 100644 index bb4c110..0000000 --- a/src/smolyak_api.jl +++ /dev/null @@ -1,208 +0,0 @@ -##### -##### Smolyak bases -##### - -export smolyak_basis - -struct SmolyakBasis{I<:SmolyakIndices,U<:UnivariateBasis} <: MultivariateBasis - smolyak_indices::I - univariate_parent::U -end - -function Base.show(io::IO, smolyak_basis::SmolyakBasis{<:SmolyakIndices{N}}) where N - (; smolyak_indices, univariate_parent) = smolyak_basis - print(io, "Sparse multivariate basis on ℝ", SuperScript(N), "\n ", smolyak_indices, - "\n using ", univariate_parent) -end - -Base.length(basis::SmolyakBasis{<:SmolyakIndices{N}}) where N = N - -function Base.getindex(basis::SmolyakBasis, i::Int) - @argcheck 1 ≤ i ≤ length(basis) BoundsError(basis, i) - basis.univariate_parent -end - -""" -$(SIGNATURES) - -Create a sparse Smolyak basis. - -# Arguments - -- `univariate_family`: should be a callable that takes a `grid_kind` and a `dimension` - parameter, eg `Chebyshev`. - -- `grid_kind`: the grid kind, eg `InteriorGrid()` etc. - -- `smolyak_parameters`: the Smolyak grid specification parameters, see - [`SmolyakParameters`](@ref). - -- `N`: the dimension. wrapped in a `Val` for type stability, a convenience constructor also - takes integers. - -## Example - -```jldoctest -julia> basis = smolyak_basis(Chebyshev, InteriorGrid(), SmolyakParameters(3), 2) -Sparse multivariate basis on ℝ² - Smolyak indexing, ∑bᵢ ≤ 3, all bᵢ ≤ 3, dimension 81 - using Chebyshev polynomials (1st kind), InteriorGrid(), dimension: 27 - -julia> dimension(basis) -81 - -julia> domain(basis) -[-1,1]² -``` - -## Properties - -*Grids nest*: increasing arguments of `SmolyakParameters` result in a refined grid that -contains points of the cruder grid. -""" -function smolyak_basis(univariate_family, grid_kind::AbstractGrid, - smolyak_parameters::SmolyakParameters, ::Val{N}) where {N} - @argcheck N ≥ 1 - smolyak_indices = SmolyakIndices{N}(grid_kind, smolyak_parameters) - univariate_parent = univariate_family(grid_kind, highest_visited_index(smolyak_indices)) - SmolyakBasis(smolyak_indices, univariate_parent) -end - -# convenience constructor -@inline function smolyak_basis(univariate_family, grid_kind::AbstractGrid, - smolyak_parameters::SmolyakParameters, N::Integer) - smolyak_basis(univariate_family, grid_kind, smolyak_parameters, Val(N)) -end - -function domain(smolyak_basis::SmolyakBasis{<:SmolyakIndices{N}}) where {N} - D = domain(smolyak_basis.univariate_parent) - coordinate_domains(Val(N), D) -end - -dimension(smolyak_basis::SmolyakBasis) = length(smolyak_basis.smolyak_indices) - -""" -$(SIGNATURES) - -Helper function to make univariate bases for a Smolyak basis. -""" -function _univariate_bases_at(smolyak_basis::SmolyakBasis{<:SmolyakIndices{N,H}}, - x) where {N,H} - (; univariate_parent) = smolyak_basis - map(x -> sacollect(SVector{H}, basis_at(univariate_parent, x)), x) -end - -function basis_at(smolyak_basis::SmolyakBasis{<:SmolyakIndices{N}}, - x::Union{Tuple,AbstractVector}) where {N} - @argcheck length(x) == N - SmolyakProduct(smolyak_basis.smolyak_indices, - _univariate_bases_at(smolyak_basis, NTuple{N}(x)), - nothing) -end - -function basis_at(smolyak_basis::SmolyakBasis{<:SmolyakIndices{N}}, - Dx::∂CoordinateExpansion) where {N} - (; ∂D, x) = Dx - @argcheck length(x) == N - SmolyakProduct(smolyak_basis.smolyak_indices, - _univariate_bases_at(smolyak_basis, x), - ∂D) -end - -struct SmolyakGridIterator{T,I,S} - smolyak_indices::I - sources::S -end - -Base.eltype(::Type{<:SmolyakGridIterator{T}}) where {T} = T - -Base.length(itr::SmolyakGridIterator) = length(itr.smolyak_indices) - -function grid(::Type{T}, - smolyak_basis::SmolyakBasis{<:SmolyakIndices{N,H}}) where {T<:Real,N,H} - (; smolyak_indices, univariate_parent) = smolyak_basis - sources = sacollect(SVector{H}, gridpoint(T, univariate_parent, i) - for i in SmolyakGridShuffle(univariate_parent.grid_kind, H)) - SmolyakGridIterator{NTuple{N,T},typeof(smolyak_indices),typeof(sources)}(smolyak_indices, sources) -end - -function Base.iterate(itr::SmolyakGridIterator, state...) - (; smolyak_indices, sources) = itr - result = iterate(smolyak_indices, state...) - result ≡ nothing && return nothing - ι, state′ = result - map(i -> sources[i], ι), state′ -end - -""" -$(SIGNATURES) - -Utility function to check is `basis1` is a subset of `basis2` with shared indices. -""" -function _is_shared_index_subset(basis1::Chebyshev{K1}, basis2::Chebyshev{K2}) where {K1,K2} - K1 == K2 && basis1.N ≤ basis2.N -end - -function is_subset_basis(basis1::SmolyakBasis{<:SmolyakIndices{N1,H1,B1,M1}}, - basis2::SmolyakBasis{<:SmolyakIndices{N2,H2,B2,M2}}) where {N1,H1,B1,M1,N2,H2,B2,M2} - (N1 == N2 && B2 ≥ B1 && M2 ≥ M1 && - # NOTE: traversal relies on the same (column major) ordering of indices in both - # bases. Testing for this is currently innocuous, as Chebyshev has this property. - # If some basis is added to the code which doesn't this should be tested for in - # `augment_coefficients` which should then use a different code path. - _is_shared_index_subset(basis1.univariate_parent, basis2.univariate_parent)) -end - -""" -$(TYPEDEF) - -Given two iterations `ι1 ∈ itr1` and `ι2 ∈ itr2`, and a vector `θ1` such that `length(θ1) == -length(itr1)`, return an iterator that returns elements of `θ1` when `ι1 == ι2` and zero -otherwise. - -# Internals - -state is a tuple of: - -- index for the next upcoming element of `θ1`, -- the next item in `itr1`, set to `(0, 0, …)` after all of them are used -- the corresponding iterator state (ignore for sentinel value `(0, 0, …)` -- state of `itr2` (only after the first call to `iterate`) -""" -struct PaddingIterator{V1,I1,I2} - θ1::V1 - itr1::I1 - itr2::I2 -end - -Base.length(itr::PaddingIterator) = length(itr.itr2) - -Base.eltype(itr::PaddingIterator) = eltype(itr.θ1) - -function Base.iterate(itr::PaddingIterator, state = (firstindex(itr.θ1), - iterate(itr.itr1)...)) - (; θ1, itr1, itr2) = itr - i, ι1, state1, state2... = state - res2 = iterate(itr2, state2...) - res2 ≡ nothing && return nothing - ι2, state2 = res2 - if ι1 == ι2 - x = itr.θ1[i] - res1 = iterate(itr.itr1, state1) - if res1 ≡ nothing - ι1 = map(_ -> 0, ι1) # sentinel ensures never visiting ι1 == ι2 branch again - else - ι1, state1 = res1 - i += 1 - end - else - x = zero(eltype(itr.θ1)) - end - x, (i, ι1, state1, state2) -end - -function augment_coefficients(basis1::SmolyakBasis, basis2::SmolyakBasis, θ1::AbstractVector) - @argcheck is_subset_basis(basis1, basis2) - @argcheck dimension(basis1) == length(θ1) - collect(PaddingIterator(θ1, basis1.smolyak_indices, basis2.smolyak_indices)) -end diff --git a/src/smolyak_traversal.jl b/src/smolyak_traversal.jl index 13b665a..0fe95bd 100644 --- a/src/smolyak_traversal.jl +++ b/src/smolyak_traversal.jl @@ -2,230 +2,138 @@ ##### 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*. +$(SIGNATURES) → accum, slack, remainders, states, levels -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. +Initialize the state for [`__smolyak_step`](@ref), states are documented there. """ -@inline function SmolyakParameters(B::Integer, M::Integer = B) - SmolyakParameters{Int(B),Int(M)}() +function __smolyak_init(family, kind, total::Int, f::F, itrs::NTuple{N,Any}) where {F,N} + slack = total + r = block_length(family, kind, 0) + remainders = ntuple(_ -> r - 1, Val(N)) + itr_results = map(_start, itrs) + states = map(last, itr_results) + levels = ntuple(_ -> 0, Val(N)) + accum = foldr(f, map(first, itr_results); init = ()) + accum, slack, remainders, states, levels end -#### -#### Nesting sizes and shuffling -#### -#### NOTE: This is not exported as we have no API for nesting univariate bases, only -#### Smolyak. When refactoring, consider exporting with a unified API. - """ -$(TYPEDEF) - -An iterator of indices for picking elements from a grid of length `len`, which should be a -valid cumulative block length. +Sentinel value for having finished iteration in `__smolyak_step`. """ -struct SmolyakGridShuffle{K} - grid_kind::K - len::Int -end - -Base.length(ι::SmolyakGridShuffle) = ι.len - -Base.eltype(::Type{<:SmolyakGridShuffle}) = Int - -### -### endpoint grid: 1, 1, 2, 4, 8, … -### +const Δ_DONE = typemax(Int) """ -$(SIGNATURES) +$(SIGNATURES) → accum, Δ, remainders′, states′, levels′ -Cumulative block length at block `b`. -""" -@inline function nesting_total_length(::Type{Chebyshev}, ::EndpointGrid, b::Int) - b == 0 ? 1 : ((1 << b) + 1) -end +Step through Smolyak indices of iterators. -""" -$(SIGNATURES) +`itrs` yield the `xs`. **The iterators should be stateless.** -Length of each block `b`. +`accum` is equivalent to `foldr(f, xs)`. `f` should map a value and a tuple to a tuple +of one more element. -!!! note - Smolyak grids use “blocks” of polynomials, each indexed by ``b == 0, …, B`, with an - increasing number of points in each. -""" -@inline function nesting_block_length(::Type{Chebyshev}, ::EndpointGrid, b::Int) - b ≤ 1 ? b + 1 : 1 << (b - 1) -end +`slack` is the `total - sum(levels)`. `Δ` is the change in `slack`. Design note: easier +to apply recursively than `slack`, caller should make the adjustment. -function Base.iterate(ι::SmolyakGridShuffle{EndpointGrid}) - i = (ι.len + 1) ÷ 2 - i, (0, 0) # step = 0 is special-cased -end +`remainders` contains the count of elements left in each level before we move to a +different combination. -function Base.iterate(ι::SmolyakGridShuffle{EndpointGrid}, (i, step)) - (; len) = ι - i == 0 && return len > 1 ? (1, (1, len - 1)) : nothing - i′ = i + step - if i′ ≤ len - i′, (i′, step) - else - step′ = step ÷ 2 - if step′ ≥ 2 - i′ = step′ ÷ 2 + 1 - i′, (i′, step′) +`states` are states of iterators. `levels` are the levels currently visited. +""" +function __smolyak_step(family, kind, each::Int, f::F, itrs::NTuple{N,Any}, + accum::NTuple{N,Any}, slack::Int, remainders::NTuple{N,Any}, + states::NTuple{N,Any}, levels::NTuple{N,Any}) where {F,N} + I1, Iτ... = itrs + a1, aτ... = accum + r1, rτ... = remainders + s1, sτ... = states + l1, lτ... = levels + if r1 > 0 # step within block + x1, s1′ = _next(I1, s1) + (f(x1, aτ), + 0, # no change in slack + (r1 - 1, rτ...), # one less element in 1 + (s1′, sτ...), # step iterator + levels) + elseif l1 < each && slack > 0 # next block, same tail + x1, s1′ = _next(I1, s1) + (f(x1, aτ), + -1, # decrease slack + (block_length(family, kind, l1 + 1) - 1, rτ...), # remaining elements: all in block + (s1′, sτ...), # step state 1 + (l1+1, lτ...)) # next level + elseif N == 1 + accum, Δ_DONE, remainders, states, levels # done with iteration + else # go into tail + aτ′, Δτ, rτ′, sτ′, lτ′ = __smolyak_step(family, kind, each, f, Iτ, aτ, slack + l1, rτ, sτ, lτ) + if Δτ == Δ_DONE + accum, Δ_DONE, remainders, states, levels # tail is done with iteration else - nothing + x1, s1 = _start(I1) + (f(x1, aτ′), + l1 + Δτ, # more slack as we reset 1 + (block_length(family, kind, 0) - 1, rτ′...), # all remaining in block 0 + (s1, sτ′...), # states with tail + (0, lτ′...)) # back to level 0 here end end end -### -### interior grid: 1, 3, 9, 27, … -### - -@inline nesting_total_length(::Type{Chebyshev}, ::InteriorGrid, b::Int) = 3^b +struct NonIncreasingSmolyakLevels{N} + total::Int + each::Int + @doc """ + $(SIGNATURES) → itr -@inline function nesting_block_length(::Type{Chebyshev}, ::InteriorGrid, b::Int) - b == 0 ? 1 : 2 * 3^(b - 1) -end + An iterable which yields `ℓ::Ntuple{N,Int}`, with the following properties: -function Base.iterate(ι::SmolyakGridShuffle{InteriorGrid}) - (; len) = ι - i0 = (len + 1) ÷ 2 # first index at this level - Δ = len # basis for step size - a = 2 # alternating as 2Δa and Δa - i0, (i0, i0, Δ, a) -end - -function Base.iterate(ι::SmolyakGridShuffle{InteriorGrid}, (i, i0, Δ, a)) - (; len) = ι - i′ = i + a * Δ - if i′ ≤ len - i′, (i′, i0, Δ, 3 - a) - else - if Δ == 1 - nothing - else - Δ = Δ ÷ 3 - i0 -= Δ - i0, (i0, i0, Δ, 2) - end + 1. each tuple `ℓ` is non-increasing (weakly decreasing), + 2. each element of each tuple is between `0` and `each` (inclusive), + 3. the sum of all elements in each tuple is not larger than total. + """ + function NonIncreasingSmolyakLevels{N}(total, each) where N + @argcheck N isa Integer && N ≥ 1 + @argcheck each ≤ total + new{N}(total, each) end end -### -### interior grid type 2: 1, 3, 7, … -### +Base.eltype(::Type{NonIncreasingSmolyakLevels{N}}) where N = NTuple{N,Int} -@inline nesting_total_length(::Type{Chebyshev}, ::InteriorGrid2, b::Int) = (1 << (b + 1)) - 1 - -@inline nesting_block_length(::Type{Chebyshev}, ::InteriorGrid2, b::Int) = 1 << b - -function Base.iterate(ι::SmolyakGridShuffle{InteriorGrid2}) - i = (ι.len + 1) ÷ 2 - i, (i, 2 * i) -end +Base.IteratorSize(::Type{<:NonIncreasingSmolyakLevels}) = Base.SizeUnknown() -function Base.iterate(ι::SmolyakGridShuffle{InteriorGrid2}, (i, step)) - i′ = i + step - if i′ ≤ ι.len - i′, (i′, step) +function __step_noninc(total::Int, each::Int, Σ_and_indices::Int...) + Σ, i1, iτ... = Σ_and_indices + if Σ < total && isempty(iτ) + # single index, increment if possible + (i1 + 1, i1 + 1) + elseif Σ < total && i1 > last(iτ) + # room to increment indices in the tail + Σ′, iτ′... = __step_noninc(total - i1, i1, Σ - i1, iτ...) + (i1 + Σ′, i1, iτ′...) + elseif !isempty(iτ) && first(iτ) < min(i1, each) && first(iτ) + i1 < total + # can increment next index, zero out tail of tail + i2 = first(iτ) + (i1 + i2 + 1, i1, i2 + 1, ntuple(_ -> 0, Val(length(iτ) - 1))...) else - step′ = step ÷ 2 - if step′ ≥ 2 - i′ = step′ ÷ 2 - i′, (i′, step′) - else - nothing - end + # increment first index; stop iteration if this is > each + (i1 + 1, i1 + 1, ntuple(_ -> 0, Val(length(iτ)))...) end end -#### -#### index traversal -#### - -function __inc_init(nesting_total_lengths, ::Val{N}, ::Val{B}) where {N,B} - indices = ntuple(_ -> 1, Val(N)) - blocks = ntuple(_ -> 0, Val(N)) - l = first(nesting_total_lengths) - limits = ntuple(_ -> l, Val(N)) - slack = B - slack, indices, blocks, limits -end - -""" -$(SIGNATURES) - -Internal implementation of the Smolyak indexing iterator. - -# Arguments - -- `nesting_total_lengths`: precalculated nesting total lengths, constant during iteration, - indexes with an offset of `1` - -- `slack`: `B - sum(blocks)`, cached - -- `indices`: current indices - -- `blocks`: block indexes - -- `limits`: limit for each index (for column-major reset) - -# Return values - -- `valid::Bool`: `false` iff there is no next element, in which case the following values - should be ignored - -- `Δ::Int`: change in `slack` - -- `indices′`, `blocks′, `limits′`: next values for corresponding arguments above, each an - `::NTuple{N,Int}` -""" -@inline function __inc(nesting_total_lengths::NTuple{Mp1,Int}, slack::Int, - indices::NTuple{N,Int}, blocks::NTuple{N,Int}, - limits::NTuple{N,Int}) where {Mp1,N} - i1, iτ... = indices - b1, bτ... = blocks - l1, lτ... = limits - if i1 < l1 # increment i1, same block - true, 0, (i1 + 1, iτ...), blocks, limits - elseif b1 < (Mp1 - 1) && slack > 0 # increment i1, next block - b1′ = b1 + 1 - true, -1, (i1 + 1, iτ...), (b1′, bτ...), (nesting_total_lengths[b1′ + 1], lτ...) +function Base.iterate(itr::NonIncreasingSmolyakLevels{N}, state = nothing) where N + (; total, each) = itr + if state ≡ nothing + indices = ntuple(_ -> 0, N) + indices, (0, indices...) else - if N == 1 # end of iteration, arbitrary value since !valid - false, 0, indices, blocks, limits - else # i1 = 1, increment tail if applicable - Δ1 = b1 - valid, Δτ, iτ′, bτ′, lτ′ = __inc(nesting_total_lengths, slack + Δ1, iτ, bτ, lτ) - valid, Δ1 + Δτ, (1, iτ′...), (0, bτ′...), (nesting_total_lengths[1], lτ′...) + state′ = __step_noninc(total, each, state...) + Σ′, indices′... = state′ + if first(indices′) ≤ each + indices′, state′ + else + nothing end end end @@ -233,147 +141,27 @@ end """ $(SIGNATURES) -Calculate the length of a [`SmolyakIndices`](@ref) iterator. Argument as in the latter. +Calculate the dimension of Smolyak basis. """ -function __smolyak_length(grid_kind::AbstractGrid, ::Val{N}, ::Val{B}, M::Int) where {N,B} - # implicit assumption: M ≤ B, enforced by the SmolyakParameters constructor - _bl(b) = nesting_block_length(Chebyshev, grid_kind, b) - c = zeros(MVector{B+1,Int}) # indexed as 0, …, B - for b in 0:M - c[b + 1] = _bl(b) - end - for n in 2:N - for b in B:(-1):0 # blocks with indices that sum to b - s = 0 - for a in 0:min(b, M) - s += _bl(a) * c[b - a + 1] +function __smolyak_length(family, kind, ::Val{N}, total::Int, each::Int) where N + L = 0 + P = factorial(N) # permutations + for ℓ in NonIncreasingSmolyakLevels{N}(total, each) + p = P # combinations, accounting for repetitions + C = 1 # will contain product of block lengths after loop below + l_prev = -1 # previous value (sentinel) + r = 1 # run counter for repeated values + for l in ℓ + C *= block_length(family, kind, l) + if l == l_prev + r += 1 + p ÷= r # account for repetitions + else + l_prev = l + r = 1 # reset counter end - # can safely overwrite since they will not be used again for n + 1 - c[b + 1] = s end + L += p * C 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′ + L end diff --git a/src/transformations.jl b/src/transformations.jl index e535d97..52d1863 100644 --- a/src/transformations.jl +++ b/src/transformations.jl @@ -2,19 +2,14 @@ ##### transformations ##### -export domain, domain_kind, transform_to, transform_from, coordinate_transformations, +export domain, domain_kind, ErrorOutsideDomain, NearestInDomain, + transform_to, transform_from, coordinate_transformations, BoundedLinear, InfRational, SemiInfRational #### #### generic api #### -""" - -!!! note - Abstract type used for code organization, not exported. -""" - """ $(TYPEDEF) @@ -22,8 +17,6 @@ An abstract type for univariate transformations. """ abstract type AbstractUnivariateTransformation end -Broadcast.broadcastable(transformation::AbstractUnivariateTransformation) = Ref(transformation) - domain_kind(::Type{<:AbstractUnivariateTransformation}) = :univariate """ @@ -51,92 +44,23 @@ values and derivatives. function transform_from end #### -#### coordinate transformations +#### handling out of bounds values #### -struct CoordinateTransformations{T<:Tuple} - transformations::T -end - -domain_kind(::Type{<:CoordinateTransformations}) = :multivariate - -function domain(coordinate_transformations::CoordinateTransformations) - coordinate_domains(map(domain, coordinate_transformations.transformations)) -end - -function Base.Tuple(coordinate_transformations::CoordinateTransformations) - coordinate_transformations.transformations -end - -function Base.show(io::IO, ct::CoordinateTransformations) - print(io, "coordinate transformations") - for t in ct.transformations - print(io, "\n ", t) - end -end - -Broadcast.broadcastable(ct::CoordinateTransformations) = Ref(ct) - """ -$(SIGNATURES) - -Wrapper for coordinate-wise transformations. To extract components, convert to Tuple. - -```jldoctest -julia> using StaticArrays - -julia> ct = coordinate_transformations(BoundedLinear(0, 2), SemiInfRational(2, 3)) -coordinate transformations - (0.0,2.0) ↔ domain [linear transformation] - (2,∞) ↔ domain [rational transformation with scale 3] - -julia> d1 = domain(Chebyshev(InteriorGrid(), 5)) -[-1,1] - -julia> dom = coordinate_domains(d1, d1) -[-1,1]² - -julia> x = transform_from(dom, ct, (0.4, 0.5)) -(1.4, 11.0) - -julia> y = transform_to(dom, ct, x) -(0.3999999999999999, 0.5) -``` +Error when the value is outside the domain. """ -function coordinate_transformations(transformations::Tuple) - CoordinateTransformations(transformations) -end - -coordinate_transformations(transformations...) = coordinate_transformations(transformations) - -function transform_to(domain::CoordinateDomains, ct::CoordinateTransformations, x::Tuple) - (; domains) = domain - (; transformations) = ct - @argcheck length(domains) == length(transformations) == length(x) - map((d, t, x) -> transform_to(d, t, x), domains, transformations, x) -end - -function transform_to(domain::CoordinateDomains{T}, ct::CoordinateTransformations, - x::AbstractVector) where T - SVector(transform_to(domain, ct, _ntuple_like(T, x))) -end - -function transform_to(domain::CoordinateDomains, ct::CoordinateTransformations, - Dx::∂CoordinateExpansion) - ∂CoordinateExpansion(Dx.∂D, transform_to(domain, ct, Dx.x)) -end +struct ErrorOutsideDomain end -function transform_from(domain::CoordinateDomains, ct::CoordinateTransformations, x::Tuple) - (; domains) = domain - (; transformations) = ct - @argcheck length(domains) == length(transformations) == length(x) - map((d, t, x) -> transform_from(d, t, x), domains, transformations, x) -end +""" +Use the nearest value in the domain. +""" +struct NearestInDomain end -function transform_from(domain::CoordinateDomains{T}, ct::CoordinateTransformations, - x::AbstractVector) where {T} - SVector(transform_from(domain, ct, _ntuple_like(T, x))) -end +""" +Default handling of values outside domain. Not part of the API. +""" +const OUTSIDE_DOMAIN = ErrorOutsideDomain() #### #### specific transformations @@ -146,57 +70,80 @@ end ### bounded linear ### -struct BoundedLinear{T <: Real} <: AbstractUnivariateTransformation - "Midpoint `m`." - m::T - "Scale `s`." - s::T - function BoundedLinear(a::T, b::T) where {T <: Real} - @argcheck isfinite(a) && isfinite(b) DomainError - s = (b - a) / 2 - m = (a + b) / 2 - @argcheck s > 0 DomainError((; a, b), "Need `a < b`.") - m, s = promote(m, s) - new{typeof(m)}(m, s) +struct BoundedLinear{T <: Real,O} <: AbstractUnivariateTransformation + "Lower limit." + lower::T + "Upper limit." + upper::T + "Handling values outside the domain." + outside_domain::O + @doc """ + $(SIGNATURES) + + Transform the domain to `y ∈ (lower, upper)`, using a linear mapping. + + `lower < upper` is enforced. + + `outside_domain` determines how values outside `[lower, upper]` are handled. See + [`ErrorOutsideDomain`](@ref) and [`NearestInDomain`](@ref). + """ + function BoundedLinear(lower::T, upper::T, + outside_domain::O = OUTSIDE_DOMAIN) where {T<:Real,O} + @argcheck isfinite(lower) && isfinite(upper) DomainError + @argcheck upper > lower DomainError((; lower, upper), "Need `lower < upper`.") + new{typeof(lower),O}(lower, upper, outside_domain) end end -function Base.show(io::IO, transformation::BoundedLinear) - (; m, s) = transformation - print(io, "(", m - s, ",", m + s, ") ↔ domain [linear transformation]") +function BoundedLinear(lower::Real, upper::Real, outside_domain = OUTSIDE_DOMAIN) + BoundedLinear(promote(lower, upper)..., outside_domain) end -""" -$(TYPEDEF) +function BoundedLinear(; lower, upper, outside_domain = OUTSIDE_DOMAIN()) + BoundedLinear(lower, upper, outside_domain) +end -Transform the domain to `y ∈ (a, b)`, using ``y = x ⋅ s + m``. +function Base.show(io::IO, transformation::BoundedLinear) + (; lower, upper, outside_domain) = transformation + print(io, "BoundedLinear(", lower, ", ", upper) + outside_domain ≠ OUTSIDE_DOMAIN && print(io, ", ", outside_domain) + print(io, ")") +end -`m` and `s` are calculated and checked by the constructor; `a < b` is enforced. -""" -BoundedLinear(a::Real, b::Real) = BoundedLinear(promote(a, b)...) +function _handle_domain(t::BoundedLinear{<:Real,ErrorOutsideDomain}, x::Real) + (; lower, upper) = t + @argcheck lower ≤ x ≤ upper DomainError(x, "outside domain") + x +end + +function _handle_domain(t::BoundedLinear{<:Real,NearestInDomain}, x::Real) + (; lower, upper) = t + clamp(x, lower, upper) +end function transform_from(::PM1, t::BoundedLinear, x::Scalar) - (; m, s) = t - x * s + m + (; lower, upper) = t + (x + 1) / 2 * (upper-lower) + lower end -function transform_to(::PM1, t::BoundedLinear, y::Real) - (; m, s) = t - (y - m) / s +function transform_to(domain::PM1, t::BoundedLinear, y::Real) + (; lower, upper, outside_domain) = t + (_handle_domain(t, y) - lower) / (upper-lower) * 2 - 1 end function transform_to(domain::PM1, t::BoundedLinear, y::𝑑Expansion{Dp1}) where Dp1 - (; m, s) = t + (; lower, upper) = t (; coefficients) = y y0, yD... = coefficients x0 = transform_to(domain, t, y0) + s = (upper - lower) / 2 xD = map(y -> y / s, yD) 𝑑Expansion(SVector(x0, xD...)) end function domain(t::BoundedLinear) - (; m, s) = t - UnivariateDomain(m - s, m + s) + (; lower, upper) = t + UnivariateDomain(lower, upper) end ### @@ -204,49 +151,43 @@ end ### struct SemiInfRational{T<:Real} <: AbstractUnivariateTransformation - "The finite endpoint `A`." - A::T - "Scale factor `L ≠ 0`." - L::T - function SemiInfRational(A::T, L::T) where {T <: Real} - @argcheck isfinite(A) DomainError - @argcheck isfinite(L) && L ≠ 0 DomainError - new{T}(A, L) + "The finite endpoint." + endpoint::T + "Scale factor." + scale::T + @doc """ + $(SIGNATURES) + + The domian transformed to `[endpoint, Inf)` (when `scale > 0`) or `(-Inf,endpoint]` + (when `scale < 0`) using ``y = endpoint + scale ⋅ (1 + x) / (1 - x)``. + + When used with Chebyshev polynomials, also known as a “rational Chebyshev” basis. + + # Example mappings for the domain ``(-1,1)`` + + - ``-1/2 ↦ endpoint + scale / 3`` + - ``0 ↦ endpoint + scale`` + - ``1/2 ↦ endpoint + 3 ⋅ scale`` + """ + function SemiInfRational(; endpoint::Real = 0, scale::Real = 1) + @argcheck isfinite(endpoint) DomainError + @argcheck isfinite(scale) && scale ≠ 0 DomainError + endpoint, scale = promote(endpoint, scale) + new{typeof(endpoint)}(endpoint, scale) end end function Base.show(io::IO, transformation::SemiInfRational) - (; A, L) = transformation - if L > 0 - D = "($A,∞)" - else - D = "(-∞,A)" - end - print(io, D, " ↔ domain [rational transformation with scale ", L, "]") + (; endpoint, scale) = transformation + print(io, "SemiInfRational(endpoint = ", endpoint, ", scale = ", scale, ")") end -""" -$(SIGNATURES) - -The domian transformed to `[A, Inf)` (when `L > 0`) or `(-Inf,A]` -(when `L < 0`) using ``y = A + L ⋅ (1 + x) / (1 - x)``. - -When used with Chebyshev polynomials, also known as a “rational Chebyshev” basis. - -# Example mappings for the domain ``(-1,1)`` - -- ``-1/2 ↦ A + L / 3`` -- ``0 ↦ A + L`` -- ``1/2 ↦ A + 3 ⋅ L`` -""" -SemiInfRational(A::Real, L::Real) = SemiInfRational(promote(A, L)...) - -transform_from(::PM1, t::SemiInfRational, x) = t.A + t.L * (1 + x) / (1 - x) +transform_from(::PM1, t::SemiInfRational, x) = t.endpoint + t.scale * (1 + x) / (1 - x) function transform_to(::PM1, t::SemiInfRational, y::Real) - (; A, L) = t - z = y - A - x = (z - L) / (z + L) + (; endpoint, scale) = t + z = y - endpoint + x = (z - scale) / (z + scale) if y == Inf || y == -Inf one(x) else @@ -255,22 +196,22 @@ function transform_to(::PM1, t::SemiInfRational, y::Real) end function transform_to(domain::PM1, t::SemiInfRational, y::𝑑Expansion{Dp1}) where Dp1 - (; A, L) = t + (; scale) = t (; 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 = (coefficients[2] * Q) / (2*L) + x1 = (coefficients[2] * Q) / (2*scale) Dp1 == 2 && return 𝑑Expansion(SVector(x0, x1)) error("$(Dp1-1)th derivative not implemented yet, open an issue.") end function domain(t::SemiInfRational) - (; L, A) = t - A = float(A) - ∞ = oftype(A, Inf) - L > 0 ? UnivariateDomain(A, ∞) : UnivariateDomain(-∞, A) + (; endpoint, scale) = t + endpoint = float(endpoint) + ∞ = oftype(endpoint, Inf) + scale > 0 ? UnivariateDomain(endpoint, ∞) : UnivariateDomain(-∞, endpoint) end ### @@ -278,40 +219,41 @@ end ### struct InfRational{T <: Real} <: AbstractUnivariateTransformation - "The center `A`." - A::T - "Scale factor `L > 0`." - L::T - function InfRational(A::T, L::T) where {T <: Real} - @argcheck isfinite(A) DomainError - @argcheck isfinite(L) && L > 0 DomainError - new{T}(A, L) + "The center" + center::T + "Scale factor" + scale::T + @doc """ + $(SIGNATURES) + + The domain transformed to `(-Inf, Inf)` using + ``y = center + scale ⋅ x / √(1 - x^2)``, with `scale > 0`. + + # Example mappings (for domain ``(-1,1)``) + + - ``0 ↦ center`` + - ``±0.5 ↦ center ± scale / √3`` + """ + function InfRational(center::T, scale::T) where {T <: Real} + @argcheck isfinite(center) DomainError + @argcheck isfinite(scale) && scale > 0 DomainError + new{T}(center, scale) end end function Base.show(io::IO, transformation::InfRational) - (; A, L) = transformation - print(io, "(-∞,∞) ↔ domain [rational transformation with center ", A, ", scale ", L, "]") + (; center, scale) = transformation + print(io, "InfRational(; center = ", center, ", scale = ", scale, ")") end -""" -$(SIGNATURES) - -The domain transformed to `(-Inf, Inf)` using ``y = A + L ⋅ x / √(1 - x^2)``, with `L > 0`. - -# Example mappings (for domain ``(-1,1)``) - -- ``0 ↦ A`` -- ``±0.5 ↦ A ± L / √3`` -""" -InfRational(A::Real, L::Real) = InfRational(promote(A, L)...) +InfRational(; center::Real = 0.0, scale::Real = 1.0) = InfRational(promote(center, scale)...) -transform_from(::PM1, T::InfRational, x::Real) = T.A + T.L * x / √(1 - abs2(x)) +transform_from(::PM1, T::InfRational, x::Real) = T.center + T.scale * x / √(1 - abs2(x)) function transform_to(::PM1, t::InfRational, y::Real) - (; A, L) = t - z = y - A - x = z / hypot(z, L) + (; center, scale) = t + z = y - center + x = z / hypot(z, scale) if isinf(y) y > 0 ? one(x) : -one(x) else @@ -320,14 +262,14 @@ function transform_to(::PM1, t::InfRational, y::Real) end function transform_to(domain::PM1, t::InfRational, y::𝑑Expansion{Dp1}) where Dp1 - (; A, L) = t + (; scale) = t (; coefficients) = y x0 = transform_to(domain, t, coefficients[1]) Dp1 == 1 && return SVector(x0) # based on Boyd (2001), Table E.5 Q = 1 - abs2(x0) sQ = √Q - x1 = (coefficients[2] * Q * sQ) / L + x1 = (coefficients[2] * Q * sQ) / scale Dp1 == 2 && return 𝑑Expansion(SVector(x0, x1)) error("$(Dp1-1)th derivative not implemented yet, open an issue.") end diff --git a/src/univariate.jl b/src/univariate.jl new file mode 100644 index 0000000..8d7dc0b --- /dev/null +++ b/src/univariate.jl @@ -0,0 +1,252 @@ +#### +#### Chebyshev polynomials on [-1,1] +#### + +export Chebyshev, Endpoints, Interior, UnivariateBasis + +#### +#### generic building blocks +#### + +""" +$(TYPEDEF) + +The first `N` Chebyhev polynomials of the first kind, defined on `[-1,1]`. +""" +struct Chebyshev end + +#### +#### basis function iterator +#### + +struct ChebyshevIterator{T} + x::T +end + +function _start(itr::ChebyshevIterator{T}) where T + (; x) = itr + _one(T), (_one(T), x) +end + +function _next(itr::ChebyshevIterator{T}, (fp, fpp)) where T + (; x) = itr + f = _sub(_mul(2, x, fp), fpp) + f::T, (f, fp) +end + +_eltype(::Type{<:ChebyshevIterator{T}}) where {T} = T + +""" +$(SIGNATURES) + +Helper function to calculate the extrema of the `N`th Chebyshev polynomial, indexed by +`1 ≤ i ≤ N` (not checked). Results are in `[-1,1]`. `N==1` (the constant) is special +cased to zero, for nesting. +""" +function _chebyshev_extremum(::Type{T}, i::Int, N::Int) where {T <: Real} + if N == 1 + cospi(1/T(2))::T # 0.0 as a practical fallback + else + cospi((N - i) ./ T(N - 1))::T + end +end + +#### +#### kinds +#### + +""" +$(TYPEDEF) + +Like [`Endpoints`](@ref), but with endpoints dropped. +""" +struct Interior end + +""" +$(SIGNATURES) + +Length of a univariate grid. +""" +function grid_length(::Chebyshev, ::Interior, level::Int) + @argcheck level ≥ 0 + (1 << (level + 1)) - 1 +end + +""" +$(SIGNATURES) + +Length of a single block, these are concatenated to form the grid. +""" +function block_length(::Chebyshev, ::Interior, level::Int) + @argcheck level ≥ 0 + 1 << level +end + +""" +$(SIGNATURES) + +Map `i` to an integer for calling [`_chebyshev_extremum_shuffle`](@ref). Interior +indices start from `2`, endpoint from `1`. Caller is responsible for making sure that +`i` is in the valid range `1:grid_length(Chebyshev(), kind, level)`, this is not checked. +""" +function _chebyshev_extremum_shuffle(kind::Interior, i::Int, level::Int) + p = ndigits(i, base = 2) - 1 # trust constant folding fast path to top_set_bit + remainder = i - (1 << p) + start = 1 << (level - p) + start + remainder * (start << 1) + 1 +end + +""" +$(TYPEDEF) + +Chebyshev-Lobatto grid. The extrema of Chebyshev polynomials, including endpoints of `[-1,1]`. + +!!! note + For small dimensions may fall back to a grid that does not contain endpoints. +""" +struct Endpoints end + +function grid_length(::Chebyshev, ::Endpoints, level::Int) + @argcheck level ≥ 0 + (level ≤ 1 ? level * 2 : (1 << level)) + 1 +end + +function block_length(::Chebyshev, ::Endpoints, level::Int) + @argcheck level ≥ 0 + level ≤ 1 ? (level + 1) : (1 << (level - 1)) +end + +function _chebyshev_extremum_shuffle(::Endpoints, i::Int, level::Int) + if i > 3 + _chebyshev_extremum_shuffle(Interior(), i - 2, level - 1) + elseif i == 1 + 1 + 1 << (level - 1) + else + 1 + (i - 2) * (1 << level) + end +end + +#### +#### univariate bases +#### + +""" +Implementation of univariate bases. Not part of the API. +""" +struct UnivariateBasis{F,K,D} <: FunctionBasis + family::F + kind::K + domain_transformation::D + level::Int + grid_level::Int + @doc """ + $(SIGNATURES) + + Univariate basis from `family`, using the given `kind`. + + `domain_transformation` determines the domain and the relevant transformation to the + domain of `family`. + + `level` is an integer, starting from `0`, specifying the number of *blocks* used to + build the grid. + """ + function UnivariateBasis(family::F, kind::K, domain_transformation::D, + level::Int, grid_level::Int = level) where {F,K,D} + @argcheck level ≥ 0 + new{F,K,D}(family, kind, domain_transformation, level, grid_level) + end +end + +function Base.show(io::IO, basis::UnivariateBasis) + (; family, kind, domain_transformation, level, grid_level) = basis + _print_dimensions(io, basis) + print(io, "UnivariateBasis(", family, ", ", kind, ", ", domain_transformation, ", ", + level) + level ≠ grid_level && print(io, ", ", grid_level) + print(io, ")") +end + +domain(U::UnivariateBasis) = domain(U.domain_transformation) + +dimension(U::UnivariateBasis) = grid_length(U.family, U.kind, U.level) + +struct UnivariateBasisAt{I} + infinite_itr::I + N::Int +end + +Base.eltype(::Type{UnivariateBasisAt{I}}) where I = _eltype(I) + +Base.length(itr::UnivariateBasisAt) = itr.N + +function Base.iterate(itr::UnivariateBasisAt, state = nothing) + (; infinite_itr, N) = itr + if state ≡ nothing + x, inner_state = _start(infinite_itr) + x, (1, inner_state) + else + i, inner_state = state + if i < N + x, inner_state′ = _next(infinite_itr, inner_state) + x, (i + 1, inner_state′) + else + nothing + end + end +end + +""" +$(SIGNATURES) + +Return an infinite iterator for the univariate basis functions, using the protocol with +[`_start`](@ref), [`_next`](@ref), etc. +""" +function _univariate_basis_itr(family::Chebyshev, domain_transformation, x::Scalar) + ChebyshevIterator(transform_to(PM1(), domain_transformation, x)) +end + +function basis_at(U::UnivariateBasis, x::Scalar) + UnivariateBasisAt(_univariate_basis_itr(U.family, U.domain_transformation, x), + dimension(U)) +end + +@concrete struct ChebyshevGrid{T} <: AbstractVector{T} + kind + domain_transformation + level::Int + N::Int + N̂::Int +end + +Base.size(g::ChebyshevGrid) = (g.N, ) + +function Base.getindex(g::ChebyshevGrid{T}, i::Int) where T + (; kind, domain_transformation, level, N̂) = g + transform_from(PM1(), domain_transformation, + _chebyshev_extremum(T, _chebyshev_extremum_shuffle(kind, i, level), N̂))::T +end + +function grid(::Type{T}, + U::UnivariateBasis{Chebyshev, + <:Union{Endpoints,Interior}}) where {T <: AbstractFloat} + (; family, kind, domain_transformation, grid_level) = U + N = grid_length(family, kind, grid_level) + ChebyshevGrid{T}(kind, domain_transformation, grid_level, N, + kind ≡ Interior() ? N̂ = N + 2 : N) +end + +function adjust_coefficients(θ1::AbstractVector, + U1::UnivariateBasis{F}, U2::UnivariateBasis{F}) where F + d1 = dimension(U1) + d2 = dimension(U2) + @argcheck length(θ1) == d1 "coefficients are not compatible with the first basis" + @argcheck(U1.domain_transformation == U2.domain_transformation, + "incompatible domain transformations") + if d2 ≤ d1 # truncate + θ1[1:d2] + else # pad with zeros + θ2 = zeros(d2) + θ2[1:d1] .= θ1 + θ2 + end +end diff --git a/src/utilities.jl b/src/utilities.jl index 80d3dbc..c1e82ec 100644 --- a/src/utilities.jl +++ b/src/utilities.jl @@ -1,3 +1,70 @@ +##### +##### utilities +##### + +#### +#### infinite iteration +#### + +""" +$(SIGNATURES) + +This package introduces its own API for infinite iteration. `_start(itr)` is not unlike +`iterate(itr)`, while `_next(itr, state) → x, state` is not unlike `iterate(itr, +state)`. Also see [`_eltype`](@ref). + +The rationale is to ease the compilation burden by ruling out `Union` types (`nothing`) +and combinatorial explosition. Julia can cope with it `iterate`, but it causes problems +with `Enzyme`. +""" +_start(itr) = iterate(itr)::Tuple + +""" +$(SIGNATURES) → x, state + +Internal API for infinite iteration. See [`_start`](@ref). +""" +_next(itr, state) = iterate(itr, state)::Tuple + +""" +$(SIGNATURES) → Type + +Internal API for infinite iteration. See [`_start`](@ref). +""" +_eltype(::Type{T}) where T = eltype(T) + +"Counting integers from 1." +struct Counting end + +_start(::Counting) = 1, 1 + +_next(::Counting, state) = state + 1, state + 1 + +_eltype(::Type{Counting}) = Int + +#### +#### printing +#### + +"Color for printing additional information." +const INFO_COLOR = :blue + +""" +$(SIGNATURES) + +Helper function to print the dimensions nicely, as a comment. Grid dimensions are only +printed when different from basis dimensions. +""" +function _print_dimensions(io::IO, basis) + d = dimension(basis) + gl = length(grid(basis)) + if !get(io, :compact, false) + printstyled(io, "# dimension: ", d; color = INFO_COLOR) + d ≠ gl && printstyled(io, ", grid length: ", gl; color = INFO_COLOR) + println(io) + end +end + const _SUPERSCRIPT_DIGITS = ['⁰', '¹', '²', '³', '⁴', '⁵', '⁶', '⁷', '⁸', '⁹'] const _SUBSCRIPT_DIGITS = ['₀', '₁', '₂', '₃', '₄', '₅', '₆', '₇', '₈', '₉'] @@ -53,6 +120,10 @@ Base.print(io::IO, s::SubScript) = print_number(io, _SUBSCRIPT_DIGITS, s.i) # end # end +#### +#### conversions +#### + """ $(SIGNATURES) diff --git a/test/Project.toml b/test/Project.toml index 700d5c6..fc7cc57 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -1,7 +1,10 @@ [deps] +Accessors = "7d9f7c33-5ae7-4f3b-8dc6-eff91059b697" Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595" BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf" +ConstructionBase = "187b0558-2788-49d3-abe0-74a17ed4e7c9" DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" +Enzyme = "7da242da-08ed-463a-9acd-ee780be4f1d9" FiniteDifferences = "26cc04aa-876d-5657-8c51-4c34ba976000" InverseFunctions = "3587e190-3f89-42d0-90ee-14403ec27112" JET = "c3a54625-cd67-489e-a8e7-0a5a0ff4e31b" @@ -14,7 +17,6 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [compat] Aqua = "0.8" -BenchmarkTools = "1" DocStringExtensions = "0.9" FiniteDifferences = "0.12" InverseFunctions = "0.1" diff --git a/test/integration_enzyme.jl b/test/integration_enzyme.jl new file mode 100644 index 0000000..40bc82e --- /dev/null +++ b/test/integration_enzyme.jl @@ -0,0 +1,16 @@ +##### +##### integration test for Enzyme +##### + +using Enzyme + +@testset "Enzyme integration" begin + basis = SmolyakBasis(Chebyshev(), Interior(), + ntuple(_ -> BoundedLinear(1.0, 2.0), Val(2)), + SmolyakLevel(total = 2, each = 2)) + g(basis, θ) = linear_combination(basis, θ, (1.5, 1.5)) + d = dimension(basis) + x = zeros(d) + dx = zeros(d) + autodiff(Reverse, g, Const(basis), Duplicated(x, dx)) +end diff --git a/test/runtests.jl b/test/runtests.jl index d119afe..80a735f 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,18 +1,20 @@ using SpectralKit, Test +import SpectralKit as SK -using JET -@testset "static analysis with JET.jl" begin - @test isempty(JET.get_reports(report_package(SpectralKit, - target_modules=(SpectralKit,), - ignored_modules = (SpectralKit.Experimental,)))) -end +# using JET +# @testset "static analysis with JET.jl" begin +# @test isempty(JET.get_reports(report_package(SpectralKit, +# target_modules=(SpectralKit,), +# #= ignored_modules = (SpectralKit.Experimental,) =#) +# )) +# end -@testset "QA with Aqua" begin - import Aqua - Aqua.test_all(SpectralKit) -end +# @testset "QA with Aqua" begin +# import Aqua +# Aqua.test_all(SpectralKit) +# end -using DocStringExtensions, StaticArrays, BenchmarkTools, FiniteDifferences +using Accessors, DocStringExtensions, StaticArrays, BenchmarkTools, FiniteDifferences include("utilities.jl") @@ -20,8 +22,9 @@ include("test_utilities.jl") include("test_derivatives.jl") include("test_domains.jl") include("test_transformations.jl") -include("test_chebyshev.jl") +include("test_univariate.jl") include("test_smolyak_traversal.jl") include("test_smolyak.jl") -include("test_generic_api.jl") # NOTE moved last as it used constructs from above -include("test_experimental.jl") # NOTE experimental code is not public API +include("integration_enzyme.jl") +# include("test_generic_api.jl") # NOTE moved last as it used constructs from above +# include("test_experimental.jl") # NOTE experimental code is not public API diff --git a/test/test_chebyshev.jl b/test/test_chebyshev.jl deleted file mode 100644 index 5c7b7a9..0000000 --- a/test/test_chebyshev.jl +++ /dev/null @@ -1,146 +0,0 @@ -#### -#### Chebyshev on [-1,1] -#### - -@testset "Chebyshev" begin - @test_throws ArgumentError Chebyshev(InteriorGrid(), 0) - @test_throws ArgumentError Chebyshev(EndpointGrid(), 0) - @test_throws TypeError Chebyshev(:invalid_grid, 10) - - for grid_kind in (InteriorGrid(), EndpointGrid()) - for N in (grid_kind ≡ InteriorGrid() ? 1 : 2):10 - - basis = Chebyshev(grid_kind, N) - @test is_function_basis(basis) - @test is_function_basis(typeof(basis)) - @test dimension(basis) == N - - # check linear combinations - for _ in 1:100 - x = rand_in_domain(basis) - bx = @inferred basis_at(basis, x) - - @test length(bx) == N - @test eltype(bx) == Float64 - @test collect(bx) ≈ [chebyshev_cos(x, i) for i in 1:N] - - θ = rand(N) - @test linear_combination(basis, θ, x) ≈ - sum(chebyshev_cos(x, i) * θ for (i,θ) in enumerate(θ)) - @test linear_combination(basis, θ, 𝑑(x))[1] ≈ - sum(chebyshev_cos_deriv(x, i) * θ for (i,θ) in enumerate(θ)) - end - - # check grid - g = @inferred collect(grid(basis)) - @test length(g) == N - if grid_kind ≡ InteriorGrid() - @test all(x -> isapprox(chebyshev_cos(x, N + 1), 0, atol = 1e-14), g) - else - @test all(x -> isapprox(chebyshev_cos_deriv(x, N), 0, atol = 1e-13), - g[2:(end-1)]) - @test g[1] == -1 - @test g[end] == 1 - end - - # augmented coefficients - for _ in 1:100 - x = rand_in_domain(basis) - θ = rand(N) - destination_basis = Chebyshev(grid_kind, N + 5) - destination_θ = augment_coefficients(basis, destination_basis, θ) - @test linear_combination(basis, θ, x) ≈ - linear_combination(destination_basis, destination_θ, x) - end - - end - end - - # compatible and incompatible grids - @testset "augment Chebyshev coefficients — errors" begin - basis = Chebyshev(InteriorGrid(), 5) - θ = randn(5) - # different grids are compatible - basis2_G = Chebyshev(EndpointGrid(), 6) - @test is_subset_basis(basis, basis2_G) - # fewer dimensions are not compatible - basis2_N = Chebyshev(InteriorGrid(), 4) - @test !is_subset_basis(basis, basis2_N) - @test_throws ArgumentError augment_coefficients(basis, basis2_N, θ) - # too few coefficients - @test_throws ArgumentError augment_coefficients(basis, basis, randn(4)) - end -end - -@testset "augmentation of transformed basis" begin - N = 5 - M = N + 3 - t = SemiInfRational(0.3, 0.9) - grid_kind = InteriorGrid() - basis = Chebyshev(grid_kind, N) ∘ t - basis′ = Chebyshev(grid_kind, M) ∘ t - @test is_subset_basis(basis, basis′) - for _ in 1:100 - x = rand_in_domain(basis) - θ = rand(N) - θ′ = augment_coefficients(basis, basis′, θ) - @test linear_combination(basis, θ, x) ≈ linear_combination(basis′, θ′, x) - end -end - -@testset "univariate derivatives" begin - basis = Chebyshev(InteriorGrid(), 5) - 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:50 - 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 - 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_smolyak.jl b/test/test_smolyak.jl index c1863a1..17f06f6 100644 --- a/test/test_smolyak.jl +++ b/test/test_smolyak.jl @@ -1,19 +1,48 @@ -using SpectralKit, Test -using SpectralKit: PM1, ∂Expansion - #### #### api #### -@testset "Smolyak API checks" begin - @test_throws MethodError smolyak_basis(Chebyshev, :invalid_grid, SmolyakParameters(3), 2) - @test_logs (:warn, "M > B replaced with M = B") SmolyakParameters(2, 4) +using ConstructionBase + +@testset "SmolyakLevels printing" begin + @test repr(SmolyakLevel(; total = 3, each = 2)) == + "SmolyakLevel(total = 3, each = 2) #= 0 ≤ ∑ℓᵢ ≤ 3, all 0 ≤ ℓᵢ ≤ 2 =#" +end + +@testset "SmolyakLevel normalization" begin + @test_logs (:warn, "‘each’ normalized to ‘total’") SmolyakLevel(total = 2, each = 4) +end + +@testset "SmolyakLevel ConstructionBase integration" begin + level = SmolyakLevel(; total = 4, each = 3) + @test ConstructionBase.setproperties(level, (total = 5,)) == + SmolyakLevel(; total = 5, each = 3) + @test ConstructionBase.setproperties(level, (each = 2, total = 5)) == + SmolyakLevel(; total = 5, each = 2) + @test_throws MethodError ConstructionBase.setproperties(level, (a_fish = 3,)) +end + +@testset "Smolyak printing" begin + basis = SmolyakBasis(Chebyshev(), Interior(), (BoundedLinear(0, 1),), SmolyakLevel(total = 3, each = 2)) + @test repr(basis) == + """# dimension: 7 +SmolyakBasis(Chebyshev(), Interior(), + (BoundedLinear(0, 1),), + SmolyakLevel(total = 3, each = 2) #= 0 ≤ ∑ℓᵢ ≤ 3, all 0 ≤ ℓᵢ ≤ 2 =#)""" + @test repr(@set basis.grid_level = SmolyakLevel(total = 4, each = 3)) == + """# dimension: 7, grid length: 15 +SmolyakBasis(Chebyshev(), Interior(), + (BoundedLinear(0, 1),), + SmolyakLevel(total = 3, each = 2) #= 0 ≤ ∑ℓᵢ ≤ 3, all 0 ≤ ℓᵢ ≤ 2 =#, + SmolyakLevel(total = 4, each = 3) #= 0 ≤ ∑ℓᵢ ≤ 4, all 0 ≤ ℓᵢ ≤ 3 =#)""" end @testset "Smolyak API sanity checks" begin f(x) = (x[1] - 3) * (x[2] + 5) # linear function, just a sanity check - basis = smolyak_basis(Chebyshev, InteriorGrid(), SmolyakParameters(3), 2) - @test @inferred(domain(basis)) ≡ coordinate_domains(PM1(), PM1()) + transformations = (BoundedLinear(2, 3), # approximation should be exact + BoundedLinear(3.0, 4.5)) + basis = SmolyakBasis(Chebyshev(), Interior(), transformations, SmolyakLevel(total = 2)) + @test @inferred(domain(basis)) ≡ domain.(transformations) g = grid(Float64, basis) iterator_sanity_checks(g) x = @inferred collect(g) @@ -25,72 +54,97 @@ end for y1 in y1 for y2 in y2 y = SVector(y1, y2) - @test linear_combination(basis, θ, y) ≈ f(y) + @test linear_combination(basis, θ, y) ≈ f(y) atol = 1e-14 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 + @test basis[1] == UnivariateBasis(Chebyshev(), Interior(), transformations[1], 2) + + # FIXME re-enable once we have derivatives + # @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 - basis = smolyak_basis(Chebyshev, InteriorGrid(), SmolyakParameters(3), 2) + t = SemiInfRational() + basis = SmolyakBasis(Chebyshev(), Interior(), (t, t), + SmolyakLevel(total = 3)) y = SVector(0.4, 0.7) θ = randn(dimension(basis)) @inferred linear_combination(basis, θ, y) @test @ballocated(linear_combination($basis, $θ, $y)) == 0 end -### -### augment coefficients -### - -@testset "Smolyak augment coefficients" begin - basis1 = smolyak_basis(Chebyshev, InteriorGrid(), SmolyakParameters(2, 2), 2) - θ1 = randn(dimension(basis1)) - - # grid ≠ - basis2_G = smolyak_basis(Chebyshev, EndpointGrid(), SmolyakParameters(2, 3), 2) - @test !is_subset_basis(basis1, basis2_G) - @test_throws ArgumentError augment_coefficients(basis1, basis2_G, θ1) - - # smolyak_parameters < - basis2_P = smolyak_basis(Chebyshev, InteriorGrid(), SmolyakParameters(2, 1), 2) - @test !is_subset_basis(basis1, basis2_P) - @test_throws ArgumentError augment_coefficients(basis1, basis2_P, θ1) - - # compatible basis - basis2 = smolyak_basis(Chebyshev, InteriorGrid(), SmolyakParameters(3, 2), 2) - θ2 = @inferred augment_coefficients(basis1, basis2, θ1) - @test length(θ2) == dimension(basis2) - @test eltype(θ2) == eltype(θ1) - for _ in 1:100 - x = (rand(), rand()) .* 4 - @test linear_combination(basis1, θ1, x) ≈ linear_combination(basis2, θ2, x) - end -end - -@testset "Smolyak nesting" begin - for grid_kind in GRIDS - for M1 in 0:5 - for M2 in (M1 + 1):5 - for B1 in 0:M1 - for B2 in (B1 + 1):M2 - basis1 = smolyak_basis(Chebyshev, grid_kind, SmolyakParameters(B1, M1), 2) - basis2 = smolyak_basis(Chebyshev, grid_kind, SmolyakParameters(B2, M2), 2) - @test is_approximate_subset(collect(grid(basis1)), collect(grid(basis2))) - end +@testset "smolyak indices" begin + for kind in KINDS + for N in 1:5 + for total in 0:4 + for each in 0:total + expected = naive_smolyak_indices(Chebyshev(), kind, Val(N), total, each) + basis = SmolyakBasis(Chebyshev(), kind, ntuple(_ -> nothing, Val(N)), + SmolyakLevel(; total, each)) + @test collect(SpectralKit.SmolyakIndices(basis)) == expected end end end end end + +@testset "Smolyak adjusted coefficients" begin + t = (BoundedLinear(2.0, 3.7), SemiInfRational()) + basis1 = SmolyakBasis(Chebyshev(), Interior(), t, SmolyakLevel(; total = 2, each = 1)) + basis2 = @set basis1.level = SmolyakLevel(; total = 3, each = 2) + θ1 = randn(dimension(basis1)) + θ2 = adjust_coefficients(θ1, basis1, basis2) + @test sum(θ1 .≠ 0) == sum(θ2 .≠ 0) +end + +# @testset "Smolyak augment coefficients" begin +# basis1 = SmolyakBasis(Chebyshev, InteriorGrid(), SmolyakParameters(2, 2), 2) +# θ1 = randn(dimension(basis1)) + +# # grid ≠ +# basis2_G = SmolyakBasis(Chebyshev, EndpointGrid(), SmolyakParameters(2, 3), 2) +# @test !is_subset_basis(basis1, basis2_G) +# @test_throws ArgumentError augment_coefficients(basis1, basis2_G, θ1) + +# # smolyak_parameters < +# basis2_P = SmolyakBasis(Chebyshev, InteriorGrid(), SmolyakParameters(2, 1), 2) +# @test !is_subset_basis(basis1, basis2_P) +# @test_throws ArgumentError augment_coefficients(basis1, basis2_P, θ1) + +# # compatible basis +# basis2 = SmolyakBasis(Chebyshev, InteriorGrid(), SmolyakParameters(3, 2), 2) +# θ2 = @inferred augment_coefficients(basis1, basis2, θ1) +# @test length(θ2) == dimension(basis2) +# @test eltype(θ2) == eltype(θ1) +# for _ in 1:100 +# x = (rand(), rand()) .* 4 +# @test linear_combination(basis1, θ1, x) ≈ linear_combination(basis2, θ2, x) +# end +# end + +# @testset "Smolyak nesting" begin +# for grid_kind in GRIDS +# for M1 in 0:5 +# for M2 in (M1 + 1):5 +# for B1 in 0:M1 +# for B2 in (B1 + 1):M2 +# basis1 = SmolyakBasis(Chebyshev, grid_kind, SmolyakParameters(B1, M1), 2) +# basis2 = SmolyakBasis(Chebyshev, grid_kind, SmolyakParameters(B2, M2), 2) +# @test is_approximate_subset(collect(grid(basis1)), collect(grid(basis2))) +# end +# end +# end +# end +# end +# end diff --git a/test/test_smolyak_traversal.jl b/test/test_smolyak_traversal.jl index ef28c0e..b262290 100644 --- a/test/test_smolyak_traversal.jl +++ b/test/test_smolyak_traversal.jl @@ -1,125 +1,38 @@ -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)) +@testset "smolyak step" begin + XS = (2, 3, 5, 7, 11) + for kind in KINDS + for N in 1:5 + for total in 1:4 + for each in 1:4 + xs = XS[1:N] + test_smolyak_step(Chebyshev(), Interior(), total, each, + (a, b) -> isempty(b) ? (a,) : (a * first(b), b...), + map(SpectralKit.ChebyshevIterator, xs)) + end 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 +@testset "nonincreasing smolyak levels" begin + for N in 1:3 + for total in 0:4 + for each in 0:total + N, total, each = 2, 3, 3 + itr = SpectralKit.NonIncreasingSmolyakLevels{N}(total, each) + @test sort(collect(itr)) == sort(nonincreasing_smolyak_levels(Val(N), total, each)) 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 +@testset "__smolyak_length" begin + for kind in KINDS + for N in 1:5 + for total in 0:4 + for each in 0:total + @test @inferred(SpectralKit.__smolyak_length(Chebyshev(), kind, Val(N), total, each)) == + length(naive_smolyak_indices(Chebyshev(), kind, Val(N), total, each)) end end end diff --git a/test/test_transformations.jl b/test/test_transformations.jl index b37d22d..0ee70e2 100644 --- a/test/test_transformations.jl +++ b/test/test_transformations.jl @@ -24,25 +24,25 @@ using SpectralKit: PM1 end @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) + @test_throws DomainError SemiInfRational(; scale = Inf) + @test_throws DomainError SemiInfRational(; scale = 0.0) + @test_throws DomainError SemiInfRational(; endpoint = NaN) - A = 3.0 - L = 4.0 - trans = SemiInfRational(A, L) + endpoint = 3.0 + scale = 4.0 + trans = SemiInfRational(; endpoint, scale) - @test extrema(domain(trans)) == (A, Inf) + @test extrema(domain(trans)) == (endpoint, Inf) for _ in 1:100 x = rand_pm1() y = transform_from(PM1(), trans, x) if x == -1 - @test y ≈ A + @test y ≈ endpoint elseif x == 1 @test y ≈ Inf else - @test A < y < Inf + @test endpoint < y < Inf end @test transform_to(PM1(), trans, y) ≈ x end @@ -58,14 +58,12 @@ end end @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) - @test_throws DomainError InfRational(NaN, 1) + @test_throws DomainError InfRational(; scale = Inf) + @test_throws DomainError InfRational(; scale = 0.0) + @test_throws DomainError InfRational(; scale = -2.0) + @test_throws DomainError InfRational(; center = NaN) - A = 0.0 - L = 1.0 - trans = InfRational(A, L) + trans = InfRational() @test extrema(domain(trans)) == (-Inf, Inf) @@ -92,36 +90,11 @@ end @test y_minf[1] == 0 end -@testset "coordinate transformations" begin - t1 = BoundedLinear(2.0, 3.0) - t2 = SemiInfRational(7.0, 1.0) - ct = coordinate_transformations(t1, t2) - md = coordinate_domains(Val(2), PM1()) - x = SVector(rand_pm1(), rand_pm1()) - y = @inferred transform_from(md, ct, x) - @test y isa SVector{2,Float64} - @test y == transform_from.(PM1(), Tuple(ct), x) - - # handle generic inputs - y2 = @inferred transform_from(md, ct, Vector(x)) - @test y2 isa SVector{2,Float64} && y2 == y - - x2 = @inferred transform_to(md, ct, [y...]) - @test x2 isa SVector{2,Float64} && all(x2 .≈ x) -end - -@testset "printing, promotion, broadcasting" begin - v = [1.0, 2.0] +@testset "printing, promotion" begin t1 = BoundedLinear(2.0, 3) - @test repr(t1) == "(2.0,3.0) ↔ domain [linear transformation]" - @test transform_to.(PM1(), t1, v) isa Vector - t2 = SemiInfRational(7.0, 1) - @test repr(t2) == "(7.0,∞) ↔ domain [rational transformation with scale 1.0]" - @test transform_to.(PM1(), t2, v) isa Vector - t3 = InfRational(0.5, 1) - @test repr(t3) == "(-∞,∞) ↔ domain [rational transformation with center 0.5, scale 1.0]" - @test transform_to.(PM1(), t3, v) isa Vector - ct = coordinate_transformations(t1, t2, t3) - @test repr(ct) == - "coordinate transformations\n " * repr(t1) * "\n " * repr(t2) * "\n " * repr(t3) + @test repr(t1) == "BoundedLinear(2.0, 3.0)" + t2 = SemiInfRational(; endpoint = 7.0) + @test repr(t2) == "SemiInfRational(endpoint = 7.0, scale = 1.0)" + t3 = InfRational(; center = 0.5, scale = 1) + @test repr(t3) == "InfRational(; center = 0.5, scale = 1.0)" end diff --git a/test/test_univariate.jl b/test/test_univariate.jl new file mode 100644 index 0000000..f2f2626 --- /dev/null +++ b/test/test_univariate.jl @@ -0,0 +1,201 @@ +#### +#### Chebyshev building blocks and univariate bases +#### + +""" +For testing Chebyshev blocks and shuffle. +""" +function _shuffle(len::Int; endpoints = true) + if !endpoints + len += 2 # for Interior, we start the algorithm with two extra points + end + if len == 1 + return [[1]] + end + shuffle = [[len ÷ 2 + 1], [1, len]] + while sum(length, shuffle) < len + used = sort!(reduce(vcat, shuffle)) + s = Int[] + for i in 1:(length(used)-1) + a = used[i] + b = used[i + 1] + c = (a + b) ÷ 2 + @assert a < c < b + push!(s, c) + end + push!(shuffle, s) + end + if !endpoints + deleteat!(shuffle, 2) + end + shuffle +end + +@testset "blocks and shuffle" begin + for kind in KINDS + for level in 0:5 + N = SK.grid_length(Chebyshev(), kind, level) + block_lengths = map(level -> SK.block_length(Chebyshev(), kind, level), 0:level) + @test sum(block_lengths) == N + S = _shuffle(N; endpoints = kind == Endpoints()) + @test block_lengths == map(length, S) + @test map(i -> SK._chebyshev_extremum_shuffle(kind, i, level), 1:N) == foldl(vcat, S) + end + end +end + +@testset "Chebyshev grid" begin + for kind in KINDS + previous = Float64[] + for level in 0:7 + b = UnivariateBasis(Chebyshev(), kind, BoundedLinear(1.0, 3.0), level) + g = collect(grid(b)) + @test is_approximate_subset(previous, g) + g = previous + end + end +end + +@testset "Chebyshev basics" begin + transformation = BoundedLinear(1.0, 3.0) + @test_throws ArgumentError UnivariateBasis(Chebyshev(), Endpoints(), transformation, -1) + for kind in KINDS + for level in 0:5 + basis = UnivariateBasis(Chebyshev(), kind, transformation, level) + @test is_function_basis(basis) + @test is_function_basis(typeof(basis)) + N = @inferred dimension(basis) + @test N ≥ 1 + + # check linear combinations + for _ in 1:100 + x = rand_in_domain(basis) + bx = @inferred basis_at(basis, x) + + @test length(bx) == N + @test eltype(bx) == Float64 + y = transform_to(PM1(), transformation, x) + @test collect(bx) ≈ [chebyshev_cos(y, i) for i in 1:N] + + θ = rand(N) + @test linear_combination(basis, θ, x) ≈ + sum(chebyshev_cos(y, i) * θ for (i,θ) in enumerate(θ)) + @test linear_combination(basis, θ, 𝑑(x))[1] ≈ + sum(chebyshev_cos_deriv(y, i) * θ for (i,θ) in enumerate(θ)) + end + + # check grid + g = @inferred collect(grid(basis)) + @test length(g) == N + a, b = extrema(domain(basis)) + @test all(a .≤ g .≤ b) + if kind ≡ Interior() + @test all(a .< g .< b) + end + @test all(x -> is_chebyshev_extrema(transform_to(PM1(), transformation, x), + kind ≡ Interior() ? N + 2 : N), g) + end + end +end + +@testset "Chebysev adjusted coefficients" begin + transformation = SemiInfRational(; endpoint = 3.0, scale = 7.0) + for kind in KINDS + level0 = 3 + basis0 = UnivariateBasis(Chebyshev(), kind, transformation, level0) + θ0 = randn(dimension(basis0)) + for Δ in 1:4 + basis = @set basis0.level = level0 + Δ + θ = adjust_coefficients(θ0, basis0, basis) + @test basis.level == basis0.level + Δ + for _ in 1:10 + x = rand_in_domain(basis0) + @test linear_combination(basis0, θ0, x) ≈ linear_combination(basis, θ, x) + end + end + end +end + +@testset "printing" begin + basis = UnivariateBasis(Chebyshev(), Interior(), SemiInfRational(), 3) + @test repr(basis) == + """# dimension: 15 +UnivariateBasis(Chebyshev(), Interior(), SemiInfRational(endpoint = 0, scale = 1), 3)""" + @test repr(@set basis.grid_level = 4) == + """# dimension: 15, grid length: 31 +UnivariateBasis(Chebyshev(), Interior(), SemiInfRational(endpoint = 0, scale = 1), 3, 4)""" +end + +@testset "different grid_level" begin + basis = UnivariateBasis(Chebyshev(), Interior(), BoundedLinear(0.0, 1.0), 2, 3) + @test length(grid(basis)) == + SpectralKit.grid_length(Chebyshev(), Interior(), 3) == + length(collect(grid(basis))) +end + +@testset "outside domain handling" begin + basis = UnivariateBasis(Chebyshev(), Interior(), BoundedLinear(0.0, 1.0), 4) + θ = randn(dimension(basis)) + @test_throws DomainError linear_combination(basis, θ, -0.5) + basis = UnivariateBasis(Chebyshev(), Interior(), BoundedLinear(0.0, 1.0, NearestInDomain()), 4) + @test linear_combination(basis, θ, -0.5) == linear_combination(basis, θ, 0.0) +end + +@testset "univariate derivatives" begin + for (transformation, N) in ((BoundedLinear(-2, 3), 5), + (SemiInfRational(endpoint = 0.7, scale = 0.3), 1), + (InfRational(center = 0.4, scale = 0.9), 1)) + basis = UnivariateBasis(Chebyshev(), Interior(), transformation, 3) + dom = domain(basis) + D = 𝑑^Val(N) + f = linear_combination(basis, randn(dimension(basis))) + for _ in 1:50 + x = rand_in_domain(basis) + y = f(D(x)) + for i in 0:N + @test y[i] ≈ DD(f, x, i; domain = dom) rtol = 1e-2 atol = 1e-2 + end + end + end +end + +@testset "endpoint continuity for derivatives" begin + # 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 + basis = UnivariateBasis(Chebyshev(), Interior(), + SemiInfRational(; endpoint = 2.3, scale = 0.7), 3) + for i in 1:dimension(basis) + θ = e_i(basis, i) + y_pinf = @inferred linear_combination(basis, θ, x_pinf) + @test y_pinf[0] == 1 + @test y_pinf[1] == 0 + end + end + + @testset "InfRational endpoints continuity" begin + basis = UnivariateBasis(Chebyshev(), Interior(), + InfRational(; center = 2.3, scale = 0.7), 3) + for i in 1:dimension(basis) + θ = e_i(basis, i) + y_pinf = @inferred linear_combination(basis, θ, x_pinf) + @test y_pinf[0] == 1 + @test y_pinf[1] == 0 + y_minf = @inferred linear_combination(basis, θ, x_minf) + @test y_minf[0] == (-1)^(i+1) + @test y_minf[1] == 0 + end + end +end + +@testset "Chebyshev allocations" begin + basis = UnivariateBasis(Chebyshev(), Interior(), SemiInfRational(), 3) + y = 0.5 + θ = randn(dimension(basis)) + @inferred linear_combination(basis, θ, y) + @test @ballocated(linear_combination($basis, $θ, $y)) == 0 +end diff --git a/test/utilities.jl b/test/utilities.jl index babf0d4..d8c8714 100644 --- a/test/utilities.jl +++ b/test/utilities.jl @@ -2,10 +2,20 @@ ##### utility functions for tests ##### -using SpectralKit: TransformedBasis, SmolyakBasis, SmolyakIndices # dispatch for rand_in_domain +"kinds we test on" +const KINDS = [Endpoints(), Interior()] + +# FIXME reenable what is needed +# using SpectralKit: TransformedBasis, SmolyakBasis, SmolyakIndices # dispatch for rand_in_domain chebyshev_cos(x, n) = cos((n - 1) * acos(x)) +""" +$(SIGNATURES) + +Derivative of the `n`th Chebyshev polynomial at `x`, using the cosine formula (special +cased at ±1). +""" function chebyshev_cos_deriv(x, n) z = cos(zero(x)) * abs2(n - 1) if x == -1 @@ -21,6 +31,20 @@ end """ $(SIGNATURES) +Test if `x` is an extrema of the `n`th Chebyshev polynomial, by checking derivatives to +be within tolerance (±1 special cased). +""" +function is_chebyshev_extrema(x, n; tol = 1e-8) + if abs(abs(x) - 1) ≤ tol + true + else + abs(chebyshev_cos_deriv(x, n)) ≤ tol + end +end + +""" +$(SIGNATURES) + Return a random value in [-1,1], putting an atomic mass on endpoints. The intention is to provide comprehensive testing for endpoints. @@ -34,27 +58,38 @@ Return a random value in the domain of the given basis, putting an atomic mass o 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))) +function rand_in_domain(basis::SpectralKit.UnivariateBasis{Chebyshev}) + transform_from(PM1(), basis.domain_transformation, rand_pm1()) end -function rand_in_domain(basis::TransformedBasis) - (; parent, transformation) = basis - transform_from(parent, transformation, rand_in_domain(parent)) -end +# function rand_in_domain(basis::SmolyakBasis{<:SmolyakIndices{N}}) where N +# (; univariate_parent) = basis +# SVector(ntuple(_ -> rand_in_domain(univariate_parent), Val(N))) +# end -"Flags (`true`) for elements in `a` that are within `atol` of some element in `b`." +# FIXME remove +# function rand_in_domain(basis::TransformedBasis) +# (; parent, transformation) = basis +# transform_from(parent, transformation, rand_in_domain(parent)) +# end + +""" +$(SIGNATURES) + +Flags (`true`) for elements in `a` that are within `atol` of some element in `b`. +""" function is_approximately_in(a, b; atol = √eps()) _same(a::Real, b::Real) = a == b || abs(a - b) ≤ atol # Inf = Inf, etc _same(a::AbstractVector, b::AbstractVector) = all(_same.(a, b)) - _same(a::Tuple, b::Tuple) = mapreduce((x, y) -> abs(x - y), max, a, b) ≤ atol + _same(a::Tuple, b::Tuple) = mapreduce((x, y) -> abs(x - y), max, a, b; init = 0.0) ≤ atol map(a -> any(b -> _same(a, b), b), a) end -"Are elements in `a` in `b`, approximately." +""" +$(SIGNATURES) + +Are elements in `a` in `b`, approximately. +""" function is_approximate_subset(a, b; atol = √eps()) sum(is_approximately_in(a, b; atol = atol)) == length(a) end @@ -80,11 +115,110 @@ function iterator_sanity_checks(itr) @test count(_ -> true, itr) == length(itr) end -"nth derivative of f at x." -function DD(f, x, n = 1; p = 10) +""" +$(SIGNATURES) + +`n`th derivative of f at x. + +`p` gives the degree of the finite difference rule. `domain` +""" +function DD(f, x, n = 1; p = 10, domain = (-Inf,Inf)) + A, B = extrema(domain) + a = x - A + b = B - x if n == 0 f(x) + elseif a == 0 + forward_fdm(p, n)(f, x) + elseif b == 0 + backward_fdm(p, n)(f, x) else - central_fdm(p, n)(f, x) + r = min(a, b) + if isfinite(r) + central_fdm(p, n, max_range = r)(f, x) + else + central_fdm(p, n)(f, x) + end + end +end + +""" +$(SIGNATURES) + +A vector of Smolyak indices traversed in column-major order. For testing. +""" +function naive_smolyak_indices(family, kind, ::Val{N}, total, each) where N + # map indices to levels + I = SpectralKit.grid_length(family, kind, total) + index_to_level = zeros(I) + i = 0 + for level in 0:total + b = SpectralKit.block_length(family, kind, level) + index_to_level[(i+1):(i+b)] .= level + i += b + end + # collect indices + indices = Vector{NTuple{N,Int}}() + for ι in Iterators.product(ntuple(_ -> 1:I, Val(N))...) + levels = map(i -> index_to_level[i], ι) + if all(levels .≤ each) && sum(levels) ≤ total + push!(indices, ι) + end + end + indices +end + +""" +$(SIGNATURES) + +Get the `n`th element from an infinite iterator. For unit testing. +""" +function _nth(infinite_itr, i) + x, state = SpectralKit._start(infinite_itr) + for _ in 2:i + x, state = SpectralKit._next(infinite_itr, state) + end + x +end + +""" +$(SIGNATURES) + +Test the Smolyak iterator implementation building blocks directly using a parallel, +naive calculation and checking invariants. +""" +function test_smolyak_step(family, kind, total, each, f, itrs::NTuple{N,Any}) where N + i = 1 + reference = naive_smolyak_indices(family, kind, Val(N), total, each) + (accum, slack, remainders, states, levels) = SpectralKit.__smolyak_init(family, kind, + total, f, itrs) + while true + @test slack + sum(levels) == total # simple sanity check + cs = map(_nth, itrs, reference[i]) + @test accum == foldr(f, cs; init = ()) + next = SpectralKit.__smolyak_step(family, kind, each, f, itrs, + accum, slack, remainders, states, levels) + (accum, Δ, remainders, states, levels) = next + Δ == SpectralKit.Δ_DONE && break + slack += Δ + i += 1 + end + @test i == length(reference) # we used up all elements + nothing +end + +""" +$(SIGNATURES) + +Implementation of non-increasing Smolyak levels, for unit testing. +""" +function nonincreasing_smolyak_levels(::Val{N}, total, each) where N + I = 0:each + c = Vector{NTuple{N,Int}}() + for i in Iterators.product(ntuple(_ -> I, Val(N))...) + if sum(i) ≤ total && issorted(i, rev = true) + push!(c, i) + end end + c end