Add evaluation method for arrays of TaylorN (#364)

* WIP: add `evaluate!` method for `TaylorN`s

* Update

* Small fix

* pow! fix for Julia 1.6

* Remove unused `@inbounds`

* Bump minor version

* Update README.md

* Small fix and add comments

* Update tests

* Update tests again

Author: PerezHz

Project.toml

@@ -1,7 +1,7 @@
 name = "TaylorSeries"
 uuid = "6aa5eb33-94cf-58f4-a9d0-e4b2c4fc25ea"
 repo = "https://github.com/JuliaDiff/TaylorSeries.jl.git"
-version = "0.17.8"
+version = "0.18.0"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
@@ -12,14 +12,14 @@ SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 [weakdeps]
 IntervalArithmetic = "d1acc4aa-44c8-5952-acd4-ba5d80a2a253"
 JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
-StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd"
+StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 
 [extensions]
 TaylorSeriesIAExt = "IntervalArithmetic"
 TaylorSeriesJLD2Ext = "JLD2"
-TaylorSeriesSAExt = "StaticArrays"
 TaylorSeriesRATExt = "RecursiveArrayTools"
+TaylorSeriesSAExt = "StaticArrays"
 
 [compat]
 Aqua = "0.8"

README.md

@@ -30,7 +30,7 @@ julia> t = Taylor1(Float64, 5)
 
 julia> exp(t)
  1.0 + 1.0 t + 0.5 t² + 0.16666666666666666 t³ + 0.041666666666666664 t⁴ + 0.008333333333333333 t⁵ + 𝒪(t⁶)
- 
+
  julia> log(1 + t)
  1.0 t - 0.5 t² + 0.3333333333333333 t³ - 0.25 t⁴ + 0.2 t⁵ + 𝒪(t⁶)
  ```
@@ -40,7 +40,7 @@ julia> x, y = set_variables("x y", order=2);
 
 julia> exp(x + y)
  1.0 + 1.0 x + 1.0 y + 0.5 x² + 1.0 x y + 0.5 y² + 𝒪(‖x‖³)
- 
+
 ```
 Differential and integral calculus on Taylor series:
 ```julia
@@ -95,6 +95,6 @@ We thank the participants of the course for putting up with the half-baked
 material and contributing energy and ideas.
 
 We acknowledge financial support from DGAPA-UNAM PAPIME grants
-PE-105911 and PE-107114, and DGAPA-PAPIIT grants IG-101113, 
-IG-100616, and IG-100819.
+PE-105911 and PE-107114, and DGAPA-PAPIIT grants IG-101113,
+IG-100616, IG-100819 and IG-101122.
 LB acknowledges support through a *Cátedra Moshinsky* (2013).

src/arithmetic.jl

@@ -162,9 +162,7 @@ for T in (:Taylor1, :HomogeneousPolynomial, :TaylorN)
 end
 
 for T in (:Taylor1, :TaylorN)
-    @eval function zero(a::$T)
-        return $T(zero.(a.coeffs))
-    end
+    @eval zero(a::$T) = $T(zero.(a.coeffs))
     @eval function one(a::$T)
         b = zero(a)
         b[0] = one(b[0])
@@ -539,25 +537,50 @@ for T in (:Taylor1, :TaylorN)
         return nothing
     end
 
-    @eval @inline function mul!(v::$T, a::$T, b::NumberNotSeries, k::Int)
+    @eval begin
         if $T == Taylor1
-            @inbounds v[k] = a[k] * b
-        else
-            @inbounds for i in eachindex(v[k])
-                v[k][i] = a[k][i] * b
+            @inline function mul!(v::$T, a::$T, b::NumberNotSeries, k::Int)
+                @inbounds v[k] = a[k] * b
+                return nothing
+            end
+            @inline function mul!(v::$T, a::NumberNotSeries, b::$T, k::Int)
+                @inbounds v[k] = a * b[k]
+                return nothing
+            end
+            @inline function muladd!(v::$T, a::$T, b::NumberNotSeries, k::Int)
+                @inbounds v[k] += a[k] * b
+                return nothing
+            end
+            @inline function muladd!(v::$T, a::NumberNotSeries, b::$T, k::Int)
+                @inbounds v[k] += a * b[k]
+                return nothing
             end
-        end
-        return nothing
-    end
-    @eval @inline function mul!(v::$T, a::NumberNotSeries, b::$T, k::Int)
-        if $T == Taylor1
-            @inbounds v[k] = a * b[k]
         else
-            @inbounds for i in eachindex(v[k])
-                v[k][i] = a * b[k][i]
+            @inline function mul!(v::$T, a::$T, b::NumberNotSeries, k::Int)
+                @inbounds for i in eachindex(v[k])
+                    v[k][i] = a[k][i] * b
+                end
+                return nothing
+            end
+            @inline function mul!(v::$T, a::NumberNotSeries, b::$T, k::Int)
+                @inbounds for i in eachindex(v[k])
+                    v[k][i] = a * b[k][i]
+                end
+                return nothing
+            end
+            @inline function muladd!(v::$T, a::$T, b::NumberNotSeries, k::Int)
+                @inbounds for i in eachindex(v[k])
+                    v[k][i] += a[k][i] * b
+                end
+                return nothing
+            end
+            @inline function muladd!(v::$T, a::NumberNotSeries, b::$T, k::Int)
+                @inbounds for i in eachindex(v[k])
+                    v[k][i] += a * b[k][i]
+                end
+                return nothing
             end
         end
-        return nothing
     end
 
     @eval @inline function mul!(v::$T, a::$T, b::NumberNotSeries)
@@ -951,7 +974,8 @@ end
 
 @inline function div!(c::Taylor1{T}, a::NumberNotSeries,
         b::Taylor1{T}, k::Int) where {T<:Number}
-    iszero(a) && !iszero(b) && zero!(c, k)
+    zero!(c, k)
+    iszero(a) && !iszero(b) && return nothing
     # order and coefficient of first factorized term
     # In this case, since a[k]=0 for k>0, we can simplify to:
     # ordfact, cdivfact = 0, a/b[0]
@@ -970,7 +994,8 @@ end
 
 @inline function div!(c::Taylor1{TaylorN{T}}, a::NumberNotSeries,
         b::Taylor1{TaylorN{T}}, k::Int) where {T<:NumberNotSeries}
-    iszero(a) && !iszero(b) && zero!(c, k)
+    zero!(c, k)
+    iszero(a) && !iszero(b) && return nothing
     # order and coefficient of first factorized term
     # In this case, since a[k]=0 for k>0, we can simplify to:
     # ordfact, cdivfact = 0, a/b[0]
@@ -1141,14 +1166,14 @@ end
         """Division does not define a Taylor1 polynomial;
         order k=$(ordfact) => coeff[$(ordfact)]=$(cdivfact).""") )
 
+    zero!(c, k)
+
     if k == 0
         # @inbounds c[0] = a[ordfact]/b[ordfact]
         @inbounds div!(c[0], a[ordfact], b[ordfact])
         return nothing
     end
 
-    @inbounds zero!(c, k)
-
     imin = max(0, k+ordfact-b.order)
     @inbounds mul!(c[k], c[imin], b[k+ordfact-imin])
     @inbounds for i = imin+1:k-1

src/evaluate.jl

@@ -187,43 +187,42 @@ function _evaluate(a::HomogeneousPolynomial{T}, vals::NTuple) where {T}
     return suma
 end
 
-function _evaluate(a::HomogeneousPolynomial{T}, vals::NTuple{N,<:TaylorN{T}}) where
-        {N,T<:NumberNotSeries}
-    # @assert length(vals) == get_numvars()
-    a.order == 0 && return a[1]*one(vals[1])
+function _evaluate!(res::TaylorN{T}, a::HomogeneousPolynomial{T}, vals::NTuple{N,<:TaylorN{T}},
+            valscache::Vector{TaylorN{T}}, aux::TaylorN{T}) where {N,T<:NumberNotSeries}
     ct = coeff_table[a.order+1]
-    suma = TaylorN(zero(T), vals[1].order)
-    #
-    # vv = power_by_squaring.(vals, ct[1])
-    vv = vals .^ ct[1]
-    tmp = zero(suma)
-    aux = one(suma)
+    for el in eachindex(valscache)
+        power_by_squaring!(valscache[el], vals[el], aux, ct[1][el])
+    end
     for (i, a_coeff) in enumerate(a.coeffs)
         iszero(a_coeff) && continue
-        # @inbounds vv .= power_by_squaring.(vals, ct[i])
-        vv .= vals .^ ct[i]
-        # tmp = prod( vv )
-        for ord in eachindex(tmp)
-            @inbounds one!(aux, vv[1], ord)
+        # valscache .= vals .^ ct[i]
+        @inbounds for el in eachindex(valscache)
+            power_by_squaring!(valscache[el], vals[el], aux, ct[i][el])
         end
-        for j in eachindex(vv)
-            for ord in eachindex(tmp)
-                zero!(tmp, ord)
-                @inbounds mul!(tmp, aux, vv[j], ord)
-            end
-            for ord in eachindex(tmp)
-                identity!(aux, tmp, ord)
-            end
+        # aux = one(valscache[1])
+        for ord in eachindex(aux)
+            @inbounds one!(aux, valscache[1], ord)
         end
-        # suma += a_coeff * tmp
-        for ord in eachindex(tmp)
-            for ordQ in eachindex(tmp[ord])
-                zero!(aux[ord], ordQ)
-                aux[ord][ordQ] = a_coeff * tmp[ord][ordQ]
-                suma[ord][ordQ] += aux[ord][ordQ]
-            end
+        for j in eachindex(valscache)
+            # aux *= valscache[j]
+            mul!(aux, valscache[j])
+        end
+        # res += a_coeff * aux
+        for ord in eachindex(aux)
+            muladd!(res, a_coeff, aux, ord)
         end
     end
+    return nothing
+end
+
+function _evaluate(a::HomogeneousPolynomial{T}, vals::NTuple{N,<:TaylorN{T}}) where
+        {N,T<:NumberNotSeries}
+    # @assert length(vals) == get_numvars()
+    a.order == 0 && return a[1]*one(vals[1])
+    suma = TaylorN(zero(T), vals[1].order)
+    valscache = [zero(val) for val in vals]
+    aux = zero(suma)
+    _evaluate!(suma, a, vals, valscache, aux)
     return suma
 end
 
@@ -319,6 +318,17 @@ _evaluate(a::TaylorN{T}, vals::NTuple, ::Val{true}) where {T<:NumberNotSeries} =
 _evaluate(a::TaylorN{T}, vals::NTuple, ::Val{false}) where {T<:Number} =
     sum( _evaluate(a, vals) )
 
+function _evaluate(a::TaylorN{T}, vals::NTuple{N,<:TaylorN}, ::Val{false}) where {N,T<:Number}
+    R = promote_type(T, TS.numtype(vals[1]))
+    res = TaylorN(zero(R), vals[1].order)
+    valscache = [zero(val) for val in vals]
+    aux = zero(res)
+    @inbounds for homPol in eachindex(a)
+        _evaluate!(res, a[homPol], vals, valscache, aux)
+    end
+    return res
+end
+
 function _evaluate(a::TaylorN{T}, vals::NTuple{N,<:Number}) where {N,T<:Number}
     R = promote_type(T, typeof(vals[1]))
     a_length = length(a)
@@ -329,13 +339,20 @@ function _evaluate(a::TaylorN{T}, vals::NTuple{N,<:Number}) where {N,T<:Number}
     return suma
 end
 
-function _evaluate(a::TaylorN{T}, vals::NTuple{N,<:TaylorN}) where {N,T<:Number}
-    R = TaylorN{promote_type(T, TS.numtype(vals[1]))}
-    a_length = length(a)
-    suma = zeros(R, a_length)
+function _evaluate!(res::Vector{TaylorN{T}}, a::TaylorN{T}, vals::NTuple{N,<:TaylorN},
+        valscache::Vector{TaylorN{T}}, aux::TaylorN{T}) where {N,T<:Number}
     @inbounds for homPol in eachindex(a)
-        suma[homPol+1] = _evaluate(a[homPol], vals)
+        _evaluate!(res[homPol+1], a[homPol], vals, valscache, aux)
     end
+    return nothing
+end
+
+function _evaluate(a::TaylorN{T}, vals::NTuple{N,<:TaylorN}) where {N,T<:Number}
+    R = promote_type(T, TS.numtype(vals[1]))
+    suma = [TaylorN(zero(R), vals[1].order) for _ in eachindex(a)]
+    valscache = [zero(val) for val in vals]
+    aux = zero(suma[1])
+    _evaluate!(suma, a, vals, valscache, aux)
     return suma
 end
 
@@ -486,6 +503,24 @@ function evaluate!(x::AbstractArray{TaylorN{T}}, δx::Array{TaylorN{T},1},
     return nothing
 end
 
+function evaluate!(a::TaylorN{T}, vals::NTuple{N,TaylorN{T}}, dest::TaylorN{T}, valscache::Vector{TaylorN{T}}, aux::TaylorN{T}) where {N,T<:Number}
+    @inbounds for homPol in eachindex(a)
+        _evaluate!(dest, a[homPol], vals, valscache, aux)
+    end
+    return nothing
+end
+
+function evaluate!(a::AbstractArray{TaylorN{T}}, vals::NTuple{N,TaylorN{T}}, dest::AbstractArray{TaylorN{T}}) where {N,T<:Number}
+    # initialize evaluation cache
+    valscache = [zero(val) for val in vals]
+    aux = zero(dest[1])
+    # loop over elements of `a`
+    for i in eachindex(a)
+        (!iszero(dest[i])) && zero!(dest[i])
+        evaluate!(a[i], vals, dest[i], valscache, aux)
+    end
+    return nothing
+end
 
 # In-place Horner methods, used when the result of an evaluation (substitution)
 # is Taylor1{}

src/functions.jl

@@ -354,12 +354,18 @@ for T in (:Taylor1, :TaylorN)
             return nothing
         end
 
-        @inline function one!(c::$T{T}, a::$T{T}, k::Int) where {T<:Number}
-            zero!(c, k)
-            if k == 0
-                @inbounds c[0] = one(a[0])
+        if $T == Taylor1
+            @inline function one!(c::$T{T}, a::$T{T}, k::Int) where {T<:Number}
+                zero!(c, k)
+                (k == 0) && (@inbounds c[0] = one(constant_term(a)))
+                return nothing
+            end
+        else
+            @inline function one!(c::$T{T}, a::$T{T}, k::Int) where {T<:Number}
+                zero!(c, k)
+                (k == 0) && (@inbounds c[0][1] = one(constant_term(a)))
+                return nothing
             end
-            return nothing
         end
 
         @inline function abs!(c::$T{T}, a::$T{T}, k::Int) where {T<:Number}

src/power.jl

@@ -62,8 +62,15 @@ end
 # this method assumes `y`, `x` and `aux` are of same order
 # TODO: add power_by_squaring! method for HomogeneousPolynomial and mixtures
 for T in (:Taylor1, :TaylorN)
-    @eval function power_by_squaring!(y::$T{T}, x::$T{T}, aux::$T{T},
+    @eval @inline function power_by_squaring_0!(y::$T{T}, x::$T{T}) where {T<:NumberNotSeries}
+        for k in eachindex(y)
+            one!(y, x, k)
+        end
+        return nothing
+    end
+    @eval @inline function power_by_squaring!(y::$T{T}, x::$T{T}, aux::$T{T},
             p::Integer) where {T<:NumberNotSeries}
+        (p == 0) && return power_by_squaring_0!(y, x)
         t = trailing_zeros(p) + 1
         p >>= t
         # aux = x
@@ -272,7 +279,7 @@ exploits `k_0`, the order of the first non-zero coefficient of `a`.
     isinteger(r) && r > 0 && (k > r*findlast(a)) && return nothing
 
     if k == lnull
-        @inbounds c[k] = (a[l0])^r
+        @inbounds c[k] = (a[l0])^float(r)
         return nothing
     end
 

test/manyvariables.jl

@@ -821,6 +821,33 @@ end
         @test float(TaylorN{Complex{Rational}}) == float(TaylorN{Complex{Float64}})
     end
 
+    @testset "Test evaluate! method for arrays of TaylorN" begin
+        x = set_variables("x", order=2, numvars=4)
+        function radntn!(y)
+            for k in eachindex(y)
+                for l in eachindex(y[k])
+                    y[k][l] = randn()
+                end
+            end
+            nothing
+        end
+        y = zero(x[1])
+        radntn!(y)
+        n = 10
+        v = [zero(x[1]) for _ in 1:n]
+        r = [zero(x[1]) for _ in 1:n] # output vector
+        radntn!.(v)
+        x1 = randn(4) .+ x
+        # warmup
+        TaylorSeries.evaluate!(v, (x1...,), r)
+        # call twice to make sure `r` is reset on second call
+        TaylorSeries.evaluate!(v, (x1...,), r)
+        r2 = evaluate.(v, Ref(x1))
+        @test r == r2
+        @test iszero(norm(r-r2, Inf))
+
+    end
+
 end
 
 @testset "Integrate for several variables" begin