Preliminary work on static analogs for Taylor1

benchmark/REQUIRE

@@ -0,0 +1,2 @@
+BenchmarkTools 0.3.1
+PkgBenchmark 0.1.1

benchmark/benchmarks.jl

@@ -0,0 +1,67 @@
+using BenchmarkTools, TaylorSeries
+
+const SUITE = BenchmarkGroup()
+
+S = SUITE["Constructors"] = BenchmarkGroup()
+for i in (05,10,20,50)
+    x = rand(i)
+    t = tuple(x...)
+    S["STaylor1(x::Array[len = $i])"] = @benchmarkable STaylor1($x)
+    S["STaylor1(x::Float64, Val($i))"] = @benchmarkable STaylor1($t)
+end
+
+function f(g, x, y, n)
+    t = x
+    for i=1:n
+        t = g(x,y)
+    end
+    t
+end
+function f(g, x, n)
+    t = x
+    for i=1:n
+        t = g(t)
+    end
+    t
+end
+n = 100
+dims = (5,20)
+
+S = SUITE["arithmetic"] = BenchmarkGroup()
+for i in dims
+    r = rand(i)
+    x = Taylor1(r)
+    x2 = Taylor1(r)
+    q = STaylor1(r)
+    q2 = STaylor1(r)
+    y = rand()
+    S["Taylor1{$i,Float64} + Float64"] = @benchmarkable f(+, $x, $y, $n)
+    S["STaylor1{$i,Float64} + Float64"] = @benchmarkable f(+, $q, $y, $n)
+    S["Taylor1{$i,Float64} - Float64"] = @benchmarkable f(-, $x, $y, $n)
+    S["STaylor1{$i,Float64} - Float64"] = @benchmarkable f(-, $q, $y, $n)
+    S["Taylor1{$i,Float64} / Float64"] = @benchmarkable f(/, $x, $y, $n)
+    S["STaylor1{$i,Float64} / Float64"] = @benchmarkable f(/, $q, $y, $n)
+    S["Taylor1{$i,Float64} * Float64"] = @benchmarkable f(*, $x, $y, $n)
+    S["STaylor1{$i,Float64} * Float64"] = @benchmarkable f(*, $q, $y, $n)
+    S["Taylor1{$i,Float64} + Taylor1{$i,Float64}"] = @benchmarkable f(+, $x, $x2, $n)
+    S["STaylor1{$i,Float64} + STaylor1{$i,Float64}"] = @benchmarkable f(+, $q, $q2, $n)
+    S["Taylor1{$i,Float64} - Taylor1{$i,Float64}"] = @benchmarkable f(-, $x, $x2, $n)
+    S["STaylor1{$i,Float64} - STaylor1{$i,Float64}"] = @benchmarkable f(-, $q, $q2, $n)
+    S["Taylor1{$i,Float64} * Taylor1{$i,Float64}"] = @benchmarkable f(*, $x, $x2, $n)
+    S["STaylor1{$i,Float64} * STaylor1{$i,Float64}"] = @benchmarkable f(*, $q, $q2, $n)
+    S["Taylor1{$i,Float64} / Taylor1{$i,Float64}"] = @benchmarkable f(/, $x, $x2, $n)
+    S["STaylor1{$i,Float64} / STaylor1{$i,Float64}"] = @benchmarkable f(/, $q, $q2, $n)
+end
+
+S = SUITE["functions"] = BenchmarkGroup()
+for i in dims
+    r = rand(i)
+    x = Taylor1(r)
+    q = STaylor1(r)
+    y = rand()
+    for g in (exp, abs, zero, one, real, imag, conj, adjoint, iszero, isnan,
+              isinf, deg2rad, rad2deg)
+        S["$(Symbol(g))(Taylor1{Float64}), len = $i"] = @benchmarkable f($g, $x, $n)
+        S["$(Symbol(g))(STaylor1{$i,Float64}),"] = @benchmarkable f($g, $q, $n)
+    end
+end

benchmark/run_benchmark.jl

@@ -0,0 +1,12 @@
+using PkgBenchmark
+results = benchmarkpkg("TaylorSeries")
+show(results)
+
+#=
+# specify tag and uncommit to benchmark versus prior tagged version
+tag =
+results = judge("TaylorSeries", tag)
+show(results)
+=#
+
+export_markdown("results.md", results)

benchmark/tune.json

@@ -0,0 +1 @@
+[{"Julia":"1.3.1","BenchmarkTools":"0.4.3"},[["BenchmarkGroup",{"data":{"Constructors":["BenchmarkGroup",{"data":{"STaylor1(x::Float64, Val(5))":["BenchmarkTools.Parameters",{"gctrial":true,"time_tolerance":0.05,"samples":10000,"evals":12,"gcsample":false,"seconds":5.0,"overhead":0.0,"memory_tolerance":0.01}],"STaylor1(x::Array[length = 10])":["BenchmarkTools.Parameters",{"gctrial":true,"time_tolerance":0.05,"samples":10000,"evals":15,"gcsample":false,"seconds":5.0,"overhead":0.0,"memory_tolerance":0.01}],"STaylor1(x::Float64, Val(10))":["BenchmarkTools.Parameters",{"gctrial":true,"time_tolerance":0.05,"samples":10000,"evals":24,"gcsample":false,"seconds":5.0,"overhead":0.0,"memory_tolerance":0.01}],"STaylor1(x::Float64, Val(50))":["BenchmarkTools.Parameters",{"gctrial":true,"time_tolerance":0.05,"samples":10000,"evals":22,"gcsample":false,"seconds":5.0,"overhead":0.0,"memory_tolerance":0.01}],"STaylor1(x::Array[length = 5])":["BenchmarkTools.Parameters",{"gctrial":true,"time_tolerance":0.05,"samples":10000,"evals":10,"gcsample":false,"seconds":5.0,"overhead":0.0,"memory_tolerance":0.01}],"STaylor1(x::Array[length = 50])":["BenchmarkTools.Parameters",{"gctrial":true,"time_tolerance":0.05,"samples":10000,"evals":10,"gcsample":false,"seconds":5.0,"overhead":0.0,"memory_tolerance":0.01}]},"tags":[]}]},"tags":[]}]]]

src/TaylorSeries.jl

@@ -42,7 +42,7 @@ import Base: zero, one, zeros, ones, isinf, isnan, iszero,
     power_by_squaring,
     rtoldefault, isfinite, isapprox, rad2deg, deg2rad
 
-export Taylor1, TaylorN, HomogeneousPolynomial, AbstractSeries
+export Taylor1, TaylorN, STaylor1, HomogeneousPolynomial, AbstractSeries
 
 export getcoeff, derivative, integrate, differentiate,
     evaluate, evaluate!, inverse, set_taylor1_varname,

src/arithmetic.jl

@@ -23,6 +23,8 @@ for T in (:Taylor1, :TaylorN)
     end
 end
 
+==(a::STaylor1, b::STaylor1) = (a.coeffs == b.coeffs)
+
 function ==(a::HomogeneousPolynomial, b::HomogeneousPolynomial)
     a.order == b.order && return a.coeffs == b.coeffs
     return iszero(a.coeffs) && iszero(b.coeffs)
@@ -31,11 +33,23 @@ end
 for T in (:Taylor1, :HomogeneousPolynomial, :TaylorN)
     @eval iszero(a::$T) = iszero(a.coeffs)
 end
+function iszero(a::STaylor1)
+    flag = true
+    for i in eachindex(a)
+        flag &= iszero(a[i])
+        if ~flag
+            break
+        end
+    end
+    flag
+end
 
 ## zero and one ##
 for T in (:Taylor1, :TaylorN), f in (:zero, :one)
     @eval ($f)(a::$T) = $T(($f)(a[0]), a.order)
 end
+zero(::STaylor1{N,T}) where {N, T<:Number} = STaylor1(zero(T),Val(N-1))
+one(::STaylor1{N,T}) where {N, T<:Number} = STaylor1(one(T),Val(N-1))
 
 function zero(a::HomogeneousPolynomial{T}) where {T<:Number}
     v = zero.(a.coeffs)
@@ -205,6 +219,48 @@ end
 -(a::Taylor1{T}, b::TaylorN{S}) where {T<:NumberNotSeries,S<:NumberNotSeries} =
     -(promote(a,b)...)
 
+@inline +(a::STaylor1{N,T}, b::STaylor1{N,T}) where {N, T<:Number} = STaylor1(add_tuples(a.coeffs, b.coeffs))
+@inline -(a::STaylor1{N,T}, b::STaylor1{N,T}) where {N, T<:Number} = STaylor1(sub_tuples(a.coeffs, b.coeffs))
+@inline +(a::STaylor1) = a
+@inline -(a::STaylor1) = STaylor1(minus_tuple(a.coeffs))
+
+@generated function add_first_tuple(tup::NTuple{N,T}, x::T) where {N, T<:Number}
+    ex = :((tup[1] + x),)
+    append!(ex.args, Any[:(tup[$i]) for i in 2:N])
+    return quote
+        Base.@_inline_meta
+        @inbounds return $ex
+    end
+end
+
+@generated function sub_first_tuple(tup::NTuple{N,T}, x::T) where {N, T<:Number}
+    ex = :((tup[1] - x),)
+    append!(ex.args, Any[:(tup[$i]) for i in 2:N])
+    return quote
+        Base.@_inline_meta
+        @inbounds return $ex
+    end
+end
+
+function +(a::STaylor1{N,T}, b::T) where {N, T<:Number}
+    STaylor1{N,T}(add_first_tuple(a.coeffs, b))
+end
+function +(b::T, a::STaylor1{N,T}) where {N, T<:Number}
+    STaylor1{N,T}(add_first_tuple(a.coeffs, b))
+end
+function -(a::STaylor1{N,T}, b::T) where {N, T<:Number}
+    STaylor1{N,T}(sub_first_tuple(a.coeffs, b))
+end
+-(b::T, a::STaylor1{N,T}) where {N, T<:Number}  = b + (-a)
+
+#+(a::STaylor1{N,T}, b::STaylor1{N,S}) where {N, T<:NumberNotSeries, S<:NumberNotSeries} = +(promote(a,b)...)
+#+(a::STaylor1{N,T}, b::S) where {N, T<:NumberNotSeries, S<:NumberNotSeries} = +(promote(a,b)...)
+#+(b::S, a::STaylor1{N,T}) where {N, T<:NumberNotSeries, S<:NumberNotSeries} = +(promote(b,a)...)
+
+#-(a::STaylor1{N,T}, b::STaylor1{N,S}) where {N, T<:NumberNotSeries, S<:NumberNotSeries} = -(promote(a,b)...)
+#-(a::STaylor1{N,T}, b::S) where {N, T<:NumberNotSeries, S<:NumberNotSeries} = -(promote(a,b)...)
+#-(b::S, a::STaylor1{N,T}) where {N, T<:NumberNotSeries, S<:NumberNotSeries} = -(promote(b,a)...)
+
 
 ## Multiplication ##
 for T in (:Taylor1, :HomogeneousPolynomial, :TaylorN)
@@ -225,6 +281,38 @@ for T in (:Taylor1, :HomogeneousPolynomial, :TaylorN)
     end
 end
 
+@generated function *(x::STaylor1{N,T}, y::STaylor1{N,T}) where {T<:Number,N}
+     ex_calc = quote end
+     append!(ex_calc.args, Any[nothing for i in 1:N])
+     syms = Symbol[Symbol("c$i") for i in 1:N]
+     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 = syms[j]
+         ex_line = :($sym = $ex_line)
+         ex_calc.args[j] = ex_line
+     end
+     exout = :(($(syms[1]),))
+     for i = 2:N
+         push!(exout.args, syms[i])
+     end
+     return quote
+                Base.@_inline_meta
+                $ex_calc
+                return STaylor1{N,T}($exout)
+             end
+end
+
+
+function *(a::STaylor1{N,T}, b::T) where {N, T<:Number}
+    STaylor1{N,T}(scale_tuple(a.coeffs, b))
+end
+function *(b::T, a::STaylor1{N,T}) where {N, T<:Number}
+    STaylor1{N,T}(scale_tuple(a.coeffs, b))
+end
+
 for T in (:HomogeneousPolynomial, :TaylorN)
 
     @eval begin
@@ -428,6 +516,49 @@ function /(a::Taylor1{T}, b::Taylor1{T}) where {T<:Number}
     return c
 end
 
+function /(a::STaylor1{N,T}, b::T) where {N, T<:Number}
+    STaylor1{N,T}(div_tuple_by_scalar(a.coeffs, b))
+end
+
+@generated function /(a::STaylor1{N,T}, b::STaylor1{N,T}) where {N,T<:Number}
+
+    ex_calc = quote end
+    append!(ex_calc.args, Any[nothing for i in 1:N])
+    syms = Symbol[Symbol("c$i") for i in 1:N]
+
+    # add error
+    ex_line = quote
+                  if iszero(b[0])
+                      throw(ArgumentError("""The 0th order STaylor1 coefficient
+                                             must be non-zero for b, (a/b)(x) is
+                                             not differentiable at x=0)."""))
+                  end
+              end
+    ex_calc.args[1] = ex_line
+
+    # add recursion relation
+    for j = 0:(N-1)
+        ex_line = :(a[$(j)])
+        for k = 1:j
+            sym = syms[j-k+1]
+            ex_line = :($ex_line - $sym*b[$k])
+        end
+        sym = syms[j+1]
+        ex_line = :($sym = ($ex_line)/b[0])
+        ex_calc.args[j+2] = ex_line
+    end
+
+    exout = :(($(syms[1]),))
+    for i = 2:N
+        push!(exout.args, syms[i])
+    end
+    return quote
+               Base.@_inline_meta
+               $ex_calc
+               return STaylor1{N,T}($exout)
+            end
+end
+
 /(a::TaylorN{T}, b::TaylorN{S}) where
     {T<:NumberNotSeriesN, S<:NumberNotSeriesN} = /(promote(a,b)...)
 

src/auxiliary.jl

@@ -104,6 +104,12 @@ function setindex!(a::Taylor1{T}, x::Array{T,1}, u::StepRange{Int,Int}) where {T
 end
 
 
+# getindex STaylor1
+getindex(a::STaylor1, n::Int) = a.coeffs[n+1]
+getindex(a::STaylor1, u::UnitRange{Int}) = a.coeffs[u .+ 1]
+getindex(a::STaylor1, c::Colon) = a.coeffs[c]
+getindex(a::STaylor1, u::StepRange{Int,Int}) = a.coeffs[u .+ 1]
+
 """
     getcoeff(a, v)
 
@@ -236,6 +242,15 @@ for T in (:Taylor1, :HomogeneousPolynomial, :TaylorN)
         @inline axes(a::$T) = ()
     end
 end
+@inline iterate(a::STaylor1{N,T}, state=0) where {N, T<:Number} = state > N-1 ? nothing : (a.coeffs[state+1], state+1)
+@inline firstindex(a::STaylor1) = 0
+@inline lastindex(a::STaylor1{N,T}) where {N, T<:Number} = N-1
+@inline eachindex(s::STaylor1{N,T}) where {N, T<:Number} = UnitRange(0, N-1)
+@inline size(s::STaylor1{N,T}) where {N, T<:Number} = N
+@inline length(s::STaylor1{N,T}) where {N, T<:Number} = N
+@inline get_order(s::STaylor1{N,T}) where {N, T<:Number} = N - 1
+@inline eltype(s::STaylor1{N,T}) where {N, T<:Number} = T
+@inline axes(a::STaylor1) = ()
 
 
 ## fixorder ##
@@ -271,6 +286,18 @@ function Base.findlast(a::Taylor1{T}) where {T<:Number}
     return last-1
 end
 
+function Base.findfirst(a::STaylor1{N,T}) where {N, T<:Number}
+    first = findfirst(x->!iszero(x), a.coeffs)
+    isa(first, Nothing) && return -1
+    return first-1
+end
+# Finds the last non-zero entry; extended to Taylor1
+function Base.findlast(a::STaylor1{N,T}) where {N, T<:Number}
+    last = findlast(x->!iszero(x), a.coeffs)
+    isa(last, Nothing) && return -1
+    return last-1
+end
+
 
 ## copyto! ##
 # Inspired from base/abstractarray.jl, line 665
@@ -297,6 +324,8 @@ and `TaylorN`. The fallback behavior is to return `a` itself if
 """
 constant_term(a::Taylor1) = a[0]
 
+constant_term(a::STaylor1) = a[0]
+
 constant_term(a::TaylorN) = a[0][1]
 
 constant_term(a::Vector{T}) where {T<:Number} = constant_term.(a)
@@ -317,3 +346,46 @@ linear_polynomial(a::TaylorN) = TaylorN([a[1]])
 linear_polynomial(a::Vector{T}) where {T<:Number} = linear_polynomial.(a)
 
 linear_polynomial(a::Number) = a
+
+
+#=
+Tuple operations taken from ForwardDiff.jl
+
+Copyright (c) 2015: Jarrett Revels, Theodore Papamarkou, Miles Lubin, and other contributors
+
+Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:
+The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
+
+=#
+function tupexpr(f, N)
+    ex = Expr(:tuple, [f(i) for i=1:N]...)
+    return quote
+        $(Expr(:meta, :inline))
+        @inbounds return $ex
+    end
+end
+
+@generated function scale_tuple(tup::NTuple{N}, x) where N
+    return tupexpr(i -> :(tup[$i] * x), N)
+end
+
+@generated function div_tuple_by_scalar(tup::NTuple{N}, x) where N
+    return tupexpr(i -> :(tup[$i] / x), N)
+end
+
+@generated function add_tuples(a::NTuple{N}, b::NTuple{N})  where N
+    return tupexpr(i -> :(a[$i] + b[$i]), N)
+end
+
+@generated function sub_tuples(a::NTuple{N}, b::NTuple{N})  where N
+    return tupexpr(i -> :(a[$i] - b[$i]), N)
+end
+
+@generated function minus_tuple(tup::NTuple{N}) where N
+    return tupexpr(i -> :(-tup[$i]), N)
+end
+
+@generated function mul_tuples(a::NTuple{N}, b::NTuple{N}, afactor, bfactor) where N
+    return tupexpr(i -> :((afactor * a[$i]) + (bfactor * b[$i])), N)
+end

src/broadcasting.jl

@@ -12,7 +12,7 @@ import .Broadcast: BroadcastStyle, Broadcasted, broadcasted
 
 # BroadcastStyle definitions and basic precedence rules
 struct Taylor1Style{T} <: Base.Broadcast.AbstractArrayStyle{0} end
-Taylor1Style{T}(::Val{N}) where {T, N}= Base.Broadcast.DefaultArrayStyle{N}()
+Taylor1Style{T}(::Val{N}) where {T, N} = Base.Broadcast.DefaultArrayStyle{N}()
 BroadcastStyle(::Type{<:Taylor1{T}}) where {T} = Taylor1Style{T}()
 BroadcastStyle(::Taylor1Style{T}, ::Base.Broadcast.DefaultArrayStyle{0}) where {T} = Taylor1Style{T}()
 BroadcastStyle(::Taylor1Style{T}, ::Base.Broadcast.DefaultArrayStyle{1}) where {T} = Base.Broadcast.DefaultArrayStyle{1}()