Merge 8b144de into 5e79642

src/TaylorSeries.jl

@@ -32,7 +32,7 @@ end
 import Base: ==, +, -, *, /, ^
 
 import Base: iterate, size, eachindex, firstindex, lastindex,
-    eltype, length, getindex, setindex!
+    eltype, length, getindex, setindex!, axes, copyto!
 
 import Base: zero, one, zeros, ones, isinf, isnan, iszero,
     convert, promote_rule, promote, show,
@@ -67,6 +67,7 @@ include("other_functions.jl")
 include("evaluate.jl")
 include("calculus.jl")
 include("dictmutfunct.jl")
+include("broadcasting.jl")
 include("printing.jl")
 
 function __init__()

src/auxiliary.jl

@@ -19,7 +19,9 @@ function resize_coeffs1!(coeffs::Array{T,1}, order::Int) where {T<:Number}
     resize!(coeffs, order+1)
     if order > lencoef-1
         z = zero(coeffs[1])
-        @__dot__ coeffs[lencoef+1:order+1] = z
+        @simd for ord in lencoef+1:order+1
+            @inbounds coeffs[ord] = z
+        end
     end
     return nothing
 end
@@ -38,7 +40,9 @@ function resize_coeffsHP!(coeffs::Array{T,1}, order::Int) where {T<:Number}
     num_coeffs == lencoef && return nothing
     resize!(coeffs, num_coeffs)
     z = zero(coeffs[1])
-    @__dot__ coeffs[lencoef+1:num_coeffs] = z
+    @simd for ord in lencoef+1:num_coeffs
+        @inbounds coeffs[ord] = z
+    end
     return nothing
 end
 
@@ -175,31 +179,30 @@ setindex!(a::TaylorN{T}, x::Array{T,1}, ::Colon) where {T<:Number} =
     (a[0:end] = x; a[:])
 
 
-## eltype, length, get_order ##
-for T in (:Taylor1, :TaylorN)
+## eltype, length, get_order, etc ##
+for T in (:Taylor1, :HomogeneousPolynomial, :TaylorN)
     @eval begin
-        iterate(a::$T, state=0) = state > a.order ? nothing : (a.coeffs[state+1], state+1)
-        eachindex(a::$T) = firstindex(a):lastindex(a)
-        # Base.iterate(rS::Iterators.Reverse{$T}, state=rS.itr.order) = state < 0 ? nothing : (a.coeffs[state], state-1)
-        eltype(::$T{S}) where {S<:Number} = S
-        length(a::$T) = length(a.coeffs)
-        size(a::$T) = (length(a),)
-        firstindex(a::$T) = 0
-        lastindex(a::$T) = a.order
-        get_order(a::$T) = a.order
+        if $T == HomogeneousPolynomial
+            @inline iterate(a::$T, state=1) = state > length(a) ? nothing : (a.coeffs[state], state+1)
+            # Base.iterate(rS::Iterators.Reverse{$T}, state=rS.itr.order) = state < 0 ? nothing : (a.coeffs[state], state-1)
+            @inline length(a::$T) = size_table[a.order+1]
+            @inline firstindex(a::$T) = 1
+            @inline lastindex(a::$T) = length(a)
+        else
+            @inline iterate(a::$T, state=0) = state > a.order ? nothing : (a.coeffs[state+1], state+1)
+            # Base.iterate(rS::Iterators.Reverse{$T}, state=rS.itr.order) = state < 0 ? nothing : (a.coeffs[state], state-1)
+            @inline length(a::$T) = length(a.coeffs)
+            @inline firstindex(a::$T) = 0
+            @inline lastindex(a::$T) = a.order
+        end
+        @inline eachindex(a::$T) = firstindex(a):lastindex(a)
+        @inline eltype(::$T{S}) where {S<:Number} = S
+        @inline size(a::$T) = size(a.coeffs)
+        @inline get_order(a::$T) = a.order
+        @inline axes(a::$T) = ()
     end
 end
 
-# Base.iterate(a::HomogeneousPolynomial, state=1) = state > a.order, nothing : (a.coeffs[state+1], state+1)
-# Base.iterate(rS::Iterators.Reverse{$T}, state=rS.itr.order) = state < 0 ? nothing : (a.coeffs[state], state-1)
-# Base.eachindex(a::HomogeneousPolynomial) = firstindex(a):lastindex(a)
-eltype(::HomogeneousPolynomial{S}) where {S<:Number} = S
-length(a::HomogeneousPolynomial) = size_table[a.order+1]#length(a.coeffs)
-size(a::HomogeneousPolynomial) = (length(a),)
-firstindex(a::HomogeneousPolynomial) = 1
-lastindex(a::HomogeneousPolynomial) = length(a)
-get_order(a::HomogeneousPolynomial) = a.order
-
 
 ## fixorder ##
 for T in (:Taylor1, :TaylorN)
@@ -232,6 +235,23 @@ function Base.findlast(a::Taylor1{T}) where {T<:Number}
     return last-1
 end
 
+
+## copyto! ##
+# Inspired from base/abstractarray.jl, line 665
+for T in (:Taylor1, :HomogeneousPolynomial, :TaylorN)
+    @eval function copyto!(dst::$T{T}, src::$T{T}) where {T}
+        length(dst) < length(src) && throw(ArgumentError(string("Destination has fewer elements than required; no copy performed")))
+        destiter = eachindex(dst)
+        y = iterate(destiter)
+        for x in src
+            dst[y[1]] = x
+            y = iterate(destiter, y[2])
+        end
+        return dst
+    end
+end
+
+
 """
     constant_term(a)
 
@@ -243,7 +263,7 @@ constant_term(a::Taylor1) = a[0]
 
 constant_term(a::TaylorN) = a[0][1]
 
-constant_term(a::Vector{T}) where {T<:Number}= a[1]
+constant_term(a::Vector{T}) where {T<:Number} = constant_term.(a)
 
 constant_term(a::Number) = a
 
@@ -258,6 +278,6 @@ linear_polynomial(a::Taylor1) = Taylor1([zero(a[1]), a[1]])
 
 linear_polynomial(a::TaylorN) = TaylorN([a[1]])
 
-linear_polynomial(a::Vector{T}) where {T<:Number} = a[1]
+linear_polynomial(a::Vector{T}) where {T<:Number} = linear_polynomial.(a)
 
 linear_polynomial(a::Number) = a

src/broadcasting.jl

@@ -0,0 +1,141 @@
+# This file is part of the TaylorSeries.jl Julia package, MIT license
+#
+# Luis Benet & David P. Sanders
+# UNAM
+#
+# MIT Expat license
+#
+
+## Broadcast for Taylor1 and TaylorN
+
+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}()
+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}()
+#
+struct HomogeneousPolynomialStyle{T} <: Base.Broadcast.AbstractArrayStyle{0} end
+HomogeneousPolynomialStyle{T}(::Val{N}) where {T, N}= Base.Broadcast.DefaultArrayStyle{N}()
+BroadcastStyle(::Type{<:HomogeneousPolynomial{T}}) where {T} = HomogeneousPolynomialStyle{T}()
+BroadcastStyle(::HomogeneousPolynomialStyle{T}, ::Base.Broadcast.DefaultArrayStyle{0}) where {T} = HomogeneousPolynomialStyle{T}()
+BroadcastStyle(::HomogeneousPolynomialStyle{T}, ::Base.Broadcast.DefaultArrayStyle{1}) where {T} = Base.Broadcast.DefaultArrayStyle{1}()
+#
+struct TaylorNStyle{T} <: Base.Broadcast.AbstractArrayStyle{0} end
+TaylorNStyle{T}(::Val{N}) where {T, N}= Base.Broadcast.DefaultArrayStyle{N}()
+BroadcastStyle(::Type{<:TaylorN{T}}) where {T} = TaylorNStyle{T}()
+BroadcastStyle(::TaylorNStyle{T}, ::Base.Broadcast.DefaultArrayStyle{0}) where {T} = TaylorNStyle{T}()
+BroadcastStyle(::TaylorNStyle{T}, ::Base.Broadcast.DefaultArrayStyle{1}) where {T} = Base.Broadcast.DefaultArrayStyle{1}()
+
+# Precedence rules for mixtures
+BroadcastStyle(::TaylorNStyle{Taylor1{T}}, ::Taylor1Style{T}) where {T} = TaylorNStyle{Taylor1{T}}()
+BroadcastStyle(::Taylor1Style{TaylorN{T}}, ::TaylorNStyle{T}) where {T} = Taylor1Style{TaylorN{T}}()
+
+# Extend eltypes so things like [1.0] .+ t work
+Base.Broadcast.eltypes(t::Tuple{Taylor1,AbstractArray}) =
+    Tuple{Base.Broadcast._broadcast_getindex_eltype([t[1]]), Base.Broadcast._broadcast_getindex_eltype(t[2])}
+Base.Broadcast.eltypes(t::Tuple{AbstractArray,Taylor1}) =
+    Tuple{Base.Broadcast._broadcast_getindex_eltype(t[1]), Base.Broadcast._broadcast_getindex_eltype([t[2]])}
+Base.Broadcast.eltypes(t::Tuple{HomogeneousPolynomial,AbstractArray}) =
+    Tuple{Base.Broadcast._broadcast_getindex_eltype([t[1]]), Base.Broadcast._broadcast_getindex_eltype(t[2])}
+Base.Broadcast.eltypes(t::Tuple{AbstractArray,HomogeneousPolynomial}) =
+    Tuple{Base.Broadcast._broadcast_getindex_eltype(t[1]), Base.Broadcast._broadcast_getindex_eltype([t[2]])}
+Base.Broadcast.eltypes(t::Tuple{TaylorN,AbstractArray}) =
+    Tuple{Base.Broadcast._broadcast_getindex_eltype([t[1]]), Base.Broadcast._broadcast_getindex_eltype(t[2])}
+Base.Broadcast.eltypes(t::Tuple{AbstractArray,TaylorN}) =
+    Tuple{Base.Broadcast._broadcast_getindex_eltype(t[1]), Base.Broadcast._broadcast_getindex_eltype([t[2]])}
+
+
+# # We follow https://docs.julialang.org/en/v1/manual/interfaces/#man-interface-iteration-1
+# "`A = find_taylor(As)` returns the first Taylor1 among the arguments."
+# find_taylor(bc::Broadcasted) = find_taylor(bc.args)
+# find_taylor(args::Tuple) = find_taylor(find_taylor(args[1]), Base.tail(args))
+# find_taylor(x) = x
+# find_taylor(a::Taylor1, rest) = a
+# find_taylor(a::HomogeneousPolynomial, rest) = a
+# find_taylor(a::TaylorN, rest) = a
+# find_taylor(::AbstractArray, rest) = find_taylor(rest)
+#
+# # Extend similar
+# function similar(bc::Broadcasted{Taylor1Style{S}}, ::Type{T}) where {S, T}
+#     # Proper promotion
+#     R = Base.Broadcast.combine_eltypes(bc.f, bc.args)
+#     # Scan the inputs for the Taylor1:
+#     A = find_taylor(bc)
+#     # Create the output
+#     return Taylor1(similar(A.coeffs, R), A.order)
+# end
+#
+# function similar(bc::Broadcasted{HomogeneousPolynomialStyle{S}}, ::Type{T}) where {S, T}
+#     # Proper promotion
+#     # combine_eltypes(f, args::Tuple) = Base._return_type(f, eltypes(args))
+#     R = Base.Broadcast.combine_eltypes(bc.f, bc.args)
+#     # Scan the inputs for the HomogeneousPolynomial:
+#     A = find_taylor(bc)
+#     # Create the output
+#     return HomogeneousPolynomial(similar(A.coeffs, R), A.order)
+# end
+#
+# function similar(bc::Broadcasted{TaylorNStyle{S}}, ::Type{T}) where {S, T}
+#     # Proper promotion
+#     R = Base.Broadcast.combine_eltypes(bc.f, bc.args)
+#     # Scan the inputs for the TaylorN:
+#     A = find_taylor(bc)
+#     # Create the output
+#     return TaylorN(similar(A.coeffs, R), A.order)
+# end
+
+
+# Adapted from Base.Broadcast.copyto!, base/broadcasting.jl, line 832
+for T in (:Taylor1, :HomogeneousPolynomial, :TaylorN)
+    @eval begin
+        @inline function copyto!(dest::$T{T}, bc::Broadcasted) where {T}
+            axes(dest) == axes(bc) || Base.Broadcast.throwdm(axes(dest), axes(bc))
+            # Performance optimization: broadcast!(identity, dest, A) is equivalent to copyto!(dest, A) if indices match
+            if bc.f === identity && bc.args isa Tuple{$T{T}} # only a single input argument to broadcast!
+                A = bc.args[1]
+                if axes(dest) == axes(A)
+                    return copyto!(dest, A)
+                end
+            end
+            bc′ = Base.Broadcast.preprocess(dest, bc)
+            copyto!(dest, bc′[1])
+            return dest
+        end
+    end
+end
+
+
+# Broadcasted extensions
+@inline broadcasted(::Taylor1Style{T}, ::Type{Float32}, a::Taylor1{T}) where {T} =
+    Taylor1(Float32.(a.coeffs), a.order)
+@inline broadcasted(::TaylorNStyle{T}, ::Type{Float32}, a::TaylorN{T}) where {T} =
+    convert(TaylorN{Float32}, a)
+
+# # This prevents broadcasting being applied to the Taylor1/TaylorN params
+# # for the mutating functions, and to act only in `k`
+# for (T, TS) in ((:Taylor1, :Taylor1Style), (:TaylorN, :TaylorNStyle))
+#     for f in (add!, subst!, sqr!, sqrt!, exp!, log!, identity!, zero!,
+#         one!, abs!, abs2!, deg2rad!, rad2deg!)
+#         @eval begin
+#             @inline function broadcasted(::$TS{T}, fn::typeof($f), r::$T{T}, a::$T{T}, k) where {T}
+#                 @inbounds for i in eachindex(k)
+#                     fn(r, a, k[i])
+#                 end
+#                 nothing
+#             end
+#         end
+#     end
+#     for f in (sincos!, tan!, asin!, acos!, atan!, sinhcosh!, tanh!)
+#         @eval begin
+#             @inline function broadcasted(::$TS{T}, fn::typeof($f), r::$T{T}, a::$T{T}, b::Taylor1{T}, k) where {T}
+#                 @inbounds for i in eachindex(k)
+#                     fn(r, a, b, k[i])
+#                 end
+#                 nothing
+#             end
+#         end
+#     end
+# end

src/evaluate.jl

@@ -126,7 +126,7 @@ function evaluate(a::Taylor1{Taylor1{T}}, x::Taylor1{T}) where {T<:Number}
 end
 
 evaluate(p::Taylor1{T}, x::Array{S}) where {T<:Number, S<:Number} =
-    evaluate.(p, x)
+    evaluate.([p], x)
 
 #function-like behavior for Taylor1
 (p::Taylor1)(x) = evaluate(p, x)

src/other_functions.jl

@@ -13,7 +13,7 @@ for T in (:Taylor1, :HomogeneousPolynomial, :TaylorN)
         @eval ($f)(a::$T) = $T(($f).(a.coeffs), a.order)
     end
 
-    @eval adjoint(a::$T) = conj.(a)
+    @eval adjoint(a::$T) = conj(a)
 
     ## isinf and isnan ##
     @eval isinf(a::$T) = any( isinf.(a.coeffs) )