Add plotgrid and fix grid (#38)

* Add plotgrid and fix grid

* Fix grammatrix for ContinuousPolynomial{0}

* fix test

* Update test_dirichlet.jl

* Update test_dirichlet.jl

* update transform

* Update arrowhead.jl

* Support singularities

* Update Project.toml

* Update Project.toml

* simplify arrowhead show

* InfiniteLinearAlgebra v0.7

* Update continuouspolynomial.jl

* Update Project.toml

* add tests

* Update test_piecewisepolynomial.jl

* DirichletPolynomial can preserve vanishing

Author: dlfivefifty

Project.toml

@@ -1,7 +1,7 @@
 name = "PiecewiseOrthogonalPolynomials"
 uuid = "4461d12d-4663-4550-8580-cb764c85e20f"
 authors = ["Sheehan Olver <[email protected]>"]
-version = "0.2.0"
+version = "0.2.1"
 
 [deps]
 ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
@@ -25,12 +25,12 @@ BandedMatrices = "0.17.33"
 BlockArrays = "0.16.25"
 BlockBandedMatrices = "0.12.2"
 ClassicalOrthogonalPolynomials = "0.11"
-ContinuumArrays = "0.14"
+ContinuumArrays = "0.15"
 FillArrays = "1.0"
 InfiniteArrays = "0.12.6, 0.13"
-InfiniteLinearAlgebra = "0.6.16"
+InfiniteLinearAlgebra = "0.6.16, 0.7"
 LazyArrays = "1.0"
-LazyBandedMatrices = "0.8.15"
+LazyBandedMatrices = "0.9.1"
 MatrixFactorizations = "2.0"
 QuasiArrays = "0.11"
 julia = "1.9"

examples/coulomb.jl

@@ -1,12 +1,83 @@
-using PiecewiseOrthogonalPolynomials, Plots
+using PiecewiseOrthogonalPolynomials, ClassicalOrthogonalPolynomials, BandedMatrices, Plots
+using Base: oneto
 
-p = reverse(2.0 .^ (-(0:10)))
+let M = 20, R = 100
+    global r = reverse(2.0 .^ (-(0:M))); r = R*[-reverse(r); 0; r]
+    global V = x -> max(-1/abs(x), -2.0 ^(M)/R)
+end
 
-C = ContinuousPolynomial{1}(p)
-x = axes(C,1)
-a = C / C \ log.(x)
+Q = DirichletPolynomial(r)
+P = ContinuousPolynomial{0}(Q)
 
+a = expand(P, V)
+
+
+
+p = 20; KR = Block.(oneto(p))
+Δ = (diff(Q)'diff(Q))[KR,KR];
+M = grammatrix(Q)[KR,KR];
+
+R = P\Q
+@time A = (R'* P' * (a .* P) *R)[KR,KR];
+
+
+λ,V = eigen((Matrix(Δ) + 10A), (Matrix(M)))
+
+
+
+g = range(-100,100; length=10000)
+p = plot(g, 10a[g]; legend=false, ylims=(-150,5), linestyle=:dash, linewidth=2)
+scatter!(r, 10a[r])
+savefig("coulombpotential.pdf")
+
+p = plot(g, 10a[g]; legend=false, ylims=(-30,5), linestyle=:dash, linewidth=2)
+Pl = Q[g,KR];
+for j = 1:100; plot!(g, Pl * real(V[:,j]) .+ λ[j], linewidth=2) end
+p
+
+savefig("coulomb.pdf")
+
+g = range(-5,5; length=1000)
+p = plot(g, 10a[g]; legend=false, ylims=(-40,0), linestyle=:dash, linewidth=2)
+Pl = Q[g,KR];
+for j = 1:10; plot!(g, Pl * real(V[:,j]) .+ λ[j], linewidth=2) end
+p
+
+savefig("coulomb.pdf")
+
+p2 = plot(g, 10a[g]; legend=false, ylims=(-2800,-2700), linestyle=:dash, linewidth=2)
+
+nanabs = x -> iszero(x) ? NaN : abs(x)
+scatter(nanabs.(10coefficients(a)[1:size(V,1)]); yscale=:log10, legend=false)
+
+KR = Block.(oneto(30)); scatter(10nanabs.(coefficients(expand(P[:,KR], V))[KR]); yscale=:log10, legend=false, yticks=10.0 .^(-20:5:5))
+savefig("coulombcoeffs.pdf")
+scatter!(nanabs.(V[:,2]); yscale=:log10)
+
+
+a = expand(P, x -> x^2)
+Δ = -weaklaplacian(Q)[KR,KR]
+M = grammatrix(Q)[KR,KR]
+
+R = P\Q
+@time A = (R'* P' * (a .* P) *R)[KR,KR];
+λ,V = eigen((Matrix(Δ) + 1000A), (Matrix(M)))
+
+
+plot(Q[:,KR] * ((Δ + 1000A) \ [1; zeros(size(Δ,1)-1)]))
+p
+
+
+r = range(-1,1; length=10)
+Q = DirichletPolynomial(r)
+λ,V = eigen(Matrix((-weaklaplacian(Q))[KR,KR]), Matrix(grammatrix(Q)[KR,KR]))
+
+plot(Q[:,KR] * V[:,4])
+
+
+
+
+
+ContinuousPolynomial{0}(r) \ ContinuousPolynomial{1}(r)
 
-a = C[:,1]
 
-C \ (a .* C[:,1])

src/PiecewiseOrthogonalPolynomials.jl

@@ -5,11 +5,11 @@ import ArrayLayouts: sublayout, sub_materialize, symmetriclayout, transposelayou
 import BandedMatrices: _BandedMatrix
 import BlockArrays: BlockSlice, block, blockindex, blockvec
 import BlockBandedMatrices: _BandedBlockBandedMatrix, AbstractBandedBlockBandedMatrix, subblockbandwidths, blockbandwidths, AbstractBandedBlockBandedLayout, layout_replace_in_print_matrix
-import ClassicalOrthogonalPolynomials: grid, ldiv, pad, adaptivetransform_ldiv, grammatrix
+import ClassicalOrthogonalPolynomials: grid, plotgrid, ldiv, pad, adaptivetransform_ldiv, grammatrix, Plan, singularities, basis_singularities, singularitiesbroadcast
 import ContinuumArrays: @simplify, factorize, TransformFactorization, AbstractBasisLayout, MemoryLayout, layout_broadcasted, ExpansionLayout, basis, plan_grid_transform, grammatrix
 import LazyArrays: paddeddata
 import LazyBandedMatrices: BlockBroadcastMatrix, BlockVec, BandedLazyLayouts, AbstractLazyBandedBlockBandedLayout, UpperOrLowerTriangular
-import Base: axes, getindex, +, -, *, /, ==, \, OneTo, oneto, replace_in_print_matrix, copy, diff, getproperty, adjoint, transpose, tail
+import Base: axes, getindex, +, -, *, /, ==, \, OneTo, oneto, replace_in_print_matrix, copy, diff, getproperty, adjoint, transpose, tail, _sum
 import LinearAlgebra: BlasInt
 import MatrixFactorizations: reversecholcopy
 

src/arrowhead.jl

@@ -25,6 +25,8 @@ const ArrowheadMatrices = Union{ArrowheadMatrix,Symmetric{<:Any,<:ArrowheadMatri
 
 subblockbandwidths(A::ArrowheadMatrices) = (1,1)
 
+BlockArrays._show_typeof(io::IO, B::ArrowheadMatrix{T}) where T = print(io, "ArrowheadMatrix{$T}")
+
 function blockbandwidths(A::ArrowheadMatrix)
     l,u = bandwidths(A.D[1])
     max(l,length(A.C)),max(u,length(A.B))
@@ -193,10 +195,10 @@ function materialize!(M::MatMulMatAdd{<:ArrowheadLayouts,<:AbstractColumnMajor,<
     _fill_lmul!(β, Y)
     for J = blockaxes(X,2)
         mul!(view(Y, Block(1), J), A.A, view(X, Block(1), J), α, one(α))
-        for k = 1:length(A.B)
+        for k = 1:min(length(A.B), blocksize(X,1)-1)
             mul!(view(Y, Block(1), J), A.B[k], view(X, Block(k+1), J), α, one(α))
         end
-        for k = 1:length(A.C)
+        for k = 1:min(length(A.C), blocksize(Y,1)-1)
             mul!(view(Y, Block(k+1), J), A.C[k], view(X, Block(1), J), α, one(α))
         end
     end

src/continuouspolynomial.jl

@@ -8,7 +8,7 @@ ContinuousPolynomial{o}(pts) where {o} = ContinuousPolynomial{o,Float64}(pts)
 ContinuousPolynomial{o,T}(P::ContinuousPolynomial) where {o,T} = ContinuousPolynomial{o,T}(P.points)
 ContinuousPolynomial{o}(P::ContinuousPolynomial) where {o} = ContinuousPolynomial{o,eltype(P)}(P)
 
-PiecewisePolynomial(P::ContinuousPolynomial{0,T}) where {T} = PiecewisePolynomial(Legendre{T}(), P.points)
+PiecewisePolynomial(P::ContinuousPolynomial{o,T}) where {o,T} = PiecewisePolynomial(Legendre{T}(), P.points)
 
 axes(B::ContinuousPolynomial{0}) = axes(PiecewisePolynomial(B))
 axes(B::ContinuousPolynomial{1}) =
@@ -44,10 +44,12 @@ end
 
 factorize(V::SubQuasiArray{T,N,<:ContinuousPolynomial{0},<:Tuple{Inclusion,BlockSlice}}, dims...) where {T,N} =
     factorize(view(PiecewisePolynomial(parent(V)), parentindices(V)...), dims...)
-grid(V::SubQuasiArray{T,N,<:ContinuousPolynomial{0},<:Tuple{Inclusion,Any}}) where {T,N} =
-    grid(view(PiecewisePolynomial(parent(V)), parentindices(V)...))
-grid(V::SubQuasiArray{T,N,<:ContinuousPolynomial,<:Tuple{Inclusion,Any}}) where {T,N} =
-    grid(view(ContinuousPolynomial{0,T}(parent(V)), parentindices(V)...))    
+
+plan_grid_transform(P::ContinuousPolynomial{0}, args...) = plan_grid_transform(PiecewisePolynomial(P), args...)
+
+for grd in (:grid, :plotgrid)
+    @eval $grd(C::ContinuousPolynomial, n...) = $grd(PiecewisePolynomial(C), n...)
+end
 
 function adaptivetransform_ldiv(Q::ContinuousPolynomial{1,V}, f::AbstractQuasiVector) where V
     T = promote_type(V, eltype(f))
@@ -94,23 +96,6 @@ end
 #######
 
 function \(P::ContinuousPolynomial{0}, C::ContinuousPolynomial{1})
-    T = promote_type(eltype(P), eltype(C))
-    # diag blocks based on
-    # L = Legendre{T}() \ Weighted(Jacobi{T}(1,1))
-    @assert P.points == C.points
-    N = length(P.points)
-    v = mortar(Fill.((convert(T, 2):2:∞) ./ (3:2:∞), N - 1))
-    z = Zeros{T}(axes(v))
-    H1 = BlockBroadcastArray(hcat, z, v)
-    M1 = BlockVcat(Zeros{T}(N, 2), H1)
-    M2 = BlockVcat(Ones{T}(N, 2) / 2, Zeros{T}((axes(v, 1), Base.OneTo(2))))
-    H3 = BlockBroadcastArray(hcat, z, -v)
-    M3 = BlockVcat(Hcat(Ones{T}(N) / 2, -Ones{T}(N) / 2), H3)
-    dat = BlockHcat(M1, M2, M3)'
-    _BandedBlockBandedMatrix(dat, axes(P, 2), (1, 1), (0, 1))
-end
-
-function \(P::ContinuousPolynomial{0, <:Any, <:AbstractRange}, C::ContinuousPolynomial{1, <:Any, <:AbstractRange})
     T = promote_type(eltype(P), eltype(C))
     @assert P.points == C.points
     v = (convert(T, 2):2:∞) ./ (3:2:∞)
@@ -132,7 +117,7 @@ function grammatrix(A::ContinuousPolynomial{0,T}) where T
     N = length(r) - 1
     hs = diff(r)
     M = grammatrix(Legendre{T}())
-    ArrowheadMatrix{T}(Diagonal(Fill(hs[1], N)), (), (), [Diagonal(M.diag[2:end] * h/2) for h in hs])
+    ArrowheadMatrix{T}(Diagonal(hs), (), (), [Diagonal(M.diag[2:end] * h/2) for h in hs])
 end
 
 function grammatrix(A::ContinuousPolynomial{0,T, <:AbstractRange}) where T
@@ -196,7 +181,7 @@ function diff(C::ContinuousPolynomial{1,T}; dims=1) where T
     H = BlockBroadcastArray(hcat, z, v)
     M = BlockVcat(Hcat(Ones{T}(N) .* [zero(T); s] , -Ones{T}(N) .* [s; zero(T)] ), H)
     P = ContinuousPolynomial{0}(C)
-    ApplyQuasiMatrix(*, P, _BandedBlockBandedMatrix(M', (axes(P, 2), axes(C, 2)), (0, 0), (0, 1)))
+    ApplyQuasiMatrix(*, P, _BandedBlockBandedMatrix(M', axes(P, 2), (0, 0), (0, 1)))
 end
 
 function weaklaplacian(C::ContinuousPolynomial{1,T,<:AbstractRange}) where T
@@ -210,3 +195,12 @@ function weaklaplacian(C::ContinuousPolynomial{1,T,<:AbstractRange}) where T
         Fill(Diagonal(convert(T, -16) .* (1:∞) .^ 2 ./ (s .* ((2:2:∞) .+ 1))), N-1)))
 end
 
+
+
+###
+# singularities
+###
+
+singularities(C::ContinuousPolynomial{λ}) where λ = C
+basis_singularities(C::ContinuousPolynomial) = C
+singularitiesbroadcast(_, C::ContinuousPolynomial) = C # Assume we stay smooth

src/dirichlet.jl

@@ -43,6 +43,11 @@ function \(Q::DirichletPolynomial, C::ContinuousPolynomial{1})
     ArrowheadMatrix(layout_getindex(Eye{T}(m), 2:m-1, :), (), (), Fill(Eye{T}(∞), m-1))
 end
 
+function \(P::ContinuousPolynomial, Q::DirichletPolynomial)
+    C = ContinuousPolynomial{1}(Q)
+    (P \ C) * (C \ Q)
+end
+
 
 function adaptivetransform_ldiv(Q::DirichletPolynomial{V}, f::AbstractQuasiVector) where V
     C = ContinuousPolynomial{1}(Q)
@@ -78,6 +83,12 @@ function grammatrix(Q::DirichletPolynomial{T, <:AbstractRange}) where T
                 (h*b/2)'), ∞, 0, 2), N)))
 end
 
+function grammatrix(Q::DirichletPolynomial)
+    C = ContinuousPolynomial{1}(Q)
+    R = C \ Q
+    R' * grammatrix(C) * R
+end
+
 @simplify function *(Ac::QuasiAdjoint{<:Any,<:DirichletPolynomial}, B::ContinuousPolynomial)
     A = Ac'
     C = ContinuousPolynomial{1}(A)
@@ -87,4 +98,17 @@ end
 function diff(Q::DirichletPolynomial{T}; dims=1) where T
     C = ContinuousPolynomial{1}(Q)
     diff(C; dims=dims) * (C \ Q)
-end
+end
+
+for grd in (:grid, :plotgrid)
+    @eval $grd(C::DirichletPolynomial, n...) = $grd(PiecewisePolynomial(C), n...)
+end
+
+###
+# singularities
+###
+
+singularities(C::DirichletPolynomial) = C
+basis_singularities(C::DirichletPolynomial) = C
+singularitiesbroadcast(_, Q::DirichletPolynomial) = ContinuousPolynomial{1}(Q) # Assume we stay smooth but might not vanish
+singularitiesbroadcast(::typeof(sin), Q::DirichletPolynomial) = Q # Smooth functions such that f(0) == 0 preserve behaviour

src/piecewisepolynomial.jl

@@ -25,28 +25,31 @@ function repeatgrid(ax, g, pts)
     ret
 end
 
-function grid(V::SubQuasiArray{T,2,<:PiecewisePolynomial,<:Tuple{Inclusion,BlockSlice}}) where {T}
-    P = parent(V)
-    _, JR = parentindices(V)
-    N = Int(last(JR))
-    g = grid(P.basis[:, OneTo(N)])
+
+function grid(P::PiecewisePolynomial, N::Block{1})
+    g = grid(P.basis, Int(N))
     repeatgrid(axes(P.basis, 1), g, P.points)
 end
 
-function grid(V::SubQuasiArray{T,N,<:PiecewisePolynomial,<:Tuple{Inclusion,Any}}) where {T,N}
-    P = parent(V)
-    kr,jr = parentindices(V)
-    J = findblock(axes(P,2), last(jr))
-    grid(view(P, kr, Block(1):J))
+function plotgrid(P::PiecewisePolynomial, N::Block{1})
+    g = plotgrid(P.basis, Int(N))
+    vec(repeatgrid(axes(P.basis, 1), g, P.points)[end:-1:1,:]) # sort
 end
 
 
-struct ApplyFactorization{T, FF, FAC<:Factorization{T}} <: Factorization{T}
+
+grid(P::PiecewisePolynomial, n::Int) = grid(P, findblock(axes(P,2),n))
+plotgrid(P::PiecewisePolynomial, n::Int) = plotgrid(P, findblock(axes(P,2),n))
+
+
+
+
+struct ApplyPlan{T, FF, FAC<:Plan{T}} <: Plan{T}
     f::FF
     F::FAC
 end
 
-\(P::ApplyFactorization, f) = P.f(P.F \ f)
+*(P::ApplyPlan, f::AbstractArray) = P.f(P.F * f)
 
 
 _perm_blockvec(X::AbstractMatrix) = BlockVec(transpose(X))
@@ -60,12 +63,15 @@ function _perm_blockvec(X::AbstractArray{T,3}) where T
     ret
 end
 
+function plan_grid_transform(P::PiecewisePolynomial, N::Block{1}, dims...)
+    x,F = plan_grid_transform(P.basis, (Int(N), length(P.points)-1, dims...), 1)
+    repeatgrid(axes(P.basis, 1), x, P.points), ApplyPlan(_perm_blockvec, F)
+end
+
 function factorize(V::SubQuasiArray{<:Any,2,<:PiecewisePolynomial,<:Tuple{Inclusion,BlockSlice}}, dims...)
     P = parent(V)
     _,JR = parentindices(V)
-    N = Int(last(JR.block))
-    x,F = plan_grid_transform(P.basis, Array{eltype(P)}(undef, N, length(P.points)-1, dims...), 1)
-    ApplyFactorization(_perm_blockvec, TransformFactorization(repeatgrid(axes(P.basis, 1), x, P.points), F))
+    TransformFactorization(plan_grid_transform(P, last(JR.block), dims...)...)
 end
 
 
@@ -105,3 +111,12 @@ function layout_broadcasted(::Tuple{ExpansionLayout{PiecewisePolynomialLayout{0}
     P = ContinuousPolynomial{0}(C)
     (a .* P) * (P \ C)
 end
+
+
+###
+# singularities
+###
+
+singularities(C::PiecewisePolynomial) = C
+basis_singularities(C::PiecewisePolynomial) = C
+singularitiesbroadcast(_, C::PiecewisePolynomial) = C # Assume we stay piecewise smooth

test/test_arrowhead.jl

@@ -103,7 +103,7 @@ import Base: oneto, OneTo
         @test blockbandwidths(L) == (1,1)
         @test subblockbandwidths(L) == (1,1)
 
-        @test stringmime("text/plain", L[Block.(OneTo(3)), Block.(OneTo(3))]) == "3×3-blocked 9×10 ArrowheadMatrix{Float64, BandedMatrix{Float64, Fill{Float64, 2, Tuple{OneTo{$Int}, OneTo{$Int}}}, OneTo{$Int}}, Tuple{BandedMatrix{Float64, Fill{Float64, 2, Tuple{OneTo{$Int}, OneTo{$Int}}}, OneTo{$Int}}}, Tuple{BandedMatrix{Float64, LazyArrays.ApplyArray{Float64, 2, typeof(vcat), Tuple{Fill{Float64, 2, Tuple{OneTo{$Int}, OneTo{$Int}}}, Fill{Float64, 2, Tuple{OneTo{$Int}, OneTo{$Int}}}}}, OneTo{$Int}}}, Vector{BandedMatrix{Float64, Matrix{Float64}, OneTo{$Int}}}}:\n  0.5   0.5    ⋅    ⋅   │   0.666667    ⋅          ⋅        │   ⋅    ⋅    ⋅ \n   ⋅    0.5   0.5   ⋅   │    ⋅         0.666667    ⋅        │   ⋅    ⋅    ⋅ \n   ⋅     ⋅    0.5  0.5  │    ⋅          ⋅         0.666667  │   ⋅    ⋅    ⋅ \n ───────────────────────┼───────────────────────────────────┼───────────────\n -0.5   0.5    ⋅    ⋅   │   0.0         ⋅          ⋅        │  0.8   ⋅    ⋅ \n   ⋅   -0.5   0.5   ⋅   │    ⋅         0.0         ⋅        │   ⋅   0.8   ⋅ \n   ⋅     ⋅   -0.5  0.5  │    ⋅          ⋅         0.0       │   ⋅    ⋅   0.8\n ───────────────────────┼───────────────────────────────────┼───────────────\n   ⋅     ⋅     ⋅    ⋅   │  -0.666667    ⋅          ⋅        │  0.0   ⋅    ⋅ \n   ⋅     ⋅     ⋅    ⋅   │    ⋅        -0.666667    ⋅        │   ⋅   0.0   ⋅ \n   ⋅     ⋅     ⋅    ⋅   │    ⋅          ⋅        -0.666667  │   ⋅    ⋅   0.0"
+        @test stringmime("text/plain", L[Block.(OneTo(3)), Block.(OneTo(3))]) == "3×3-blocked 9×10 ArrowheadMatrix{Float64}:\n  0.5   0.5    ⋅    ⋅   │   0.666667    ⋅          ⋅        │   ⋅    ⋅    ⋅ \n   ⋅    0.5   0.5   ⋅   │    ⋅         0.666667    ⋅        │   ⋅    ⋅    ⋅ \n   ⋅     ⋅    0.5  0.5  │    ⋅          ⋅         0.666667  │   ⋅    ⋅    ⋅ \n ───────────────────────┼───────────────────────────────────┼───────────────\n -0.5   0.5    ⋅    ⋅   │   0.0         ⋅          ⋅        │  0.8   ⋅    ⋅ \n   ⋅   -0.5   0.5   ⋅   │    ⋅         0.0         ⋅        │   ⋅   0.8   ⋅ \n   ⋅     ⋅   -0.5  0.5  │    ⋅          ⋅         0.0       │   ⋅    ⋅   0.8\n ───────────────────────┼───────────────────────────────────┼───────────────\n   ⋅     ⋅     ⋅    ⋅   │  -0.666667    ⋅          ⋅        │  0.0   ⋅    ⋅ \n   ⋅     ⋅     ⋅    ⋅   │    ⋅        -0.666667    ⋅        │   ⋅   0.0   ⋅ \n   ⋅     ⋅     ⋅    ⋅   │    ⋅          ⋅        -0.666667  │   ⋅    ⋅   0.0"
     end
 
     @testset "Helmholtz solve" begin