From a96e34ba55c6a58674453f27858e5410cd448b85 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?El=C3=ADas=20Snorrason?= Date: Mon, 20 Aug 2018 15:19:46 +0000 Subject: [PATCH 1/3] Implementation for vector valued Weeks{T}.func --- src/weeks.jl | 38 ++++++++++++++++++++++++++++++-------- 1 file changed, 30 insertions(+), 8 deletions(-) diff --git a/src/weeks.jl b/src/weeks.jl index 3aae2c6..95f29b5 100644 --- a/src/weeks.jl +++ b/src/weeks.jl @@ -15,12 +15,12 @@ mutable struct Weeks{T} <: AbstractWeeks Nterms::Int sigma::Float64 b::Float64 - coefficients::Array{T,1} + coefficients::Array{T,2} end function _get_coefficients(func, Nterms, sigma, b, ::Type{T}) where T <:Number a0 = _wcoeff(func, Nterms, sigma, b, T) - return a0[Nterms+1:2*Nterms] + return a0[Nterms+1:2*Nterms,:] end _get_coefficients(w::Weeks{T}) where T <: Number = _get_coefficients(w.func, w.Nterms, w.sigma, w.b, T) @@ -229,10 +229,16 @@ function _wcoeff(F, N, sig, b) s = sig .+ imaginary_unit * y FF0 = map(F, s) FF = [FF1 * (b + imaginary_unit * y1) for (FF1,y1) in zip(FF0,y)] - a = (FF |> AbstractFFTs.fftshift |> FFTW.fft |> AbstractFFTs.fftshift) / (2 * N) - return exp.(Complex(zero(h),-one(h)) * n*h/2) .* a + FFTranspose = Array{Complex{Float64},2}(undef,(length(y),2)) + transpose!(FFTranspose,reshape(vcat(FF...),(2,length(y)))) + FFp = FFTW.plan_fft!(FFTranspose,1,flags=FFTW.MEASURE) + a = colFFTwshift(FFp,FFTranspose) / (2 * N) + return exp.(Complex(zero(h),-one(h)) * n*h/2) .* a end + + + function _laguerre(a::AbstractVector,x::AbstractArray) N = length(a) - 1 unp1 = zero(x) @@ -247,13 +253,15 @@ function _laguerre(a::AbstractVector,x::AbstractArray) return unm1 end -function _laguerre(a::AbstractVector,x) - n = length(a) - 1 + + +function _laguerre(a::AbstractArray,x) + n = size(a)[1] - 1 unp1 = zero(x) - un = a[n+1]*one(x) + un = a[n+1,:] .* one(x) local unm1 for n in n:-1:1 - unm1 = (1//n)*(2*n-1-x) * un - n/(n+1)*unp1 + a[n] + unm1 = (1//n)*(2*n-1-x) .* un .- n/(n+1)*unp1 + a[n,:] unp1 = un un = unm1 end @@ -283,3 +291,17 @@ function _werr2t(b, F, N, sig) sa2 = sum(abs.(@view a[3*N+1:4*N])) return log(sa2) end + + +function colFFTwshift(plan::N, F::Array{T,2}) where {T<:Number,N<:FFTW.cFFTWPlan} + + # Performs a columnwise FFT along with two shifts when computing the Laguerre coefficients + # of a vector valued function "weeks.func". + # + # plan is a predetermined transform plan of the columns of F. + # + # Samples of each vector element of "w.func" are stored in the columns of the + # array F. + + return AbstractFFTs.fftshift(plan*AbstractFFTs.fftshift(F)) +end From 985a90e4014696b651d47a637b55b88e83906a67 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?El=C3=ADas=20Snorrason?= Date: Mon, 20 Aug 2018 15:55:47 +0000 Subject: [PATCH 2/3] Arbitrary length of vector Weeks.func --- src/weeks.jl | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/src/weeks.jl b/src/weeks.jl index 95f29b5..1b48073 100644 --- a/src/weeks.jl +++ b/src/weeks.jl @@ -1,6 +1,7 @@ # import Optim # broken import AbstractFFTs import FFTW +using LinearAlgebra: transpose! abstract type AbstractWeeks <: AbstractILt end @@ -221,6 +222,7 @@ _wcoeff(F, N, sig, b, ::Type{T}) where T <: Real = real(_wcoeff(F, N, sig, b)) _wcoeff(F, N, sig, b, ::Type{T}) where T <: Complex = _wcoeff(F, N, sig, b) function _wcoeff(F, N, sig, b) + FN = length(F(complex(1.01))) n = -N:N-1 # FIXME: remove 1 and test h = pi / N # FIXME: what data type ? th = h .* (n .+ 1//2) @@ -229,8 +231,8 @@ function _wcoeff(F, N, sig, b) s = sig .+ imaginary_unit * y FF0 = map(F, s) FF = [FF1 * (b + imaginary_unit * y1) for (FF1,y1) in zip(FF0,y)] - FFTranspose = Array{Complex{Float64},2}(undef,(length(y),2)) - transpose!(FFTranspose,reshape(vcat(FF...),(2,length(y)))) + FFTranspose = Array{Complex{Float64},2}(undef,(length(y),FN)) + transpose!(FFTranspose,reshape(vcat(FF...),(FN,length(y)))) FFp = FFTW.plan_fft!(FFTranspose,1,flags=FFTW.MEASURE) a = colFFTwshift(FFp,FFTranspose) / (2 * N) return exp.(Complex(zero(h),-one(h)) * n*h/2) .* a From eec0dabb2db739400a2c33fd12466988ae408c89 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?El=C3=ADas=20Snorrason?= Date: Mon, 20 Aug 2018 17:25:59 +0000 Subject: [PATCH 3/3] Test vector valued functions --- test/weeks_test.jl | 18 ++++++++++++++++++ 1 file changed, 18 insertions(+) diff --git a/test/weeks_test.jl b/test/weeks_test.jl index e0f7032..add4812 100644 --- a/test/weeks_test.jl +++ b/test/weeks_test.jl @@ -59,3 +59,21 @@ end let Fc = Weeks(Fcomplex, 1024, datatype=Complex), trange = range(0.0, stop=30.0, length=1000) @test isapprox(Fc.(trange), fcomplex.(trange), atol=0.001) end + + +### Complex Vector +Nr = collect(1:64) +function fcomplex(t) + α = complex(-0.3,6.0) + return Nr .*exp(α*t) +end + +function Fcomplex(s::Complex{Float64})::Vector{Complex{Float64}} + # Laplace domain + α = complex(-0.3,6.0) + return Nr ./(s-α) +end + +let Ft = Weeks(Fcomplex,512,10.0,1.0,datatype=Complex) + @test isapprox(Ft(0.1), fcomplex(0.1),atol=1.0E-8) +end