Added mapping of points to deal with general Rays (#12)

* Added mapping of points to deal with general Rays

* Added unit test for general rays

* Working on scaled rays.

* Incorporated SRay into Ray

Author: tomtrogdon

.gitignore

@@ -5,3 +5,4 @@ deps/deps.jl
 test/.DS_Store
 .DS_Store
 Manifest.toml
+dev_code.jl

src/Domains/Domains.jl

@@ -13,11 +13,11 @@ isless(d2::Ray{true,T2},d1::IntervalOrSegment{T1}) where {T1<:Real,T2<:Real} = d
 Base.setdiff(d::Union{AbstractInterval,Segment,Ray,Line}, ptsin::UnionDomain{AS}) where {AS <: AbstractVector{P}} where {P <: Point} =
     affine_setdiff(d, ptsin)
 
-Base.setdiff(d::Union{AbstractInterval,Segment,Ray,Line}, ptsin::WrappedDomain{<:AbstractVector}) = 
+Base.setdiff(d::Union{AbstractInterval,Segment,Ray,Line}, ptsin::WrappedDomain{<:AbstractVector}) =
     affine_setdiff(d, ptsin)
 
-Base.setdiff(d::Union{AbstractInterval,Segment,Ray,Line}, ptsin::AbstractVector{<:Number}) = 
+Base.setdiff(d::Union{AbstractInterval,Segment,Ray,Line}, ptsin::AbstractVector{<:Number}) =
     ApproxFunBase._affine_setdiff(d, ptsin)
 
-Base.setdiff(d::Union{AbstractInterval,Segment,Ray,Line}, ptsin::Number) = 
-    ApproxFunBase._affine_setdiff(d, ptsin)    
+Base.setdiff(d::Union{AbstractInterval,Segment,Ray,Line}, ptsin::Number) =
+    ApproxFunBase._affine_setdiff(d, ptsin)

src/Domains/Ray.jl

@@ -1,40 +1,33 @@
-
-
 export Ray
 
-
-
-## Standard interval
-
-# angle is (false==0) and π (true==1)
-# or ranges from (-1,1].  We use 1 as 1==true.
-#orientation true means oriented out
 """
-    Ray{a}(c,o)
+    Ray{a}(c,L,o)
 
-represents a ray at angle `a` starting at `c`, with orientation out to
+represents a scaled ray (with scale factor L) at angle `a` starting at `c`, with orientation out to
 infinity (`o = true`) or back from infinity (`o = false`).
 """
 struct Ray{angle,T<:Number} <: SegmentDomain{T}
     center::T
+    L::T
     orientation::Bool
-    Ray{angle,T}(c,o) where {angle,T} = new{angle,T}(c,o)
-    Ray{angle,T}(c) where {angle,T} = new{angle,T}(c,true)
-    Ray{angle,T}() where {angle,T} = new{angle,T}(zero(T),true)
+    Ray{angle,T}(c,L,o) where {angle,T} = new{angle,T}(c,L,o)
+    Ray{angle,T}(c,o) where {angle,T} = new{angle,T}(c,one(T),o)
+    Ray{angle,T}(c) where {angle,T} = new{angle,T}(c,one(T),true)
+    Ray{angle,T}() where {angle,T} = new{angle,T}(zero(T),one(T),true)
     Ray{angle,T}(r::Ray{angle,T}) where {angle,T} = r
 end
 
-const RealRay{T} = Union{Ray{false,T},Ray{true,T}}
-
-Ray{a}(c,o) where {a} = Ray{a,typeof(c)}(c,o)
+Ray{a}(c,L,o) where {a} = Ray{a,typeof(c)}(c,L,o)
+Ray{a}(c,o) where {a} = Ray{a,typeof(c)}(c,one(typeof(c)),o)
 Ray{a}(c::Number) where {a} = Ray{a,typeof(c)}(c)
 Ray{a}() where {a} = Ray{a,Float64}()
 
 angle(d::Ray{a}) where {a} = a*π
 
 # ensure the angle is always in (-1,1]
-Ray(c,a,o) = Ray{a==0 ? false : (abs(a)≈(1.0π) ? true : mod(a/π-1,-2)+1),typeof(c)}(c,o)
-Ray(c,a) = Ray(c,a,true)
+Ray(c,a,L,o) = Ray{a==0 ? false : (abs(a)≈(1.0π) ? true : mod(a/π-1,-2)+1),typeof(c)}(c,L,o)
+Ray(c,a,o) = Ray{a==0 ? false : (abs(a)≈(1.0π) ? true : mod(a/π-1,-2)+1),typeof(c)}(c,one(typeof(c)),o)
+Ray(c,a) = Ray(c,a,one(typeof(c)),true)
 
 Ray() = Ray{false}()
 
@@ -47,9 +40,9 @@ function convert(::Type{Ray}, d::AbstractInterval)
     @assert abs(a)==Inf || abs(b)==Inf
 
     if abs(b)==Inf
-        Ray(a,angle(b),true)
+        Ray(a,angle(b),one(typeof(a)),true)
     else #abs(a)==Inf
-        Ray(b,angle(a),false)
+        Ray(b,angle(a),one(typeof(a)),false)
     end
 end
 Ray(d::AbstractInterval) = convert(Ray, d)
@@ -66,7 +59,7 @@ isempty(::Ray) = false
 
 function mobiuspars(d::Ray)
     s=(d.orientation ? 1 : -1)
-    α=conj(cisangle(d))
+    α=conj(cisangle(d))/d.L
     c=d.center
     s*α,-s*(1+α*c),α,1-α*c
 end
@@ -84,40 +77,37 @@ ray_invfromcanonicalD(x) = (x-1)^2/2
 
 
 for op in (:ray_tocanonical,:ray_tocanonicalD)
-    @eval $op(o,x)=(o ? 1 : -1)*$op(x)
+    @eval $op(L,o,x)= (o ? 1 : -1)*$op(x/L)
 end
-ray_fromcanonical(o,x)=ray_fromcanonical((o ? 1 : -1)*x)
-ray_fromcanonicalD(o,x)=(o ? 1 : -1)*ray_fromcanonicalD((o ? 1 : -1)*x)
-ray_invfromcanonicalD(o,x)=(o ? 1 : -1)*ray_invfromcanonicalD((o ? 1 : -1)*x)
+
+ray_fromcanonical(L,o,x)=L*ray_fromcanonical((o ? 1 : -1)*x)
+ray_fromcanonicalD(L,o,x)=L*(o ? 1 : -1)*ray_fromcanonicalD((o ? 1 : -1)*x)
+ray_invfromcanonicalD(L,o,x)=L*(o ? 1 : -1)*ray_invfromcanonicalD((o ? 1 : -1)*x)
 
 cisangle(::Ray{a}) where {a}=cis(a*π)
 cisangle(::Ray{false})=1
 cisangle(::Ray{true})=-1
 
 tocanonical(d::Ray,x) =
-    ray_tocanonical(d.orientation,conj(cisangle(d)).*(x-d.center))
+    ray_tocanonical(d.L,d.orientation,conj(cisangle(d)).*(x-d.center))
 tocanonicalD(d::Ray,x) =
-    conj(cisangle(d)).*ray_tocanonicalD(d.orientation,conj(cisangle(d)).*(x-d.center))
-fromcanonical(d::Ray,x) = cisangle(d)*ray_fromcanonical(d.orientation,x)+d.center
-fromcanonical(d::Ray{false},x) = ray_fromcanonical(d.orientation,x)+d.center
-fromcanonical(d::Ray{true},x) = -ray_fromcanonical(d.orientation,x)+d.center
-fromcanonicalD(d::Ray,x) = cisangle(d)*ray_fromcanonicalD(d.orientation,x)
-invfromcanonicalD(d::Ray,x) = conj(cisangle(d))*ray_invfromcanonicalD(d.orientation,x)
-
-
-
-
+    conj(cisangle(d)).*ray_tocanonicalD(d.L,d.orientation,conj(cisangle(d)).*(x-d.center))
+fromcanonical(d::Ray,x) = cisangle(d)*ray_fromcanonical(d.L,d.orientation,x)+d.center
+fromcanonical(d::Ray{false},x) = ray_fromcanonical(d.L,d.orientation,x)+d.center
+fromcanonical(d::Ray{true},x) = -ray_fromcanonical(d.L,d.orientation,x)+d.center
+fromcanonicalD(d::Ray,x) = cisangle(d)*ray_fromcanonicalD(d.L,d.orientation,x)
+invfromcanonicalD(d::Ray,x) = conj(cisangle(d))*ray_invfromcanonicalD(d.L,d.orientation,x)
 arclength(d::Ray) = Inf
 
-==(d::Ray{a},m::Ray{a}) where {a} = d.center == m.center
+==(d::Ray{a},m::Ray{a}) where {a} = d.center == m.center && d.L == m.L
 
 
 mappoint(a::Ray{false}, b::Ray{false}, x::Number) =
-    x - a.center + b.center
+    b.center + b.L/a.L*(x - a.center)
 
 
 function mappoint(a::Ray, b::Ray, x::Number)
     d = x - a.center;
-    k = d * exp((angle(b)-angle(d))*im)
+    k = b.L/a.L * d * exp((angle(b)-angle(d))*im)
     k + b.center
 end

src/Spaces/Laguerre/Laguerre.jl

@@ -14,6 +14,7 @@ export Laguerre, NormalizedLaguerre, LaguerreWeight, WeightedLaguerre
 # p_{n+1} = (A_n x + B_n)p_n - C_n p_{n-1}
 #####
 
+
 """
 `Laguerre(α)` is a space spanned by generalized Laguerre polynomials `Lₙᵅ(x)` 's
 on `(0, Inf)`, which satisfy the differential equations
@@ -33,7 +34,7 @@ Laguerre() = Laguerre(0)
 NormalizedLaguerre(α) = NormalizedPolynomialSpace(Laguerre(α))
 NormalizedLaguerre() = NormalizedLaguerre(0)
 
-spacescompatible(A::Laguerre,B::Laguerre) = A.α ≈ B.α
+spacescompatible(A::Laguerre,B::Laguerre) = A.α ≈ B.α && B.domain == A.domain
 
 canonicaldomain(::Laguerre) = Ray()
 domain(d::Laguerre) = d.domain
@@ -101,31 +102,62 @@ laguerrel(n::AbstractRange,S::Laguerre,v) = laguerrel(n,S.a,S.b,v)
 laguerrel(n,S::Laguerre,v) = laguerrel(n,S.a,S.b,v)
 
 
+#struct LaguerreTransformPlan{T,TT}
+#    space::Laguerre{TT}
+#    points::Vector{T}
+#    weights::Vector{T}
+#end
+
 struct LaguerreTransformPlan{T,TT}
     space::Laguerre{TT}
     points::Vector{T}
-    weights::Vector{T}
+    transform::Array{T}
+    c1::Vector{T}
+    c2::Vector{T}
 end
 
-plan_transform(S::Laguerre,v::AbstractVector) = LaguerreTransformPlan(S,gausslaguerre(length(v),1.0S.α)...)
-function *(plan::LaguerreTransformPlan,vals)
-#    @assert S==plan.space
-    x,w = plan.points, plan.weights
-    V=laguerrel(0:length(vals)-1,plan.space.α,x)'
-    #w2=w.*x.^(S.α-plan.space.α)   # need to weight if plan is different
-    w2=w
-    nrm=(V.^2)*w2
-    V*(w2.*vals)./nrm
+function Lag_conv(n::Int64,α::Float64)
+    out = ones(Float64,n)
+    out[1] = 1/gamma(α+1)
+    for i in 2:n
+        out[i] = out[i-1]*(i-1.0)/(i-1.0+α)
+    end
+    return sqrt.(out)
+end
+
+function lag_transform(n::Int64,α::Float64)
+    egn = eigen(SymTridiagonal([(2i + α + 1.) for i = 0:n-1],[-sqrt(i*(i+α)) for i = 1:n-1]))
+    U = egn.vectors
+    ns = sqrt(gamma(α+1))*U[1,:]
+    egn.values, U, Lag_conv(n,α), ns
 end
 
+# The old transform plan.  Is this faster for n < 150?
+#plan_transform(S::Laguerre,v::AbstractVector) = LaguerreTransformPlan(S,gausslaguerre(length(v),1.0S.α)...)
+#function *(plan::LaguerreTransformPlan,vals)
+#    @assert S==plan.space
+#    x,w = plan.points, plan.weights
+#    V=laguerrel(0:length(vals)-1,plan.space.α,x)'
+#    #w2=w.*x.^(S.α-plan.space.α)   # need to weight if plan is different
+#    w2=w
+#    nrm=(V.^2)*w2
+#    V*(w2.*vals)./nrm
+#end
+
+plan_transform(S::Laguerre,v::AbstractVector) = LaguerreTransformPlan(S,lag_transform(length(v),1.0S.α)...)
+function *(plan::LaguerreTransformPlan,vals)
+    x,w,c1,c2 = plan.points, plan.transform, plan.c1, plan.c2
+    c1.*(w*(c2.*vals))
+end
 
-points(L::Laguerre,n) = gausslaguerre(n,1.0L.α)[1]
 
+points(L::Laguerre,n) = map(x -> mappoint(Ray(),L.domain,x), gausslaguerre(n,1.0L.α)[1])
 
 
+# TODO: Separate off a real case for the real axis.
 Derivative(L::Laguerre,k) =
-    DerivativeWrapper(SpaceOperator(ToeplitzOperator(Float64[],[zeros(k);(-1.)^k]),
-                                    L,Laguerre(L.α+k)))
+    DerivativeWrapper(SpaceOperator(ToeplitzOperator(Complex{Float64}[],[zeros(k);(-1.)^k]*conj(cisangle(L.domain))/L.domain.L),
+                                    L,Laguerre(L.α+k,L.domain)))
 
 
 union_rule(A::Laguerre,B::Laguerre) = Laguerre(min(A.α,B.α))
@@ -164,15 +196,12 @@ struct LaguerreWeight{S,T} <: WeightSpace{S,Ray{false,Float64},Float64}
     L::T
     space::S
 end
-
-
 """
     LaguerreWeight(α, space)
 
 weights `space` by `x^α * exp(-x)`.
 """
 LaguerreWeight(α, space::Space) = LaguerreWeight(α, one(α),space)
-
 """
     WeightedLaguerre(α)
 
@@ -214,8 +243,9 @@ last(f::Fun{<:LaguerreWeight,T}) where T = zero(T)
 function Derivative(sp::LaguerreWeight,k)
     @assert sp.α == 0
     if k==1
+        c = conj(cisangle(sp.space.domain))/sp.space.domain.L
         D=Derivative(sp.space)
-        D2=D-sp.L*I
+        D2=c*D-(c*sp.L)*I
         DerivativeWrapper(SpaceOperator(D2,sp,LaguerreWeight(sp.α,sp.L,rangespace(D2))),1)
     else
         D=Derivative(sp)
@@ -300,6 +330,3 @@ function Multiplication(f::Fun{LaguerreWeight{H,T}},S::Laguerre) where {H<:Lague
     rs=rangespace(M)
     MultiplicationWrapper(f,SpaceOperator(M,S,LaguerreWeight(space(f).α, space(f).L, rs)))
 end
-
-
-

test/LaguerreTest.jl

@@ -1,7 +1,42 @@
 using ApproxFunOrthogonalPolynomials, ApproxFunBase, SpecialFunctions, Test
     import ApproxFunBase: testbandedoperator
 
-@testset "Laguerre and WeightedLagueree" begin
+
+@testset "Laguerre and WeightedLaguerre" begin
+
+    @testset "General scaled rays" begin
+        r = Ray(-1.0,0.0,2.0,true)
+        L = Laguerre(1.0,r)
+        f = x -> exp.(-x)
+        F = Fun(f, L)
+        @test F(.3) ≈ -F'(.3)
+
+        r = Ray(1.0,0.0,2.0,false)
+        L = Laguerre(1.0,r)
+        f = x -> exp.(-x)
+        F = Fun(f, L, 100)
+        @test F(1.3) ≈ -F'(1.3)
+
+        r = Ray(1.0,π,2.0,true)
+        L = Laguerre(1.0,r)
+        f = x -> exp.(x)
+        F = Fun(f, L, 100)
+        @test F(.3) ≈ F'(.3)
+
+        r = Ray(-3.0,0.0,2.0,true)
+        L = Laguerre(1.0,r)
+        f = x -> exp.(-x)
+        F = Fun(f, L)  # overflow/underflow issues beyond 190ish
+        @test abs(F(0.0) - f(0.)) < 1e-5
+        @test abs(F(-2.0) - f(-2.)) < 1e-5
+
+        r = Ray(3.0,π,.1,false)
+        L = Laguerre(1.0,r)
+        f = x -> exp.(-x^2)
+        F = Fun(f, L)  # overflow/underflow issues beyond 190ish
+        @test abs(F(-1.0) - f(-1.0)) < 1e-5
+    end
+
     @testset "Evaluation" begin
         f=Fun(Laguerre(0.), [1,2,3])
         @test f(0.1) ≈ 5.215