current state

src/SemiclassicalOrthogonalPolynomials.jl

@@ -420,11 +420,23 @@ end
 
 function \(A::SemiclassicalJacobi, B::SemiclassicalJacobi{T}) where {T} 
     if A.b == -1 && B.b ≠ -1
-        return ApplyArray(inv, B \ A)
+        #=
+        Too many memory layout issues - colsupport(inv(B \ A), 1) == 1:∞ without wrapping UpperTriangulation since 
+            LazyArrays.applylayout(typeof(*), 
+                TriangularLayout{'U', 'N', LazyArrays.LazyLayout}(), 
+                TriangularLayout{'U', 'N', LazyArrays.LazyLayout}(), 
+                BidiagonalLayout{LazyArrays.LazyLayout, LazyArrays.LazyLayout}()
+            ) = LazyArrays.ApplyLayout{typeof(*)}().
+        Additionally, since there's a method ambiguity (that could be fixed using 
+            ext = Base.get_extension(LazyArrays, :LazyArraysBandedMatricesExt)
+            Base.copy(L::Ldiv{<:ext.BandedLazyLayouts,<:LazyArrays.AbstractLazyLayout}) = LazyArrays.lazymaterialize(\, L.A, L.B)
+        ), we use ApplyArray instead of inv(B \ A).
+        =#
+        return UpperTriangular(ApplyArray(inv, B \ A)) 
     elseif B.b == -1 && A.b ≠ -1
         # First convert Bᵗᵃ⁻¹ᶜ into Bᵗᵃ⁰ᶜ
-        Bᵗᵃ⁰ᶜ = SemiclassicalJacobi(B.t, B.a, zero(B.b), B.c) 
-        Bᵗᵃ¹ᶜ = SemiclassicalJacobi(B.t, B.a, one(B.a), B.c)
+        Bᵗᵃ⁰ᶜ = SemiclassicalJacobi(B.t, B.a, zero(B.b), B.c, A) 
+        Bᵗᵃ¹ᶜ = SemiclassicalJacobi(B.t, B.a, one(B.a), B.c, A)
         Rᵦₐ₁ᵪᵗᵃ⁰ᶜ = Weighted(Bᵗᵃ⁰ᶜ) \ Weighted(Bᵗᵃ¹ᶜ)
         b1 = Rᵦₐ₁ᵪᵗᵃ⁰ᶜ[band(0)]
         b0 = Vcat(one(T), Rᵦₐ₁ᵪᵗᵃ⁰ᶜ[band(-1)])

src/derivatives.jl

@@ -55,7 +55,8 @@ function diff(P::SemiclassicalJacobi{T}; dims=1) where {T}
         Dₐ₀ᵪᵃ⁺¹¹ᶜ⁺¹ = diff(Pᵗᵃ⁰ᶜ)
         Pᵗᵃ⁺¹¹ᶜ⁺¹ = Dₐ₀ᵪᵃ⁺¹¹ᶜ⁺¹.args[1]
         Pᵗᵃ⁺¹⁰ᶜ⁺¹ = SemiclassicalJacobi(P.t, P.a + 1, zero(P.b), P.c + 1, Pᵗᵃ⁰ᶜ)
-        Rₐ₊₁₁ᵪ₊₁ᵗᵃ⁺¹⁰ᶜ⁺¹ = ApplyArray(inv, Pᵗᵃ⁺¹¹ᶜ⁺¹ \ Pᵗᵃ⁺¹⁰ᶜ⁺¹) #  
+        # Rₐ₊₁₁ᵪ₊₁ᵗᵃ⁺¹⁰ᶜ⁺¹ = ApplyArray(inv, Pᵗᵃ⁺¹¹ᶜ⁺¹ \ Pᵗᵃ⁺¹⁰ᶜ⁺¹) #  
+        Rₐ₊₁₁ᵪ₊₁ᵗᵃ⁺¹⁰ᶜ⁺¹ = Pᵗᵃ⁺¹⁰ᶜ⁺¹ \ Pᵗᵃ⁺¹¹ᶜ⁺¹ 
         Dₐ₋₁ᵪᵃ⁺¹⁰ᶜ⁺¹ = Rₐ₊₁₁ᵪ₊₁ᵗᵃ⁺¹⁰ᶜ⁺¹ * coefficients(Dₐ₀ᵪᵃ⁺¹¹ᶜ⁺¹) * Rᵦₐ₁ᵪᵗᵃ⁰ᶜ # bidigaonalconjugation data and avoid inv
         b2 = Vcat(zero(T), zero(T), Dₐ₋₁ᵪᵃ⁺¹⁰ᶜ⁺¹[band(1)])
         b1 = Vcat(zero(T), Dₐ₋₁ᵪᵃ⁺¹⁰ᶜ⁺¹[band(0)])
@@ -65,9 +66,6 @@ function diff(P::SemiclassicalJacobi{T}; dims=1) where {T}
     end
 end
 
-
-
-
 function divdiff(wP::Weighted{<:Any,<:SemiclassicalJacobi}, wQ::Weighted{<:Any,<:SemiclassicalJacobi})
     Q,P = wQ.P,wP.P
     ((-sum(orthogonalityweight(Q))/sum(orthogonalityweight(P))) * (Q \ diff(P))')

test/runtests.jl

@@ -161,7 +161,7 @@ end
         @testset "Mass matrix" begin
             @test (T'*(w_T .* T))[1:10,1:10] ≈ sum(w_T)I
             M = W'*(w_W .* W);
-            @test_broken (sum(w_T)*inv(R)'L)[1:10,1:10] ≈ M[1:10,1:10]
+            @test (sum(w_T)*inv(R)'L)[1:10,1:10] ≈ M[1:10,1:10]
             @test T[0.1,1:10]' ≈ W[0.1,1:10]' * R[1:10,1:10]
         end
     end

test/test_neg1b.jl

@@ -39,57 +39,116 @@ end
     end
 end
 
-@testset "Families" begin 
+@testset "Families" begin
     # Not implemented efficiently currently. 
     # TODO: Fix cholesky_jacobimatrix when Δb == 0 and b == -1 so that building families works (currently, it would return a "Polynomials must be orthonormal" error)
-    t = 2.0 
+    t = 2.0
     P = SemiclassicalJacobi.(t, -1//2:13//2, -1.0, -1//2:13//2) # // so that we get a UnitRange, else broadcasting into a Family fails
     for (i, p) in enumerate(P)
         @test jacobimatrix(p)[1:100, 1:100] == jacobimatrix(SemiclassicalJacobi(t, (-1//2:13//2)[i], -1.0, (-1//2:13//2)[i]))[1:100, 1:100]
     end
 end
 
+#=
 @testset "Expansions" begin
     Ps = SemiclassicalJacobi.(2, -1//2:5//2, -1.0, -1//2:5//2)
     for P in Ps
-        @show 1
         for (idx, g) in enumerate((x -> exp(x) + sin(x), x -> (1 - x) * cos(x^3), x -> 5.0 + (1 - x)))
+            @show idx
             f = expand(P, g)
             for x in LinRange(0, 1, 100)
-                @test f[x] ≈ g(x) atol=1e-9
+                @test f[x] ≈ g(x) atol = 1e-9
             end
             x = axes(P, 1)
             @test P[:, 1:20] \ g.(x) ≈ coefficients(f)[1:20]
-            if idx == 2 
+            if idx == 2
                 @test coefficients(f)[1] ≈ 0 atol = 1e-9
-            elseif idx == 3 
+            elseif idx == 3
                 @test coefficients(f)[1:2] ≈ [5.0, 1.0]
-                @test coefficients(f)[3:1000] ≈ zeros(1000-3+1)
+                @test coefficients(f)[3:1000] ≈ zeros(1000 - 3 + 1)
             end
         end
-    end 
+    end
 end
+=#
+
+Ps = SemiclassicalJacobi.(2, -1//2:5//2, -1.0, -1//2:5//2)
+g = x -> exp(x) + sin(x)
+expand(Ps[1], g) # works 
+expand(Ps[2], g) # works 
+expand(Ps[3], g) 
+
+P = SemiclassicalJacobi(2.0, 1.5, -1.0, 1.5)
+g = x -> exp(x) + sin(x)
+f = g.(axes(P, 1))
+F = P * (P \ f)
+
+R = SemiclassicalJacobi(P.t, mod(P.a,-1), mod(P.b,-1), mod(P.c,-1))
+R̃ = SemiclassicalOrthogonalPolynomials.toclassical(R)
+P \ 
+
+using ArrayLayouts
+A, B = Q, Q \ f
+# M = Mul(A, B)[0.2]
+
+R̃ = SemiclassicalOrthogonalPolynomials.toclassical(SemiclassicalJacobi(Q.t, mod(Q.a, -1), mod(Q.b, -1), mod(Q.c, -1)))
+P, Q = Q, R̃
+
+R = SemiclassicalJacobi(P.t, mod(P.a, -1), mod(P.b, -1), mod(P.c, -1))
+R̃ = SemiclassicalOrthogonalPolynomials.toclassical(R)
+
+P \ R
+
+
+A, b = M.A, M.B
+
+
+Ba = reverse(B.args)
+LazyA
+LazyArrays._mul_colsupport(1, reverse(B.args)...)
+
+ext = Base.get_extension(LazyArrays, :LazyArraysBandedMatricesExt)
+Base.copy(L::Ldiv{<:ext.BandedLazyLayouts,<:LazyArrays.AbstractLazyLayout}) = LazyArrays.lazymaterialize(\, L.A, L.B)
+
+using ArrayLayouts
+
+@inline LazyArrays.simplifiable(M::LazyArrays.Mul{<:LazyArrays.ApplyLayout{typeof(\)},<:LazyArrays.DiagonalLayout{<:LazyArrays.AbstractFillLayout}}) = Val(false)
+@inline LazyArrays.simplifiable(M::LazyArrays.Mul{LazyArrays.ApplyLayout{typeof(\)},<:LazyArrays.DiagonalLayout{<:LazyArrays.AbstractFillLayout}}) = Val(false) # need both this and the above
+@inline LazyArrays.simplifiable(M::LazyArrays.Mul{<:LazyArrays.AbstractInvLayout,<:LazyArrays.DiagonalLayout{<:LazyArrays.AbstractFillLayout}}) = Val(false)
+@inline function Base.copy(M::Mul{LazyArrays.ApplyLayout{typeof(\)},StyleB} where {StyleB<:(ArrayLayouts.DiagonalLayout{<:ArrayLayouts.AbstractFillLayout})})
+    A, B = LazyArrays.arguments(\, M.A)
+    A \ (B * M.B)
+end
+@inline Base.copy(M::Mul{LazyArrays.ApplyLayout{typeof(\)},<:ext.BandedLazyLayouts}) = LazyArrays.simplify(M)
+
+
+SemiclassicalOrthogonalPolynomials.arguments(F)
+
+
+
+toclassical(SemiclassicalJacobi(Q.t, mod(Q.a, -1), mod(Q.b, -1), mod(Q.c, -1)))
+
 
 A = SemiclassicalJacobi(2.0, -0.5, -0.0, -0.5)
 B = SemiclassicalJacobi(2.0, 0.5, 0.0, 0.5)
 A \ B
-Q, P = A, B 
-Qt, Qa, Qb, Qc = Q.t, Q.a, Q.b, Q.c 
-Δa = Qa-P.a
-Δb = Qb-P.b
-Δc = Qc-P.c
+Q, P = A, B
+Qt, Qa, Qb, Qc = Q.t, Q.a, Q.b, Q.c
+Δa = Qa - P.a
+Δb = Qb - P.b
+Δc = Qc - P.c
 
 @testset "Differentiation" begin
     t, a, b, c = 2.0, 1.0, -1.0, 1.0
     Rᵦₐ₁ᵪᵗᵃ⁰ᶜ = Weighted(SemiclassicalJacobi(t, a, 0.0, c)) \ Weighted(SemiclassicalJacobi(t, a, 1.0, c))
     Dₐ₀ᵪᵃ⁺¹¹ᶜ⁺¹ = diff(SemiclassicalJacobi(t, a, 0.0, c))
-    Rₐ₊₁₁ᵪ₊₁ᵗᵃ⁺¹⁰ᶜ⁺¹ = ApplyArray(inv, SemiclassicalJacobi(t, a+1, 1.0, c+1) \ SemiclassicalJacobi(t, a+1, 0.0, c + 1))
+    Rₐ₊₁₁ᵪ₊₁ᵗᵃ⁺¹⁰ᶜ⁺¹ = ApplyArray(inv, SemiclassicalJacobi(t, a + 1, 1.0, c + 1) \ SemiclassicalJacobi(t, a + 1, 0.0, c + 1))
     Dₐ₋₁ᵪᵃ⁺¹⁰ᶜ⁺¹ = Rₐ₊₁₁ᵪ₊₁ᵗᵃ⁺¹⁰ᶜ⁺¹ * Dₐ₀ᵪᵃ⁺¹¹ᶜ⁺¹.args[2] * Rᵦₐ₁ᵪᵗᵃ⁰ᶜ
     b2 = Vcat(0.0, 0.0, Dₐ₋₁ᵪᵃ⁺¹⁰ᶜ⁺¹[band(1)])
     b1 = Vcat(0.0, Dₐ₋₁ᵪᵃ⁺¹⁰ᶜ⁺¹[band(0)])
     data = Hcat(b2, b1)'
     D = _BandedMatrix(data, ∞, -1, 2)
-    @test Hcat(Zeros(∞), Dₐ₋₁ᵪᵃ⁺¹⁰ᶜ⁺¹)[1:100, 1:100] ≈ D[1:100, 1:100] 
+    @test Hcat(Zeros(∞), Dₐ₋₁ᵪᵃ⁺¹⁰ᶜ⁺¹)[1:100, 1:100] ≈ D[1:100, 1:100]
     P = SemiclassicalJacobi(t, a, b, c)
     DP = diff(P)
     @test DP.args[2][1:100, 1:100] ≈ D[1:100, 1:100]
@@ -104,4 +163,4 @@ Qt, Qa, Qb, Qc = Q.t, Q.a, Q.b, Q.c
             @test df[x] ≈ dg[x]
         end
     end
-end
+end