Skip to content

Commit

Permalink
Rank of laguerre coefficients in constructor, fixes JuliaMath#4
Browse files Browse the repository at this point in the history 
The coefficient-rank in the Weeks{T}-type is no longer declared
explicitly. Instead, it is supplied as an argument to the outer
constructor.

Array- and complex-tests have been reformatted and grouped into
testsets.
  • Loading branch information
@elisno
elisno committed Dec 13, 2018
1 parent 005feac commit bd27190
Show file tree
Hide file tree
Showing 3 changed files with 101 additions and 66 deletions.
22 changes: 17 additions & 5 deletions src/weeks.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -15,7 +15,7 @@ mutable struct Weeks{T} <: AbstractWeeks
Nterms::Int
sigma::Float64
b::Float64
coefficients::Array{T,1}
coefficients::Array{T}
end

function _get_coefficients(func, Nterms, sigma, b, ::Type{T}) where T <:Number
Expand All @@ -36,6 +36,7 @@ _get_array_coefficients(w::Weeks{T}) where T <: Number = _get_array_coefficients
#_set_array_coefficients(w::Weeks) = (w.coefficients = _get_coefficients(w))

const weeks_default_num_terms = 64
const weeks_default_ndims = 0

"""
w::Weeks{datatype} = Weeks(func::Function, Nterms::Integer=64, sigma=1.0, b=1.0; datatype=Float64)
Expand All @@ -59,14 +60,25 @@ julia> ft(pi/2)
```
"""
function Weeks(func::Function, Nterms::Integer=weeks_default_num_terms,
sigma=1.0, b=1.0; datatype=Float64)
sigma=1.0, b=1.0; datatype=Float64,rank::Integer=weeks_default_ndims)
outdatatype = datatype == Complex ? Complex{Float64} : datatype # allow `Complex` as abbrev for Complex{Float64}
return Weeks{outdatatype}(func, Nterms, sigma, b, _get_coefficients(func, Nterms, sigma, b, outdatatype))
if rank == 0
# Use scalar-functionality
return Weeks{outdatatype}(func, Nterms, sigma, b, _get_coefficients(func, Nterms, sigma, b, outdatatype))
else
# Use array-functionality
return Weeks{outdatatype}(func, Nterms, sigma, b, _get_array_coefficients(func, Nterms, sigma, b, outdatatype))
end
end

function eval_ilt(w::Weeks, t)
L = _laguerre(w.coefficients, 2 * w.b * t)
return L * exp((w.sigma - w.b) * t)
coeffs = w.coefficients
if ndims(coeffs)==1
L = _laguerre(coeffs, 2 * w.b * t)
else
L = reshape(mapslices(i -> _laguerre(i,2 * w.b * t),coeffs,dims=(1)),size(coeffs)[2:end])
end
return L .* exp((w.sigma - w.b) * t)
end

function eval_ilt(w::Weeks, t::AbstractArray)
Expand Down
116 changes: 68 additions & 48 deletions test/arrayweeks_test.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11,60 +11,80 @@ function Fcomplex(s)
return 1 / (s - α)
end

# Test scalar-functionality of arrcoeff, Fcomplex=(s-a)⁻¹
let arr = try InverseLaplace._arrcoeff(Fcomplex,4,1.0,1.0);
InverseLaplace._arrcoeff(Fcomplex,4,1.0,1.0)
catch
InverseLaplace._arrcoeff(Fcomplex,4,1.0,1.0)
end, scal = InverseLaplace._wcoeff(Fcomplex,4,1.0,1.0)
@testset "Array functionality" begin
@testset "Internal functions" begin
# Test scalar-functionality of arrcoeff, Fcomplex=(s-a)⁻¹
let arr = try InverseLaplace._arrcoeff(Fcomplex,4,1.0,1.0);
InverseLaplace._arrcoeff(Fcomplex,4,1.0,1.0)
catch
InverseLaplace._arrcoeff(Fcomplex,4,1.0,1.0)
end, scal = InverseLaplace._wcoeff(Fcomplex,4,1.0,1.0)

@test arr == scal
end

# Test for ordering of array valued coefficients (along first dimension)
let arr = try InverseLaplace._arrcoeff(arrFcomplex,4,1.0,1.0);
InverseLaplace._arrcoeff(arrFcomplex,4,1.0,1.0)
catch
InverseLaplace._arrcoeff(arrFcomplex,4,1.0,1.0)
end, scal = InverseLaplace._wcoeff(Fcomplex,4,1.0,1.0)
@test arr[:,1,1] == scal
@test arr[:,1,2] == 2 .* arr[:,1,1]
@test isapprox(arr[:,2,1] , 3 .* arr[:,1,1], atol=1E-15)
@test arr[:,2,2] == 4 .* arr[:,1,1]
end
@test arr == scal
end

# Compare _get_coefficients and _laguerre for array functions and scalar functions
let arr = try InverseLaplace._get_array_coefficients(arrFcomplex,4,1.0,1.0,Complex);
InverseLaplace._get_array_coefficients(arrFcomplex,4,1.0,1.0,Complex)
catch
InverseLaplace._get_array_coefficients(arrFcomplex,4,1.0,1.0,Complex)
end, scal = InverseLaplace._get_coefficients(Fcomplex,4,1.0,1.0,Complex)
@test arr[:,1,1] == scal
@test arr[:,1,2] == 2 .* scal
@test isapprox(arr[:,2,1] , 3 .* scal, atol=1E-15)
@test arr[:,2,2] == 4 .* scal

laguerreeval = reshape(mapslices(i -> InverseLaplace._laguerre(i,1.0),arr,dims=(1)),(2,2))
@test isapprox(laguerreeval, [1 2; 3 4] .* InverseLaplace._laguerre(scal,1.0), atol = 1E-15)
end
# Test for ordering of array valued coefficients (along first dimension)
let arr = try InverseLaplace._arrcoeff(arrFcomplex,4,1.0,1.0);
InverseLaplace._arrcoeff(arrFcomplex,4,1.0,1.0)
catch
InverseLaplace._arrcoeff(arrFcomplex,4,1.0,1.0)
end, scal = InverseLaplace._wcoeff(Fcomplex,4,1.0,1.0)
@test arr[:,1,1] == scal
@test arr[:,1,2] == 2 .* arr[:,1,1]
@test isapprox(arr[:,2,1] , 3 .* arr[:,1,1], atol=1E-15)
@test arr[:,2,2] == 4 .* arr[:,1,1]
end

# Compare _get_coefficients and _laguerre for array functions and scalar functions
let arr = try InverseLaplace._get_array_coefficients(arrFcomplex,4,1.0,1.0,Complex);
InverseLaplace._get_array_coefficients(arrFcomplex,4,1.0,1.0,Complex)
catch
InverseLaplace._get_array_coefficients(arrFcomplex,4,1.0,1.0,Complex)
end, scal = InverseLaplace._get_coefficients(Fcomplex,4,1.0,1.0,Complex)
@test arr[:,1,1] == scal
@test arr[:,1,2] == 2 .* scal
@test isapprox(arr[:,2,1] , 3 .* scal, atol=1E-15)
@test arr[:,2,2] == 4 .* scal

laguerreeval = reshape(mapslices(i -> InverseLaplace._laguerre(i,1.0),arr,dims=(1)),(2,2))
@test isapprox(laguerreeval, [1 2; 3 4] .* InverseLaplace._laguerre(scal,1.0), atol = 1E-15)
end
end
@testset "ILT for arrays" begin
# Test the ILT calculation for arrays and compare with the scalar Weeks method
let
# Weeks parameters
N = 4
σ = 1.0
b = 0.5
t = 1.0
testingrank = 2
funcdims = size(arrFcomplex(rand(Complex{Float64}))) # (2,2)

# Test the ILT calculation for arrays and compare with the scalar Weeks method
let
# Weeks parameters
σ = 1.0
b = 0.5
t = 1.0
# Evaluating ILT for array valued function
coef = InverseLaplace._get_array_coefficients(arrFcomplex,4,σ,b,Complex)
lag = reshape(mapslices(i -> InverseLaplace._laguerre(i,2 * b * t),coef,dims=(1)),(2,2))
inverse = lag .* exp((σ - b) * t)

# Evaluating ILT for array valued function
coef = InverseLaplace._get_array_coefficients(arrFcomplex,4,σ,b,Complex)
lag = reshape(mapslices(i -> InverseLaplace._laguerre(i,2 * b * t),coef,dims=(1)),(2,2))
inverse = lag .* exp((σ - b) * t)
# Evaluating ILT for scalar function, then multiplying by an array
Ft = Weeks(Fcomplex,4,σ,b,datatype=Complex)
Ft_array_eval = [1 2;3 4] .* Ft(1.0)

# Evaluating ILT for scalar function, then multiplying by an array
Ft = Weeks(Fcomplex,4,σ,b,datatype=Complex)
Ft_array_eval = [1 2;3 4] .* Ft(1.0)
@test isapprox(Ft_array_eval , inverse, atol=1E-15)

@test isapprox(Ft_array_eval , inverse, atol=1E-15)
# Test the array constructor for Weeks{T}, both real and complex
Arr_c = Weeks(arrFcomplex,N,σ,b,datatype=Complex,rank=testingrank)
Scal_c = Weeks(Fcomplex,N,σ,b,datatype=Complex)
Arr_r = Weeks(arrFcomplex,N,σ,b,rank=testingrank)
Scal_r = Weeks(Fcomplex,N,σ,b)
@test ndims(Arr_c.coefficients) == testingrank + 1
@test size(Arr_c.coefficients) == (Arr_c.Nterms, 2,2)
@test ndims(Arr_r.coefficients) == testingrank + 1
@test size(Arr_r.coefficients) == (Arr_r.Nterms,2,2)
@test real(Arr_c.coefficients) == Arr_r.coefficients
@test isapprox(Arr_c(t),[1 2; 3 4] .* Scal_c(t),atol=1E-15)
@test isapprox(Arr_r(t),[1 2; 3 4] .* Scal_r(t),atol=1E-15)
@test isapprox(real(Arr_c(t)),Arr_r(t),atol=1E-15)
end
end
end
29 changes: 16 additions & 13 deletions test/weeks_test.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -43,19 +43,22 @@ setparameters(fle,2.0,2.0,80)
@test c1 != fle.coefficients

### Complex
@testset "Complex-functionality" begin
function Fcomplex(s)
# Laplace domain
α = complex(-0.3, 6.0)
return 1 / (s - α)
end

function Fcomplex(s)
# Laplace domain
α = complex(-0.3, 6.0)
return 1 / (s - α)
end

function fcomplex(t)
# Time domain
α = complex(-0.3, 6.0)
return exp(α * t)
end
function fcomplex(t)
# Time domain
α = complex(-0.3, 6.0)
return exp(α * t)
end

let Fc = Weeks(Fcomplex, 256, 1.0, 10.0, datatype=Complex), trange = range(0.0, stop=15.0, length=5)
@test isapprox(Fc.(trange), fcomplex.(trange), atol=0.0001)
let Fc = Weeks(Fcomplex, 256, 1.0, 10.0, datatype=Complex), trange = range(0.0, stop=15.0, length=5)
@test isapprox(Fc.(trange), fcomplex.(trange), atol=0.0001)
Fr = Weeks(Fcomplex, 256, 1.0, 10.0)
@test real(Fc.(trange)) == Fr.(trange)
end
end

0 comments on commit bd27190

Please sign in to comment.