From cff012802fab30a7cebd0847588fd9cc98361769 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tam=C3=A1s=20K=2E=20Papp?= Date: Fri, 14 Aug 2026 12:48:08 +0200 Subject: [PATCH 01/34] add comments/docstrings --- src/utilities.jl | 12 ++++++++++++ test/utilities.jl | 12 ++++++++++-- 2 files changed, 22 insertions(+), 2 deletions(-) diff --git a/src/utilities.jl b/src/utilities.jl index 80d3dbc..5b08d56 100644 --- a/src/utilities.jl +++ b/src/utilities.jl @@ -1,3 +1,11 @@ +##### +##### utilities +##### + +#### +#### printing +#### + const _SUPERSCRIPT_DIGITS = ['⁰', '¹', '²', '³', '⁴', '⁵', '⁶', '⁷', '⁸', '⁹'] const _SUBSCRIPT_DIGITS = ['₀', '₁', '₂', '₃', '₄', '₅', '₆', '₇', '₈', '₉'] @@ -53,6 +61,10 @@ Base.print(io::IO, s::SubScript) = print_number(io, _SUBSCRIPT_DIGITS, s.i) # end # end +#### +#### conversions +#### + """ $(SIGNATURES) diff --git a/test/utilities.jl b/test/utilities.jl index babf0d4..be71a92 100644 --- a/test/utilities.jl +++ b/test/utilities.jl @@ -46,7 +46,11 @@ function rand_in_domain(basis::TransformedBasis) transform_from(parent, transformation, rand_in_domain(parent)) end -"Flags (`true`) for elements in `a` that are within `atol` of some element in `b`." +""" +$(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)) @@ -54,7 +58,11 @@ function is_approximately_in(a, b; atol = √eps()) 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 From 31a2747464e760f7c70b289f0dc3081a4165ac38 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tam=C3=A1s=20K=2E=20Papp?= Date: Mon, 17 Aug 2026 16:55:13 +0200 Subject: [PATCH 02/34] rework transformations API --- src/transformations.jl | 165 +++++++++++++++++------------------ test/test_transformations.jl | 46 +++++----- test/utilities.jl | 23 ++--- 3 files changed, 114 insertions(+), 120 deletions(-) diff --git a/src/transformations.jl b/src/transformations.jl index e535d97..530e8da 100644 --- a/src/transformations.jl +++ b/src/transformations.jl @@ -9,12 +9,6 @@ export domain, domain_kind, transform_to, transform_from, coordinate_transformat #### generic api #### -""" - -!!! note - Abstract type used for code organization, not exported. -""" - """ $(TYPEDEF) @@ -151,11 +145,19 @@ struct BoundedLinear{T <: Real} <: AbstractUnivariateTransformation 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`.") + @doc """ + $(SIGNATURES) + + Transform the domain to `y ∈ (lower, upper)`, using a linear mapping. + + `lower < upper` is enforced. + """ + function BoundedLinear(; lower::Real, upper::Real) + @argcheck isfinite(lower) && isfinite(upper) DomainError + lower, upper = promote(lower, upper) + s = (upper - lower) / 2 + m = (lower + upper) / 2 + @argcheck s > 0 DomainError((; lower, upper), "Need `lower < upper`.") m, s = promote(m, s) new{typeof(m)}(m, s) end @@ -166,15 +168,6 @@ function Base.show(io::IO, transformation::BoundedLinear) print(io, "(", m - s, ",", m + s, ") ↔ domain [linear transformation]") end -""" -$(TYPEDEF) - -Transform the domain to `y ∈ (a, b)`, using ``y = x ⋅ s + m``. - -`m` and `s` are calculated and checked by the constructor; `a < b` is enforced. -""" -BoundedLinear(a::Real, b::Real) = BoundedLinear(promote(a, b)...) - function transform_from(::PM1, t::BoundedLinear, x::Scalar) (; m, s) = t x * s + m @@ -204,49 +197,48 @@ 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,∞)" + (; endpoint, scale) = transformation + if scale > 0 + D = "($(endpoint),∞)" else - D = "(-∞,A)" + D = "(-∞,$(endpoint))" end - print(io, D, " ↔ domain [rational transformation with scale ", L, "]") + print(io, D, " ↔ domain [rational transformation with 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 +247,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 +270,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, "(-∞,∞) ↔ domain [rational transformation with 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 +313,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/test/test_transformations.jl b/test/test_transformations.jl index b37d22d..1c326e1 100644 --- a/test/test_transformations.jl +++ b/test/test_transformations.jl @@ -1,11 +1,11 @@ using SpectralKit: PM1 @testset "bounded linear domain transformations" begin - @test_throws DomainError BoundedLinear(-1.0, Inf) - @test_throws DomainError BoundedLinear(-1.0, -2.0) + @test_throws DomainError BoundedLinear(; lower = -1.0, upper = Inf) + @test_throws DomainError BoundedLinear(; lower = -1.0, upper = -2.0) A, B = 1, 5 - trans = BoundedLinear(A, B) + trans = BoundedLinear(; lower = A, upper = B) @test extrema(domain(trans)) == (A, B) @@ -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) @@ -93,8 +91,8 @@ end end @testset "coordinate transformations" begin - t1 = BoundedLinear(2.0, 3.0) - t2 = SemiInfRational(7.0, 1.0) + t1 = BoundedLinear(; lower = 2.0, upper = 3.0) + t2 = SemiInfRational(; endpoint = 7.0, scale = 1.0) ct = coordinate_transformations(t1, t2) md = coordinate_domains(Val(2), PM1()) x = SVector(rand_pm1(), rand_pm1()) @@ -112,13 +110,13 @@ end @testset "printing, promotion, broadcasting" begin v = [1.0, 2.0] - t1 = BoundedLinear(2.0, 3) + t1 = BoundedLinear(; lower = 2.0, upper = 3) @test repr(t1) == "(2.0,3.0) ↔ domain [linear transformation]" @test transform_to.(PM1(), t1, v) isa Vector - t2 = SemiInfRational(7.0, 1) + t2 = SemiInfRational(; endpoint = 7.0) @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) + t3 = InfRational(; center = 0.5, scale = 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) diff --git a/test/utilities.jl b/test/utilities.jl index be71a92..fcf56a5 100644 --- a/test/utilities.jl +++ b/test/utilities.jl @@ -34,17 +34,20 @@ 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 + +# FIXME remove +# function rand_in_domain(basis::TransformedBasis) +# (; parent, transformation) = basis +# transform_from(parent, transformation, rand_in_domain(parent)) +# end """ $(SIGNATURES) @@ -54,7 +57,7 @@ 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 From b9b29e618d16981f2a61b26f93dda7a1f0c32049 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tam=C3=A1s=20K=2E=20Papp?= Date: Mon, 17 Aug 2026 16:55:24 +0200 Subject: [PATCH 03/34] start a changelog --- CHANGELOG.md | 20 ++++++++++++++++++++ 1 file changed, 20 insertions(+) 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)] From a420ad96ea77a60b7ce99b6dce57972a0d0c78a9 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tam=C3=A1s=20K=2E=20Papp?= Date: Mon, 17 Aug 2026 16:56:04 +0200 Subject: [PATCH 04/34] rework chebyshev and levels --- Project.toml | 2 + src/SpectralKit.jl | 9 +- src/chebyshev.jl | 142 ---------------------------- src/generic_api.jl | 191 +++++++++++++++---------------------- src/univariate.jl | 205 ++++++++++++++++++++++++++++++++++++++++ test/runtests.jl | 11 ++- test/test_chebyshev.jl | 146 ---------------------------- test/test_univariate.jl | 196 ++++++++++++++++++++++++++++++++++++++ 8 files changed, 490 insertions(+), 412 deletions(-) delete mode 100644 src/chebyshev.jl create mode 100644 src/univariate.jl delete mode 100644 test/test_chebyshev.jl create mode 100644 test/test_univariate.jl diff --git a/Project.toml b/Project.toml index a646fac..6c0e3e5 100644 --- a/Project.toml +++ b/Project.toml @@ -9,6 +9,7 @@ projects = ["test", "docs"] [deps] ArgCheck = "dce04be8-c92d-5529-be00-80e4d2c0e197" Compat = "34da2185-b29b-5c13-b0c7-acf172513d20" +ConcreteStructs = "2569d6c7-a4a2-43d3-a901-331e8e4be471" DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" InverseFunctions = "3587e190-3f89-42d0-90ee-14403ec27112" OrderedCollections = "bac558e1-5e72-5ebc-8fee-abe8a469f55d" @@ -17,6 +18,7 @@ StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" [compat] ArgCheck = "1, 2" Compat = "4.18.0" +ConcreteStructs = "0.2.7" DocStringExtensions = "0.8, 0.9" InverseFunctions = "0.1" OrderedCollections = "1, 2" diff --git a/src/SpectralKit.jl b/src/SpectralKit.jl index 4757c2c..1221609 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("smolyak_traversal.jl") -include("smolyak_api.jl") -include("experimental.jl") # experimental code is not part of the API, see its module docstring +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 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/generic_api.jl b/src/generic_api.jl index 1b00d76..7a76948 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_basis, 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,100 @@ function collocation_matrix(basis, x = grid(basis)) end """ -`$(FUNCTIONNAME)(basis1, basis2, θ1)` +`$(FUNCTIONNAME)(basis, Δ)` -Return a set of coefficients `θ2` for `basis2` such that -```julia -linear_combination(basis1, θ1, x) == linear_combination(basis2, θ2, x) -``` -for any `x` in the domain. In practice this means padding with zeros. - -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). +Make `basis` thicker (`Δ > 0`) or thinner (`Δ < 0`). Return `nothing` if this is not +possible. """ -function augment_coefficients end +function adjust_basis end """ -$(SIGNATURES) +`$(FUNCTIONNAME)(θ1, basis1, basis2)` -Return a `Bool` indicating whether coefficients in `basis1` can be augmented to `basis2` -with [`augment_coefficients`](@ref). +Return a set of coefficients `θ2` for `basis2` such that is either augmentes +```julia +linear_combination(basis1, θ1, x) ≈ linear_combination(basis2, θ2, x) +``` +for any `x` in the domain. In practice this means either padding with zeros, or +truncating the coefficients. -!!! 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). +Throw a `ArgumentError` if this is not possible. See [`adjust_basis`](@ref). """ -is_subset_basis(basis1::FunctionBasis, basis2::FunctionBasis) = false +function adjust_coefficients end -#### -#### transformed basis -#### +# FIXME remove this +# #### +# #### transformed basis +# #### -""" -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 +# """ +# 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 -function Base.:(∘)(parent::FunctionBasis, transformation) - TransformedBasis(parent, transformation) -end +# function Base.:(∘)(parent::FunctionBasis, transformation) +# TransformedBasis(parent, transformation) +# end -Base.parent(basis::TransformedBasis) = basis.parent +# Base.parent(basis::TransformedBasis) = basis.parent -""" -$(SIGNATURES) +# """ +# $(SIGNATURES) -Return the transformation of transformed bases, or `nothing` it not applicable. -""" -transformation(basis::TransformedBasis) = basis.transformation +# Return the transformation of transformed bases, or `nothing` it not applicable. +# """ +# transformation(basis::TransformedBasis) = basis.transformation -transformation(::FunctionBasis) = nothing +# transformation(::FunctionBasis) = nothing -domain(basis::TransformedBasis) = domain(basis.transformation) +# domain(basis::TransformedBasis) = domain(basis.transformation) -dimension(basis::TransformedBasis) = dimension(basis.parent) +# dimension(basis::TransformedBasis) = dimension(basis.parent) -function basis_at(basis::TransformedBasis, x) - (; parent, transformation) = basis - basis_at(parent, transform_to(domain(parent), transformation, x)) -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 +# 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 -function Base.:(∘)(linear_combination::LinearCombination, transformation) - (; basis, θ) = linear_combination - LinearCombination(basis ∘ transformation, θ) -end +# function Base.:(∘)(linear_combination::LinearCombination, transformation) +# (; basis, θ) = linear_combination +# LinearCombination(basis ∘ transformation, θ) +# end -Base.length(basis::TransformedBasis{<:MultivariateBasis}) = length(basis.parent) +# Base.length(basis::TransformedBasis{<:MultivariateBasis}) = length(basis.parent) -function Base.getindex(basis::TransformedBasis{<:MultivariateBasis}, i::Int) - (; parent, transformation) = basis - TransformedBasis(parent[i], Tuple(transformation)[i]) -end +# function Base.getindex(basis::TransformedBasis{<:MultivariateBasis}, i::Int) +# (; parent, transformation) = basis +# TransformedBasis(parent[i], Tuple(transformation)[i]) +# end -function is_subset_basis(basis1::TransformedBasis, basis2::TransformedBasis) - basis1.transformation ≡ basis2.transformation && - is_subset_basis(basis1.parent, basis2.parent) -end +# function is_subset_basis(basis1::TransformedBasis, basis2::TransformedBasis) +# basis1.transformation ≡ basis2.transformation && +# is_subset_basis(basis1.parent, basis2.parent) +# end -function augment_coefficients(basis1::TransformedBasis, basis2::TransformedBasis, θ1) - @argcheck is_subset_basis(basis1, basis2) - augment_coefficients(basis1.parent, basis2.parent, θ1) -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_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/univariate.jl b/src/univariate.jl new file mode 100644 index 0000000..83448ef --- /dev/null +++ b/src/univariate.jl @@ -0,0 +1,205 @@ +#### +#### Chebyshev polynomials on [-1,1] +#### + +export Chebyshev, Endpoints, Interior, univariate_basis + +#### +#### 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 + N::Int +end + +Base.eltype(::Type{<:ChebyshevIterator{T}}) where {T} = T + +Base.length(itr::ChebyshevIterator) = itr.N + +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 + +""" +$(SIGNATURES) + +Helper function to calculate the extrema of the `N`th Chebyshev polynomial, indexed by +`1 ≤ i ≤ N` (not checked). `N==1` (the constant) is special cases 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 + +struct ChebyshevShuffle + N::Int + endpoints::Bool + ChebyshevShuffle(N::Int; endpoints) = new(N, endpoints) +end + +function Base.length(shuffle::ChebyshevShuffle) + (; N, endpoints) = shuffle + endpoints ? N : max(1, N - 2) +end + +Base.eltype(::Type{ChebyshevShuffle}) = Int + +function Base.iterate(shuffle::ChebyshevShuffle, state = (0, -1)) + (; N, endpoints) = shuffle + i, step = state + if step == -1 # sentinel for first element + i = (N + 1) ÷ 2 + if endpoints # go to 1 + i′ = 1 + step′ = N - 1 + else # skip endpoints, go to next layer + step′ = (N - 1) ÷ 2 + i′ = step′ ÷ 2 + 1 + end + i, (i′, step′) + elseif step ≤ 1 # N = 1, iteration is done + nothing + else + i′ = i + step + if i′ > N # overrun, halve step and back + step = step ÷ 2 + i′ = step ÷ 2 + 1 + end + i, (i′, step) + end +end + +#### +#### grids +#### + +""" +$(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 ≥ 1 + (level ≤ 2 ? (level - 1) * 2 : (1 << (level - 1))) + 1 +end + +function block_length(::Chebyshev, ::Endpoints, level::Int) + @argcheck level ≥ 1 + level ≤ 2 ? level : 1 << (level - 2) +end + +""" +$(TYPEDEF) + +Like [`Endpoints`](@ref), but with endpoints dropped. +""" +struct Interior end + +function grid_length(::Chebyshev, ::Interior, level::Int) + @argcheck level ≥ 1 + (1 << level) - 1 +end + +function block_length(::Chebyshev, ::Interior, level::Int) + @argcheck level ≥ 1 + 1 << (level - 1) +end + +""" +Implementation of univariate bases. Not part of the API. +""" +@concrete struct UnivariateBasis <: FunctionBasis + family + grid_kind + domain_transformation + level +end + +""" +$(SIGNATURES) + +Univariate basis from `family`, using the given `grid_kind`.. + +`level` is an integer, starting from `1`, specifying the number of *blocks* used to +build the grid, which in turn determine + +""" +function univariate_basis(family, grid_kind, domain_transformation, level) + @argcheck level ≥ 1 + UnivariateBasis(family, grid_kind, domain_transformation, level) +end + +domain(U::UnivariateBasis) = domain(U.domain_transformation) + +dimension(U::UnivariateBasis) = grid_length(U.family, U.grid_kind, U.level) + +function basis_at(U::UnivariateBasis{Chebyshev}, x::Scalar) + ChebyshevIterator(transform_to(PM1(), U.domain_transformation, x), dimension(U)) +end + +function grid(::Type{T}, + U::UnivariateBasis{Chebyshev,K}) where {T <: AbstractFloat, + K <: Union{Interior,Endpoints}} + (; family, grid_kind, domain_transformation, level) = U + N = grid_length(family, grid_kind, level) + if K ≡ Interior + N += 2 # account for dropped endpoints + endpoints = false + else + endpoints = true + end + (transform_from(PM1(), U.domain_transformation, _chebyshev_extremum(T, i, N)) + for i in ChebyshevShuffle(N; endpoints)) +end + +function adjust_basis(U::UnivariateBasis, Δ::Int) + level′ = U.level + Δ + if level′ > 0 + UnivariateBasis(U.family, U.grid_kind, U.domain_transformation, U.level - Δ) + else + nothing + end +end + +function adjust_coefficients(θ1::AbstractVector, U1::UnivariateBasis{Chebyshev}, U2::UnivariateBasis{Chebyshev}) + 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 + θ1[1:d1] + else + θ2 = zeros(d2) + θ2[1:d1] .= θ1 + θ2 + end +end diff --git a/test/runtests.jl b/test/runtests.jl index d119afe..3c28b3b 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,4 +1,5 @@ using SpectralKit, Test +import SpectralKit as SK using JET @testset "static analysis with JET.jl" begin @@ -20,8 +21,8 @@ include("test_utilities.jl") include("test_derivatives.jl") include("test_domains.jl") include("test_transformations.jl") -include("test_chebyshev.jl") -include("test_smolyak_traversal.jl") -include("test_smolyak.jl") -include("test_generic_api.jl") # NOTE moved last as it used constructs from above -include("test_experimental.jl") # NOTE experimental code is not public API +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 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_univariate.jl b/test/test_univariate.jl new file mode 100644 index 0000000..e7d2124 --- /dev/null +++ b/test/test_univariate.jl @@ -0,0 +1,196 @@ +#### +#### Chebyshev on [-1,1] +#### + +""" +For testing Chebyshev blocks and shuffle. +""" +function _shuffle(len::Int; endpoints = true) + if len == 1 + return [[1]] + end + if !endpoints + len += 2 # for Interior, we start the algorithm with two extra points + 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 [Interior(), Endpoints()] + for level in 1:5 + N = SK.grid_length(Chebyshev(), kind, level) + block_lengths = map(level -> SK.block_length(Chebyshev(), kind, level), 1:level) + @test sum(block_lengths) == N + S = _shuffle(N; endpoints = kind == Endpoints()) + @test block_lengths == map(length, S) + end + end +end + + +@testset "Chebyshev grid" begin + for kind in [Interior(), Endpoints()] + previous = Float64[] + for level in 1:7 + b = univariate_basis(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 univariate_basis(Chebyshev(), Endpoints(), transformation, 0) + for grid_kind in (Interior(), Endpoints()) + for level in 1:5 + basis = univariate_basis(Chebyshev(), grid_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) + @test all(x -> isapprox(chebyshev_cos_deriv(transform_to(PM1(), transformation, x), N + 2), + 0, atol = 1e-13), g) + end + end +end + +@testset "Chebysev adjusted basis" begin +transformation = SemiInfRational + level = 3 + for grid_kind in (Interior(), Endpoints()) + basis = univariate_basis(Chebyshev(), grid_kind, + for Δ in 1:4 + + +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 From cf8540cccd5143097619ad5f39860072e54c5066 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tam=C3=A1s=20K=2E=20Papp?= Date: Mon, 17 Aug 2026 17:05:20 +0200 Subject: [PATCH 05/34] test adjustment code --- src/univariate.jl | 8 ++++---- test/test_univariate.jl | 19 ++++++++++++++----- 2 files changed, 18 insertions(+), 9 deletions(-) diff --git a/src/univariate.jl b/src/univariate.jl index 83448ef..4e2d2e4 100644 --- a/src/univariate.jl +++ b/src/univariate.jl @@ -183,7 +183,7 @@ end function adjust_basis(U::UnivariateBasis, Δ::Int) level′ = U.level + Δ if level′ > 0 - UnivariateBasis(U.family, U.grid_kind, U.domain_transformation, U.level - Δ) + UnivariateBasis(U.family, U.grid_kind, U.domain_transformation, level′) else nothing end @@ -195,9 +195,9 @@ function adjust_coefficients(θ1::AbstractVector, U1::UnivariateBasis{Chebyshev} @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 - θ1[1:d1] - else + if d2 ≤ d1 # truncate + θ1[1:d2] + else # pad with zeros θ2 = zeros(d2) θ2[1:d1] .= θ1 θ2 diff --git a/test/test_univariate.jl b/test/test_univariate.jl index e7d2124..b3bdc1a 100644 --- a/test/test_univariate.jl +++ b/test/test_univariate.jl @@ -97,13 +97,22 @@ end end @testset "Chebysev adjusted basis" begin -transformation = SemiInfRational - level = 3 + transformation = SemiInfRational(; endpoint = 3.0, scale = 7.0) for grid_kind in (Interior(), Endpoints()) - basis = univariate_basis(Chebyshev(), grid_kind, + level0 = 3 + basis0 = univariate_basis(Chebyshev(), grid_kind, transformation, level0) + θ0 = randn(dimension(basis0)) for Δ in 1:4 - - + basis = adjust_basis(basis0, Δ) + θ = 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 + @test adjust_basis(basis0, -3) ≡ nothing end # # compatible and incompatible grids From 5e90056c057aef3fc0effc957eb54eb5015998ce Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tam=C3=A1s=20K=2E=20Papp?= Date: Mon, 17 Aug 2026 17:09:31 +0200 Subject: [PATCH 06/34] add back test for derivatives --- test/test_univariate.jl | 65 ++++++++++------------------------------- 1 file changed, 16 insertions(+), 49 deletions(-) diff --git a/test/test_univariate.jl b/test/test_univariate.jl index b3bdc1a..0d24ec1 100644 --- a/test/test_univariate.jl +++ b/test/test_univariate.jl @@ -115,55 +115,22 @@ end @test adjust_basis(basis0, -3) ≡ nothing 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 "univariate derivatives" begin + for (transformation, N) in ((BoundedLinear(lower = -2, upper = 3), 5), + (SemiInfRational(endpoint = 0.7, scale = 0.3), 1), + (InfRational(center = 0.4, scale = 0.9), 1)) + basis = univariate_basis(Chebyshev(), Interior(), transformation, 3) + 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) atol = 1e-6 + end + end + end +end # @testset "endpoint continuity for derivatives" begin # N = 10 From c0955fd583bc13be7dd272ed13544c5ad7064057 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tam=C3=A1s=20K=2E=20Papp?= Date: Tue, 18 Aug 2026 13:45:09 +0200 Subject: [PATCH 07/34] minor test fixes --- test/runtests.jl | 21 +++++----- test/test_univariate.jl | 89 +++++++++++++++++++---------------------- test/utilities.jl | 20 +++++++++ 3 files changed, 73 insertions(+), 57 deletions(-) diff --git a/test/runtests.jl b/test/runtests.jl index 3c28b3b..09feb7a 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,17 +1,18 @@ 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 diff --git a/test/test_univariate.jl b/test/test_univariate.jl index 0d24ec1..4030054 100644 --- a/test/test_univariate.jl +++ b/test/test_univariate.jl @@ -48,7 +48,8 @@ end for kind in [Interior(), Endpoints()] previous = Float64[] for level in 1:7 - b = univariate_basis(Chebyshev(), kind, BoundedLinear(1.0, 3.0), level) + b = univariate_basis(Chebyshev(), kind, + BoundedLinear(; lower = 1.0, upper = 3.0), level) g = collect(grid(b)) @test is_approximate_subset(previous, g) g = previous @@ -58,15 +59,15 @@ end @testset "Chebyshev basics" begin - transformation = BoundedLinear(1.0, 3.0) + transformation = BoundedLinear(; lower = 1.0, upper = 3.0) @test_throws ArgumentError univariate_basis(Chebyshev(), Endpoints(), transformation, 0) - for grid_kind in (Interior(), Endpoints()) + for kind in (Interior(), Endpoints()) for level in 1:5 - basis = univariate_basis(Chebyshev(), grid_kind, transformation, level) + basis = univariate_basis(Chebyshev(), kind, transformation, level) @test is_function_basis(basis) @test is_function_basis(typeof(basis)) N = @inferred dimension(basis) - @test N > 1 + @test N ≥ 1 # check linear combinations for _ in 1:100 @@ -90,17 +91,17 @@ end @test length(g) == N a, b = extrema(domain(basis)) @test all(a .≤ g .≤ b) - @test all(x -> isapprox(chebyshev_cos_deriv(transform_to(PM1(), transformation, x), N + 2), - 0, atol = 1e-13), g) + @test all(x -> is_chebyshev_extrema(transform_to(PM1(), transformation, x), + kind ≡ Interior() ? N + 2 : N), g) end end end @testset "Chebysev adjusted basis" begin transformation = SemiInfRational(; endpoint = 3.0, scale = 7.0) - for grid_kind in (Interior(), Endpoints()) + for kind in (Interior(), Endpoints()) level0 = 3 - basis0 = univariate_basis(Chebyshev(), grid_kind, transformation, level0) + basis0 = univariate_basis(Chebyshev(), kind, transformation, level0) θ0 = randn(dimension(basis0)) for Δ in 1:4 basis = adjust_basis(basis0, Δ) @@ -111,8 +112,8 @@ end @test linear_combination(basis0, θ0, x) ≈ linear_combination(basis, θ, x) end end + @test adjust_basis(basis0, -3) ≡ nothing end - @test adjust_basis(basis0, -3) ≡ nothing end @testset "univariate derivatives" begin @@ -132,41 +133,35 @@ 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) +@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 = univariate_basis(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 -# 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 + @testset "InfRational endpoints continuity" begin + basis = univariate_basis(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 diff --git a/test/utilities.jl b/test/utilities.jl index fcf56a5..c40570e 100644 --- a/test/utilities.jl +++ b/test/utilities.jl @@ -6,6 +6,12 @@ using SpectralKit: TransformedBasis, SmolyakBasis, SmolyakIndices # dispatch for 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 +27,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-10) + 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. From fab34de1b7a214f1329ad3d4fdc4cfbf289d03b8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tam=C3=A1s=20K=2E=20Papp?= Date: Tue, 18 Aug 2026 14:29:03 +0200 Subject: [PATCH 08/34] rename grid_kind to kind, reenable what is needed --- src/univariate.jl | 16 ++++++++-------- test/utilities.jl | 3 ++- 2 files changed, 10 insertions(+), 9 deletions(-) diff --git a/src/univariate.jl b/src/univariate.jl index 4e2d2e4..f8c23c2 100644 --- a/src/univariate.jl +++ b/src/univariate.jl @@ -138,7 +138,7 @@ Implementation of univariate bases. Not part of the API. """ @concrete struct UnivariateBasis <: FunctionBasis family - grid_kind + kind domain_transformation level end @@ -146,20 +146,20 @@ end """ $(SIGNATURES) -Univariate basis from `family`, using the given `grid_kind`.. +Univariate basis from `family`, using the given `kind`. `level` is an integer, starting from `1`, specifying the number of *blocks* used to build the grid, which in turn determine """ -function univariate_basis(family, grid_kind, domain_transformation, level) +function univariate_basis(family, kind, domain_transformation, level) @argcheck level ≥ 1 - UnivariateBasis(family, grid_kind, domain_transformation, level) + UnivariateBasis(family, kind, domain_transformation, level) end domain(U::UnivariateBasis) = domain(U.domain_transformation) -dimension(U::UnivariateBasis) = grid_length(U.family, U.grid_kind, U.level) +dimension(U::UnivariateBasis) = grid_length(U.family, U.kind, U.level) function basis_at(U::UnivariateBasis{Chebyshev}, x::Scalar) ChebyshevIterator(transform_to(PM1(), U.domain_transformation, x), dimension(U)) @@ -168,8 +168,8 @@ end function grid(::Type{T}, U::UnivariateBasis{Chebyshev,K}) where {T <: AbstractFloat, K <: Union{Interior,Endpoints}} - (; family, grid_kind, domain_transformation, level) = U - N = grid_length(family, grid_kind, level) + (; family, kind, domain_transformation, level) = U + N = grid_length(family, kind, level) if K ≡ Interior N += 2 # account for dropped endpoints endpoints = false @@ -183,7 +183,7 @@ end function adjust_basis(U::UnivariateBasis, Δ::Int) level′ = U.level + Δ if level′ > 0 - UnivariateBasis(U.family, U.grid_kind, U.domain_transformation, level′) + UnivariateBasis(U.family, U.kind, U.domain_transformation, level′) else nothing end diff --git a/test/utilities.jl b/test/utilities.jl index c40570e..3897dbf 100644 --- a/test/utilities.jl +++ b/test/utilities.jl @@ -2,7 +2,8 @@ ##### utility functions for tests ##### -using SpectralKit: TransformedBasis, SmolyakBasis, SmolyakIndices # dispatch for rand_in_domain +# FIXME reenable what is needed +# using SpectralKit: TransformedBasis, SmolyakBasis, SmolyakIndices # dispatch for rand_in_domain chebyshev_cos(x, n) = cos((n - 1) * acos(x)) From bf0dcf26e6b640741bb45a5da0c09364be207df6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tam=C3=A1s=20K=2E=20Papp?= Date: Tue, 18 Aug 2026 15:10:17 +0200 Subject: [PATCH 09/34] make Chebyshev iterator infinite --- src/univariate.jl | 14 ++++++-------- 1 file changed, 6 insertions(+), 8 deletions(-) diff --git a/src/univariate.jl b/src/univariate.jl index f8c23c2..593ba9e 100644 --- a/src/univariate.jl +++ b/src/univariate.jl @@ -21,23 +21,21 @@ struct Chebyshev end struct ChebyshevIterator{T} x::T - N::Int end Base.eltype(::Type{<:ChebyshevIterator{T}}) where {T} = T -Base.length(itr::ChebyshevIterator) = itr.N +Base.IteratorSize(::Type{<:ChebyshevIterator}) = Base.IsInfinite() function Base.iterate(itr::ChebyshevIterator{T}) where T (; x) = itr - _one(T), (2, _one(T), x) + _one(T), (_one(T), x) end -function Base.iterate(itr::ChebyshevIterator{T}, (i, fp, fpp)) where T - (; x, N) = itr - i > N && return nothing +function Base.iterate(itr::ChebyshevIterator{T}, (fp, fpp)) where T + (; x) = itr f = _sub(_mul(2, x, fp), fpp) - f::T, (i + 1, f, fp) + f::T, (f, fp) end """ @@ -162,7 +160,7 @@ domain(U::UnivariateBasis) = domain(U.domain_transformation) dimension(U::UnivariateBasis) = grid_length(U.family, U.kind, U.level) function basis_at(U::UnivariateBasis{Chebyshev}, x::Scalar) - ChebyshevIterator(transform_to(PM1(), U.domain_transformation, x), dimension(U)) + Iterators.take(ChebyshevIterator(transform_to(PM1(), U.domain_transformation, x)), dimension(U)) end function grid(::Type{T}, From ec123d37d64fc3978871b6facd8b56d96db9a21d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tam=C3=A1s=20K=2E=20Papp?= Date: Tue, 18 Aug 2026 17:10:07 +0200 Subject: [PATCH 10/34] Clean up grid calculations, tighter testing. --- src/univariate.jl | 118 ++++++++++++++++++++-------------------- test/test_univariate.jl | 12 ++-- 2 files changed, 66 insertions(+), 64 deletions(-) diff --git a/src/univariate.jl b/src/univariate.jl index 593ba9e..1f3c804 100644 --- a/src/univariate.jl +++ b/src/univariate.jl @@ -42,7 +42,8 @@ end $(SIGNATURES) Helper function to calculate the extrema of the `N`th Chebyshev polynomial, indexed by -`1 ≤ i ≤ N` (not checked). `N==1` (the constant) is special cases to zero, for nesting. +`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 @@ -52,47 +53,50 @@ function _chebyshev_extremum(::Type{T}, i::Int, N::Int) where {T <: Real} end end -struct ChebyshevShuffle - N::Int - endpoints::Bool - ChebyshevShuffle(N::Int; endpoints) = new(N, endpoints) -end +#### +#### kinds +#### + +""" +$(TYPEDEF) + +Like [`Endpoints`](@ref), but with endpoints dropped. +""" +struct Interior end + +""" +$(SIGNATURES) -function Base.length(shuffle::ChebyshevShuffle) - (; N, endpoints) = shuffle - endpoints ? N : max(1, N - 2) +Length of a univariate grid. +""" +function grid_length(::Chebyshev, ::Interior, level::Int) + @argcheck level ≥ 1 + (1 << level) - 1 end -Base.eltype(::Type{ChebyshevShuffle}) = Int +""" +$(SIGNATURES) -function Base.iterate(shuffle::ChebyshevShuffle, state = (0, -1)) - (; N, endpoints) = shuffle - i, step = state - if step == -1 # sentinel for first element - i = (N + 1) ÷ 2 - if endpoints # go to 1 - i′ = 1 - step′ = N - 1 - else # skip endpoints, go to next layer - step′ = (N - 1) ÷ 2 - i′ = step′ ÷ 2 + 1 - end - i, (i′, step′) - elseif step ≤ 1 # N = 1, iteration is done - nothing - else - i′ = i + step - if i′ > N # overrun, halve step and back - step = step ÷ 2 - i′ = step ÷ 2 + 1 - end - i, (i′, step) - end +Length of a single block, these are concatenated to form the grid. +""" +function block_length(::Chebyshev, ::Interior, level::Int) + @argcheck level ≥ 1 + 1 << (level - 1) end -#### -#### grids -#### +""" +$(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) # trust constant folding fast path to top_set_bit + remainder = i - (1 << (p - 1)) + start = 1 << (level - p) + start + remainder * (start << 1) + 1 +end """ $(TYPEDEF) @@ -114,22 +118,19 @@ function block_length(::Chebyshev, ::Endpoints, level::Int) level ≤ 2 ? level : 1 << (level - 2) end -""" -$(TYPEDEF) - -Like [`Endpoints`](@ref), but with endpoints dropped. -""" -struct Interior end - -function grid_length(::Chebyshev, ::Interior, level::Int) - @argcheck level ≥ 1 - (1 << level) - 1 +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 - 2) + else + 1 + (i - 2) * (1 << (level - 1)) + end end -function block_length(::Chebyshev, ::Interior, level::Int) - @argcheck level ≥ 1 - 1 << (level - 1) -end +#### +#### univariate bases +#### """ Implementation of univariate bases. Not part of the API. @@ -163,19 +164,18 @@ function basis_at(U::UnivariateBasis{Chebyshev}, x::Scalar) Iterators.take(ChebyshevIterator(transform_to(PM1(), U.domain_transformation, x)), dimension(U)) end -function grid(::Type{T}, - U::UnivariateBasis{Chebyshev,K}) where {T <: AbstractFloat, - K <: Union{Interior,Endpoints}} +function grid(::Type{T}, U::UnivariateBasis{Chebyshev}) where {T <: AbstractFloat} (; family, kind, domain_transformation, level) = U N = grid_length(family, kind, level) - if K ≡ Interior - N += 2 # account for dropped endpoints - endpoints = false + if kind ≡ Interior() + N̂ = N + 2 else - endpoints = true + @assert kind ≡ Endpoints() + N̂ = N end - (transform_from(PM1(), U.domain_transformation, _chebyshev_extremum(T, i, N)) - for i in ChebyshevShuffle(N; endpoints)) + (transform_from(PM1(), U.domain_transformation, + _chebyshev_extremum(T, _chebyshev_extremum_shuffle(kind, i, level), N̂)) + for i in 1:N) end function adjust_basis(U::UnivariateBasis, Δ::Int) diff --git a/test/test_univariate.jl b/test/test_univariate.jl index 4030054..5aa071a 100644 --- a/test/test_univariate.jl +++ b/test/test_univariate.jl @@ -6,12 +6,12 @@ For testing Chebyshev blocks and shuffle. """ function _shuffle(len::Int; endpoints = true) - if len == 1 - return [[1]] - end 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)) @@ -39,11 +39,11 @@ end @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 [Interior(), Endpoints()] previous = Float64[] @@ -57,7 +57,6 @@ end end end - @testset "Chebyshev basics" begin transformation = BoundedLinear(; lower = 1.0, upper = 3.0) @test_throws ArgumentError univariate_basis(Chebyshev(), Endpoints(), transformation, 0) @@ -91,6 +90,9 @@ end @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 From 357eb5def638afe7cebaca1ea9e5544faab60364 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tam=C3=A1s=20K=2E=20Papp?= Date: Tue, 18 Aug 2026 17:35:12 +0200 Subject: [PATCH 11/34] reindex levels from 0 --- src/univariate.jl | 29 ++++++++++++++--------------- test/test_univariate.jl | 4 ++-- test/utilities.jl | 2 +- 3 files changed, 17 insertions(+), 18 deletions(-) diff --git a/src/univariate.jl b/src/univariate.jl index 1f3c804..48e7459 100644 --- a/src/univariate.jl +++ b/src/univariate.jl @@ -70,8 +70,8 @@ $(SIGNATURES) Length of a univariate grid. """ function grid_length(::Chebyshev, ::Interior, level::Int) - @argcheck level ≥ 1 - (1 << level) - 1 + @argcheck level ≥ 0 + (1 << (level + 1)) - 1 end """ @@ -80,8 +80,8 @@ $(SIGNATURES) Length of a single block, these are concatenated to form the grid. """ function block_length(::Chebyshev, ::Interior, level::Int) - @argcheck level ≥ 1 - 1 << (level - 1) + @argcheck level ≥ 0 + 1 << level end """ @@ -92,8 +92,8 @@ indices start from `2`, endpoint from `1`. Caller is responsible for making sure `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) # trust constant folding fast path to top_set_bit - remainder = i - (1 << (p - 1)) + 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 @@ -109,22 +109,22 @@ Chebyshev-Lobatto grid. The extrema of Chebyshev polynomials, including endpoint struct Endpoints end function grid_length(::Chebyshev, ::Endpoints, level::Int) - @argcheck level ≥ 1 - (level ≤ 2 ? (level - 1) * 2 : (1 << (level - 1))) + 1 + @argcheck level ≥ 0 + (level ≤ 1 ? level * 2 : (1 << level)) + 1 end function block_length(::Chebyshev, ::Endpoints, level::Int) - @argcheck level ≥ 1 - level ≤ 2 ? level : 1 << (level - 2) + @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 - 2) + 1 + 1 << (level - 1) else - 1 + (i - 2) * (1 << (level - 1)) + 1 + (i - 2) * (1 << level) end end @@ -147,9 +147,8 @@ $(SIGNATURES) Univariate basis from `family`, using the given `kind`. -`level` is an integer, starting from `1`, specifying the number of *blocks* used to -build the grid, which in turn determine - +`level` is an integer, starting from `0`, specifying the number of *blocks* used to +build the grid. """ function univariate_basis(family, kind, domain_transformation, level) @argcheck level ≥ 1 diff --git a/test/test_univariate.jl b/test/test_univariate.jl index 5aa071a..da1f735 100644 --- a/test/test_univariate.jl +++ b/test/test_univariate.jl @@ -33,9 +33,9 @@ end @testset "blocks and shuffle" begin for kind in [Interior(), Endpoints()] - for level in 1:5 + for level in 0:5 N = SK.grid_length(Chebyshev(), kind, level) - block_lengths = map(level -> SK.block_length(Chebyshev(), kind, level), 1: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) diff --git a/test/utilities.jl b/test/utilities.jl index 3897dbf..1798dd4 100644 --- a/test/utilities.jl +++ b/test/utilities.jl @@ -31,7 +31,7 @@ $(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-10) +function is_chebyshev_extrema(x, n; tol = 1e-8) if abs(abs(x) - 1) ≤ tol true else From 620df901d9bbed93deed1a696f2e7991ce93c587 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tam=C3=A1s=20K=2E=20Papp?= Date: Wed, 19 Aug 2026 13:02:25 +0200 Subject: [PATCH 12/34] add KINDS variable to tests --- src/smolyak_traversal.jl | 175 ++------------------------------------- test/test_univariate.jl | 8 +- test/utilities.jl | 3 + 3 files changed, 12 insertions(+), 174 deletions(-) diff --git a/src/smolyak_traversal.jl b/src/smolyak_traversal.jl index 13b665a..51881e6 100644 --- a/src/smolyak_traversal.jl +++ b/src/smolyak_traversal.jl @@ -2,181 +2,16 @@ ##### Smolyak implementation details ##### -### -### Blocking parameters -### - -export SmolyakParameters - -struct SmolyakParameters{B,M} - function SmolyakParameters{B,M}() where {B,M} - @argcheck B isa Int && B ≥ 0 - @argcheck M isa Int && M ≥ 0 - M > B && @warn "M > B replaced with M = B" M B - new{B,min(B,M)}() # maintain M ≤ B - end -end - -function Base.show(io::IO, ::SmolyakParameters{B,M}) where {B,M} - print(io, "Smolyak parameters, ∑bᵢ ≤ $(B), all bᵢ ≤ $(M)") -end - -""" -$(SIGNATURES) - -Parameters for Smolyak grids that are *independent of the dimension of the domain*. - -Polynomials are organized into blocks (of eg `1, 2, 2, 4, 8, 16, …`) polynomials (and -corresponding gridpoints), indexed with a *block index* `b` that starts at `0`. `B ≥ ∑ -bᵢ` and `0 ≤ bᵢ ≤ M` constrain the number of blocks along each dimension `i`. - -`M > B` is not an error, but will be normalized to `M = B` with a warning. -""" -@inline function SmolyakParameters(B::Integer, M::Integer = B) - SmolyakParameters{Int(B),Int(M)}() -end - -#### -#### Nesting sizes and shuffling -#### -#### 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. -""" -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, … -### - -""" -$(SIGNATURES) - -Cumulative block length at block `b`. -""" -@inline function nesting_total_length(::Type{Chebyshev}, ::EndpointGrid, b::Int) - b == 0 ? 1 : ((1 << b) + 1) -end - -""" -$(SIGNATURES) - -Length of each block `b`. - -!!! 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 - -function Base.iterate(ι::SmolyakGridShuffle{EndpointGrid}) - i = (ι.len + 1) ÷ 2 - i, (0, 0) # step = 0 is special-cased -end - -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′) - else - nothing - end - end -end - -### -### interior grid: 1, 3, 9, 27, … -### - -@inline nesting_total_length(::Type{Chebyshev}, ::InteriorGrid, b::Int) = 3^b - -@inline function nesting_block_length(::Type{Chebyshev}, ::InteriorGrid, b::Int) - b == 0 ? 1 : 2 * 3^(b - 1) -end - -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 - end -end - -### -### interior grid type 2: 1, 3, 7, … -### - -@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 - -function Base.iterate(ι::SmolyakGridShuffle{InteriorGrid2}, (i, step)) - i′ = i + step - if i′ ≤ ι.len - i′, (i′, step) - else - step′ = step ÷ 2 - if step′ ≥ 2 - i′ = step′ ÷ 2 - i′, (i′, step′) - else - nothing - end - end -end - #### #### index traversal #### -function __inc_init(nesting_total_lengths, ::Val{N}, ::Val{B}) where {N,B} +function __inc_init(family, kind, ::Val{N}, total) where N indices = ntuple(_ -> 1, Val(N)) - blocks = ntuple(_ -> 0, Val(N)) - l = first(nesting_total_lengths) - limits = ntuple(_ -> l, Val(N)) - slack = B + levels = ntuple(_ -> 0, Val(N)) + b0 = block_length(family, kind, 0) + limits = ntuple(_ -> b0, Val(N)) + slack = total slack, indices, blocks, limits end diff --git a/test/test_univariate.jl b/test/test_univariate.jl index da1f735..416cbac 100644 --- a/test/test_univariate.jl +++ b/test/test_univariate.jl @@ -32,7 +32,7 @@ function _shuffle(len::Int; endpoints = true) end @testset "blocks and shuffle" begin - for kind in [Interior(), Endpoints()] + 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) @@ -45,7 +45,7 @@ end end @testset "Chebyshev grid" begin - for kind in [Interior(), Endpoints()] + for kind in KINDS previous = Float64[] for level in 1:7 b = univariate_basis(Chebyshev(), kind, @@ -60,7 +60,7 @@ end @testset "Chebyshev basics" begin transformation = BoundedLinear(; lower = 1.0, upper = 3.0) @test_throws ArgumentError univariate_basis(Chebyshev(), Endpoints(), transformation, 0) - for kind in (Interior(), Endpoints()) + for kind in KINDS for level in 1:5 basis = univariate_basis(Chebyshev(), kind, transformation, level) @test is_function_basis(basis) @@ -101,7 +101,7 @@ end @testset "Chebysev adjusted basis" begin transformation = SemiInfRational(; endpoint = 3.0, scale = 7.0) - for kind in (Interior(), Endpoints()) + for kind in KINDS level0 = 3 basis0 = univariate_basis(Chebyshev(), kind, transformation, level0) θ0 = randn(dimension(basis0)) diff --git a/test/utilities.jl b/test/utilities.jl index 1798dd4..411324e 100644 --- a/test/utilities.jl +++ b/test/utilities.jl @@ -2,6 +2,9 @@ ##### utility functions for tests ##### +"kinds we test on" +const KINDS = [Endpoints(), Interior()] + # FIXME reenable what is needed # using SpectralKit: TransformedBasis, SmolyakBasis, SmolyakIndices # dispatch for rand_in_domain From 98d254bce4137cf7b8a551a7c5f0d268a7ad731b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tam=C3=A1s=20K=2E=20Papp?= Date: Fri, 21 Aug 2026 17:40:10 +0200 Subject: [PATCH 13/34] reimplement Smolyak traversal more flexibly --- src/SpectralKit.jl | 2 +- src/smolyak_traversal.jl | 373 +++++++++++++++++++-------------- src/univariate.jl | 3 +- test/runtests.jl | 2 +- test/test_smolyak_traversal.jl | 133 ++---------- test/utilities.jl | 54 +++++ 6 files changed, 291 insertions(+), 276 deletions(-) diff --git a/src/SpectralKit.jl b/src/SpectralKit.jl index 1221609..271fa3f 100644 --- a/src/SpectralKit.jl +++ b/src/SpectralKit.jl @@ -12,7 +12,7 @@ include("domains.jl") include("transformations.jl") include("generic_api.jl") include("univariate.jl") -# include("smolyak_traversal.jl") +include("smolyak_traversal.jl") # include("smolyak_api.jl") # include("experimental.jl") # experimental code is not part of the API, see its module docstring diff --git a/src/smolyak_traversal.jl b/src/smolyak_traversal.jl index 51881e6..2d5d07c 100644 --- a/src/smolyak_traversal.jl +++ b/src/smolyak_traversal.jl @@ -2,85 +2,145 @@ ##### Smolyak implementation details ##### -#### -#### index traversal -#### - -function __inc_init(family, kind, ::Val{N}, total) where N - indices = ntuple(_ -> 1, Val(N)) - levels = ntuple(_ -> 0, Val(N)) - b0 = block_length(family, kind, 0) - limits = ntuple(_ -> b0, Val(N)) +function __smolyak_init(family, kind, total::Int, f::F, itrs::NTuple{N}) where {F,N} slack = total - slack, indices, blocks, limits + r = block_length(family, kind, 0) + remainders = ntuple(_ -> r - 1, Val(N)) + itr_results = map(iterate, itrs) + states = map(last, itr_results) + levels = ntuple(_ -> 0, Val(N)) + x, cached... = reverse(accumulate(f, reverse(map(first, itr_results)))) + x, slack, remainders, states, cached, levels end """ -$(SIGNATURES) - -Internal implementation of the Smolyak indexing iterator. - -# Arguments +$(SIGNATURES) → x, Δ, remainders′, states′, cached′, levels′ -- `nesting_total_lengths`: precalculated nesting total lengths, constant during iteration, - indexes with an offset of `1` +Step through Smolyak indices of iterators. -- `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}` +`itrs` yield the `xs`, `(x, cached...)` is equivalent to `reverse(accumulate(reverse(...)))` applied to t """ -@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τ...) - 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τ′...) - end +function __smolyak_step(family, kind, each::Int, f::F, itrs::NTuple{N}, + slack::Int, remainders::NTuple{N}, states::NTuple{N}, + cached, levels::NTuple{N}) where {F,N} + @assert length(cached) == N - 1 + I1, Iτ... = itrs + s1, sτ... = states + r1, rτ... = remainders + l1, lτ... = levels + if r1 > 0 # step within block + x1, s1′ = iterate(I1, s1) + (f(x1, cached), # product + 0, # no change in slack + (r1 - 1, rτ...), # one less element in 1 + (s1′, sτ...), # step iterator + cached, + levels) + elseif l1 < each && slack > 0 # next block, same tail + x1, s1′ = iterate(I1, s1) + (f(x1, cached), # step to next block + -1, # decrease slack + (block_length(family, kind, l1 + 1) - 1, rτ...), # remaining elements: all in block + + (s1′, sτ...), # step state 1 + cached, + (l1+1, lτ...)) # next level + elseif N == 1 + nothing # done with iteration + else # go into tail + _, cτ... = cached + next = __smolyak_step(family, kind, each, f, Iτ, + slack + l1, rτ, sτ, cτ, lτ) + next ≡ nothing && return nothing + xτ, Δτ, rτ′, sτ′, cτ′, lτ′ = next + x1, s1 = iterate(I1) + cached′ = (xτ, cτ′...) + (f(x1, cached′), + l1 + Δτ, # more slack as we reset 1 + (block_length(family, kind, 0) - 1, rτ′...), # all remaining in block 0 + (s1, sτ′...), # states with tail + cached′, # cache updated as above + (0, lτ′...)) # back to level 0 here end end +#### +#### index traversal +#### + +# function __inc_init(family, kind, ::Val{N}, total) where N +# indices = ntuple(_ -> 1, Val(N)) +# levels = ntuple(_ -> 0, Val(N)) +# b0 = block_length(family, kind, 0) +# limits = ntuple(_ -> b0, Val(N)) +# slack = total +# slack, indices, levels, limits +# end + +# """ +# $(SIGNATURES) → valid, Δ, indices, levels, limits + +# Internal implementation of the Smolyak indexing iterator. + +# # Arguments + +# - `slack`: `total - sum(blocks)`, cached + +# - `indices`: current indices + +# - `levels`: level along each dimension + +# - `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′`, `levels′, `limits′`: next values for corresponding arguments above, each an +# `::NTuple{N,Int}` +# """ +# @inline function __inc(family, kind, each::Int, # unchanged +# slack::Int, indices::NTuple{N,Int}, levels::NTuple{N,Int}, +# limits::NTuple{N,Int}) where N +# i1, iτ... = indices +# l1, lτ... = levels +# L1, Lτ... = limits +# if i1 < L1 # increment i1, same block +# true, 0, (i1 + 1, iτ...), levels, limits +# elseif l1 < each && slack > 0 # increment i1, next block +# l1′ = l1 + 1 +# true, -1, (i1 + 1, iτ...), (l1′, lτ...), (grid_length(family, kind, l1′), Lτ...) +# else +# if N == 1 # end of iteration, arbitrary value since !valid +# false, 0, indices, levels, limits +# else # i1 = 1, increment tail if applicable +# Δ1 = l1 +# valid, Δτ, iτ′, lτ′, Lτ′ = __inc(family, kind, each, slack + Δ1, iτ, lτ, Lτ) +# valid, Δ1 + Δτ, (1, iτ′...), (0, lτ′...), (grid_length(family, kind, 0), Lτ′...) +# end +# end +# end + """ $(SIGNATURES) Calculate the length of a [`SmolyakIndices`](@ref) iterator. Argument as in the latter. """ -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 +function __smolyak_length(family, kind, N::Int, total::Int, each::Int) + # implicit assumption: each ≤ total, enforced by the SmolyakParameters constructor + _bl(b) = block_length(family, kind, b) + c = zeros(Int, total + 1) # indexed as 0, …, total + each = min(each, total) + for b in 0:each c[b + 1] = _bl(b) end for n in 2:N - for b in B:(-1):0 # blocks with indices that sum to b + for b in total:(-1):0 # blocks with indices that sum to b s = 0 - for a in 0:min(b, M) + for a in 0:min(b, each) s += _bl(a) * c[b - a + 1] end # can safely overwrite since they will not be used again for n + 1 @@ -90,125 +150,124 @@ function __smolyak_length(grid_kind::AbstractGrid, ::Val{N}, ::Val{B}, M::Int) w sum(c) end -""" -$(TYPEDEF) +# """ +# $(TYPEDEF) -Indexing specification in a Smolyak basis/interpolation. +# Indexing specification in a Smolyak basis/interpolation. -# Type parameters +# # Type parameters -- `N`: the dimension of indices +# - `N`: the dimension of indices -- `H`: highest index visited for all dimensions +# - `H`: highest index visited for all dimensions -- `B ≥ 0`: sum of block indices, starting from `0` (ie `B = 0` has just one element), +# - `B ≥ 0`: sum of block indices, starting from `0` (ie `B = 0` has just one element), -- `M`: upper bound on each block index +# - `M`: upper bound on each block index -# Constructor +# # Constructor -Takes the dimension `N` as a parameter, `grid_kind`, and a `SmolyakParameters` object, -calculating everything else. +# Takes the dimension `N` as a parameter, `kind`, and a `SmolyakParameters` object, +# calculating everything else. -# Details +# # Details -Consider positive integer indices `(i1, …, iN)`, each starting at one. +# 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. +# Let `ℓ(b) = nesting_total_length(Chebyshev, grid_knid, kind, b)`, and `l1` denote the +# smallest integer such that `i1 ≤ ℓ(l1)`, 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: +# 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` +# 1. `1 ≤ i1 ≤ ℓ(M)`, …, `1 ≤ iN ≤ ℓ(M)`, +# 2. `0 ≤ l1 ≤ M`, …, `1 ≤ bN ≤ M`, +# 3. `l1 + … + 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 +# 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}(kind, smolyak_parameters::SmolyakParameters{B,M}) where {N,B,M} +# @argcheck N ≥ 1 +# Mp1 = M + 1 +# len = __smolyak_length(kind, Val(N), Val(B), M) +# first_block_length = nesting_total_length(Chebyshev, kind, 0) +# nesting_total_lengths = ntuple(bp1 -> nesting_total_length(Chebyshev, 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 +# 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 +# @inline highest_visited_index(::SmolyakIndices{N,H}) where {N,H} = H -Base.eltype(::Type{<:SmolyakIndices{N}}) where N = NTuple{N,Int} +# Base.eltype(::Type{<:SmolyakIndices{N}}) where N = NTuple{N,Int} -@inline Base.length(ι::SmolyakIndices) = ι.len +# @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{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 +# @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 -#### +# #### +# #### 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 +# struct SmolyakProduct{I<:SmolyakIndices,S<:Tuple,P} +# smolyak_indices::I +# sources::S +# product_kind::P +# @doc """ +# $(SIGNATURES) -Base.length(smolyak_product::SmolyakProduct) = length(smolyak_product.smolyak_indices) +# An iterator conceptually equivalent to -function Base.eltype(::Type{SmolyakProduct{I,S,P}}) where {I,S,P} - _product_type(P, fieldtypes(S)) -end +# ``` +# [prod(getindex.(sources, indices)) for indices in smolyak_indices] +# ``` -@inline function Base.iterate(smolyak_product::SmolyakProduct, state...) - (; smolyak_indices, sources, product_kind) = smolyak_product - itr = iterate(smolyak_indices, state...) - itr ≡ nothing && return nothing - indices, state′ = itr - _product(product_kind, map(getindex, sources, indices)), state′ -end +# using [`_product`](@ref) instead to account for derivatives. Detailed docs of the +# arguments are in [`SmolyakIndices`](@ref). + +# Caller should arrange the elements of `sources` in the correct order, see +# [`nested_extrema_indices`](@ref). Each element in `sources` should have at least +# `H` elements (cf type parameters of [`SmolyakIndices`](@ref)), this is not checked. +# """ +# function SmolyakProduct(smolyak_indices::I, sources::S, +# product_kind::P) where {N,I<:SmolyakIndices{N},S,P} +# @argcheck length(sources) == N +# new{I,S,P}(smolyak_indices, sources, product_kind) +# end +# end + +# Base.length(smolyak_product::SmolyakProduct) = length(smolyak_product.smolyak_indices) + +# function Base.eltype(::Type{SmolyakProduct{I,S,P}}) where {I,S,P} +# _product_type(P, fieldtypes(S)) +# end + +# @inline function Base.iterate(smolyak_product::SmolyakProduct, state...) +# (; smolyak_indices, sources, product_kind) = smolyak_product +# itr = iterate(smolyak_indices, state...) +# itr ≡ nothing && return nothing +# indices, state′ = itr +# _product(product_kind, map(getindex, sources, indices)), state′ +# end diff --git a/src/univariate.jl b/src/univariate.jl index 48e7459..87ec4d5 100644 --- a/src/univariate.jl +++ b/src/univariate.jl @@ -186,7 +186,8 @@ function adjust_basis(U::UnivariateBasis, Δ::Int) end end -function adjust_coefficients(θ1::AbstractVector, U1::UnivariateBasis{Chebyshev}, U2::UnivariateBasis{Chebyshev}) +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" diff --git a/test/runtests.jl b/test/runtests.jl index 09feb7a..42682f4 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -23,7 +23,7 @@ include("test_derivatives.jl") include("test_domains.jl") include("test_transformations.jl") include("test_univariate.jl") -# include("test_smolyak_traversal.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 diff --git a/test/test_smolyak_traversal.jl b/test/test_smolyak_traversal.jl index ef28c0e..ccf1f51 100644 --- a/test/test_smolyak_traversal.jl +++ b/test/test_smolyak_traversal.jl @@ -1,125 +1,26 @@ -using SpectralKit: nesting_total_length, nesting_block_length, SmolyakIndices, - __smolyak_length, SmolyakGridShuffle, SmolyakProduct - -"grids we test on" -GRIDS = (EndpointGrid(), InteriorGrid(), InteriorGrid2()) - -### -### block parameters -### - - -@testset "printing SmolyakParameters" begin - @test repr(SmolyakParameters(3, 2)) == "Smolyak parameters, ∑bᵢ ≤ 3, all bᵢ ≤ 2" -end - - -#### -#### blocks -#### - -@testset "block length" begin - for grid_kind in GRIDS - for b in 0:5 - nA = nesting_total_length(Chebyshev, grid_kind, b) - gA = grid(Chebyshev(grid_kind, nA)) - gB = grid(Chebyshev(grid_kind, nesting_total_length(Chebyshev, grid_kind, b + 1))) - - @test is_approximate_subset(gA, gB) - @test sum(b -> nesting_block_length(Chebyshev, grid_kind, b), 0:b) == nA - end - end -end - -""" -Collect shuffled indices for the given `grid_kind` for block indices `0, …, b`, returned as -a `Vector{Vector{Int}}`. For testing. -""" -function shuffled_indices_upto_b(grid_kind, b) - _grid(b) = grid(Chebyshev(grid_kind, nesting_total_length(Chebyshev, grid_kind, b))) - g0 = _grid(b) - indices = Vector{Vector{Int}}() - for b in b:(-1):1 - in_b = is_approximately_in(g0, _grid(b)) - notin_bm1 = .!is_approximately_in(g0, _grid(b - 1)) - mask = in_b .& notin_bm1 - push!(indices, findall(mask)) - end - push!(indices, [(length(g0) + 1) ÷ 2]) - reverse(indices) -end - -@testset "block shuffle" begin - @testset "endpoint" begin - for grid_kind in GRIDS - for b in 0:6 - len = nesting_total_length(Chebyshev, grid_kind, b) - ι = SmolyakGridShuffle(grid_kind, len) - @test length(ι) == len - @test @inferred eltype(ι) == Int - @test collect(ι) == reduce(vcat, shuffled_indices_upto_b(grid_kind, b)) - end - end - end -end - -#### -#### traversal -#### - -""" -Naive implementation of Smolyan index iteration, traversing a `CartesianIndices` and keeping -valid indexes. For testing/comparison. Returns a vector of `indexes => blocks` pairs. -""" -function smolyak_indices_check(grid_kind, N, B, M) - m = nesting_total_length(Chebyshev, grid_kind, M) - b_table = fill(M, m) - for b in (M-1):(-1):0 - b_table[1:nesting_total_length(Chebyshev, grid_kind, b)] .= b - end - T = NTuple{N,Int} - result = Vector{Pair{T,T}}() - for ι in CartesianIndices(ntuple(_ -> 1:m, N)) - ix = Tuple(ι) - blocks = map(i -> b_table[i], ix) - if sum(blocks) ≤ B - push!(result, ix => blocks) - end - end - result -end - -@testset "Smolyak indices" begin - for grid_kind in GRIDS - for B in 0:3 - for M in 0:B - for N in 1:4 - ι = SmolyakIndices{N}(grid_kind, SmolyakParameters(B, M)) - x1 = @inferred collect(ι) - x2 = first.(smolyak_indices_check(grid_kind, N, B, M)) - len = @inferred __smolyak_length(grid_kind, Val(N), Val(B), M) - @test x1 == x2 - @test len == length(x1) == length(ι) == length(x2) +@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), + map(SpectralKit.ChebyshevIterator, xs)) end end end end end -@testset "Smolyak product primitives" begin - for grid_kind in GRIDS - for B in 0:3 - for M in 0:B - for N in 1:4 - ι = SmolyakIndices{N}(grid_kind, SmolyakParameters(B, M)) - ℓ = nesting_total_length(Chebyshev, grid_kind, min(B,M)) - sources = ntuple(_ -> rand(SVector{ℓ, Float64}), Val(N)) - P = SmolyakProduct(ι, sources, nothing) - @test length(ι) == length(P) - @test eltype(P) == Float64 - for (i, p) in zip(ι, P) - @test prod(getindex.(sources, i)) ≈ p - end +@testset "__smolyak_length" begin + for kind in KINDS + for N in 1:5 + for total in 1:4 + for each in 1:4 + @test SpectralKit.__smolyak_length(Chebyshev(), kind, N, total, each) == + length(naive_smolyak_indices(Chebyshev(), kind, Val(N), total, each)) end end end diff --git a/test/utilities.jl b/test/utilities.jl index 411324e..b9d5ea0 100644 --- a/test/utilities.jl +++ b/test/utilities.jl @@ -123,3 +123,57 @@ function DD(f, x, n = 1; p = 10) central_fdm(p, n)(f, x) 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) + +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}) where N + i = 1 + reference = naive_smolyak_indices(family, kind, Val(N), total, each) + x, slack, remainders, states, cached, levels = SpectralKit.__smolyak_init(family, kind, total, f, itrs) + while true + @test slack + sum(levels) == total # simple sanity check + cs = map((itr, i) -> first(Iterators.drop(itr, i - 1)), + itrs, reference[i]) + expected_x, expected_cached... = reverse(accumulate(f, reverse(cs))) + @test x == expected_x + @test cached == expected_cached + next = SpectralKit.__smolyak_step(family, kind, each, f, itrs, slack, remainders, + states, cached, levels) + next ≡ nothing && break + (x, Δ, remainders, states, cached, levels) = next + slack += Δ + i += 1 + end + @test i == length(reference) # we used up all elements + nothing +end From f4c72831b88ca79cb63bd82c5da85c0aed4f0fcf Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tam=C3=A1s=20K=2E=20Papp?= Date: Sun, 23 Aug 2026 10:08:05 +0200 Subject: [PATCH 14/34] typo corrections --- src/derivatives.jl | 2 +- src/smolyak_traversal.jl | 7 ++++++- 2 files changed, 7 insertions(+), 2 deletions(-) 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/smolyak_traversal.jl b/src/smolyak_traversal.jl index 2d5d07c..c598d18 100644 --- a/src/smolyak_traversal.jl +++ b/src/smolyak_traversal.jl @@ -18,7 +18,12 @@ $(SIGNATURES) → x, Δ, remainders′, states′, cached′, levels′ Step through Smolyak indices of iterators. -`itrs` yield the `xs`, `(x, cached...)` is equivalent to `reverse(accumulate(reverse(...)))` applied to t +`itrs` yield the `xs`. + +`(x, cached...)` is equivalent to `reverse(accumulate(reverse(...)))` applied to the `xs`. + +`Δ` is the change in `slack`. Design note: easier to apply recursively than `slack`, +caller should make the adjustment. """ function __smolyak_step(family, kind, each::Int, f::F, itrs::NTuple{N}, slack::Int, remainders::NTuple{N}, states::NTuple{N}, From c2894feacf22fa8266241dbcbe7473af468499b6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tam=C3=A1s=20K=2E=20Papp?= Date: Mon, 24 Aug 2026 17:27:33 +0200 Subject: [PATCH 15/34] rework basics of Smolyak API --- src/SpectralKit.jl | 2 +- src/smolyak_api.jl | 348 ++++++++++++++++++++------------- src/smolyak_traversal.jl | 54 ++--- src/univariate.jl | 31 ++- test/runtests.jl | 2 +- test/test_smolyak.jl | 117 +++++------ test/test_smolyak_traversal.jl | 2 +- test/test_univariate.jl | 2 +- test/utilities.jl | 15 +- 9 files changed, 333 insertions(+), 240 deletions(-) diff --git a/src/SpectralKit.jl b/src/SpectralKit.jl index 271fa3f..074d0ee 100644 --- a/src/SpectralKit.jl +++ b/src/SpectralKit.jl @@ -13,7 +13,7 @@ include("transformations.jl") include("generic_api.jl") include("univariate.jl") include("smolyak_traversal.jl") -# include("smolyak_api.jl") +include("smolyak_api.jl") # include("experimental.jl") # experimental code is not part of the API, see its module docstring end # module diff --git a/src/smolyak_api.jl b/src/smolyak_api.jl index bb4c110..6fced9a 100644 --- a/src/smolyak_api.jl +++ b/src/smolyak_api.jl @@ -2,26 +2,86 @@ ##### Smolyak bases ##### -export smolyak_basis +export SmolyakLevel, smolyak_basis -struct SmolyakBasis{I<:SmolyakIndices,U<:UnivariateBasis} <: MultivariateBasis - smolyak_indices::I - univariate_parent::U +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 ≥ 1 + @argcheck each ≥ 1 + if each > total + @warn "‘each’ normalized to ‘total’" each total + each = total + end + new(total, each) + end 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) +""" +$(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, itrs, g, + state = nothing) + if state ≡ nothing + accum, state... = __smolyak_init(family, kind, level.total, f, itrs) + g(accum), (accum, state...) + else + result = __smolyak_step(family, kind, level.each, f, itrs, + state...) + if result ≡ nothing + nothing + else + slack = state[2] + accum, Δ, rest... = result + g(accum), (accum, slack + Δ, rest...) + end + end end -Base.length(basis::SmolyakBasis{<:SmolyakIndices{N}}) where N = N +function Base.show(io::IO, level::SmolyakLevel) + (; total, each) = level + print(io, "Smolyak parameters, ∑ℓᵢ ≤ $(total), all ℓᵢ ≤ $(each)") +end -function Base.getindex(basis::SmolyakBasis, i::Int) - @argcheck 1 ≤ i ≤ length(basis) BoundsError(basis, i) - basis.univariate_parent +@concrete struct SmolyakBasis <: MultivariateBasis + family + kind + domain_transformations + level 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) @@ -29,19 +89,20 @@ Create a sparse Smolyak basis. # Arguments -- `univariate_family`: should be a callable that takes a `grid_kind` and a `dimension` - parameter, eg `Chebyshev`. +- `family`: univariate function family, eg `Chebyshev`. -- `grid_kind`: the grid kind, eg `InteriorGrid()` etc. +- `kind`: the grid kind, eg `Interior()` or `Endpoints()`. -- `smolyak_parameters`: the Smolyak grid specification parameters, see - [`SmolyakParameters`](@ref). +- `domain_transformations` + +- `smolyak_level`: the Smolyak level specificaion, see [`SmolyakLevel`](@ref). - `N`: the dimension. wrapped in a `Val` for type stability, a convenience constructor also takes integers. ## Example +FIXME these examples need to be updated ```jldoctest julia> basis = smolyak_basis(Chebyshev, InteriorGrid(), SmolyakParameters(3), 2) Sparse multivariate basis on ℝ² @@ -60,149 +121,162 @@ julia> domain(basis) *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) +function smolyak_basis(univariate_family, kind, domain_transformations::Tuple, level::SmolyakLevel) + SmolyakBasis(univariate_family, kind, domain_transformations, level) 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)) +function domain(smolyak_basis::SmolyakBasis) + map(domain, smolyak_basis.domain_transformations) end -function domain(smolyak_basis::SmolyakBasis{<:SmolyakIndices{N}}) where {N} - D = domain(smolyak_basis.univariate_parent) - coordinate_domains(Val(N), D) +function dimension(smolyak_basis::SmolyakBasis) + (; family, kind, domain_transformations, level) = smolyak_basis + N = length(domain_transformations) + __smolyak_length(family, kind, N, level.total, level.each) 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) +struct SmolyakBasisAt{I,P,F,K,L<:SmolyakLevel} + family::F + kind::K + level::L + itrs::I + product_kind::P 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) +function Base.eltype(::Type{<:SmolyakBasisAt{I,P}}) where {I,P} + _product_type(P, map(eltype, fieldtypes(I))) 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) +function Base.length(itr::SmolyakBasisAt) + (; family, kind, level, itrs, product_kind) = itr + N = length(itrs) + __smolyak_length(family, kind, N, level.total, level.each) end -struct SmolyakGridIterator{T,I,S} - smolyak_indices::I - sources::S +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, + (a, b) -> isempty(b) ? (a,) : (a * first(b), b...), + itrs, first, state) 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) +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) -> basis_at(univariate_basis(family, kind, d, level.each), x), + x, domain_transformations) + SmolyakBasisAt(family, kind, level, itrs, nothing) 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′ +function basis_at(smolyak_basis::SmolyakBasis, x::SVector) + basis_at(smolyak_basis, Tuple(x)) 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 +# function basis_at(smolyak_basis::SmolyakBasis, Dx::∂CoordinateExpansion) +# (; family, kind, domain_transformations, level) = smolyak_basis +# (; ∂D, x) = Dx +# itrs = map((x, d) -> basis_at(univariate_basis(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 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)) +function Base.eltype(::Type{<:SmolyakGrid{I}}) where {I} + Tuple{map(eltype, fieldtypes(I))...} 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: +function Base.length(itr::SmolyakGrid) + __smolyak_length(itr.family, itr.kind, length(itr.itrs), itr.level.total, itr.level.each) +end -- 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) +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 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)) +function grid(::Type{T}, smolyak_basis::SmolyakBasis) where {T<:AbstractFloat} + (; family, kind, domain_transformations, level) = smolyak_basis + (; each) = level + itrs = map(d -> grid(T, univariate_basis(family, kind, d, each)), domain_transformations) + SmolyakGrid(family, kind, level, itrs) 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 c598d18..60efd84 100644 --- a/src/smolyak_traversal.jl +++ b/src/smolyak_traversal.jl @@ -2,69 +2,73 @@ ##### Smolyak implementation details ##### -function __smolyak_init(family, kind, total::Int, f::F, itrs::NTuple{N}) where {F,N} +""" +$(SIGNATURES) → accum, slack, remainders, states, levels + +Initialize the state for [`__smolyak_step`](@ref), states are documented there. +""" +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(iterate, itrs) states = map(last, itr_results) levels = ntuple(_ -> 0, Val(N)) - x, cached... = reverse(accumulate(f, reverse(map(first, itr_results)))) - x, slack, remainders, states, cached, levels + accum = foldr(f, map(first, itr_results); init = ()) + accum, slack, remainders, states, levels end """ -$(SIGNATURES) → x, Δ, remainders′, states′, cached′, levels′ +$(SIGNATURES) → accum, Δ, remainders′, states′, levels′ Step through Smolyak indices of iterators. -`itrs` yield the `xs`. +`itrs` yield the `xs`. **The iterators should be stateless.** + +`accum` is equivalent to `foldr(f, xs)`. `f` should map a value and a tuple to a tuple +of one more element. -`(x, cached...)` is equivalent to `reverse(accumulate(reverse(...)))` applied to the `xs`. +`slack` is the `total - sum(levels)`. `Δ` is the change in `slack`. Design note: easier +to apply recursively than `slack`, caller should make the adjustment. -`Δ` is the change in `slack`. Design note: easier to apply recursively than `slack`, -caller should make the adjustment. +`remainders` contains the count of elements left in each level before we move to a +different combination. + +`states` are states of iterators. `levels` are the levels currently visited. """ -function __smolyak_step(family, kind, each::Int, f::F, itrs::NTuple{N}, - slack::Int, remainders::NTuple{N}, states::NTuple{N}, - cached, levels::NTuple{N}) where {F,N} - @assert length(cached) == N - 1 +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 - s1, sτ... = states + a1, aτ... = accum r1, rτ... = remainders + s1, sτ... = states l1, lτ... = levels if r1 > 0 # step within block x1, s1′ = iterate(I1, s1) - (f(x1, cached), # product + (f(x1, aτ), 0, # no change in slack (r1 - 1, rτ...), # one less element in 1 (s1′, sτ...), # step iterator - cached, levels) elseif l1 < each && slack > 0 # next block, same tail x1, s1′ = iterate(I1, s1) - (f(x1, cached), # step to next block + (f(x1, aτ), -1, # decrease slack (block_length(family, kind, l1 + 1) - 1, rτ...), # remaining elements: all in block - (s1′, sτ...), # step state 1 - cached, (l1+1, lτ...)) # next level elseif N == 1 nothing # done with iteration else # go into tail - _, cτ... = cached - next = __smolyak_step(family, kind, each, f, Iτ, - slack + l1, rτ, sτ, cτ, lτ) + next = __smolyak_step(family, kind, each, f, Iτ, aτ, slack + l1, rτ, sτ, lτ) next ≡ nothing && return nothing - xτ, Δτ, rτ′, sτ′, cτ′, lτ′ = next + aτ′, Δτ, rτ′, sτ′, lτ′ = next x1, s1 = iterate(I1) - cached′ = (xτ, cτ′...) - (f(x1, cached′), + (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 - cached′, # cache updated as above (0, lτ′...)) # back to level 0 here end end diff --git a/src/univariate.jl b/src/univariate.jl index 87ec4d5..a8b1e5e 100644 --- a/src/univariate.jl +++ b/src/univariate.jl @@ -163,18 +163,29 @@ function basis_at(U::UnivariateBasis{Chebyshev}, x::Scalar) Iterators.take(ChebyshevIterator(transform_to(PM1(), U.domain_transformation, x)), dimension(U)) end -function grid(::Type{T}, U::UnivariateBasis{Chebyshev}) where {T <: AbstractFloat} +@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, level) = U N = grid_length(family, kind, level) - if kind ≡ Interior() - N̂ = N + 2 - else - @assert kind ≡ Endpoints() - N̂ = N - end - (transform_from(PM1(), U.domain_transformation, - _chebyshev_extremum(T, _chebyshev_extremum_shuffle(kind, i, level), N̂)) - for i in 1:N) + ChebyshevGrid{T}(kind, domain_transformation, level, N, + kind ≡ Interior() ? N̂ = N + 2 : N) end function adjust_basis(U::UnivariateBasis, Δ::Int) diff --git a/test/runtests.jl b/test/runtests.jl index 42682f4..6b4fab4 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -24,6 +24,6 @@ include("test_domains.jl") include("test_transformations.jl") include("test_univariate.jl") include("test_smolyak_traversal.jl") -# include("test_smolyak.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 diff --git a/test/test_smolyak.jl b/test/test_smolyak.jl index c1863a1..18a5612 100644 --- a/test/test_smolyak.jl +++ b/test/test_smolyak.jl @@ -1,19 +1,21 @@ -using SpectralKit, Test -using SpectralKit: PM1, ∂Expansion - #### #### api #### +@testset "printing SmolyakLevels" begin + @test repr(SmolyakLevel(; total = 3, each = 2)) == "Smolyak parameters, ∑ℓᵢ ≤ 3, all ℓᵢ ≤ 2" +end + @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) + @test_logs (:warn, "‘each’ normalized to ‘total’") SmolyakLevel(total = 2, each = 4) 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(lower = 2, upper = 3), # approximation should be exact + BoundedLinear(lower = 3.0, upper = 4.5)) + basis = smolyak_basis(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,25 +27,28 @@ 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 + # 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 = smolyak_basis(Chebyshev(), Interior(), (t, t), + SmolyakLevel(total = 3)) y = SVector(0.4, 0.7) θ = randn(dimension(basis)) @inferred linear_combination(basis, θ, y) @@ -54,43 +59,43 @@ end ### augment coefficients ### -@testset "Smolyak augment coefficients" begin - basis1 = smolyak_basis(Chebyshev, InteriorGrid(), SmolyakParameters(2, 2), 2) - θ1 = randn(dimension(basis1)) +# @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) +# # 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) +# # 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 +# # 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 - end - end - end - 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 +# end +# end +# end +# end +# end diff --git a/test/test_smolyak_traversal.jl b/test/test_smolyak_traversal.jl index ccf1f51..5ee1ac0 100644 --- a/test/test_smolyak_traversal.jl +++ b/test/test_smolyak_traversal.jl @@ -6,7 +6,7 @@ for each in 1:4 xs = XS[1:N] test_smolyak_step(Chebyshev(), Interior(), total, each, - (a, b) -> isempty(b) ? a : a * first(b), + (a, b) -> isempty(b) ? @show((a,)) : @show((a * first(b), b...)), map(SpectralKit.ChebyshevIterator, xs)) end end diff --git a/test/test_univariate.jl b/test/test_univariate.jl index 416cbac..c13eb66 100644 --- a/test/test_univariate.jl +++ b/test/test_univariate.jl @@ -129,7 +129,7 @@ end x = rand_in_domain(basis) y = f(D(x)) for i in 0:N - @test y[i] ≈ DD(f, x, i) atol = 1e-6 + @test y[i] ≈ DD(f, x, i) rtol = 1e-4 atol = 1e-4 end end end diff --git a/test/utilities.jl b/test/utilities.jl index b9d5ea0..f32ef0e 100644 --- a/test/utilities.jl +++ b/test/utilities.jl @@ -156,21 +156,20 @@ $(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}) where N +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) - x, slack, remainders, states, cached, levels = SpectralKit.__smolyak_init(family, kind, total, f, itrs) + (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((itr, i) -> first(Iterators.drop(itr, i - 1)), itrs, reference[i]) - expected_x, expected_cached... = reverse(accumulate(f, reverse(cs))) - @test x == expected_x - @test cached == expected_cached - next = SpectralKit.__smolyak_step(family, kind, each, f, itrs, slack, remainders, - states, cached, levels) + @test accum == foldr(f, cs; init = ()) + next = SpectralKit.__smolyak_step(family, kind, each, f, itrs, + accum, slack, remainders, states, levels) next ≡ nothing && break - (x, Δ, remainders, states, cached, levels) = next + (accum, Δ, remainders, states, levels) = next slack += Δ i += 1 end From 5ba39f5b01165d698ab3e8cd47158ec7dc40d69c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tam=C3=A1s=20K=2E=20Papp?= Date: Tue, 25 Aug 2026 17:06:23 +0200 Subject: [PATCH 16/34] make iterator allocation-free --- src/smolyak_api.jl | 6 +- src/smolyak_traversal.jl | 113 +++++++++++++++++++++++++-------- test/test_smolyak_traversal.jl | 19 +++++- test/utilities.jl | 16 +++++ 4 files changed, 122 insertions(+), 32 deletions(-) diff --git a/src/smolyak_api.jl b/src/smolyak_api.jl index 6fced9a..2325f77 100644 --- a/src/smolyak_api.jl +++ b/src/smolyak_api.jl @@ -132,7 +132,7 @@ end function dimension(smolyak_basis::SmolyakBasis) (; family, kind, domain_transformations, level) = smolyak_basis N = length(domain_transformations) - __smolyak_length(family, kind, N, level.total, level.each) + __smolyak_length(family, kind, Val(N), level.total, level.each) end struct SmolyakBasisAt{I,P,F,K,L<:SmolyakLevel} @@ -150,7 +150,7 @@ end function Base.length(itr::SmolyakBasisAt) (; family, kind, level, itrs, product_kind) = itr N = length(itrs) - __smolyak_length(family, kind, N, level.total, level.each) + __smolyak_length(family, kind, Val(N), level.total, level.each) end function Base.iterate(itr::SmolyakBasisAt, state = nothing) @@ -193,7 +193,7 @@ function Base.eltype(::Type{<:SmolyakGrid{I}}) where {I} end function Base.length(itr::SmolyakGrid) - __smolyak_length(itr.family, itr.kind, length(itr.itrs), itr.level.total, itr.level.each) + __smolyak_length(itr.family, itr.kind, Val(length(itr.itrs)), itr.level.total, itr.level.each) end function Base.iterate(itr::SmolyakGrid, state = nothing) diff --git a/src/smolyak_traversal.jl b/src/smolyak_traversal.jl index 60efd84..f78a1aa 100644 --- a/src/smolyak_traversal.jl +++ b/src/smolyak_traversal.jl @@ -73,6 +73,93 @@ function __smolyak_step(family, kind, each::Int, f::F, itrs::NTuple{N,Any}, end end +struct NonIncreasingSmolyakLevels{N} + total::Int + each::Int + @doc """ + $(SIGNATURES) → itr + + An iterable which yields `ℓ::Ntuple{N,Int}`, with the following properties: + + 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 + +Base.eltype(::Type{NonIncreasingSmolyakLevels{N}}) where N = NTuple{N,Int} + +Base.IteratorSize(::Type{<:NonIncreasingSmolyakLevels}) = Base.SizeUnknown() + +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 + # increment first index; stop iteration if this is > each + (i1 + 1, i1 + 1, ntuple(_ -> 0, Val(length(iτ)))...) + end +end + +function Base.iterate(itr::NonIncreasingSmolyakLevels{N}, state = nothing) where N + (; total, each) = itr + if state ≡ nothing + indices = ntuple(_ -> 0, N) + indices, (0, indices...) + else + state′ = __step_noninc(total, each, state...) + Σ′, indices′... = state′ + if first(indices′) ≤ each + indices′, state′ + else + nothing + end + end +end + +""" +$(SIGNATURES) + +Calculate the length of a [`SmolyakIndices`](@ref) iterator. Argument as in the latter. +""" +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 + end + L += p * C + end + L +end + + #### #### index traversal #### @@ -133,32 +220,6 @@ end # end # end -""" -$(SIGNATURES) - -Calculate the length of a [`SmolyakIndices`](@ref) iterator. Argument as in the latter. -""" -function __smolyak_length(family, kind, N::Int, total::Int, each::Int) - # implicit assumption: each ≤ total, enforced by the SmolyakParameters constructor - _bl(b) = block_length(family, kind, b) - c = zeros(Int, total + 1) # indexed as 0, …, total - each = min(each, total) - for b in 0:each - c[b + 1] = _bl(b) - end - for n in 2:N - for b in total:(-1):0 # blocks with indices that sum to b - s = 0 - for a in 0:min(b, each) - s += _bl(a) * c[b - a + 1] - end - # can safely overwrite since they will not be used again for n + 1 - c[b + 1] = s - end - end - sum(c) -end - # """ # $(TYPEDEF) diff --git a/test/test_smolyak_traversal.jl b/test/test_smolyak_traversal.jl index 5ee1ac0..c977ccd 100644 --- a/test/test_smolyak_traversal.jl +++ b/test/test_smolyak_traversal.jl @@ -6,7 +6,7 @@ for each in 1:4 xs = XS[1:N] test_smolyak_step(Chebyshev(), Interior(), total, each, - (a, b) -> isempty(b) ? @show((a,)) : @show((a * first(b), b...)), + (a, b) -> isempty(b) ? (a,) : (a * first(b), b...), map(SpectralKit.ChebyshevIterator, xs)) end end @@ -14,12 +14,25 @@ end 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_length" begin for kind in KINDS for N in 1:5 for total in 1:4 - for each in 1:4 - @test SpectralKit.__smolyak_length(Chebyshev(), kind, N, total, each) == + for each in 1:total + # kind, N, total, each = Endpoints(), 4, 2, 1 + @test @inferred(SpectralKit.__smolyak_length(Chebyshev(), kind, Val(N), total, each)) == length(naive_smolyak_indices(Chebyshev(), kind, Val(N), total, each)) end end diff --git a/test/utilities.jl b/test/utilities.jl index f32ef0e..20cda26 100644 --- a/test/utilities.jl +++ b/test/utilities.jl @@ -176,3 +176,19 @@ function test_smolyak_step(family, kind, total, each, f, itrs::NTuple{N,Any}) wh @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 From 6b7c0df86fe2fed67286a796a067d2acdb53a9b1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tam=C3=A1s=20K=2E=20Papp?= Date: Wed, 26 Aug 2026 15:20:04 +0200 Subject: [PATCH 17/34] remove adjust_basis, rely on fields --- src/generic_api.jl | 12 ++---------- src/univariate.jl | 9 --------- test/test_univariate.jl | 7 +++---- 3 files changed, 5 insertions(+), 23 deletions(-) diff --git a/src/generic_api.jl b/src/generic_api.jl index 7a76948..64dcfce 100644 --- a/src/generic_api.jl +++ b/src/generic_api.jl @@ -3,7 +3,7 @@ ##### export is_function_basis, dimension, basis_at, linear_combination, grid, collocation_matrix, - adjust_basis, adjust_coefficients + adjust_coefficients """ $(TYPEDEF) @@ -191,14 +191,6 @@ function collocation_matrix(basis, x = grid(basis)) C end -""" -`$(FUNCTIONNAME)(basis, Δ)` - -Make `basis` thicker (`Δ > 0`) or thinner (`Δ < 0`). Return `nothing` if this is not -possible. -""" -function adjust_basis end - """ `$(FUNCTIONNAME)(θ1, basis1, basis2)` @@ -209,7 +201,7 @@ linear_combination(basis1, θ1, x) ≈ linear_combination(basis2, θ2, x) for any `x` in the domain. In practice this means either padding with zeros, or truncating the coefficients. -Throw a `ArgumentError` if this is not possible. See [`adjust_basis`](@ref). +Throw an error if this is not possible. """ function adjust_coefficients end diff --git a/src/univariate.jl b/src/univariate.jl index a8b1e5e..72d57ac 100644 --- a/src/univariate.jl +++ b/src/univariate.jl @@ -188,15 +188,6 @@ function grid(::Type{T}, kind ≡ Interior() ? N̂ = N + 2 : N) end -function adjust_basis(U::UnivariateBasis, Δ::Int) - level′ = U.level + Δ - if level′ > 0 - UnivariateBasis(U.family, U.kind, U.domain_transformation, level′) - else - nothing - end -end - function adjust_coefficients(θ1::AbstractVector, U1::UnivariateBasis{F}, U2::UnivariateBasis{F}) where F d1 = dimension(U1) diff --git a/test/test_univariate.jl b/test/test_univariate.jl index c13eb66..b44509d 100644 --- a/test/test_univariate.jl +++ b/test/test_univariate.jl @@ -1,5 +1,5 @@ #### -#### Chebyshev on [-1,1] +#### Chebyshev building blocks and univariate bases #### """ @@ -99,14 +99,14 @@ end end end -@testset "Chebysev adjusted basis" begin +@testset "Chebysev adjusted coefficients" begin transformation = SemiInfRational(; endpoint = 3.0, scale = 7.0) for kind in KINDS level0 = 3 basis0 = univariate_basis(Chebyshev(), kind, transformation, level0) θ0 = randn(dimension(basis0)) for Δ in 1:4 - basis = adjust_basis(basis0, Δ) + basis = @set basis0.level = level0 + Δ θ = adjust_coefficients(θ0, basis0, basis) @test basis.level == basis0.level + Δ for _ in 1:10 @@ -114,7 +114,6 @@ end @test linear_combination(basis0, θ0, x) ≈ linear_combination(basis, θ, x) end end - @test adjust_basis(basis0, -3) ≡ nothing end end From 86029c6dc770a6ff3bf72feb3f144a8fc33b43ed Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tam=C3=A1s=20K=2E=20Papp?= Date: Wed, 26 Aug 2026 16:16:23 +0200 Subject: [PATCH 18/34] use type constructors instead of functions --- src/smolyak_api.jl | 141 ++++++++++++++++++++++++++-------------- src/univariate.jl | 38 +++++------ test/Project.toml | 1 + test/runtests.jl | 2 +- test/test_smolyak.jl | 18 ++--- test/test_univariate.jl | 24 +++---- 6 files changed, 134 insertions(+), 90 deletions(-) diff --git a/src/smolyak_api.jl b/src/smolyak_api.jl index 2325f77..80bcb27 100644 --- a/src/smolyak_api.jl +++ b/src/smolyak_api.jl @@ -2,7 +2,7 @@ ##### Smolyak bases ##### -export SmolyakLevel, smolyak_basis +export SmolyakLevel, SmolyakBasis struct SmolyakLevel total::Int @@ -62,68 +62,70 @@ function Base.show(io::IO, level::SmolyakLevel) print(io, "Smolyak parameters, ∑ℓᵢ ≤ $(total), all ℓᵢ ≤ $(each)") end -@concrete struct SmolyakBasis <: MultivariateBasis - family - kind - domain_transformations - level +struct SmolyakBasis{F,K,D} <: MultivariateBasis + family::F + kind::K + domain_transformations::D + level::SmolyakLevel + @doc """ + $(SIGNATURES) + """ + function SmolyakBasis(family::F, kind::K, domain_transformations::D, + level::SmolyakLevel) where {F,K,D<:Tuple} + new{F,K,D}(family, kind, domain_transformations, level) + end 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.show(io::IO, basis::SmolyakBasis) where N + (; family, kind, domain_transformations, level) = basis + lead = "SmolyakBasis(" + next = ",\n" * ' '^length(lead) + print(io, "SmolyakBasis(", family, next, kind, next, domain_transformations, next, + level, + ") # dimension: ", + dimension(basis)) +end -# function Base.getindex(basis::SmolyakBasis, i::Int) -# @argcheck 1 ≤ i ≤ length(basis) BoundsError(basis, i) -# basis.univariate_parent -# end -""" -$(SIGNATURES) +# """ +# $(SIGNATURES) -Create a sparse Smolyak basis. +# Create a sparse Smolyak basis. -# Arguments +# # Arguments -- `family`: univariate function family, eg `Chebyshev`. +# - `family`: univariate function family, eg `Chebyshev`. -- `kind`: the grid kind, eg `Interior()` or `Endpoints()`. +# - `kind`: the grid kind, eg `Interior()` or `Endpoints()`. -- `domain_transformations` +# - `domain_transformations` -- `smolyak_level`: the Smolyak level specificaion, see [`SmolyakLevel`](@ref). +# - `smolyak_level`: the Smolyak level specificaion, see [`SmolyakLevel`](@ref). -- `N`: the dimension. wrapped in a `Val` for type stability, a convenience constructor also - takes integers. +# - `N`: the dimension. wrapped in a `Val` for type stability, a convenience constructor also +# takes integers. -## Example +# ## 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 +# 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> dimension(basis) +# 81 -julia> domain(basis) -[-1,1]² -``` +# julia> domain(basis) +# [-1,1]² +# ``` -## Properties +# ## Properties -*Grids nest*: increasing arguments of `SmolyakParameters` result in a refined grid that -contains points of the cruder grid. -""" -function smolyak_basis(univariate_family, kind, domain_transformations::Tuple, level::SmolyakLevel) - SmolyakBasis(univariate_family, kind, domain_transformations, level) -end +# *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) @@ -164,7 +166,7 @@ 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) -> basis_at(univariate_basis(family, kind, d, level.each), x), + itrs = map((x, d) -> basis_at(UnivariateBasis(family, kind, d, level.each), x), x, domain_transformations) SmolyakBasisAt(family, kind, level, itrs, nothing) end @@ -176,7 +178,7 @@ 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(univariate_basis(family, kind, d, each), x), +# itrs = map((x, d) -> basis_at(UnivariateBasis(family, kind, d, each), x), # x, domain_transformations) # BasisAt(family, kind, total, each, itrs, ∂D) # end @@ -204,10 +206,51 @@ end function grid(::Type{T}, smolyak_basis::SmolyakBasis) where {T<:AbstractFloat} (; family, kind, domain_transformations, level) = smolyak_basis (; each) = level - itrs = map(d -> grid(T, univariate_basis(family, kind, d, each)), domain_transformations) + itrs = map(d -> grid(T, UnivariateBasis(family, kind, d, each)), domain_transformations) SmolyakGrid(family, kind, level, itrs) end +@concrete struct SmolyakIndices + family + kind + level + itrs +end + +""" +$(SIGNATURES) + +Iterate through indices of +""" +function smolyak_indices(basis::SmolyakBasis) + (; family, kind, domain_transformations, level) = basis + itr1 = 1:level.each + itrs = ntuple(_ -> itr1, 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, smolyak_indices(basis1)) + θ[ι] = x + end + z = zero(T) + [get(θ, ι, z) for ι in smolyak_indices(basis2)] +end + # """ # $(SIGNATURES) diff --git a/src/univariate.jl b/src/univariate.jl index 72d57ac..9c72712 100644 --- a/src/univariate.jl +++ b/src/univariate.jl @@ -2,7 +2,7 @@ #### Chebyshev polynomials on [-1,1] #### -export Chebyshev, Endpoints, Interior, univariate_basis +export Chebyshev, Endpoints, Interior, UnivariateBasis #### #### generic building blocks @@ -135,24 +135,24 @@ end """ Implementation of univariate bases. Not part of the API. """ -@concrete struct UnivariateBasis <: FunctionBasis - family - kind - domain_transformation - level -end - -""" -$(SIGNATURES) - -Univariate basis from `family`, using the given `kind`. - -`level` is an integer, starting from `0`, specifying the number of *blocks* used to -build the grid. -""" -function univariate_basis(family, kind, domain_transformation, level) - @argcheck level ≥ 1 - UnivariateBasis(family, kind, domain_transformation, level) +struct UnivariateBasis{F,K,D} <: FunctionBasis + family::F + kind::K + domain_transformation::D + level::Int + @doc """ + $(SIGNATURES) + + Univariate basis from `family`, using the given `kind`. + + `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) where {F,K,D} + @argcheck level ≥ 0 + new{F,K,D}(family, kind, domain_transformation, level) + end end domain(U::UnivariateBasis) = domain(U.domain_transformation) diff --git a/test/Project.toml b/test/Project.toml index 700d5c6..956fdc9 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -1,4 +1,5 @@ [deps] +Accessors = "7d9f7c33-5ae7-4f3b-8dc6-eff91059b697" Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595" BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf" DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" diff --git a/test/runtests.jl b/test/runtests.jl index 6b4fab4..2a169e8 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -14,7 +14,7 @@ import SpectralKit as SK # Aqua.test_all(SpectralKit) # end -using DocStringExtensions, StaticArrays, BenchmarkTools, FiniteDifferences +using Accessors, DocStringExtensions, StaticArrays, BenchmarkTools, FiniteDifferences include("utilities.jl") diff --git a/test/test_smolyak.jl b/test/test_smolyak.jl index 18a5612..a071b21 100644 --- a/test/test_smolyak.jl +++ b/test/test_smolyak.jl @@ -14,7 +14,7 @@ end f(x) = (x[1] - 3) * (x[2] + 5) # linear function, just a sanity check transformations = (BoundedLinear(lower = 2, upper = 3), # approximation should be exact BoundedLinear(lower = 3.0, upper = 4.5)) - basis = smolyak_basis(Chebyshev(), Interior(), transformations, SmolyakLevel(total = 2)) + basis = SmolyakBasis(Chebyshev(), Interior(), transformations, SmolyakLevel(total = 2)) @test @inferred(domain(basis)) ≡ domain.(transformations) g = grid(Float64, basis) iterator_sanity_checks(g) @@ -47,8 +47,8 @@ end @testset "Smolyak API allocations" begin t = SemiInfRational() - basis = smolyak_basis(Chebyshev(), Interior(), (t, t), - SmolyakLevel(total = 3)) + basis = SmolyakBasis(Chebyshev(), Interior(), (t, t), + SmolyakLevel(total = 3)) y = SVector(0.4, 0.7) θ = randn(dimension(basis)) @inferred linear_combination(basis, θ, y) @@ -60,21 +60,21 @@ end ### # @testset "Smolyak augment coefficients" begin -# basis1 = smolyak_basis(Chebyshev, InteriorGrid(), SmolyakParameters(2, 2), 2) +# basis1 = SmolyakBasis(Chebyshev, InteriorGrid(), SmolyakParameters(2, 2), 2) # θ1 = randn(dimension(basis1)) # # grid ≠ -# basis2_G = smolyak_basis(Chebyshev, EndpointGrid(), SmolyakParameters(2, 3), 2) +# 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 = smolyak_basis(Chebyshev, InteriorGrid(), SmolyakParameters(2, 1), 2) +# 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 = smolyak_basis(Chebyshev, InteriorGrid(), SmolyakParameters(3, 2), 2) +# 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) @@ -90,8 +90,8 @@ end # 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) +# 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 diff --git a/test/test_univariate.jl b/test/test_univariate.jl index b44509d..e0d4d6e 100644 --- a/test/test_univariate.jl +++ b/test/test_univariate.jl @@ -47,9 +47,9 @@ end @testset "Chebyshev grid" begin for kind in KINDS previous = Float64[] - for level in 1:7 - b = univariate_basis(Chebyshev(), kind, - BoundedLinear(; lower = 1.0, upper = 3.0), level) + for level in 0:7 + b = UnivariateBasis(Chebyshev(), kind, + BoundedLinear(; lower = 1.0, upper = 3.0), level) g = collect(grid(b)) @test is_approximate_subset(previous, g) g = previous @@ -59,10 +59,10 @@ end @testset "Chebyshev basics" begin transformation = BoundedLinear(; lower = 1.0, upper = 3.0) - @test_throws ArgumentError univariate_basis(Chebyshev(), Endpoints(), transformation, 0) + @test_throws ArgumentError UnivariateBasis(Chebyshev(), Endpoints(), transformation, -1) for kind in KINDS - for level in 1:5 - basis = univariate_basis(Chebyshev(), kind, transformation, level) + 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) @@ -103,7 +103,7 @@ end transformation = SemiInfRational(; endpoint = 3.0, scale = 7.0) for kind in KINDS level0 = 3 - basis0 = univariate_basis(Chebyshev(), kind, transformation, level0) + basis0 = UnivariateBasis(Chebyshev(), kind, transformation, level0) θ0 = randn(dimension(basis0)) for Δ in 1:4 basis = @set basis0.level = level0 + Δ @@ -121,7 +121,7 @@ end for (transformation, N) in ((BoundedLinear(lower = -2, upper = 3), 5), (SemiInfRational(endpoint = 0.7, scale = 0.3), 1), (InfRational(center = 0.4, scale = 0.9), 1)) - basis = univariate_basis(Chebyshev(), Interior(), transformation, 3) + basis = UnivariateBasis(Chebyshev(), Interior(), transformation, 3) D = 𝑑^Val(N) f = linear_combination(basis, randn(dimension(basis))) for _ in 1:50 @@ -142,8 +142,8 @@ end x_minf = 𝑑(-Inf) @testset "SemiInfRational endpoints continuity" begin - basis = univariate_basis(Chebyshev(), Interior(), - SemiInfRational(; endpoint = 2.3, scale = 0.7), 3) + 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) @@ -153,8 +153,8 @@ end end @testset "InfRational endpoints continuity" begin - basis = univariate_basis(Chebyshev(), Interior(), - InfRational(; center = 2.3, scale = 0.7), 3) + 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) From db6f6c28aa21cb52082edd49d84d75d4462098a0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tam=C3=A1s=20K=2E=20Papp?= Date: Wed, 26 Aug 2026 16:21:50 +0200 Subject: [PATCH 19/34] show method for UnivariateBasis, fix for SmolyakBasis --- src/smolyak_api.jl | 6 ++---- src/univariate.jl | 8 ++++++++ 2 files changed, 10 insertions(+), 4 deletions(-) diff --git a/src/smolyak_api.jl b/src/smolyak_api.jl index 80bcb27..fcb0c75 100644 --- a/src/smolyak_api.jl +++ b/src/smolyak_api.jl @@ -76,14 +76,12 @@ struct SmolyakBasis{F,K,D} <: MultivariateBasis end end -function Base.show(io::IO, basis::SmolyakBasis) where N +function Base.show(io::IO, basis::SmolyakBasis) (; family, kind, domain_transformations, level) = basis lead = "SmolyakBasis(" next = ",\n" * ' '^length(lead) print(io, "SmolyakBasis(", family, next, kind, next, domain_transformations, next, - level, - ") # dimension: ", - dimension(basis)) + level, ") # dimension: ", dimension(basis)) end diff --git a/src/univariate.jl b/src/univariate.jl index 9c72712..8dcf30a 100644 --- a/src/univariate.jl +++ b/src/univariate.jl @@ -155,6 +155,14 @@ struct UnivariateBasis{F,K,D} <: FunctionBasis end end +function Base.show(io::IO, basis::UnivariateBasis) + (; family, kind, domain_transformation, level) = basis + lead = "UnivariateBasis(" + next = ",\n" * ' '^length(lead) + print(io, "UnivariateBasis(", family, next, kind, next, domain_transformation, next, + level, ") # dimension: ", dimension(basis)) +end + domain(U::UnivariateBasis) = domain(U.domain_transformation) dimension(U::UnivariateBasis) = grid_length(U.family, U.kind, U.level) From 440b151a0f320b2be497070dbb278a36dd001528 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tam=C3=A1s=20K=2E=20Papp?= Date: Wed, 26 Aug 2026 16:30:42 +0200 Subject: [PATCH 20/34] remove coordinate transformations (unneeded) --- src/transformations.jl | 88 ------------------------------------ test/test_transformations.jl | 21 --------- 2 files changed, 109 deletions(-) diff --git a/src/transformations.jl b/src/transformations.jl index 530e8da..65f62c8 100644 --- a/src/transformations.jl +++ b/src/transformations.jl @@ -44,94 +44,6 @@ values and derivatives. """ function transform_from end -#### -#### coordinate transformations -#### - -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) -``` -""" -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 - -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 - -function transform_from(domain::CoordinateDomains{T}, ct::CoordinateTransformations, - x::AbstractVector) where {T} - SVector(transform_from(domain, ct, _ntuple_like(T, x))) -end - #### #### specific transformations #### diff --git a/test/test_transformations.jl b/test/test_transformations.jl index 1c326e1..e6145d4 100644 --- a/test/test_transformations.jl +++ b/test/test_transformations.jl @@ -90,24 +90,6 @@ end @test y_minf[1] == 0 end -@testset "coordinate transformations" begin - t1 = BoundedLinear(; lower = 2.0, upper = 3.0) - t2 = SemiInfRational(; endpoint = 7.0, scale = 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] t1 = BoundedLinear(; lower = 2.0, upper = 3) @@ -119,7 +101,4 @@ end t3 = InfRational(; center = 0.5, scale = 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) end From 8e9feabd59b0a6e148c067d13ca69f9221213c73 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tam=C3=A1s=20K=2E=20Papp?= Date: Wed, 26 Aug 2026 16:54:41 +0200 Subject: [PATCH 21/34] make boundedlinear use lower, upper; rework printing --- src/transformations.jl | 43 +++++++++++++++--------------------- test/test_transformations.jl | 6 ++--- 2 files changed, 21 insertions(+), 28 deletions(-) diff --git a/src/transformations.jl b/src/transformations.jl index 65f62c8..f196eff 100644 --- a/src/transformations.jl +++ b/src/transformations.jl @@ -53,10 +53,10 @@ function transform_from end ### struct BoundedLinear{T <: Real} <: AbstractUnivariateTransformation - "Midpoint `m`." - m::T - "Scale `s`." - s::T + "Lower limit." + lower::T + "Upper limit." + upper::T @doc """ $(SIGNATURES) @@ -67,41 +67,39 @@ struct BoundedLinear{T <: Real} <: AbstractUnivariateTransformation function BoundedLinear(; lower::Real, upper::Real) @argcheck isfinite(lower) && isfinite(upper) DomainError lower, upper = promote(lower, upper) - s = (upper - lower) / 2 - m = (lower + upper) / 2 - @argcheck s > 0 DomainError((; lower, upper), "Need `lower < upper`.") - m, s = promote(m, s) - new{typeof(m)}(m, s) + @argcheck upper > lower DomainError((; lower, upper), "Need `lower < upper`.") + new{typeof(lower)}(lower, upper) end end function Base.show(io::IO, transformation::BoundedLinear) - (; m, s) = transformation - print(io, "(", m - s, ",", m + s, ") ↔ domain [linear transformation]") + (; lower, upper) = transformation + print(io, "BoundedLinear(lower = ", lower, ", upper = ", 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 + (; lower, upper) = 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 ### @@ -137,12 +135,7 @@ end function Base.show(io::IO, transformation::SemiInfRational) (; endpoint, scale) = transformation - if scale > 0 - D = "($(endpoint),∞)" - else - D = "(-∞,$(endpoint))" - end - print(io, D, " ↔ domain [rational transformation with scale ", scale, "]") + print(io, "SemiInfRational(endpoint = ", endpoint, ", scale = ", scale, ")") end transform_from(::PM1, t::SemiInfRational, x) = t.endpoint + t.scale * (1 + x) / (1 - x) @@ -206,7 +199,7 @@ end function Base.show(io::IO, transformation::InfRational) (; center, scale) = transformation - print(io, "(-∞,∞) ↔ domain [rational transformation with center ", center, ", scale ", scale, "]") + print(io, "InfRational(; center = ", center, ", scale = ", scale, ")") end InfRational(; center::Real = 0.0, scale::Real = 1.0) = InfRational(promote(center, scale)...) diff --git a/test/test_transformations.jl b/test/test_transformations.jl index e6145d4..d796d01 100644 --- a/test/test_transformations.jl +++ b/test/test_transformations.jl @@ -93,12 +93,12 @@ end @testset "printing, promotion, broadcasting" begin v = [1.0, 2.0] t1 = BoundedLinear(; lower = 2.0, upper = 3) - @test repr(t1) == "(2.0,3.0) ↔ domain [linear transformation]" + @test repr(t1) == "BoundedLinear(lower = 2.0, upper = 3.0)" @test transform_to.(PM1(), t1, v) isa Vector t2 = SemiInfRational(; endpoint = 7.0) - @test repr(t2) == "(7.0,∞) ↔ domain [rational transformation with scale 1.0]" + @test repr(t2) == "SemiInfRational(endpoint = 7.0, scale = 1.0)" @test transform_to.(PM1(), t2, v) isa Vector t3 = InfRational(; center = 0.5, scale = 1) - @test repr(t3) == "(-∞,∞) ↔ domain [rational transformation with center 0.5, scale 1.0]" + @test repr(t3) == "InfRational(; center = 0.5, scale = 1.0)" @test transform_to.(PM1(), t3, v) isa Vector end From 363a9a0bdb819b914f7e43d2474635703dc08c76 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tam=C3=A1s=20K=2E=20Papp?= Date: Wed, 26 Aug 2026 16:57:49 +0200 Subject: [PATCH 22/34] remove support for broadcasting --- src/transformations.jl | 2 -- test/test_transformations.jl | 6 +----- 2 files changed, 1 insertion(+), 7 deletions(-) diff --git a/src/transformations.jl b/src/transformations.jl index f196eff..7d83b5d 100644 --- a/src/transformations.jl +++ b/src/transformations.jl @@ -16,8 +16,6 @@ An abstract type for univariate transformations. """ abstract type AbstractUnivariateTransformation end -Broadcast.broadcastable(transformation::AbstractUnivariateTransformation) = Ref(transformation) - domain_kind(::Type{<:AbstractUnivariateTransformation}) = :univariate """ diff --git a/test/test_transformations.jl b/test/test_transformations.jl index d796d01..e495091 100644 --- a/test/test_transformations.jl +++ b/test/test_transformations.jl @@ -90,15 +90,11 @@ end @test y_minf[1] == 0 end -@testset "printing, promotion, broadcasting" begin - v = [1.0, 2.0] +@testset "printing, promotion" begin t1 = BoundedLinear(; lower = 2.0, upper = 3) @test repr(t1) == "BoundedLinear(lower = 2.0, upper = 3.0)" - @test transform_to.(PM1(), t1, v) isa Vector t2 = SemiInfRational(; endpoint = 7.0) @test repr(t2) == "SemiInfRational(endpoint = 7.0, scale = 1.0)" - @test transform_to.(PM1(), t2, v) isa Vector t3 = InfRational(; center = 0.5, scale = 1) @test repr(t3) == "InfRational(; center = 0.5, scale = 1.0)" - @test transform_to.(PM1(), t3, v) isa Vector end From 5f3d9fd2ebf63f11eac48be672bb7aac2eb49dce Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tam=C3=A1s=20K=2E=20Papp?= Date: Wed, 26 Aug 2026 17:28:19 +0200 Subject: [PATCH 23/34] implement index traversal --- src/smolyak_api.jl | 98 ++++------------- src/smolyak_traversal.jl | 185 +-------------------------------- test/test_smolyak.jl | 16 +++ test/test_smolyak_traversal.jl | 5 +- 4 files changed, 37 insertions(+), 267 deletions(-) diff --git a/src/smolyak_api.jl b/src/smolyak_api.jl index fcb0c75..227bb04 100644 --- a/src/smolyak_api.jl +++ b/src/smolyak_api.jl @@ -21,8 +21,8 @@ struct SmolyakLevel If `each > total`, the normalization `each = total` is set with a warning. """ function SmolyakLevel(; total::Int, each::Int = total) - @argcheck total ≥ 1 - @argcheck each ≥ 1 + @argcheck total ≥ 0 + @argcheck each ≥ 0 if each > total @warn "‘each’ normalized to ‘total’" each total each = total @@ -208,21 +208,32 @@ function grid(::Type{T}, smolyak_basis::SmolyakBasis) where {T<:AbstractFloat} SmolyakGrid(family, kind, level, itrs) end -@concrete struct SmolyakIndices +#### +#### index traversal +#### + +@concrete struct SmolyakIndices{I} family kind level - itrs + 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 +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 smolyak_indices(basis::SmolyakBasis) +function SmolyakIndices(basis::SmolyakBasis) (; family, kind, domain_transformations, level) = basis - itr1 = 1:level.each + itr1 = Iterators.countfrom(1, 1) itrs = ntuple(_ -> itr1, Val(length(domain_transformations))) SmolyakIndices(family, kind, level, itrs) end @@ -248,76 +259,3 @@ function adjust_coefficients(θ1::AbstractVector{T}, basis1::SmolyakBasis, z = zero(T) [get(θ, ι, z) for ι in smolyak_indices(basis2)] 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 f78a1aa..b5952c6 100644 --- a/src/smolyak_traversal.jl +++ b/src/smolyak_traversal.jl @@ -134,7 +134,7 @@ end """ $(SIGNATURES) -Calculate the length of a [`SmolyakIndices`](@ref) iterator. Argument as in the latter. +Calculate the dimension of Smolyak basis. """ function __smolyak_length(family, kind, ::Val{N}, total::Int, each::Int) where N L = 0 @@ -158,186 +158,3 @@ function __smolyak_length(family, kind, ::Val{N}, total::Int, each::Int) where N end L end - - -#### -#### index traversal -#### - -# function __inc_init(family, kind, ::Val{N}, total) where N -# indices = ntuple(_ -> 1, Val(N)) -# levels = ntuple(_ -> 0, Val(N)) -# b0 = block_length(family, kind, 0) -# limits = ntuple(_ -> b0, Val(N)) -# slack = total -# slack, indices, levels, limits -# end - -# """ -# $(SIGNATURES) → valid, Δ, indices, levels, limits - -# Internal implementation of the Smolyak indexing iterator. - -# # Arguments - -# - `slack`: `total - sum(blocks)`, cached - -# - `indices`: current indices - -# - `levels`: level along each dimension - -# - `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′`, `levels′, `limits′`: next values for corresponding arguments above, each an -# `::NTuple{N,Int}` -# """ -# @inline function __inc(family, kind, each::Int, # unchanged -# slack::Int, indices::NTuple{N,Int}, levels::NTuple{N,Int}, -# limits::NTuple{N,Int}) where N -# i1, iτ... = indices -# l1, lτ... = levels -# L1, Lτ... = limits -# if i1 < L1 # increment i1, same block -# true, 0, (i1 + 1, iτ...), levels, limits -# elseif l1 < each && slack > 0 # increment i1, next block -# l1′ = l1 + 1 -# true, -1, (i1 + 1, iτ...), (l1′, lτ...), (grid_length(family, kind, l1′), Lτ...) -# else -# if N == 1 # end of iteration, arbitrary value since !valid -# false, 0, indices, levels, limits -# else # i1 = 1, increment tail if applicable -# Δ1 = l1 -# valid, Δτ, iτ′, lτ′, Lτ′ = __inc(family, kind, each, slack + Δ1, iτ, lτ, Lτ) -# valid, Δ1 + Δτ, (1, iτ′...), (0, lτ′...), (grid_length(family, kind, 0), Lτ′...) -# end -# end -# 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, `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 `l1` denote the -# smallest integer such that `i1 ≤ ℓ(l1)`, 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 ≤ l1 ≤ M`, …, `1 ≤ bN ≤ M`, -# 3. `l1 + … + 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}(kind, smolyak_parameters::SmolyakParameters{B,M}) where {N,B,M} -# @argcheck N ≥ 1 -# Mp1 = M + 1 -# len = __smolyak_length(kind, Val(N), Val(B), M) -# first_block_length = nesting_total_length(Chebyshev, kind, 0) -# nesting_total_lengths = ntuple(bp1 -> nesting_total_length(Chebyshev, kind, bp1 - 1), -# Val(Mp1)) -# H = last(nesting_total_lengths) -# new{N,H,B,M,Mp1}(len, nesting_total_lengths) -# end -# end - -# function Base.show(io::IO, smolyak_indices::SmolyakIndices{N,H,B,M}) where {N,H,B,M} -# (; len) = smolyak_indices -# print(io, "Smolyak indexing, ∑bᵢ ≤ $(B), all bᵢ ≤ $(M), dimension $(len)") -# end - -# @inline highest_visited_index(::SmolyakIndices{N,H}) where {N,H} = H - -# Base.eltype(::Type{<:SmolyakIndices{N}}) where N = NTuple{N,Int} - -# @inline Base.length(ι::SmolyakIndices) = ι.len - -# @inline function Base.iterate(ι::SmolyakIndices{N,H,B}) where {N,H,B} -# slack, indices, blocks, limits = __inc_init(ι.nesting_total_lengths, Val(N), Val(B)) -# indices, (slack, indices, blocks, limits) -# end - -# @inline function Base.iterate(ι::SmolyakIndices, (slack, indices, blocks, limits)) -# valid, Δ, indices′, blocks′, limits′ = __inc(ι.nesting_total_lengths, slack, indices, -# blocks, limits) -# valid || return nothing -# slack′ = slack + Δ -# indices′, (slack′, indices′, blocks′, limits′) -# end - -# #### -# #### product traversal -# #### - -# struct SmolyakProduct{I<:SmolyakIndices,S<:Tuple,P} -# smolyak_indices::I -# sources::S -# product_kind::P -# @doc """ -# $(SIGNATURES) - -# An iterator conceptually equivalent to - -# ``` -# [prod(getindex.(sources, indices)) for indices in smolyak_indices] -# ``` - -# using [`_product`](@ref) instead to account for derivatives. Detailed docs of the -# arguments are in [`SmolyakIndices`](@ref). - -# Caller should arrange the elements of `sources` in the correct order, see -# [`nested_extrema_indices`](@ref). Each element in `sources` should have at least -# `H` elements (cf type parameters of [`SmolyakIndices`](@ref)), this is not checked. -# """ -# function SmolyakProduct(smolyak_indices::I, sources::S, -# product_kind::P) where {N,I<:SmolyakIndices{N},S,P} -# @argcheck length(sources) == N -# new{I,S,P}(smolyak_indices, sources, product_kind) -# end -# end - -# Base.length(smolyak_product::SmolyakProduct) = length(smolyak_product.smolyak_indices) - -# function Base.eltype(::Type{SmolyakProduct{I,S,P}}) where {I,S,P} -# _product_type(P, fieldtypes(S)) -# end - -# @inline function Base.iterate(smolyak_product::SmolyakProduct, state...) -# (; smolyak_indices, sources, product_kind) = smolyak_product -# itr = iterate(smolyak_indices, state...) -# itr ≡ nothing && return nothing -# indices, state′ = itr -# _product(product_kind, map(getindex, sources, indices)), state′ -# end diff --git a/test/test_smolyak.jl b/test/test_smolyak.jl index a071b21..ff170cd 100644 --- a/test/test_smolyak.jl +++ b/test/test_smolyak.jl @@ -55,6 +55,22 @@ end @test @ballocated(linear_combination($basis, $θ, $y)) == 0 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 + + ### ### augment coefficients ### diff --git a/test/test_smolyak_traversal.jl b/test/test_smolyak_traversal.jl index c977ccd..b262290 100644 --- a/test/test_smolyak_traversal.jl +++ b/test/test_smolyak_traversal.jl @@ -29,9 +29,8 @@ end @testset "__smolyak_length" begin for kind in KINDS for N in 1:5 - for total in 1:4 - for each in 1:total - # kind, N, total, each = Endpoints(), 4, 2, 1 + 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 From b65f7749cf3952ac99002db60b024647c653cb74 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tam=C3=A1s=20K=2E=20Papp?= Date: Wed, 26 Aug 2026 19:14:32 +0200 Subject: [PATCH 24/34] BoundedLinear API, test adjust_coefficients, traversal incidental: printing changes --- src/smolyak_api.jl | 9 +++++---- src/transformations.jl | 9 ++++++--- src/univariate.jl | 5 +++-- test/test_smolyak.jl | 19 ++++++++++++------- test/test_transformations.jl | 10 +++++----- test/test_univariate.jl | 9 ++++----- 6 files changed, 35 insertions(+), 26 deletions(-) diff --git a/src/smolyak_api.jl b/src/smolyak_api.jl index 227bb04..2394980 100644 --- a/src/smolyak_api.jl +++ b/src/smolyak_api.jl @@ -59,7 +59,8 @@ end function Base.show(io::IO, level::SmolyakLevel) (; total, each) = level - print(io, "Smolyak parameters, ∑ℓᵢ ≤ $(total), all ℓᵢ ≤ $(each)") + explanation = "∑ℓᵢ ≤ $(total), all ℓᵢ ≤ $(each)" + print(io, "SmolyakLevel(total = $(total), each = $(each)) #= $(explanation) =#") end struct SmolyakBasis{F,K,D} <: MultivariateBasis @@ -81,7 +82,7 @@ function Base.show(io::IO, basis::SmolyakBasis) lead = "SmolyakBasis(" next = ",\n" * ' '^length(lead) print(io, "SmolyakBasis(", family, next, kind, next, domain_transformations, next, - level, ") # dimension: ", dimension(basis)) + level, ") # dimension: ", dimension(basis)) end @@ -253,9 +254,9 @@ function adjust_coefficients(θ1::AbstractVector{T}, basis1::SmolyakBasis, return copy(θ1) end θ = Dict{NTuple{length(domain_transformations),Int},T}() - for (x, ι) in zip(θ1, smolyak_indices(basis1)) + for (x, ι) in zip(θ1, SmolyakIndices(basis1)) θ[ι] = x end z = zero(T) - [get(θ, ι, z) for ι in smolyak_indices(basis2)] + [get(θ, ι, z) for ι in SmolyakIndices(basis2)] end diff --git a/src/transformations.jl b/src/transformations.jl index 7d83b5d..50f2763 100644 --- a/src/transformations.jl +++ b/src/transformations.jl @@ -62,17 +62,20 @@ struct BoundedLinear{T <: Real} <: AbstractUnivariateTransformation `lower < upper` is enforced. """ - function BoundedLinear(; lower::Real, upper::Real) + function BoundedLinear(lower::T, upper::T) where {T<:Real} @argcheck isfinite(lower) && isfinite(upper) DomainError - lower, upper = promote(lower, upper) @argcheck upper > lower DomainError((; lower, upper), "Need `lower < upper`.") new{typeof(lower)}(lower, upper) end end +BoundedLinear(lower::Real, upper::Real) = BoundedLinear(promote(lower, upper)...) + +BoundedLinear(; lower, upper) = BoundedLinear(lower, upper) + function Base.show(io::IO, transformation::BoundedLinear) (; lower, upper) = transformation - print(io, "BoundedLinear(lower = ", lower, ", upper = ", upper, ")") + print(io, "BoundedLinear(", lower, ", ", upper, ")") end function transform_from(::PM1, t::BoundedLinear, x::Scalar) diff --git a/src/univariate.jl b/src/univariate.jl index 8dcf30a..1dbd36b 100644 --- a/src/univariate.jl +++ b/src/univariate.jl @@ -159,8 +159,9 @@ function Base.show(io::IO, basis::UnivariateBasis) (; family, kind, domain_transformation, level) = basis lead = "UnivariateBasis(" next = ",\n" * ' '^length(lead) - print(io, "UnivariateBasis(", family, next, kind, next, domain_transformation, next, - level, ") # dimension: ", dimension(basis)) + print(io, "UnivariateBasis(", family, + next, kind, next, domain_transformation, next, + level, ") # dimension: ", dimension(basis)) end domain(U::UnivariateBasis) = domain(U.domain_transformation) diff --git a/test/test_smolyak.jl b/test/test_smolyak.jl index ff170cd..47f9e2b 100644 --- a/test/test_smolyak.jl +++ b/test/test_smolyak.jl @@ -3,7 +3,8 @@ #### @testset "printing SmolyakLevels" begin - @test repr(SmolyakLevel(; total = 3, each = 2)) == "Smolyak parameters, ∑ℓᵢ ≤ 3, all ℓᵢ ≤ 2" + @test repr(SmolyakLevel(; total = 3, each = 2)) == + "SmolyakLevel(total = 3, each = 2) #= ∑ℓᵢ ≤ 3, all ℓᵢ ≤ 2 =#" end @testset "Smolyak API checks" begin @@ -12,8 +13,8 @@ end @testset "Smolyak API sanity checks" begin f(x) = (x[1] - 3) * (x[2] + 5) # linear function, just a sanity check - transformations = (BoundedLinear(lower = 2, upper = 3), # approximation should be exact - BoundedLinear(lower = 3.0, upper = 4.5)) + 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) @@ -70,10 +71,14 @@ end end end - -### -### augment coefficients -### +@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) diff --git a/test/test_transformations.jl b/test/test_transformations.jl index e495091..0ee70e2 100644 --- a/test/test_transformations.jl +++ b/test/test_transformations.jl @@ -1,11 +1,11 @@ using SpectralKit: PM1 @testset "bounded linear domain transformations" begin - @test_throws DomainError BoundedLinear(; lower = -1.0, upper = Inf) - @test_throws DomainError BoundedLinear(; lower = -1.0, upper = -2.0) + @test_throws DomainError BoundedLinear(-1.0, Inf) + @test_throws DomainError BoundedLinear(-1.0, -2.0) A, B = 1, 5 - trans = BoundedLinear(; lower = A, upper = B) + trans = BoundedLinear(A, B) @test extrema(domain(trans)) == (A, B) @@ -91,8 +91,8 @@ end end @testset "printing, promotion" begin - t1 = BoundedLinear(; lower = 2.0, upper = 3) - @test repr(t1) == "BoundedLinear(lower = 2.0, upper = 3.0)" + t1 = BoundedLinear(2.0, 3) + @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) diff --git a/test/test_univariate.jl b/test/test_univariate.jl index e0d4d6e..cbf4d34 100644 --- a/test/test_univariate.jl +++ b/test/test_univariate.jl @@ -48,8 +48,7 @@ end for kind in KINDS previous = Float64[] for level in 0:7 - b = UnivariateBasis(Chebyshev(), kind, - BoundedLinear(; lower = 1.0, upper = 3.0), level) + b = UnivariateBasis(Chebyshev(), kind, BoundedLinear(1.0, 3.0), level) g = collect(grid(b)) @test is_approximate_subset(previous, g) g = previous @@ -58,7 +57,7 @@ end end @testset "Chebyshev basics" begin - transformation = BoundedLinear(; lower = 1.0, upper = 3.0) + transformation = BoundedLinear(1.0, 3.0) @test_throws ArgumentError UnivariateBasis(Chebyshev(), Endpoints(), transformation, -1) for kind in KINDS for level in 0:5 @@ -118,7 +117,7 @@ end end @testset "univariate derivatives" begin - for (transformation, N) in ((BoundedLinear(lower = -2, upper = 3), 5), + 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) @@ -128,7 +127,7 @@ end x = rand_in_domain(basis) y = f(D(x)) for i in 0:N - @test y[i] ≈ DD(f, x, i) rtol = 1e-4 atol = 1e-4 + @test y[i] ≈ DD(f, x, i) rtol = 1e-3 atol = 1e-3 end end end From deadd78f35b9cf64887d769c91ed2300263305b3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tam=C3=A1s=20K=2E=20Papp?= Date: Sun, 30 Aug 2026 07:46:30 +0200 Subject: [PATCH 25/34] add Enzyme integration tests --- test/Project.toml | 2 +- test/integration_enzyme.jl | 16 ++++++++++++++++ test/runtests.jl | 1 + test/test_univariate.jl | 8 ++++++++ 4 files changed, 26 insertions(+), 1 deletion(-) create mode 100644 test/integration_enzyme.jl diff --git a/test/Project.toml b/test/Project.toml index 956fdc9..f71a66e 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -3,6 +3,7 @@ Accessors = "7d9f7c33-5ae7-4f3b-8dc6-eff91059b697" Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595" BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf" 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" @@ -15,7 +16,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 2a169e8..80a735f 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -25,5 +25,6 @@ include("test_transformations.jl") include("test_univariate.jl") include("test_smolyak_traversal.jl") include("test_smolyak.jl") +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_univariate.jl b/test/test_univariate.jl index cbf4d34..3454ca9 100644 --- a/test/test_univariate.jl +++ b/test/test_univariate.jl @@ -165,3 +165,11 @@ end 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 From 2d3b2d28390cbd4070c865616e66bb10faa35748 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tam=C3=A1s=20K=2E=20Papp?= Date: Sun, 30 Aug 2026 08:33:38 +0200 Subject: [PATCH 26/34] rework iteration to fix Enzyme inference --- src/smolyak_api.jl | 28 +++++++++++----------- src/smolyak_traversal.jl | 33 +++++++++++++++----------- src/univariate.jl | 50 +++++++++++++++++++++++++++++++++------- src/utilities.jl | 40 ++++++++++++++++++++++++++++++++ test/utilities.jl | 18 ++++++++++++--- 5 files changed, 130 insertions(+), 39 deletions(-) diff --git a/src/smolyak_api.jl b/src/smolyak_api.jl index 2394980..f3e8ca4 100644 --- a/src/smolyak_api.jl +++ b/src/smolyak_api.jl @@ -39,19 +39,17 @@ A wrapper for iterating through Smolyak indices. See [`__smolyak__init`](@ref) a `g` transforms the value. """ -function __smolyak_iterate(family, kind, level::SmolyakLevel, f, itrs, g, - state = nothing) +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 - result = __smolyak_step(family, kind, level.each, f, itrs, - state...) - if result ≡ nothing + accum, Δ, rest... = __smolyak_step(family, kind, level.each, f, itrs, state...) + if Δ == Δ_DONE nothing else slack = state[2] - accum, Δ, rest... = result g(accum), (accum, slack + Δ, rest...) end end @@ -145,7 +143,7 @@ struct SmolyakBasisAt{I,P,F,K,L<:SmolyakLevel} end function Base.eltype(::Type{<:SmolyakBasisAt{I,P}}) where {I,P} - _product_type(P, map(eltype, fieldtypes(I))) + _product_type(P, map(_eltype, fieldtypes(I))) end function Base.length(itr::SmolyakBasisAt) @@ -154,19 +152,20 @@ function Base.length(itr::SmolyakBasisAt) __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, - (a, b) -> isempty(b) ? (a,) : (a * first(b), b...), - itrs, first, state) + __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) -> basis_at(UnivariateBasis(family, kind, d, level.each), x), - x, domain_transformations) + itrs = map((x, d) -> _univariate_basis_itr(family, d, x), x, domain_transformations) SmolyakBasisAt(family, kind, level, itrs, nothing) end @@ -190,7 +189,7 @@ struct SmolyakGrid{I,F,K,S<:SmolyakLevel} end function Base.eltype(::Type{<:SmolyakGrid{I}}) where {I} - Tuple{map(eltype, fieldtypes(I))...} + Tuple{map(_eltype, fieldtypes(I))...} end function Base.length(itr::SmolyakGrid) @@ -234,8 +233,7 @@ Note: only used as a building block for [`adjust_coefficients`](@ref). """ function SmolyakIndices(basis::SmolyakBasis) (; family, kind, domain_transformations, level) = basis - itr1 = Iterators.countfrom(1, 1) - itrs = ntuple(_ -> itr1, Val(length(domain_transformations))) + itrs = ntuple(_ -> Counting(), Val(length(domain_transformations))) SmolyakIndices(family, kind, level, itrs) end diff --git a/src/smolyak_traversal.jl b/src/smolyak_traversal.jl index b5952c6..0fe95bd 100644 --- a/src/smolyak_traversal.jl +++ b/src/smolyak_traversal.jl @@ -11,13 +11,18 @@ function __smolyak_init(family, kind, total::Int, f::F, itrs::NTuple{N,Any}) whe slack = total r = block_length(family, kind, 0) remainders = ntuple(_ -> r - 1, Val(N)) - itr_results = map(iterate, itrs) + 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 +""" +Sentinel value for having finished iteration in `__smolyak_step`. +""" +const Δ_DONE = typemax(Int) + """ $(SIGNATURES) → accum, Δ, remainders′, states′, levels′ @@ -45,31 +50,33 @@ function __smolyak_step(family, kind, each::Int, f::F, itrs::NTuple{N,Any}, s1, sτ... = states l1, lτ... = levels if r1 > 0 # step within block - x1, s1′ = iterate(I1, s1) + 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′ = iterate(I1, s1) + 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 - nothing # done with iteration + accum, Δ_DONE, remainders, states, levels # done with iteration else # go into tail - next = __smolyak_step(family, kind, each, f, Iτ, aτ, slack + l1, rτ, sτ, lτ) - next ≡ nothing && return nothing - aτ′, Δτ, rτ′, sτ′, lτ′ = next - x1, s1 = iterate(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 + 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 + 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 diff --git a/src/univariate.jl b/src/univariate.jl index 1dbd36b..7ae435e 100644 --- a/src/univariate.jl +++ b/src/univariate.jl @@ -23,21 +23,19 @@ struct ChebyshevIterator{T} x::T end -Base.eltype(::Type{<:ChebyshevIterator{T}}) where {T} = T - -Base.IteratorSize(::Type{<:ChebyshevIterator}) = Base.IsInfinite() - -function Base.iterate(itr::ChebyshevIterator{T}) where T +function _start(itr::ChebyshevIterator{T}) where T (; x) = itr _one(T), (_one(T), x) end -function Base.iterate(itr::ChebyshevIterator{T}, (fp, fpp)) where T +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) @@ -168,8 +166,44 @@ domain(U::UnivariateBasis) = domain(U.domain_transformation) dimension(U::UnivariateBasis) = grid_length(U.family, U.kind, U.level) -function basis_at(U::UnivariateBasis{Chebyshev}, x::Scalar) - Iterators.take(ChebyshevIterator(transform_to(PM1(), U.domain_transformation, x)), dimension(U)) +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} diff --git a/src/utilities.jl b/src/utilities.jl index 5b08d56..86a4e74 100644 --- a/src/utilities.jl +++ b/src/utilities.jl @@ -2,6 +2,46 @@ ##### 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 #### diff --git a/test/utilities.jl b/test/utilities.jl index 20cda26..b0d2fd6 100644 --- a/test/utilities.jl +++ b/test/utilities.jl @@ -153,6 +153,19 @@ 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. """ @@ -163,13 +176,12 @@ function test_smolyak_step(family, kind, total, each, f, itrs::NTuple{N,Any}) wh total, f, itrs) while true @test slack + sum(levels) == total # simple sanity check - cs = map((itr, i) -> first(Iterators.drop(itr, i - 1)), - itrs, reference[i]) + 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) - next ≡ nothing && break (accum, Δ, remainders, states, levels) = next + Δ == SpectralKit.Δ_DONE && break slack += Δ i += 1 end From 2457e5bba021980a89b389df264b2a455cc8c506 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tam=C3=A1s=20K=2E=20Papp?= Date: Mon, 31 Aug 2026 15:06:00 +0200 Subject: [PATCH 27/34] add ConstructionBase method for SmolyakLevel --- Project.toml | 2 ++ src/smolyak_api.jl | 22 +++++++++++++++++----- test/Project.toml | 1 + test/test_smolyak.jl | 17 ++++++++++++++--- 4 files changed, 34 insertions(+), 8 deletions(-) diff --git a/Project.toml b/Project.toml index 6c0e3e5..12d2a7d 100644 --- a/Project.toml +++ b/Project.toml @@ -10,6 +10,7 @@ projects = ["test", "docs"] 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" @@ -19,6 +20,7 @@ StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" 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/smolyak_api.jl b/src/smolyak_api.jl index f3e8ca4..37b1429 100644 --- a/src/smolyak_api.jl +++ b/src/smolyak_api.jl @@ -4,6 +4,8 @@ export SmolyakLevel, SmolyakBasis +import ConstructionBase + struct SmolyakLevel total::Int each::Int @@ -31,6 +33,21 @@ struct SmolyakLevel end end +function Base.show(io::IO, level::SmolyakLevel) + (; total, each) = level + explanation = "0 ≤ ∑ℓᵢ ≤ $(total), all 0 ≤ ℓᵢ ≤ $(each)" + print(io, "SmolyakLevel(total = $(total), each = $(each)) #= $(explanation) =#") +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) @@ -55,11 +72,6 @@ function __smolyak_iterate(family, kind, level::SmolyakLevel, f::F, itrs, g::G, end end -function Base.show(io::IO, level::SmolyakLevel) - (; total, each) = level - explanation = "∑ℓᵢ ≤ $(total), all ℓᵢ ≤ $(each)" - print(io, "SmolyakLevel(total = $(total), each = $(each)) #= $(explanation) =#") -end struct SmolyakBasis{F,K,D} <: MultivariateBasis family::F diff --git a/test/Project.toml b/test/Project.toml index f71a66e..fc7cc57 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -2,6 +2,7 @@ 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" diff --git a/test/test_smolyak.jl b/test/test_smolyak.jl index 47f9e2b..f9b3d36 100644 --- a/test/test_smolyak.jl +++ b/test/test_smolyak.jl @@ -2,15 +2,26 @@ #### api #### -@testset "printing SmolyakLevels" begin +using ConstructionBase + +@testset "SmolyakLevels printing" begin @test repr(SmolyakLevel(; total = 3, each = 2)) == - "SmolyakLevel(total = 3, each = 2) #= ∑ℓᵢ ≤ 3, all ℓᵢ ≤ 2 =#" + "SmolyakLevel(total = 3, each = 2) #= 0 ≤ ∑ℓᵢ ≤ 3, all 0 ≤ ℓᵢ ≤ 2 =#" end -@testset "Smolyak API checks" begin +@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 API sanity checks" begin f(x) = (x[1] - 3) * (x[2] + 5) # linear function, just a sanity check transformations = (BoundedLinear(2, 3), # approximation should be exact From 0cb15ce886de1e43def6eeafe3c647db24378434 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tam=C3=A1s=20K=2E=20Papp?= Date: Mon, 31 Aug 2026 15:08:42 +0200 Subject: [PATCH 28/34] make file names consistent --- src/SpectralKit.jl | 2 +- src/{smolyak_api.jl => smolyak.jl} | 0 2 files changed, 1 insertion(+), 1 deletion(-) rename src/{smolyak_api.jl => smolyak.jl} (100%) diff --git a/src/SpectralKit.jl b/src/SpectralKit.jl index 074d0ee..ea2198d 100644 --- a/src/SpectralKit.jl +++ b/src/SpectralKit.jl @@ -13,7 +13,7 @@ include("transformations.jl") include("generic_api.jl") include("univariate.jl") include("smolyak_traversal.jl") -include("smolyak_api.jl") +include("smolyak.jl") # include("experimental.jl") # experimental code is not part of the API, see its module docstring end # module diff --git a/src/smolyak_api.jl b/src/smolyak.jl similarity index 100% rename from src/smolyak_api.jl rename to src/smolyak.jl From 4bd1062cbb0223f1f7b49407d5e7c264d8a82007 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tam=C3=A1s=20K=2E=20Papp?= Date: Mon, 31 Aug 2026 19:46:36 +0200 Subject: [PATCH 29/34] decouple grid_level in univariate --- src/univariate.jl | 19 ++++++++++++------- test/test_univariate.jl | 7 +++++++ 2 files changed, 19 insertions(+), 7 deletions(-) diff --git a/src/univariate.jl b/src/univariate.jl index 7ae435e..b86f2cf 100644 --- a/src/univariate.jl +++ b/src/univariate.jl @@ -138,6 +138,7 @@ struct UnivariateBasis{F,K,D} <: FunctionBasis kind::K domain_transformation::D level::Int + grid_level::Int @doc """ $(SIGNATURES) @@ -147,19 +148,23 @@ struct UnivariateBasis{F,K,D} <: FunctionBasis build the grid. """ function UnivariateBasis(family::F, kind::K, domain_transformation::D, - level::Int) where {F,K,D} + level::Int, grid_level::Int = level) where {F,K,D} @argcheck level ≥ 0 - new{F,K,D}(family, kind, domain_transformation, level) + 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) = basis + (; family, kind, domain_transformation, level, grid_level) = basis lead = "UnivariateBasis(" next = ",\n" * ' '^length(lead) print(io, "UnivariateBasis(", family, next, kind, next, domain_transformation, next, - level, ") # dimension: ", dimension(basis)) + level) + if grid_level ≠ level + print(io, next, grid_level) + end + print(io, ") # dimension: ", dimension(basis)) end domain(U::UnivariateBasis) = domain(U.domain_transformation) @@ -225,9 +230,9 @@ end function grid(::Type{T}, U::UnivariateBasis{Chebyshev, <:Union{Endpoints,Interior}}) where {T <: AbstractFloat} - (; family, kind, domain_transformation, level) = U - N = grid_length(family, kind, level) - ChebyshevGrid{T}(kind, domain_transformation, level, N, + (; 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 diff --git a/test/test_univariate.jl b/test/test_univariate.jl index 3454ca9..c327cd2 100644 --- a/test/test_univariate.jl +++ b/test/test_univariate.jl @@ -116,6 +116,13 @@ end end 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 "univariate derivatives" begin for (transformation, N) in ((BoundedLinear(-2, 3), 5), (SemiInfRational(endpoint = 0.7, scale = 0.3), 1), From a58ffd5b0332bc6a78f331904ae7d6e24752d41c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tam=C3=A1s=20K=2E=20Papp?= Date: Tue, 1 Sep 2026 09:10:21 +0200 Subject: [PATCH 30/34] decouple grid_level from level in Smolyak --- src/smolyak.jl | 15 ++++++++------- 1 file changed, 8 insertions(+), 7 deletions(-) diff --git a/src/smolyak.jl b/src/smolyak.jl index 37b1429..a446582 100644 --- a/src/smolyak.jl +++ b/src/smolyak.jl @@ -72,18 +72,19 @@ function __smolyak_iterate(family, kind, level::SmolyakLevel, f::F, itrs, g::G, end end - struct SmolyakBasis{F,K,D} <: MultivariateBasis family::F kind::K domain_transformations::D level::SmolyakLevel + grid_level::SmolyakLevel @doc """ $(SIGNATURES) """ function SmolyakBasis(family::F, kind::K, domain_transformations::D, - level::SmolyakLevel) where {F,K,D<:Tuple} - new{F,K,D}(family, kind, domain_transformations, level) + level::SmolyakLevel, + grid_level::SmolyakLevel = level) where {F,K,D<:Tuple} + new{F,K,D}(family, kind, domain_transformations, level, grid_level) end end @@ -214,14 +215,14 @@ function Base.iterate(itr::SmolyakGrid, state = nothing) end function grid(::Type{T}, smolyak_basis::SmolyakBasis) where {T<:AbstractFloat} - (; family, kind, domain_transformations, level) = smolyak_basis - (; each) = level + (; 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, level, itrs) + SmolyakGrid(family, kind, grid_level, itrs) end #### -#### index traversal +#### index traversal for adjustment #### @concrete struct SmolyakIndices{I} From 52883ae9506914022790069ea6bac5ac59c00dac Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tam=C3=A1s=20K=2E=20Papp?= Date: Tue, 1 Sep 2026 09:23:24 +0200 Subject: [PATCH 31/34] make printing a bit nicer --- src/smolyak.jl | 19 ++++++++++++++----- src/univariate.jl | 12 ++++-------- src/utilities.jl | 19 +++++++++++++++++++ test/test_smolyak.jl | 17 ++++++++++++++++- test/test_univariate.jl | 10 ++++++++++ 5 files changed, 63 insertions(+), 14 deletions(-) diff --git a/src/smolyak.jl b/src/smolyak.jl index a446582..d0933ca 100644 --- a/src/smolyak.jl +++ b/src/smolyak.jl @@ -35,8 +35,11 @@ end function Base.show(io::IO, level::SmolyakLevel) (; total, each) = level - explanation = "0 ≤ ∑ℓᵢ ≤ $(total), all 0 ≤ ℓᵢ ≤ $(each)" - print(io, "SmolyakLevel(total = $(total), each = $(each)) #= $(explanation) =#") + 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) @@ -89,11 +92,17 @@ struct SmolyakBasis{F,K,D} <: MultivariateBasis end function Base.show(io::IO, basis::SmolyakBasis) - (; family, kind, domain_transformations, level) = basis + (; family, kind, domain_transformations, level, grid_level) = basis + _print_dimensions(io, basis) lead = "SmolyakBasis(" next = ",\n" * ' '^length(lead) - print(io, "SmolyakBasis(", family, next, kind, next, domain_transformations, next, - level, ") # dimension: ", dimension(basis)) + print(io, lead, family, ", ", kind, next, + domain_transformations, next, + level) + if grid_level ≠ level + print(io, next, grid_level) + end + print(io, ")") end diff --git a/src/univariate.jl b/src/univariate.jl index b86f2cf..01a6a8b 100644 --- a/src/univariate.jl +++ b/src/univariate.jl @@ -156,15 +156,11 @@ end function Base.show(io::IO, basis::UnivariateBasis) (; family, kind, domain_transformation, level, grid_level) = basis - lead = "UnivariateBasis(" - next = ",\n" * ' '^length(lead) - print(io, "UnivariateBasis(", family, - next, kind, next, domain_transformation, next, + _print_dimensions(io, basis) + print(io, "UnivariateBasis(", family, ", ", kind, ", ", domain_transformation, ", ", level) - if grid_level ≠ level - print(io, next, grid_level) - end - print(io, ") # dimension: ", dimension(basis)) + level ≠ grid_level && print(io, ", ", grid_level) + print(io, ")") end domain(U::UnivariateBasis) = domain(U.domain_transformation) diff --git a/src/utilities.jl b/src/utilities.jl index 86a4e74..c1e82ec 100644 --- a/src/utilities.jl +++ b/src/utilities.jl @@ -46,6 +46,25 @@ _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 = ['₀', '₁', '₂', '₃', '₄', '₅', '₆', '₇', '₈', '₉'] diff --git a/test/test_smolyak.jl b/test/test_smolyak.jl index f9b3d36..04a6a26 100644 --- a/test/test_smolyak.jl +++ b/test/test_smolyak.jl @@ -6,7 +6,7 @@ using ConstructionBase @testset "SmolyakLevels printing" begin @test repr(SmolyakLevel(; total = 3, each = 2)) == - "SmolyakLevel(total = 3, each = 2) #= 0 ≤ ∑ℓᵢ ≤ 3, all 0 ≤ ℓᵢ ≤ 2 =#" + "SmolyakLevel(total = 3, each = 2) #= 0 ≤ ∑ℓᵢ ≤ 3, all 0 ≤ ℓᵢ ≤ 2 =#" end @testset "SmolyakLevel normalization" begin @@ -22,6 +22,21 @@ end @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 transformations = (BoundedLinear(2, 3), # approximation should be exact diff --git a/test/test_univariate.jl b/test/test_univariate.jl index c327cd2..3bed66e 100644 --- a/test/test_univariate.jl +++ b/test/test_univariate.jl @@ -116,6 +116,16 @@ 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)) == From 5e3de5832019daf4061e1f49ea6dc1b8fc82c9ab Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tam=C3=A1s=20K=2E=20Papp?= Date: Wed, 2 Sep 2026 09:46:03 +0200 Subject: [PATCH 32/34] add outside domain handling --- src/transformations.jl | 67 +++++++++++++++++++++++++++++++++-------- test/test_univariate.jl | 11 ++++++- test/utilities.jl | 24 +++++++++++++-- 3 files changed, 86 insertions(+), 16 deletions(-) diff --git a/src/transformations.jl b/src/transformations.jl index 50f2763..52d1863 100644 --- a/src/transformations.jl +++ b/src/transformations.jl @@ -2,7 +2,8 @@ ##### 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 #### @@ -42,6 +43,25 @@ values and derivatives. """ function transform_from end +#### +#### handling out of bounds values +#### + +""" +Error when the value is outside the domain. +""" +struct ErrorOutsideDomain end + +""" +Use the nearest value in the domain. +""" +struct NearestInDomain end + +""" +Default handling of values outside domain. Not part of the API. +""" +const OUTSIDE_DOMAIN = ErrorOutsideDomain() + #### #### specific transformations #### @@ -50,42 +70,65 @@ function transform_from end ### bounded linear ### -struct BoundedLinear{T <: Real} <: AbstractUnivariateTransformation +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) where {T<:Real} + 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)}(lower, upper) + new{typeof(lower),O}(lower, upper, outside_domain) end end -BoundedLinear(lower::Real, upper::Real) = BoundedLinear(promote(lower, upper)...) +function BoundedLinear(lower::Real, upper::Real, outside_domain = OUTSIDE_DOMAIN) + BoundedLinear(promote(lower, upper)..., outside_domain) +end -BoundedLinear(; lower, upper) = BoundedLinear(lower, upper) +function BoundedLinear(; lower, upper, outside_domain = OUTSIDE_DOMAIN()) + BoundedLinear(lower, upper, outside_domain) +end function Base.show(io::IO, transformation::BoundedLinear) - (; lower, upper) = transformation - print(io, "BoundedLinear(", lower, ", ", upper, ")") + (; lower, upper, outside_domain) = transformation + print(io, "BoundedLinear(", lower, ", ", upper) + outside_domain ≠ OUTSIDE_DOMAIN && print(io, ", ", outside_domain) + print(io, ")") end -function transform_from(::PM1, t::BoundedLinear, x::Scalar) +function _handle_domain(t::BoundedLinear{<:Real,ErrorOutsideDomain}, x::Real) (; lower, upper) = t - (x+1) / 2 * (upper-lower) + lower + @argcheck lower ≤ x ≤ upper DomainError(x, "outside domain") + x end -function transform_to(::PM1, t::BoundedLinear, y::Real) +function _handle_domain(t::BoundedLinear{<:Real,NearestInDomain}, x::Real) (; lower, upper) = t - (y-lower) / (upper-lower) * 2 - 1 + clamp(x, lower, upper) +end + +function transform_from(::PM1, t::BoundedLinear, x::Scalar) + (; lower, upper) = t + (x + 1) / 2 * (upper-lower) + lower +end + +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 diff --git a/test/test_univariate.jl b/test/test_univariate.jl index 3bed66e..f2f2626 100644 --- a/test/test_univariate.jl +++ b/test/test_univariate.jl @@ -133,18 +133,27 @@ end 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) rtol = 1e-3 atol = 1e-3 + @test y[i] ≈ DD(f, x, i; domain = dom) rtol = 1e-2 atol = 1e-2 end end end diff --git a/test/utilities.jl b/test/utilities.jl index b0d2fd6..d8c8714 100644 --- a/test/utilities.jl +++ b/test/utilities.jl @@ -115,12 +115,30 @@ 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 From 20ede742f2ccc4a6a5a3fe786bad4e6fbbc5022b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tam=C3=A1s=20K=2E=20Papp?= Date: Wed, 16 Sep 2026 11:04:40 +0200 Subject: [PATCH 33/34] add getindex for smolyak --- src/smolyak.jl | 30 ++++++++++++++++++------------ src/univariate.jl | 3 +++ test/test_smolyak.jl | 2 ++ 3 files changed, 23 insertions(+), 12 deletions(-) diff --git a/src/smolyak.jl b/src/smolyak.jl index d0933ca..9c58fcf 100644 --- a/src/smolyak.jl +++ b/src/smolyak.jl @@ -83,6 +83,18 @@ struct SmolyakBasis{F,K,D} <: MultivariateBasis 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, @@ -105,24 +117,18 @@ function Base.show(io::IO, basis::SmolyakBasis) 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 # """ # $(SIGNATURES) -# Create a sparse Smolyak basis. - -# # Arguments - -# - `family`: univariate function family, eg `Chebyshev`. - -# - `kind`: the grid kind, eg `Interior()` or `Endpoints()`. - -# - `domain_transformations` +Create a sparse Smolyak basis. -# - `smolyak_level`: the Smolyak level specificaion, see [`SmolyakLevel`](@ref). +# Arguments -# - `N`: the dimension. wrapped in a `Val` for type stability, a convenience constructor also -# takes integers. # ## Example diff --git a/src/univariate.jl b/src/univariate.jl index 01a6a8b..8d7dc0b 100644 --- a/src/univariate.jl +++ b/src/univariate.jl @@ -144,6 +144,9 @@ struct UnivariateBasis{F,K,D} <: FunctionBasis 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. """ diff --git a/test/test_smolyak.jl b/test/test_smolyak.jl index 04a6a26..17f06f6 100644 --- a/test/test_smolyak.jl +++ b/test/test_smolyak.jl @@ -58,6 +58,8 @@ end end 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 From 5820eb48c46062d25aedc7e30c52148c46d579e7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tam=C3=A1s=20K=2E=20Papp?= Date: Wed, 16 Sep 2026 11:15:23 +0200 Subject: [PATCH 34/34] fix typo --- src/smolyak.jl | 8 -------- 1 file changed, 8 deletions(-) diff --git a/src/smolyak.jl b/src/smolyak.jl index 9c58fcf..be4e219 100644 --- a/src/smolyak.jl +++ b/src/smolyak.jl @@ -122,14 +122,6 @@ function Base.getindex(basis::SmolyakBasis, i::Integer) UnivariateBasis(family, kind, domain_transformations[i], level.each, grid_level.each) end -# """ -# $(SIGNATURES) - -Create a sparse Smolyak basis. - -# Arguments - - # ## Example # FIXME these examples need to be updated