mapped Chebyshev U Off Hilbert (#61)

* mapped Chebyshev U Off Hilbert

* Update stieltjes.jl

Author: dlfivefifty

Project.toml

@@ -1,7 +1,7 @@
 name = "ClassicalOrthogonalPolynomials"
 uuid = "b30e2e7b-c4ee-47da-9d5f-2c5c27239acd"
 authors = ["Sheehan Olver <[email protected]>"]
-version = "0.4.4"
+version = "0.4.5"
 
 [deps]
 ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"

src/stieltjes.jl

@@ -103,6 +103,10 @@ end
 
 
 
+##
+# mapped
+###
+
 @simplify function *(H::Hilbert, w::SubQuasiArray{<:Any,1})
     T = promote_type(eltype(H), eltype(w))
     m = parentindices(w)[1]
@@ -113,14 +117,24 @@ end
     (inv.(x .- x') * P)[m]
 end
 
+
 @simplify function *(H::Hilbert, wP::Weighted{<:Any,<:SubQuasiArray{<:Any,2}})
     T = promote_type(eltype(H), eltype(wP))
     kr,jr = parentindices(wP.P)
-        # TODO: mapping other geometries
-    @assert axes(H,1) == axes(H,2) == axes(wP,1)
     P = parent(wP.P)
-    x = axes(P,1)
-    (inv.(x .- x') * Weighted(P))[kr,jr]
+    x = axes(H,1)
+    t = axes(H,2)
+    t̃ = axes(P,1)
+    if x == t
+        (inv.(t̃ .- t̃') * Weighted(P))[kr,jr]
+    else
+        M = affine(t,t̃)
+        @assert x isa Inclusion
+        a,b = first(x),last(x)
+        x̃ = Inclusion((M.A * a .+ M.b)..(M.A * b .+ M.b)) # map interval to new interval
+        Q̃,M = arguments(*, inv.(x̃ .- t̃') * Weighted(P))
+        parent(Q̃)[affine(x,axes(parent(Q̃),1)),:] * M
+    end
 end
 
 

test/test_stieltjes.jl

@@ -190,13 +190,25 @@ end
     end
 
     @testset "OffHilbert" begin
-        U = ChebyshevU()
-        W = Weighted(U)
-        t = axes(U,1)
-        x = Inclusion(2..3)
-        T̃ = chebyshevt(2..3)
-        H = T̃ \ inv.(x .- t') * W
-        @test T̃[2.3,1:100]' * H[1:100,1:100] ≈ (inv.(2.3 .- t') * W)[:,1:100]
+        @testset "ChebyshevU" begin
+            U = ChebyshevU()
+            W = Weighted(U)
+            t = axes(U,1)
+            x = Inclusion(2..3)
+            T = chebyshevt(2..3)
+            H = T \ inv.(x .- t') * W
+            @test T[2.3,1:100]' * H[1:100,1:100] ≈ (inv.(2.3 .- t') * W)[:,1:100]
+        end
+
+        @testset "mapped" begin
+            U = chebyshevu(-1..0)
+            W = Weighted(U)
+            t = axes(U,1)
+            x = Inclusion(2..3)
+            T = chebyshevt(2..3)
+            H = T \ inv.(x .- t') * W
+            @test T[2.3,1:100]' * H[1:100,1:100] ≈ (inv.(2.3 .- t') * W)[:,1:100]
+        end
     end
 end