Supporting a static Taylor1 or TaylorN type?

[mewilhel]

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")

[lbenet]

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.

[dpsanders]

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.

[mewilhel]

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?

[lbenet]

For me it's ok if you would like to have it here.