From 399ae4f6857c69d4ab4fe6ce9e3a763f32db7531 Mon Sep 17 00:00:00 2001 From: Sheehan Olver Date: Fri, 10 Oct 2025 15:50:43 +0100 Subject: [PATCH 1/3] Support vec with RectPolynomial and fix sum ordering --- Project.toml | 2 +- src/MultivariateOrthogonalPolynomials.jl | 2 +- src/rect.jl | 19 ++++++++++++++++--- test/test_rect.jl | 17 ++++++++++++----- 4 files changed, 30 insertions(+), 10 deletions(-) diff --git a/Project.toml b/Project.toml index c0c5e44..85766f2 100644 --- a/Project.toml +++ b/Project.toml @@ -39,7 +39,7 @@ InfiniteArrays = "0.15" InfiniteLinearAlgebra = "0.9, 0.10" LazyArrays = "2.3.1" LazyBandedMatrices = "0.11.3" -QuasiArrays = "0.12" +QuasiArrays = "0.12.2" RecurrenceRelationships = "0.2" SpecialFunctions = "1, 2" StaticArrays = "1" diff --git a/src/MultivariateOrthogonalPolynomials.jl b/src/MultivariateOrthogonalPolynomials.jl index 5d1f42c..a8684ed 100644 --- a/src/MultivariateOrthogonalPolynomials.jl +++ b/src/MultivariateOrthogonalPolynomials.jl @@ -10,7 +10,7 @@ import Base: axes, in, ==, +, -, /, *, ^, \, copy, copyto!, OneTo, getindex, siz import Base.Broadcast: Broadcasted, broadcasted, DefaultArrayStyle import DomainSets: boundary, EuclideanDomain -import QuasiArrays: LazyQuasiMatrix, LazyQuasiArrayStyle, domain +import QuasiArrays: LazyQuasiMatrix, LazyQuasiArrayStyle, domain, vec_layout import ContinuumArrays: @simplify, Weight, weight, grid, plotgrid, TransformFactorization, ExpansionLayout, plotvalues, unweighted, plan_transform, checkpoints, transform_ldiv, AbstractBasisLayout, basis_axes, Inclusion, grammatrix, weaklaplacian, layout_broadcasted, laplacian, abslaplacian, laplacian_axis, abslaplacian_axis, diff_layout, operatororder, broadcastbasis import ArrayLayouts: MemoryLayout, sublayout, sub_materialize diff --git a/src/rect.jl b/src/rect.jl index 9ee2454..8c03e1c 100644 --- a/src/rect.jl +++ b/src/rect.jl @@ -27,6 +27,8 @@ const RectPolynomial{T, PP} = KronPolynomial{2, T, PP} axes(P::KronPolynomial) = (Inclusion(Ɨ(map(domain, axes.(P.args, 1))...)), _krontrav_axes(axes.(P.args, 2)...)) + + function getindex(P::RectPolynomial{T}, xy::StaticVector{2}, Jj::BlockIndex{1})::T where T a,b = P.args J,j = Int(block(Jj)),blockindex(Jj) @@ -164,9 +166,9 @@ end ## sum -function Base._sum(P::RectPolynomial, dims) +function Base._sum(P::RectPolynomial, dims::Int) @assert dims == 1 - KronTrav(sum.(P.args; dims=1)...) + KronTrav(reverse(sum.(P.args; dims=1))...) end ## multiplication @@ -183,4 +185,15 @@ function layout_broadcasted(::Tuple{ExpansionLayout{KronOPLayout{2}},KronOPLayou end -broadcastbasis(::typeof(+), A::KronPolynomial, B::KronPolynomial) = KronPolynomial(broadcastbasis.(+, A.args, B.args)...) \ No newline at end of file +broadcastbasis(::typeof(+), A::KronPolynomial, B::KronPolynomial) = KronPolynomial(broadcastbasis.(+, A.args, B.args)...) + + + +#### +# vec +#### + +function vec_layout(::ExpansionLayout{KronOPLayout{2}}, f) + A,B = basis(f).args + A*invdiagtrav(coefficients(f))*B' +end \ No newline at end of file diff --git a/test/test_rect.jl b/test/test_rect.jl index b7cf90f..c5090e1 100644 --- a/test/test_rect.jl +++ b/test/test_rect.jl @@ -148,11 +148,12 @@ using Base: oneto end @testset "sum" begin - P = RectPolynomial(Legendre(),Legendre()) - pā‚€ = expand(P, š± -> 1) - @test sum(pā‚€) ā‰ˆ 4.0 - f = expand(P, splat((x,y) -> exp(cos(x^2*y)))) - @test sum(f) ā‰ˆ 10.546408460894801 # empirical + for P in (RectPolynomial(Legendre(),Legendre()), RectPolynomial(Legendre(),Chebyshev())) + pā‚€ = expand(P, š± -> 1) + @test sum(pā‚€) ā‰ˆ 4.0 + f = expand(P, splat((x,y) -> exp(cos(x^2*y)))) + @test sum(f) ā‰ˆ 10.546408460894801 # empirical + end end @testset "KronTrav bug" begin @@ -204,4 +205,10 @@ using Base: oneto šœ = expand(RectPolynomial(Legendre(),Jacobi(1,0)),splat((x,y) -> cos(x*sin(y)))) end + + @testset "vec" begin + P = RectPolynomial(Legendre(),Chebyshev()) + f = expand(P, splat((x,y) -> cos((x-0.1)*exp(y)))) + @test vec(f)[0.1,0.2] ā‰ˆ f[SVector(0.1,0.2)] + end end From 3f60520f2668de05541a940e576bc34529faa5e7 Mon Sep 17 00:00:00 2001 From: Sheehan Olver Date: Mon, 13 Oct 2025 15:18:49 +0100 Subject: [PATCH 2/3] reshape and vec --- Project.toml | 2 +- src/MultivariateOrthogonalPolynomials.jl | 6 ++--- src/rect.jl | 10 ++++++-- test/test_rect.jl | 30 ++++++++++++++++++++---- 4 files changed, 38 insertions(+), 10 deletions(-) diff --git a/Project.toml b/Project.toml index 85766f2..c71c254 100644 --- a/Project.toml +++ b/Project.toml @@ -30,7 +30,7 @@ BandedMatrices = "1" BlockArrays = "1.0" BlockBandedMatrices = "0.13" ClassicalOrthogonalPolynomials = "0.15" -ContinuumArrays = "0.19" +ContinuumArrays = "0.19.5" DomainSets = "0.7" FastTransforms = "0.17" FillArrays = "1.0" diff --git a/src/MultivariateOrthogonalPolynomials.jl b/src/MultivariateOrthogonalPolynomials.jl index a8684ed..4dd6aaa 100644 --- a/src/MultivariateOrthogonalPolynomials.jl +++ b/src/MultivariateOrthogonalPolynomials.jl @@ -10,8 +10,8 @@ import Base: axes, in, ==, +, -, /, *, ^, \, copy, copyto!, OneTo, getindex, siz import Base.Broadcast: Broadcasted, broadcasted, DefaultArrayStyle import DomainSets: boundary, EuclideanDomain -import QuasiArrays: LazyQuasiMatrix, LazyQuasiArrayStyle, domain, vec_layout -import ContinuumArrays: @simplify, Weight, weight, grid, plotgrid, TransformFactorization, ExpansionLayout, plotvalues, unweighted, plan_transform, checkpoints, transform_ldiv, AbstractBasisLayout, basis_axes, Inclusion, grammatrix, weaklaplacian, layout_broadcasted, laplacian, abslaplacian, laplacian_axis, abslaplacian_axis, diff_layout, operatororder, broadcastbasis +import QuasiArrays: LazyQuasiMatrix, LazyQuasiArrayStyle, domain, vec_layout, reshape_layout +import ContinuumArrays: @simplify, Weight, weight, grid, plotgrid, TransformFactorization, ExpansionLayout, plotvalues, unweighted, plan_transform, checkpoints, transform_ldiv, AbstractBasisLayout, basis_axes, Inclusion, grammatrix, weaklaplacian, layout_broadcasted, laplacian, abslaplacian, laplacian_axis, abslaplacian_axis, diff_layout, operatororder, broadcastbasis, KronExpansionLayout import ArrayLayouts: MemoryLayout, sublayout, sub_materialize import BlockArrays: block, blockindex, BlockSlice, viewblock, blockcolsupport, AbstractBlockStyle, BlockStyle @@ -21,7 +21,7 @@ import LazyArrays: arguments, paddeddata, LazyArrayStyle, LazyLayout, PaddedLayo import LazyBandedMatrices: LazyBandedBlockBandedLayout, AbstractBandedBlockBandedLayout, AbstractLazyBandedBlockBandedLayout, _krontrav_axes, DiagTravLayout, invdiagtrav, ApplyBandedBlockBandedLayout, krontrav import InfiniteArrays: InfiniteCardinal, OneToInf -import ClassicalOrthogonalPolynomials: jacobimatrix, Weighted, orthogonalityweight, HalfWeighted, WeightedBasis, pad, recurrencecoefficients, clenshaw, weightedgrammatrix, Clenshaw +import ClassicalOrthogonalPolynomials: jacobimatrix, Weighted, orthogonalityweight, HalfWeighted, WeightedBasis, pad, recurrencecoefficients, clenshaw, weightedgrammatrix, Clenshaw, OPLayout import HarmonicOrthogonalPolynomials: BivariateOrthogonalPolynomial, MultivariateOrthogonalPolynomial, Plan, AngularMomentum, angularmomentum, BlockOneTo, BlockRange1, interlace, MultivariateOPLayout, AbstractMultivariateOPLayout, MAX_PLOT_BLOCKS diff --git a/src/rect.jl b/src/rect.jl index 8c03e1c..acd17d3 100644 --- a/src/rect.jl +++ b/src/rect.jl @@ -190,10 +190,16 @@ broadcastbasis(::typeof(+), A::KronPolynomial, B::KronPolynomial) = KronPolynom #### -# vec +# reshape/vec #### -function vec_layout(::ExpansionLayout{KronOPLayout{2}}, f) +function reshape_layout(::ExpansionLayout{KronOPLayout{2}}, f) A,B = basis(f).args A*invdiagtrav(coefficients(f))*B' +end + +# TODO: generalise +function vec_layout(::KronExpansionLayout{OPLayout, OPLayout}, f) + A,C,Bc = arguments(f) + RectPolynomial(A, Bc') * DiagTrav(C) end \ No newline at end of file diff --git a/test/test_rect.jl b/test/test_rect.jl index c5090e1..43b820c 100644 --- a/test/test_rect.jl +++ b/test/test_rect.jl @@ -1,7 +1,7 @@ using MultivariateOrthogonalPolynomials, ClassicalOrthogonalPolynomials, StaticArrays, LinearAlgebra, BlockArrays, FillArrays, Base64, LazyBandedMatrices, Test using ClassicalOrthogonalPolynomials: expand, coefficients, recurrencecoefficients using MultivariateOrthogonalPolynomials: weaklaplacian, ClenshawKron -using ContinuumArrays: plotgridvalues +using ContinuumArrays: plotgridvalues, ExpansionLayout using Base: oneto @testset "RectPolynomial" begin @@ -206,9 +206,31 @@ using Base: oneto šœ = expand(RectPolynomial(Legendre(),Jacobi(1,0)),splat((x,y) -> cos(x*sin(y)))) end - @testset "vec" begin - P = RectPolynomial(Legendre(),Chebyshev()) + @testset "reshape/vec" begin + P = RectPolynomial(Legendre(),Chebyshev()) f = expand(P, splat((x,y) -> cos((x-0.1)*exp(y)))) - @test vec(f)[0.1,0.2] ā‰ˆ f[SVector(0.1,0.2)] + F = reshape(f) + @test F[0.1,0.2] ā‰ˆ f[SVector(0.1,0.2)] + @test vec(F)[SVector(0.1,0.2)] ā‰ˆ f[SVector(0.1,0.2)] + + g = F[:,0.2] + h = F[0.1,:] + @test MemoryLayout(g) isa ExpansionLayout + @test MemoryLayout(h) isa ExpansionLayout + @test g[0.1] ā‰ˆ f[SVector(0.1,0.2)] + @test h[0.2] ā‰ˆ f[SVector(0.1,0.2)] + + @test sum(F; dims=1)[1,0.2] ā‰ˆ exp(-0.2)*(sin(0.9exp(0.2)) + sin(1.1exp(0.2))) + # TODO: should be matrix but isn't because of InfiniteArrays/src/reshapedarray.jl:77 + @test_broken sum(F; dims=2)[0.1,1] ā‰ˆ 2 + @test sum(F; dims=2)[0.1] ā‰ˆ 2 + end + + @testset "sample" begin + P = RectPolynomial(Legendre(),Legendre()) + f = expand(P, splat((x,y) -> exp(x*cos(y-0.1)))) + F = reshape(f) + @test sum(F; dims=1)[1,0.2] ā‰ˆ 2.346737615950585 + @test sum(F; dims=2)[0.1,1] ā‰ˆ 2.1748993079723618 end end From efd3c180b7133140430f259c1528bcdd8a9a296b Mon Sep 17 00:00:00 2001 From: Sheehan Olver Date: Mon, 13 Oct 2025 20:51:34 +0100 Subject: [PATCH 3/3] Update test_rect.jl --- test/test_rect.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_rect.jl b/test/test_rect.jl index 43b820c..dfd68a5 100644 --- a/test/test_rect.jl +++ b/test/test_rect.jl @@ -1,4 +1,4 @@ -using MultivariateOrthogonalPolynomials, ClassicalOrthogonalPolynomials, StaticArrays, LinearAlgebra, BlockArrays, FillArrays, Base64, LazyBandedMatrices, Test +using MultivariateOrthogonalPolynomials, ClassicalOrthogonalPolynomials, StaticArrays, LinearAlgebra, BlockArrays, FillArrays, Base64, LazyBandedMatrices, ArrayLayouts, Test using ClassicalOrthogonalPolynomials: expand, coefficients, recurrencecoefficients using MultivariateOrthogonalPolynomials: weaklaplacian, ClenshawKron using ContinuumArrays: plotgridvalues, ExpansionLayout