Skip to content

Commit

Permalink
Broadcasting (#207)
Browse files Browse the repository at this point in the history
* Extend axes and enforce inlining

* Add broadcasting.jl, where broadcasting is implemented  for Taylor1

* Add broadcasting tests (Taylor1 only)

* Correct a method of evaluate, and comment some broken tests

* Add parameterized BroadcastStyle, extend similar, and tests

The extensions for similar (Taylor1 and Vector{Taylor1}) are needed for 
nested Taylor1

* Extend broadcasting to TaylorN and HomogeneousPolynomials

One test of manyvariables.jl was commented

* Avoid broadcasting in resize_coeffs1! and resize_coeffsHP!

* Fix similar for nested Taylor1

* Redesign broadcast

BroadcastStyle's are subtypes of  AbstractArrayStyle{0}

* Fix broadcasting assignment

* Further fixes and cleanup

Inner mutating functions are not extended to be broadcasted

* Use zero in similar to avoid undef's (with mixtures)

And fix the vector case of constant_term and linear_polynomial

* Fixes, cleanup, and tests for broadcasting with mixtures

* No need to overload similar

* Add broadcasted extension for Float32 with tests
  • Loading branch information
lbenet authored and dpsanders committed Jun 3, 2019
1 parent 5d04aba commit 106eb96
Show file tree
Hide file tree
Showing 9 changed files with 381 additions and 60 deletions.
3 changes: 2 additions & 1 deletion src/TaylorSeries.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand Down Expand Up @@ -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__()
Expand Down
70 changes: 45 additions & 25 deletions src/auxiliary.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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

Expand Down Expand Up @@ -209,31 +213,30 @@ function setindex!(a::TaylorN{T}, x::Array{T,1}, u::StepRange{Int,Int}) where {T
end


## 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)
Expand Down Expand Up @@ -266,6 +269,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)
Expand All @@ -277,7 +297,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

Expand All @@ -292,6 +312,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
141 changes: 141 additions & 0 deletions src/broadcasting.jl
Original file line number Diff line number Diff line change
@@ -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
2 changes: 1 addition & 1 deletion src/evaluate.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
2 changes: 1 addition & 1 deletion src/other_functions.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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) )
Expand Down
Loading

0 comments on commit 106eb96

Please sign in to comment.