TaylorSeries.jl
TaylorSeries.jl copied to clipboard
Supporting a static Taylor1 or TaylorN type?
I was just wondering if anyone has thoughts around possibly supporting statically typed versions of Taylor series in this package. I've thrown together a sample implementation and it seems like there is a reasonable speed up to be had.
addition speedup: 63.0x
subtraction speedup: 63.0x
scalar multiplication speedup: 14.0x
multiplication speed up: 10.0x
Sample Snippet:
import Base: +, -, *
struct STaylor1{N,T}
coeffs::SVector{N,T}
end
+(x::STaylor1{N,T}, y::STaylor1{N,T}) where {T<:Number,N} = STaylor1{N,T}(x.coeffs + y.coeffs)
-(x::STaylor1{N,T}, y::STaylor1{N,T}) where {T<:Number,N} = STaylor1{N,T}(x.coeffs - y.coeffs)
*(x::STaylor1{N,T}, y::T) where {T<:Number,N} = STaylor1{N,T}(y*x.coeffs)
*(y::T, x::STaylor1{N,T}) where {T<:Number,N} = STaylor1{N,T}(y*x.coeffs)
@generated function *(x::STaylor1{N,T}, y::STaylor1{N,T}) where {T<:Number,N}
ex_calc = quote end
for j = 1:N
ex_line = :(x.coeffs[1]*y.coeffs[$j])
for k = 2:j
ex_line = :($ex_line + x.coeffs[$k]*y.coeffs[$(j-k+1)])
end
sym = Symbol("c$j")
ex_line = :($sym = $ex_line)
push!(ex_calc.args, ex_line)
end
exout = :((c1,))
for i = 2:N
sym = Symbol("c$i")
push!(exout.args, sym)
end
return quote
$ex_calc
tup = $exout
svec = SVector{N,T}(tup)
return STaylor1{N,T}(svec)
end
end
Benchmarking script:
using StaticArrays, BenchmarkTools, TaylorSeries
order = 5
a = 1.1
s1 = @SVector rand(order)
s2 = @SVector rand(order)
st1 = STaylor1{order,Float64}(s1)
st2 = STaylor1{order,Float64}(s2)
out1 = Taylor1(Float64[s1[i] for i in 1:length(s1)])
t1 = Taylor1(Float64[s1[i] for i in 1:length(s1)])
t2 = Taylor1(Float64[s2[i] for i in 1:length(s2)])
function f(g, x, y, n)
t = x
for i=1:n
t = g(t,x)
t = g(t,y)
end
t
end
n = 500
@btime f($+, $st1, $st2, $n)
@btime f($+, $t1, $t2, $n)
static_time = @belapsed f($+, $st1, $st2, $n)
var_time = @belapsed f($+, $t1, $t2, $n)
println("addition speedup: $(floor(var_time/static_time))x")
@btime f($-, $st1, $st2, $n)
@btime f($-, $t1, $t2, $n)
static_time = @belapsed f($-, $st1, $st2, $n)
var_time = @belapsed f($-, $t1, $t2, $n)
println("subtraction speedup: $(floor(var_time/static_time))x")
@btime f($*, $a, $st2, $n)
@btime f($*, $a, $t2, $n)
static_time = @belapsed f($*, $a, $st2, $n)
var_time = @belapsed f($*, $a, $t2, $n)
println("scalar multiplication speedup: $(floor(var_time/static_time))x")
@btime f($*, $st1, $st2, $n)
@btime f($*, $t1, $t2, $n)
static_time = @belapsed f($*, $st1, $st2, $n)
var_time = @belapsed f($*, $t1, $t2, $n)
println("multiplication speed up: $(floor(var_time/static_time))x")
That would be truly fantastic! I have been thinking about this lately, since I learned recently about Setfield.jl, which may be handy. Would you like to open a PR so we can check what would be running and what not? I guess the speed up would also go to TaylorIntegration.jl and TaylorModels.jl, so I am really interested on this.
I have some (unreleased) work along these lines that basically does exactly what you're doing: https://github.com/dpsanders/StaticTaylorSeries.jl
We should definitely put together something there.
Excellent! I'm definitely interested in pursuing this. Would you all prefer the work to be done on a fork of TaylorSeries.jl or shall I start playing with the StaticTaylorSeries.jl repository?
For me it's ok if you would like to have it here.