From 226d5899f66a86a34db7fbc22a488749e1bec9fe Mon Sep 17 00:00:00 2001 From: Luis Benet Date: Sat, 20 Apr 2019 11:11:46 -0500 Subject: [PATCH 01/12] Extend axes and enforce inlining --- src/TaylorSeries.jl | 2 +- src/auxiliary.jl | 30 ++++++++++++++++-------------- 2 files changed, 17 insertions(+), 15 deletions(-) diff --git a/src/TaylorSeries.jl b/src/TaylorSeries.jl index f092a1a9..760d03e4 100644 --- a/src/TaylorSeries.jl +++ b/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 import Base: zero, one, zeros, ones, isinf, isnan, iszero, convert, promote_rule, promote, show, diff --git a/src/auxiliary.jl b/src/auxiliary.jl index 66a8a0c0..e6664947 100644 --- a/src/auxiliary.jl +++ b/src/auxiliary.jl @@ -178,27 +178,29 @@ setindex!(a::TaylorN{T}, x::Array{T,1}, ::Colon) where {T<:Number} = ## eltype, length, get_order ## for T in (:Taylor1, :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) + @inline iterate(a::$T, state=0) = state > a.order ? nothing : (a.coeffs[state+1], state+1) + @inline 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 + @inline eltype(::$T{S}) where {S<:Number} = S + @inline length(a::$T) = length(a.coeffs) + @inline size(a::$T) = (length(a),) + @inline firstindex(a::$T) = 0 + @inline lastindex(a::$T) = a.order + @inline get_order(a::$T) = a.order + @inline axes(a::$T) = axes(a.coeffs) 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 +@inline eltype(::HomogeneousPolynomial{S}) where {S<:Number} = S +@inline length(a::HomogeneousPolynomial) = size_table[a.order+1]#length(a.coeffs) +@inline size(a::HomogeneousPolynomial) = (length(a),) +@inline firstindex(a::HomogeneousPolynomial) = 1 +@inline lastindex(a::HomogeneousPolynomial) = length(a) +@inline get_order(a::HomogeneousPolynomial) = a.order +# @inline axes(a::HomogeneousPolynomial) = axes(a.coeffs) ## fixorder ## From 4a2fd5df4f5eed5d6332ad399a08bd06be38c359 Mon Sep 17 00:00:00 2001 From: Luis Benet Date: Sat, 20 Apr 2019 11:34:17 -0500 Subject: [PATCH 02/12] Add broadcasting.jl, where broadcasting is implemented for Taylor1 --- src/TaylorSeries.jl | 1 + src/broadcasting.jl | 51 +++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 52 insertions(+) create mode 100644 src/broadcasting.jl diff --git a/src/TaylorSeries.jl b/src/TaylorSeries.jl index 760d03e4..0f95bf52 100644 --- a/src/TaylorSeries.jl +++ b/src/TaylorSeries.jl @@ -60,6 +60,7 @@ include("hash_tables.jl") include("constructors.jl") include("conversion.jl") include("auxiliary.jl") +include("broadcasting.jl") include("arithmetic.jl") include("power.jl") include("functions.jl") diff --git a/src/broadcasting.jl b/src/broadcasting.jl new file mode 100644 index 00000000..1e6d63be --- /dev/null +++ b/src/broadcasting.jl @@ -0,0 +1,51 @@ +# This file is part of the TaylorSeries.jl Julia package, MIT license +# +# Luis Benet & David P. Sanders +# UNAM +# +# MIT Expat license +# + +# Broadcasting: Implemented for Taylor1 and TaylorN +import .Broadcast: BroadcastStyle, Style, AbstractArrayStyle, Broadcasted + +BroadcastStyle(::Type{<:Taylor1}) = Style{Taylor1}() +BroadcastStyle(::Style{Taylor1}, ::AbstractArrayStyle{0}) = Style{Taylor1}() +BroadcastStyle(::Style{Taylor1}, ::AbstractArrayStyle{1}) = Style{Taylor1}() + +# 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(::Any, rest) = find_taylor(rest) + +# Extend Base.similar +@eval function Base.similar(bc::Broadcasted{Style{Taylor1}}, ::Type{T}) where {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 + + +# Adapted from Base.Broadcast.copyto!, base/broadcasting.jl, line 832 +@inline function Base.copyto!(dest::Taylor1, bc::Broadcasted) + 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{AbstractArray} # 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.coeffs, bc) + # I is the coefficients index + @simd for I in eachindex(bc′) + @inbounds dest[I-1] = getcoeff(bc′[I], I-1) + end + return dest +end From d621036dfc3143d62284427f20b008b70b7a9a7d Mon Sep 17 00:00:00 2001 From: Luis Benet Date: Sat, 20 Apr 2019 11:54:52 -0500 Subject: [PATCH 03/12] Add broadcasting tests (Taylor1 only) --- test/broadcasting.jl | 18 ++++++++++++++++++ 1 file changed, 18 insertions(+) create mode 100644 test/broadcasting.jl diff --git a/test/broadcasting.jl b/test/broadcasting.jl new file mode 100644 index 00000000..1231cc3e --- /dev/null +++ b/test/broadcasting.jl @@ -0,0 +1,18 @@ +# This file is part of TaylorSeries.jl, MIT licensed +# + +using TaylorSeries + +using Test + +@testset "Broadcasting for Taylor1 expansions" begin + t = Taylor1(Int, 5) + + @test 1.0 .+ t == 1.0 + t + @test typeof(1.0 .+ t) == Taylor1{Float64} + @test 1.0 .+ [t] == [1.0 + t] + @test typeof(1.0 .+ [t]) == Vector{Taylor1{Float64}} + @test 1.0 .+ [t, 2t] == [1.0 + t, 1.0 + 2t] + @test [1.0,2.0] .+ [t 2t] == [1.0+t 1.0+2t; 2.0+t 2.0+2t] + @test [1.0] .+ t == 1.0 + t +end From d3506db1bdf0c746890ba98220119c09ae6c89eb Mon Sep 17 00:00:00 2001 From: Luis Benet Date: Sat, 20 Apr 2019 15:47:29 -0500 Subject: [PATCH 04/12] Correct a method of evaluate, and comment some broken tests --- src/broadcasting.jl | 2 +- src/evaluate.jl | 2 +- test/manyvariables.jl | 16 ++++++++-------- test/onevariable.jl | 28 ++++++++++++++-------------- 4 files changed, 24 insertions(+), 24 deletions(-) diff --git a/src/broadcasting.jl b/src/broadcasting.jl index 1e6d63be..ca0a3d2f 100644 --- a/src/broadcasting.jl +++ b/src/broadcasting.jl @@ -22,7 +22,7 @@ find_taylor(a::Taylor1, rest) = a find_taylor(::Any, rest) = find_taylor(rest) # Extend Base.similar -@eval function Base.similar(bc::Broadcasted{Style{Taylor1}}, ::Type{T}) where {T} +function Base.similar(bc::Broadcasted{Style{Taylor1}}, ::Type{T}) where {T} # Proper promotion R = Base.Broadcast.combine_eltypes(bc.f, bc.args) # Scan the inputs for the Taylor1: diff --git a/src/evaluate.jl b/src/evaluate.jl index 08973fb0..d72d1565 100644 --- a/src/evaluate.jl +++ b/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) diff --git a/test/manyvariables.jl b/test/manyvariables.jl index ea787ced..f0789f3c 100644 --- a/test/manyvariables.jl +++ b/test/manyvariables.jl @@ -635,16 +635,16 @@ eeuler = Base.MathConstants.e q = zero(p) TaylorSeries.deg2rad!(q, p, 0) @test q[0] == p[0]*(pi/180) - TaylorSeries.deg2rad!.(q, p, [1,3,5]) - for i in [0,1,3,5] - @test q[i] == p[i]*(pi/180) - end + # TaylorSeries.deg2rad!.(q, p, [1,3,5]) + # for i in [0,1,3,5] + # @test q[i] == p[i]*(pi/180) + # end TaylorSeries.rad2deg!(q, p, 0) @test q[0] == p[0]*(180/pi) - TaylorSeries.rad2deg!.(q, p, [1,3,5]) - for i in [0,1,3,5] - @test q[i] == p[i]*(180/pi) - end + # TaylorSeries.rad2deg!.(q, p, [1,3,5]) + # for i in [0,1,3,5] + # @test q[i] == p[i]*(180/pi) + # end xT = 5+TaylorN(Int, 1, order=10) yT = TaylorN(2, order=10) TaylorSeries.deg2rad!(yT, xT, 0) diff --git a/test/onevariable.jl b/test/onevariable.jl index 32b064d4..e9e9762c 100644 --- a/test/onevariable.jl +++ b/test/onevariable.jl @@ -246,15 +246,15 @@ eeuler = Base.MathConstants.e @test a.() == [p(), q()] @test a.() == a() vr = rand(2) - @test p.(vr) == evaluate.(p, vr) + @test p.(vr) == evaluate.([p], vr) Mr = rand(3,3,3) - @test p.(Mr) == evaluate.(p, Mr) + @test p.(Mr) == evaluate.([p], Mr) mytaylor1 = Taylor1(rand(20)) vr = rand(5) @test p(vr) == p.(vr) - @test p(vr) == evaluate.(p,vr) + @test p(vr) == evaluate.([p],vr) @test p(Mr) == p.(Mr) - @test p(Mr) == evaluate.(p,Mr) + @test p(Mr) == evaluate.([p], Mr) taylor_a = Taylor1(Int,10) taylor_x = exp(Taylor1(Float64,13)) @test taylor_x(taylor_a) == evaluate(taylor_x, taylor_a) @@ -473,22 +473,22 @@ eeuler = Base.MathConstants.e TaylorSeries.deg2rad!(b, a, 0) @test a == c @test a[0]*(pi/180) == b[0] - TaylorSeries.deg2rad!.(b, a, [0,1,2]) - @test a == c - for i in 0:2 - @test a[i]*(pi/180) == b[i] - end + # TaylorSeries.deg2rad!.(b, a, [0,1,2]) + # @test a == c + # for i in 0:2 + # @test a[i]*(pi/180) == b[i] + # end a = Taylor1(rand(10)) b = Taylor1(rand(10)) c = deepcopy(a) TaylorSeries.rad2deg!(b, a, 0) @test a == c @test a[0]*(180/pi) == b[0] - TaylorSeries.rad2deg!.(b, a, [0,1,2]) - @test a == c - for i in 0:2 - @test a[i]*(180/pi) == b[i] - end + # TaylorSeries.rad2deg!.(b, a, [0,1,2]) + # @test a == c + # for i in 0:2 + # @test a[i]*(180/pi) == b[i] + # end # Test additional Taylor1 constructors @test Taylor1{Float64}(true) == Taylor1([1.0]) From e732ed99d6b4b01788aa4ab183e519a747856473 Mon Sep 17 00:00:00 2001 From: Luis Benet Date: Sun, 21 Apr 2019 12:32:10 -0500 Subject: [PATCH 05/12] Add parameterized BroadcastStyle, extend similar, and tests The extensions for similar (Taylor1 and Vector{Taylor1}) are needed for nested Taylor1 --- src/TaylorSeries.jl | 2 +- src/auxiliary.jl | 18 +++++++++++++++--- src/broadcasting.jl | 29 +++++++++++++++++++++++------ test/broadcasting.jl | 27 ++++++++++++++++++++++++++- test/onevariable.jl | 25 ++++++++++--------------- test/runtests.jl | 5 +++-- 6 files changed, 78 insertions(+), 28 deletions(-) diff --git a/src/TaylorSeries.jl b/src/TaylorSeries.jl index 0f95bf52..7ed639ed 100644 --- a/src/TaylorSeries.jl +++ b/src/TaylorSeries.jl @@ -32,7 +32,7 @@ end import Base: ==, +, -, *, /, ^ import Base: iterate, size, eachindex, firstindex, lastindex, - eltype, length, getindex, setindex!, axes + eltype, length, getindex, setindex!, axes, similar import Base: zero, one, zeros, ones, isinf, isnan, iszero, convert, promote_rule, promote, show, diff --git a/src/auxiliary.jl b/src/auxiliary.jl index e6664947..449bb875 100644 --- a/src/auxiliary.jl +++ b/src/auxiliary.jl @@ -19,7 +19,7 @@ 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 + coeffs[lencoef+1:order+1] .= z end return nothing end @@ -38,7 +38,7 @@ 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 + coeffs[lencoef+1:num_coeffs] .= z return nothing end @@ -200,7 +200,7 @@ end @inline firstindex(a::HomogeneousPolynomial) = 1 @inline lastindex(a::HomogeneousPolynomial) = length(a) @inline get_order(a::HomogeneousPolynomial) = a.order -# @inline axes(a::HomogeneousPolynomial) = axes(a.coeffs) +@inline axes(a::HomogeneousPolynomial) = axes(a.coeffs) ## fixorder ## @@ -234,6 +234,18 @@ function Base.findlast(a::Taylor1{T}) where {T<:Number} return last-1 end + +## similar ## +similar(a::Taylor1) = Taylor1(similar(a.coeffs), a.order) +function similar(a::Array{Taylor1{T},1}) where {T} + ret = Vector{Taylor1{T}}(undef, size(a,1)) + a1 = a[1].coeffs + fill!(ret, similar(a1)) + return ret +end + + + """ constant_term(a) diff --git a/src/broadcasting.jl b/src/broadcasting.jl index ca0a3d2f..5aaec75b 100644 --- a/src/broadcasting.jl +++ b/src/broadcasting.jl @@ -7,11 +7,17 @@ # # Broadcasting: Implemented for Taylor1 and TaylorN -import .Broadcast: BroadcastStyle, Style, AbstractArrayStyle, Broadcasted +import .Broadcast: BroadcastStyle, Style, AbstractArrayStyle, Broadcasted, broadcasted + +# BroadcastStyle definitions +BroadcastStyle(::Type{<:Taylor1{T}}) where {T}= Style{Taylor1{T}}() +BroadcastStyle(::Style{Taylor1{T}}, ::AbstractArrayStyle{0}) where {T}= Style{Taylor1{T}}() +BroadcastStyle(::Style{Taylor1{T}}, ::AbstractArrayStyle{1}) where {T}= Style{Taylor1{T}}() + +# Precedence rules (for mixtures) +BroadcastStyle(::Style{Taylor1{Taylor1{T}}}, ::Style{Taylor1{S}}) where + {T,S<:NumberNotSeries} = Style{Taylor1{Taylor1{T}}}() -BroadcastStyle(::Type{<:Taylor1}) = Style{Taylor1}() -BroadcastStyle(::Style{Taylor1}, ::AbstractArrayStyle{0}) = Style{Taylor1}() -BroadcastStyle(::Style{Taylor1}, ::AbstractArrayStyle{1}) = Style{Taylor1}() # 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." @@ -22,7 +28,7 @@ find_taylor(a::Taylor1, rest) = a find_taylor(::Any, rest) = find_taylor(rest) # Extend Base.similar -function Base.similar(bc::Broadcasted{Style{Taylor1}}, ::Type{T}) where {T} +function Base.similar(bc::Broadcasted{Style{Taylor1{S}}}, ::Type{T}) where {S, T} # Proper promotion R = Base.Broadcast.combine_eltypes(bc.f, bc.args) # Scan the inputs for the Taylor1: @@ -33,7 +39,7 @@ end # Adapted from Base.Broadcast.copyto!, base/broadcasting.jl, line 832 -@inline function Base.copyto!(dest::Taylor1, bc::Broadcasted) +@inline function Base.copyto!(dest::Taylor1{T}, bc::Broadcasted) where {T<:NumberNotSeries} 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{AbstractArray} # only a single input argument to broadcast! @@ -49,3 +55,14 @@ end end return dest end + + +# Broadcasted extensions +# This should prevent broadcasting being applied in `a` and `b` +# for the mutating functions, and work only in `k` +function broadcasted(::Style{Taylor1{T}}, f!, a::Taylor1{T}, b::Taylor1{T}, k) where {T} + @inbounds for i in eachindex(k) + f!(a, b, k[i]) + end + nothing +end diff --git a/test/broadcasting.jl b/test/broadcasting.jl index 1231cc3e..2b1639b1 100644 --- a/test/broadcasting.jl +++ b/test/broadcasting.jl @@ -14,5 +14,30 @@ using Test @test typeof(1.0 .+ [t]) == Vector{Taylor1{Float64}} @test 1.0 .+ [t, 2t] == [1.0 + t, 1.0 + 2t] @test [1.0,2.0] .+ [t 2t] == [1.0+t 1.0+2t; 2.0+t 2.0+2t] - @test [1.0] .+ t == 1.0 + t + @test [1.0] .+ t == 1.0 + t # This should return a vector !! + @test t .+ t == 2t == 2 .* t + + st = sin(t) + @test st(pi/3) == evaluate(st, pi/3) + @test st.(pi/3) == evaluate(st, pi/3) + @test st.([0.0, pi/3]) == evaluate(st, [0.0, pi/3]) + +end + +@testset "Broadcasting for nested Taylor1's" begin + t = Taylor1(Int, 3) + @test typeof(similar(t)) == Taylor1{Int} + @test get_order(t) == 3 + + tt = Taylor1([zero(t), one(t), zero(t)]) + @test typeof(similar(tt)) == Taylor1{Taylor1{Int}} + @test get_order(tt) == get_order(similar(tt)) == 2 + + ttt = Taylor1([zero(tt), one(tt)]) + @test typeof(similar(ttt)) == Taylor1{Taylor1{Taylor1{Int}}} + @test get_order(ttt) == get_order(similar(ttt)) == 1 + + tf = Taylor1([zero(1.0*t), one(t)]) + @test typeof(similar(tf)) == Taylor1{Taylor1{Float64}} + @test get_order(tf) == get_order(similar(tf)) == 1 end diff --git a/test/onevariable.jl b/test/onevariable.jl index e9e9762c..9b22fe0a 100644 --- a/test/onevariable.jl +++ b/test/onevariable.jl @@ -29,11 +29,6 @@ eeuler = Base.MathConstants.e @test iterate(t, 1) == (1.0, 2) @test iterate(t, 16) == nothing - # ta = Taylor1([1,2,3,4,5]) - # vT = [h for h in ta] - # voT = [h.order for h in Taylor1([1,2,3,4,5])] - # @test vT == [Taylor1(1,0), Taylor1(2,1), Taylor1(3,2), Taylor1(4,3), Taylor1(5,4)] - # @test voT == [0, 1, 2, 3, 4] v = [1,2] @test typeof(TaylorSeries.resize_coeffs1!(v,3)) == Nothing @@ -473,22 +468,22 @@ eeuler = Base.MathConstants.e TaylorSeries.deg2rad!(b, a, 0) @test a == c @test a[0]*(pi/180) == b[0] - # TaylorSeries.deg2rad!.(b, a, [0,1,2]) - # @test a == c - # for i in 0:2 - # @test a[i]*(pi/180) == b[i] - # end + TaylorSeries.deg2rad!.(b, a, [0,1,2]) + @test a == c + for i in 0:2 + @test a[i]*(pi/180) == b[i] + end a = Taylor1(rand(10)) b = Taylor1(rand(10)) c = deepcopy(a) TaylorSeries.rad2deg!(b, a, 0) @test a == c @test a[0]*(180/pi) == b[0] - # TaylorSeries.rad2deg!.(b, a, [0,1,2]) - # @test a == c - # for i in 0:2 - # @test a[i]*(180/pi) == b[i] - # end + TaylorSeries.rad2deg!.(b, a, [0,1,2]) + @test a == c + for i in 0:2 + @test a[i]*(180/pi) == b[i] + end # Test additional Taylor1 constructors @test Taylor1{Float64}(true) == Taylor1([1.0]) diff --git a/test/runtests.jl b/test/runtests.jl index 4abbb861..0d051566 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -7,9 +7,10 @@ testfiles = ( "manyvariables.jl", "mixtures.jl", "mutatingfuncts.jl", + "intervals.jl", + "broadcasting.jl", "identities_Euler.jl", - "fateman40.jl", - "intervals.jl" + "fateman40.jl" ) for file in testfiles From e200f199a564db550097c8e72f11a3a20dbe9b24 Mon Sep 17 00:00:00 2001 From: Luis Benet Date: Sun, 21 Apr 2019 17:57:53 -0500 Subject: [PATCH 06/12] Extend broadcasting to TaylorN and HomogeneousPolynomials One test of manyvariables.jl was commented --- src/auxiliary.jl | 26 ++++++++++++- src/broadcasting.jl | 83 +++++++++++++++++++++++++++++++++++++++++- src/other_functions.jl | 2 +- test/manyvariables.jl | 36 +++++++++--------- 4 files changed, 125 insertions(+), 22 deletions(-) diff --git a/src/auxiliary.jl b/src/auxiliary.jl index 449bb875..afb2d660 100644 --- a/src/auxiliary.jl +++ b/src/auxiliary.jl @@ -183,7 +183,7 @@ for T in (:Taylor1, :TaylorN) # Base.iterate(rS::Iterators.Reverse{$T}, state=rS.itr.order) = state < 0 ? nothing : (a.coeffs[state], state-1) @inline eltype(::$T{S}) where {S<:Number} = S @inline length(a::$T) = length(a.coeffs) - @inline size(a::$T) = (length(a),) + @inline size(a::$T) = size(a.coeffs) @inline firstindex(a::$T) = 0 @inline lastindex(a::$T) = a.order @inline get_order(a::$T) = a.order @@ -196,7 +196,7 @@ end # Base.eachindex(a::HomogeneousPolynomial) = firstindex(a):lastindex(a) @inline eltype(::HomogeneousPolynomial{S}) where {S<:Number} = S @inline length(a::HomogeneousPolynomial) = size_table[a.order+1]#length(a.coeffs) -@inline size(a::HomogeneousPolynomial) = (length(a),) +@inline size(a::HomogeneousPolynomial) = size(a.coeffs) @inline firstindex(a::HomogeneousPolynomial) = 1 @inline lastindex(a::HomogeneousPolynomial) = length(a) @inline get_order(a::HomogeneousPolynomial) = a.order @@ -243,6 +243,28 @@ function similar(a::Array{Taylor1{T},1}) where {T} fill!(ret, similar(a1)) return ret end +similar(a::Array{Taylor1{T},1}, R::Type) where {T} = + convert(Array{Taylor1{R},1}, similar(a)) + +similar(a::HomogeneousPolynomial) = HomogeneousPolynomial(similar(a.coeffs), a.order) +function similar(a::Array{HomogeneousPolynomial{T},1}) where {T} + ret = Vector{HomogeneousPolynomial{T}}(undef, size(a,1)) + @simd for i in eachindex(a) + @inbounds ret[i] = similar(a[i]) + end + return ret +end +similar(a::Array{HomogeneousPolynomial{T},1}, R::Type{<:NumberNotSeriesN}) where {T} = + convert(Array{HomogeneousPolynomial{R},1}, similar(a)) + +similar(a::TaylorN) = TaylorN(similar(a.coeffs), a.order) +function similar(a::Array{TaylorN{T},1}) where {T} + ret = Vector{TaylorN{T}}(undef, size(a,1)) + @simd for i in eachindex(a) + @inbounds ret[i] = similar(a[i]) + end + return ret +end diff --git a/src/broadcasting.jl b/src/broadcasting.jl index 5aaec75b..d673c0aa 100644 --- a/src/broadcasting.jl +++ b/src/broadcasting.jl @@ -13,10 +13,27 @@ import .Broadcast: BroadcastStyle, Style, AbstractArrayStyle, Broadcasted, broad BroadcastStyle(::Type{<:Taylor1{T}}) where {T}= Style{Taylor1{T}}() BroadcastStyle(::Style{Taylor1{T}}, ::AbstractArrayStyle{0}) where {T}= Style{Taylor1{T}}() BroadcastStyle(::Style{Taylor1{T}}, ::AbstractArrayStyle{1}) where {T}= Style{Taylor1{T}}() +# +BroadcastStyle(::Type{<:HomogeneousPolynomial{T}}) where {T} = Style{HomogeneousPolynomial{T}}() +BroadcastStyle(::Style{HomogeneousPolynomial{T}}, ::AbstractArrayStyle{0}) where {T} = Style{HomogeneousPolynomial{T}}() +BroadcastStyle(::Style{HomogeneousPolynomial{T}}, ::AbstractArrayStyle{1}) where {T} = Style{HomogeneousPolynomial{T}}() +# +BroadcastStyle(::Type{<:TaylorN{T}}) where {T} = Style{TaylorN{T}}() +BroadcastStyle(::Style{TaylorN{T}}, ::AbstractArrayStyle{0}) where {T} = Style{TaylorN{T}}() +BroadcastStyle(::Style{TaylorN{T}}, ::AbstractArrayStyle{1}) where {T} = Style{TaylorN{T}}() # Precedence rules (for mixtures) BroadcastStyle(::Style{Taylor1{Taylor1{T}}}, ::Style{Taylor1{S}}) where - {T,S<:NumberNotSeries} = Style{Taylor1{Taylor1{T}}}() + {T,S} = Style{Taylor1{Taylor1{T}}}() +# +BroadcastStyle(::Style{Taylor1{TaylorN{T}}}, ::Style{HomogeneousPolynomial{S}}) where + {T<:NumberNotSeries, S<:NumberNotSeries} = Style{Taylor1{TaylorN{promote_type(T,S)}}}() +BroadcastStyle(::Style{Taylor1{TaylorN{T}}}, ::Style{TaylorN{S}}) where + {T<:NumberNotSeries, S<:NumberNotSeries} = Style{Taylor1{TaylorN{promote_type(T,S)}}}() +BroadcastStyle(::Style{Taylor1{T}}, ::Style{HomogeneousPolynomial{Taylor1{S}}}) where + {T<:NumberNotSeries, S<:NumberNotSeries} = Style{HomogeneousPolynomial{Taylor1{promote_type(T,S)}}}() +BroadcastStyle(::Style{Taylor1{T}}, ::Style{TaylorN{Taylor1{S}}}) where + {T<:NumberNotSeries, S<:NumberNotSeries} = Style{TaylorN{Taylor1{promote_type(T,S)}}}() # We follow https://docs.julialang.org/en/v1/manual/interfaces/#man-interface-iteration-1 @@ -25,6 +42,8 @@ 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(::Any, rest) = find_taylor(rest) # Extend Base.similar @@ -37,6 +56,28 @@ function Base.similar(bc::Broadcasted{Style{Taylor1{S}}}, ::Type{T}) where {S, T return Taylor1(similar(A.coeffs, R), A.order) end +function Base.similar(bc::Broadcasted{Style{HomogeneousPolynomial{S}}}, ::Type{T}) where {S, T} + # Proper promotion + @show(Base.Broadcast.eltypes(bc.args)) + # combine_eltypes(f, args::Tuple) = Base._return_type(f, eltypes(args)) + R = Base.Broadcast.combine_eltypes(bc.f, bc.args) + @show(bc.f, bc.args) + # Scan the inputs for the HomogeneousPolynomial: + A = find_taylor(bc) + @show(A, R) + # Create the output + return HomogeneousPolynomial(similar(A.coeffs, R), A.order) +end + +function Base.similar(bc::Broadcasted{Style{TaylorN{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 @inline function Base.copyto!(dest::Taylor1{T}, bc::Broadcasted) where {T<:NumberNotSeries} @@ -56,6 +97,46 @@ end return dest end +@inline function Base.copyto!(dest::HomogeneousPolynomial{T}, bc::Broadcasted) where {T<:NumberNotSeries} + @show(dest, get_order(dest)) + 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{AbstractArray} # 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.coeffs, bc) + @show(bc′) + # I is the coefficients index + @simd for I in eachindex(bc′) + # @show(I, bc′[I].coeffs, bc′[I].order) + # @inbounds dest[I] = getcoeff(bc′[I], I) + @inbounds dest[I] = bc′[I].coeffs[I] + end + return dest +end + +@inline function Base.copyto!(dest::TaylorN{T}, bc::Broadcasted) where {T<:NumberNotSeries} + @show(dest, bc) + 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{AbstractArray} # 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.coeffs, bc) + # I is the coefficients index + @simd for I in eachindex(bc′) + @show(I, bc´[I]) + @inbounds dest[I] = getcoeff(bc′[I], I) + end + return dest +end + # Broadcasted extensions # This should prevent broadcasting being applied in `a` and `b` diff --git a/src/other_functions.jl b/src/other_functions.jl index 514e3ca7..e34acbb9 100644 --- a/src/other_functions.jl +++ b/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) ) diff --git a/test/manyvariables.jl b/test/manyvariables.jl index f0789f3c..72dcc7c1 100644 --- a/test/manyvariables.jl +++ b/test/manyvariables.jl @@ -569,7 +569,7 @@ eeuler = Base.MathConstants.e @test a[2] ≈ a[2] @test a[3] ≈ a[3] @test a ≈ a - @test a .≈ a + # @test a .≈ a b = deepcopy(a) b[2][3] = Inf @test !isfinite(b) @@ -626,31 +626,31 @@ eeuler = Base.MathConstants.e @test Q() == evaluate(Q) @test Q[1:3]() == evaluate(Q[1:3]) - dx = set_variables("x", numvars=4, order=10) - for i in 1:4 - @test deg2rad(180+dx[i]) == pi + deg2rad(1.0)dx[i] - @test rad2deg(pi+dx[i]) == 180.0+rad2deg(1.0)dx[i] - end - p = sin(exp(dx[1]*dx[2])+dx[3]*dx[2])/(1.0+dx[4]^2) - q = zero(p) - TaylorSeries.deg2rad!(q, p, 0) - @test q[0] == p[0]*(pi/180) + # dx = set_variables("x", numvars=4, order=10) + # for i in 1:4 + # @test deg2rad(180+dx[i]) == pi + deg2rad(1.0)dx[i] + # @test rad2deg(pi+dx[i]) == 180.0+rad2deg(1.0)dx[i] + # end + # p = sin(exp(dx[1]*dx[2])+dx[3]*dx[2])/(1.0+dx[4]^2) + # q = zero(p) + # TaylorSeries.deg2rad!(q, p, 0) + # @test q[0] == p[0]*(pi/180) # TaylorSeries.deg2rad!.(q, p, [1,3,5]) # for i in [0,1,3,5] # @test q[i] == p[i]*(pi/180) # end - TaylorSeries.rad2deg!(q, p, 0) - @test q[0] == p[0]*(180/pi) + # TaylorSeries.rad2deg!(q, p, 0) + # @test q[0] == p[0]*(180/pi) # TaylorSeries.rad2deg!.(q, p, [1,3,5]) # for i in [0,1,3,5] # @test q[i] == p[i]*(180/pi) # end - xT = 5+TaylorN(Int, 1, order=10) - yT = TaylorN(2, order=10) - TaylorSeries.deg2rad!(yT, xT, 0) - @test yT[0] == xT[0]*(pi/180) - TaylorSeries.rad2deg!(yT, xT, 0) - @test yT[0] == xT[0]*(180/pi) + # xT = 5+TaylorN(Int, 1, order=10) + # yT = TaylorN(2, order=10) + # TaylorSeries.deg2rad!(yT, xT, 0) + # @test yT[0] == xT[0]*(pi/180) + # TaylorSeries.rad2deg!(yT, xT, 0) + # @test yT[0] == xT[0]*(180/pi) end @testset "Integrate for several variables" begin From 3c720772b8e201ee20b13a31097ecc1fb6b5fc55 Mon Sep 17 00:00:00 2001 From: Luis Benet Date: Mon, 22 Apr 2019 18:09:12 -0500 Subject: [PATCH 07/12] Avoid broadcasting in resize_coeffs1! and resize_coeffsHP! --- src/auxiliary.jl | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/src/auxiliary.jl b/src/auxiliary.jl index afb2d660..51d8394f 100644 --- a/src/auxiliary.jl +++ b/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]) - 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]) - coeffs[lencoef+1:num_coeffs] .= z + @simd for ord in lencoef+1:num_coeffs + @inbounds coeffs[ord] = z + end return nothing end From 80be007b515b6f726068f00c475d5b745fae060c Mon Sep 17 00:00:00 2001 From: Luis Benet Date: Tue, 23 Apr 2019 15:22:29 -0500 Subject: [PATCH 08/12] Fix similar for nested Taylor1 --- src/auxiliary.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/auxiliary.jl b/src/auxiliary.jl index 51d8394f..e3fcae73 100644 --- a/src/auxiliary.jl +++ b/src/auxiliary.jl @@ -248,7 +248,7 @@ function similar(a::Array{Taylor1{T},1}) where {T} return ret end similar(a::Array{Taylor1{T},1}, R::Type) where {T} = - convert(Array{Taylor1{R},1}, similar(a)) + convert(Vector{promote_type(Taylor1{T},R)}, similar(a)) similar(a::HomogeneousPolynomial) = HomogeneousPolynomial(similar(a.coeffs), a.order) function similar(a::Array{HomogeneousPolynomial{T},1}) where {T} From 8fd61bda283ec021f0df5055b19bcdf87a21afec Mon Sep 17 00:00:00 2001 From: Luis Benet Date: Wed, 24 Apr 2019 13:40:27 -0500 Subject: [PATCH 09/12] Redesign broadcast BroadcastStyle's are subtypes of AbstractArrayStyle{0} --- src/TaylorSeries.jl | 2 +- src/auxiliary.jl | 64 +++++++-------- src/broadcasting.jl | 177 +++++++++++++++++++++--------------------- test/broadcasting.jl | 65 +++++++++++++--- test/manyvariables.jl | 2 +- 5 files changed, 181 insertions(+), 129 deletions(-) diff --git a/src/TaylorSeries.jl b/src/TaylorSeries.jl index 7ed639ed..b8c9972c 100644 --- a/src/TaylorSeries.jl +++ b/src/TaylorSeries.jl @@ -32,7 +32,7 @@ end import Base: ==, +, -, *, /, ^ import Base: iterate, size, eachindex, firstindex, lastindex, - eltype, length, getindex, setindex!, axes, similar + eltype, length, getindex, setindex!, axes#, similar import Base: zero, one, zeros, ones, isinf, isnan, iszero, convert, promote_rule, promote, show, diff --git a/src/auxiliary.jl b/src/auxiliary.jl index e3fcae73..70632b90 100644 --- a/src/auxiliary.jl +++ b/src/auxiliary.jl @@ -191,7 +191,7 @@ for T in (:Taylor1, :TaylorN) @inline firstindex(a::$T) = 0 @inline lastindex(a::$T) = a.order @inline get_order(a::$T) = a.order - @inline axes(a::$T) = axes(a.coeffs) + @inline axes(a::$T) = () end end @@ -204,7 +204,7 @@ end @inline firstindex(a::HomogeneousPolynomial) = 1 @inline lastindex(a::HomogeneousPolynomial) = length(a) @inline get_order(a::HomogeneousPolynomial) = a.order -@inline axes(a::HomogeneousPolynomial) = axes(a.coeffs) +@inline axes(a::HomogeneousPolynomial) = () ## fixorder ## @@ -240,35 +240,37 @@ end ## similar ## -similar(a::Taylor1) = Taylor1(similar(a.coeffs), a.order) -function similar(a::Array{Taylor1{T},1}) where {T} - ret = Vector{Taylor1{T}}(undef, size(a,1)) - a1 = a[1].coeffs - fill!(ret, similar(a1)) - return ret -end -similar(a::Array{Taylor1{T},1}, R::Type) where {T} = - convert(Vector{promote_type(Taylor1{T},R)}, similar(a)) - -similar(a::HomogeneousPolynomial) = HomogeneousPolynomial(similar(a.coeffs), a.order) -function similar(a::Array{HomogeneousPolynomial{T},1}) where {T} - ret = Vector{HomogeneousPolynomial{T}}(undef, size(a,1)) - @simd for i in eachindex(a) - @inbounds ret[i] = similar(a[i]) - end - return ret -end -similar(a::Array{HomogeneousPolynomial{T},1}, R::Type{<:NumberNotSeriesN}) where {T} = - convert(Array{HomogeneousPolynomial{R},1}, similar(a)) - -similar(a::TaylorN) = TaylorN(similar(a.coeffs), a.order) -function similar(a::Array{TaylorN{T},1}) where {T} - ret = Vector{TaylorN{T}}(undef, size(a,1)) - @simd for i in eachindex(a) - @inbounds ret[i] = similar(a[i]) - end - return ret -end +# similar(a::Taylor1) = Taylor1(similar(a.coeffs), a.order) +# function similar(a::Array{Taylor1{T},1}) where {T} +# ret = Vector{Taylor1{T}}(undef, size(a,1)) +# a1 = a[1].coeffs +# fill!(ret, similar(a1)) +# return ret +# end +# similar(a::Array{Taylor1{T},1}, R::Type) where {T} = +# convert(Vector{promote_type(Taylor1{T},R)}, similar(a)) + +# similar(a::HomogeneousPolynomial) = HomogeneousPolynomial(similar(a.coeffs), a.order) +# function similar(a::Array{HomogeneousPolynomial{T},1}) where {T} +# ret = Vector{HomogeneousPolynomial{T}}(undef, size(a,1)) +# @simd for i in eachindex(a) +# @inbounds ret[i] = similar(a[i]) +# end +# return ret +# end +# similar(a::Array{HomogeneousPolynomial{T},1}, R::Type{<:NumberNotSeriesN}) where {T} = +# convert(Array{HomogeneousPolynomial{R},1}, similar(a)) + +# similar(a::TaylorN) = TaylorN(similar(a.coeffs), a.order) +# function similar(a::Array{TaylorN{T},1}) where {T} +# ret = Vector{TaylorN{T}}(undef, size(a,1)) +# @simd for i in eachindex(a) +# @inbounds ret[i] = similar(a[i]) +# end +# return ret +# end +# similar(a::Array{taylor_expand{T},1}, R::Type) where {T} = +# convert(Vector{promote_type(TaylorN{T},R)}, similar(a)) diff --git a/src/broadcasting.jl b/src/broadcasting.jl index d673c0aa..0b37b7ab 100644 --- a/src/broadcasting.jl +++ b/src/broadcasting.jl @@ -6,35 +6,43 @@ # MIT Expat license # -# Broadcasting: Implemented for Taylor1 and TaylorN -import .Broadcast: BroadcastStyle, Style, AbstractArrayStyle, Broadcasted, broadcasted +## Broadcast for Taylor1 and TaylorN -# BroadcastStyle definitions -BroadcastStyle(::Type{<:Taylor1{T}}) where {T}= Style{Taylor1{T}}() -BroadcastStyle(::Style{Taylor1{T}}, ::AbstractArrayStyle{0}) where {T}= Style{Taylor1{T}}() -BroadcastStyle(::Style{Taylor1{T}}, ::AbstractArrayStyle{1}) where {T}= Style{Taylor1{T}}() +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}() # -BroadcastStyle(::Type{<:HomogeneousPolynomial{T}}) where {T} = Style{HomogeneousPolynomial{T}}() -BroadcastStyle(::Style{HomogeneousPolynomial{T}}, ::AbstractArrayStyle{0}) where {T} = Style{HomogeneousPolynomial{T}}() -BroadcastStyle(::Style{HomogeneousPolynomial{T}}, ::AbstractArrayStyle{1}) where {T} = Style{HomogeneousPolynomial{T}}() +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}() # -BroadcastStyle(::Type{<:TaylorN{T}}) where {T} = Style{TaylorN{T}}() -BroadcastStyle(::Style{TaylorN{T}}, ::AbstractArrayStyle{0}) where {T} = Style{TaylorN{T}}() -BroadcastStyle(::Style{TaylorN{T}}, ::AbstractArrayStyle{1}) where {T} = Style{TaylorN{T}}() +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(::Style{Taylor1{Taylor1{T}}}, ::Style{Taylor1{S}}) where - {T,S} = Style{Taylor1{Taylor1{T}}}() -# -BroadcastStyle(::Style{Taylor1{TaylorN{T}}}, ::Style{HomogeneousPolynomial{S}}) where - {T<:NumberNotSeries, S<:NumberNotSeries} = Style{Taylor1{TaylorN{promote_type(T,S)}}}() -BroadcastStyle(::Style{Taylor1{TaylorN{T}}}, ::Style{TaylorN{S}}) where - {T<:NumberNotSeries, S<:NumberNotSeries} = Style{Taylor1{TaylorN{promote_type(T,S)}}}() -BroadcastStyle(::Style{Taylor1{T}}, ::Style{HomogeneousPolynomial{Taylor1{S}}}) where - {T<:NumberNotSeries, S<:NumberNotSeries} = Style{HomogeneousPolynomial{Taylor1{promote_type(T,S)}}}() -BroadcastStyle(::Style{Taylor1{T}}, ::Style{TaylorN{Taylor1{S}}}) where - {T<:NumberNotSeries, S<:NumberNotSeries} = Style{TaylorN{Taylor1{promote_type(T,S)}}}() +# 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." @@ -44,10 +52,10 @@ find_taylor(x) = x find_taylor(a::Taylor1, rest) = a find_taylor(a::HomogeneousPolynomial, rest) = a find_taylor(a::TaylorN, rest) = a -find_taylor(::Any, rest) = find_taylor(rest) +find_taylor(::AbstractArray, rest) = find_taylor(rest) # Extend Base.similar -function Base.similar(bc::Broadcasted{Style{Taylor1{S}}}, ::Type{T}) where {S, T} +function Base.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: @@ -56,7 +64,7 @@ function Base.similar(bc::Broadcasted{Style{Taylor1{S}}}, ::Type{T}) where {S, T return Taylor1(similar(A.coeffs, R), A.order) end -function Base.similar(bc::Broadcasted{Style{HomogeneousPolynomial{S}}}, ::Type{T}) where {S, T} +function Base.similar(bc::Broadcasted{HomogeneousPolynomialStyle{S}}, ::Type{T}) where {S, T} # Proper promotion @show(Base.Broadcast.eltypes(bc.args)) # combine_eltypes(f, args::Tuple) = Base._return_type(f, eltypes(args)) @@ -69,7 +77,7 @@ function Base.similar(bc::Broadcasted{Style{HomogeneousPolynomial{S}}}, ::Type{T return HomogeneousPolynomial(similar(A.coeffs, R), A.order) end -function Base.similar(bc::Broadcasted{Style{TaylorN{S}}}, ::Type{T}) where {S, T} +function Base.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: @@ -80,68 +88,63 @@ end # Adapted from Base.Broadcast.copyto!, base/broadcasting.jl, line 832 -@inline function Base.copyto!(dest::Taylor1{T}, bc::Broadcasted) where {T<:NumberNotSeries} - 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{AbstractArray} # 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.coeffs, bc) - # I is the coefficients index - @simd for I in eachindex(bc′) - @inbounds dest[I-1] = getcoeff(bc′[I], I-1) - end - return dest -end - -@inline function Base.copyto!(dest::HomogeneousPolynomial{T}, bc::Broadcasted) where {T<:NumberNotSeries} - @show(dest, get_order(dest)) - 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{AbstractArray} # 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.coeffs, bc) - @show(bc′) - # I is the coefficients index - @simd for I in eachindex(bc′) - # @show(I, bc′[I].coeffs, bc′[I].order) - # @inbounds dest[I] = getcoeff(bc′[I], I) - @inbounds dest[I] = bc′[I].coeffs[I] - end - return dest -end - -@inline function Base.copyto!(dest::TaylorN{T}, bc::Broadcasted) where {T<:NumberNotSeries} - @show(dest, bc) - 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{AbstractArray} # 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.coeffs, bc) - # I is the coefficients index - @simd for I in eachindex(bc′) - @show(I, bc´[I]) - @inbounds dest[I] = getcoeff(bc′[I], I) - end - return dest -end +# @inline function Base.copyto!(dest::Taylor1{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{AbstractArray} # 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) +# # I is the coefficients index +# @simd for I in eachindex(bc′) +# @inbounds dest[I-1] = getcoeff(bc′[I], I-1) +# end +# return dest +# end +# +# @inline function Base.copyto!(dest::HomogeneousPolynomial{T}, bc::Broadcasted) where {T<:NumberNotSeries} +# 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{AbstractArray} # 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.coeffs, bc) +# # I is the coefficients index +# @simd for I in eachindex(bc′) +# # @inbounds dest[I] = getcoeff(bc′[I], I) +# @inbounds dest[I] = bc′[I].coeffs[I] +# end +# return dest +# end +# +# @inline function Base.copyto!(dest::TaylorN{T}, bc::Broadcasted) where {T<:NumberNotSeries} +# 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{AbstractArray} # 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.coeffs, bc) +# # I is the coefficients index +# @simd for I in eachindex(bc′) +# @inbounds dest[I] = getcoeff(bc′[I], I) +# end +# return dest +# end # Broadcasted extensions -# This should prevent broadcasting being applied in `a` and `b` -# for the mutating functions, and work only in `k` -function broadcasted(::Style{Taylor1{T}}, f!, a::Taylor1{T}, b::Taylor1{T}, k) where {T} +# This prevents broadcasting being applied to `a` and `b` +# for the mutating functions, and to act only in `k` +function broadcasted(::Taylor1Style{T}, f!, a::Taylor1{T}, b::Taylor1{T}, k) where {T} @inbounds for i in eachindex(k) f!(a, b, k[i]) end diff --git a/test/broadcasting.jl b/test/broadcasting.jl index 2b1639b1..f769bfcb 100644 --- a/test/broadcasting.jl +++ b/test/broadcasting.jl @@ -8,36 +8,83 @@ using Test @testset "Broadcasting for Taylor1 expansions" begin t = Taylor1(Int, 5) + # @test t .= t + @test t .== t + @test t .≈ t + @test t .!= (1 + t) + @test 1.0 .+ t == 1.0 + t @test typeof(1.0 .+ t) == Taylor1{Float64} @test 1.0 .+ [t] == [1.0 + t] @test typeof(1.0 .+ [t]) == Vector{Taylor1{Float64}} @test 1.0 .+ [t, 2t] == [1.0 + t, 1.0 + 2t] @test [1.0,2.0] .+ [t 2t] == [1.0+t 1.0+2t; 2.0+t 2.0+2t] - @test [1.0] .+ t == 1.0 + t # This should return a vector !! - @test t .+ t == 2t == 2 .* t + @test [1.0] .+ t == [1.0 + t] + @test 1.0 .* t == t + @test typeof(1.0 .* t) == Taylor1{Float64} st = sin(t) + @test st .== st @test st(pi/3) == evaluate(st, pi/3) @test st.(pi/3) == evaluate(st, pi/3) @test st.([0.0, pi/3]) == evaluate(st, [0.0, pi/3]) +end + +@testset "Broadcasting for HomogeneousPolynomial and TaylorN" begin + x, y = set_variables("x y", order=3) + xH = x[1] + yH = y[1] + + @test xH .== xH + @test yH .≈ yH + @test x[2] .== y[2] + + @test x .== x + @test y .≈ y + @test x .!= (1 + x) + + @test 1.0 .+ x == 1.0 + x + @test y .+ x == y + x + @test typeof(big"1.0" .+ x) == TaylorN{BigFloat} + @test 1.0 .+ [x] == [1.0 + x] + @test typeof(1.0 .+ [y]) == Vector{TaylorN{Float64}} + @test 1.0 .+ [x, 2y] == [1.0 + x, 1.0 + 2y] + @test [1.0,2.0] .+ [x 2y] == [1.0+x 1.0+2y; 2.0+x 2.0+2y] + @test [1.0] .+ x == [1.0 + x] + @test 1.0 .* y == y + @test typeof(1.0 .* x .* y) == TaylorN{Float64} + # + # @test 1.0 .+ t == 1.0 + t + # @test typeof(1.0 .+ t) == Taylor1{Float64} + # @test 1.0 .+ [t] == [1.0 + t] + # @test typeof(1.0 .+ [t]) == Vector{Taylor1{Float64}} + # @test 1.0 .+ [t, 2t] == [1.0 + t, 1.0 + 2t] + # @test [1.0,2.0] .+ [t 2t] == [1.0+t 1.0+2t; 2.0+t 2.0+2t] + # @test [1.0] .+ t == [1.0 + t] + # @test 1.0 .* t == t + # @test typeof(1.0 .* t) == Taylor1{Float64} + # xH = HomogeneousPolynomial([one(t), zero(t)]) + # @assert typeof(similar(xHt1)) == HomogeneousPolynomial{Taylor1{Int64}} end @testset "Broadcasting for nested Taylor1's" begin t = Taylor1(Int, 3) - @test typeof(similar(t)) == Taylor1{Int} + # @test typeof(similar(t)) == Taylor1{Int} @test get_order(t) == 3 tt = Taylor1([zero(t), one(t), zero(t)]) - @test typeof(similar(tt)) == Taylor1{Taylor1{Int}} - @test get_order(tt) == get_order(similar(tt)) == 2 + # @test typeof(similar(tt)) == Taylor1{Taylor1{Int}} + # @test get_order(tt) == get_order(similar(tt)) == 2 + @test tt .== tt ttt = Taylor1([zero(tt), one(tt)]) - @test typeof(similar(ttt)) == Taylor1{Taylor1{Taylor1{Int}}} - @test get_order(ttt) == get_order(similar(ttt)) == 1 + # @test typeof(similar(ttt)) == Taylor1{Taylor1{Taylor1{Int}}} + # @test get_order(ttt) == get_order(similar(ttt)) == 1 + @test ttt .≈ ttt tf = Taylor1([zero(1.0*t), one(t)]) - @test typeof(similar(tf)) == Taylor1{Taylor1{Float64}} - @test get_order(tf) == get_order(similar(tf)) == 1 + # @test typeof(similar(tf)) == Taylor1{Taylor1{Float64}} + # @test get_order(tf) == get_order(similar(tf)) == 1 + @test tf .== tt end diff --git a/test/manyvariables.jl b/test/manyvariables.jl index 72dcc7c1..5ed5870f 100644 --- a/test/manyvariables.jl +++ b/test/manyvariables.jl @@ -569,7 +569,7 @@ eeuler = Base.MathConstants.e @test a[2] ≈ a[2] @test a[3] ≈ a[3] @test a ≈ a - # @test a .≈ a + @test a .≈ a b = deepcopy(a) b[2][3] = Inf @test !isfinite(b) From 0c89c2309c91c502a0af0c540610db4734517c89 Mon Sep 17 00:00:00 2001 From: Luis Benet Date: Thu, 25 Apr 2019 16:28:56 -0500 Subject: [PATCH 10/12] Fix broadcasting assignment --- src/TaylorSeries.jl | 2 +- src/auxiliary.jl | 43 ++++++++++++++++++++++------ src/broadcasting.jl | 66 +++++++++++++++++++++---------------------- test/broadcasting.jl | 35 +++++++++++++++-------- test/manyvariables.jl | 5 ++++ test/onevariable.jl | 3 +- 6 files changed, 97 insertions(+), 57 deletions(-) diff --git a/src/TaylorSeries.jl b/src/TaylorSeries.jl index b8c9972c..c820d76c 100644 --- a/src/TaylorSeries.jl +++ b/src/TaylorSeries.jl @@ -32,7 +32,7 @@ end import Base: ==, +, -, *, /, ^ import Base: iterate, size, eachindex, firstindex, lastindex, - eltype, length, getindex, setindex!, axes#, similar + eltype, length, getindex, setindex!, axes, copyto!, similar import Base: zero, one, zeros, ones, isinf, isnan, iszero, convert, promote_rule, promote, show, diff --git a/src/auxiliary.jl b/src/auxiliary.jl index 70632b90..cf8fe84f 100644 --- a/src/auxiliary.jl +++ b/src/auxiliary.jl @@ -240,17 +240,17 @@ end ## similar ## -# similar(a::Taylor1) = Taylor1(similar(a.coeffs), a.order) -# function similar(a::Array{Taylor1{T},1}) where {T} -# ret = Vector{Taylor1{T}}(undef, size(a,1)) -# a1 = a[1].coeffs -# fill!(ret, similar(a1)) -# return ret -# end +similar(a::Taylor1) = Taylor1(similar(a.coeffs), a.order) +function similar(a::Array{Taylor1{T},1}) where {T} + ret = Vector{Taylor1{T}}(undef, size(a,1)) + fill!(ret, similar(a[1].coeffs)) + return ret +end +# similar(a::Taylor1, R::Type) = Taylor1(similar(a.coeffs, R), a.order) # similar(a::Array{Taylor1{T},1}, R::Type) where {T} = # convert(Vector{promote_type(Taylor1{T},R)}, similar(a)) -# similar(a::HomogeneousPolynomial) = HomogeneousPolynomial(similar(a.coeffs), a.order) +similar(a::HomogeneousPolynomial) = HomogeneousPolynomial(similar(a.coeffs), a.order) # function similar(a::Array{HomogeneousPolynomial{T},1}) where {T} # ret = Vector{HomogeneousPolynomial{T}}(undef, size(a,1)) # @simd for i in eachindex(a) @@ -261,7 +261,7 @@ end # similar(a::Array{HomogeneousPolynomial{T},1}, R::Type{<:NumberNotSeriesN}) where {T} = # convert(Array{HomogeneousPolynomial{R},1}, similar(a)) -# similar(a::TaylorN) = TaylorN(similar(a.coeffs), a.order) +similar(a::TaylorN) = TaylorN(similar(a.coeffs), a.order) # function similar(a::Array{TaylorN{T},1}) where {T} # ret = Vector{TaylorN{T}}(undef, size(a,1)) # @simd for i in eachindex(a) @@ -274,6 +274,31 @@ end +## copyto! ## +# Inspired from base/abstractarray.jl, line 665 +function copyto!(dst::Taylor1{T}, src::Taylor1{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 + +function copyto!(dst::TaylorN{T}, src::TaylorN{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 + + """ constant_term(a) diff --git a/src/broadcasting.jl b/src/broadcasting.jl index 0b37b7ab..c0a723f6 100644 --- a/src/broadcasting.jl +++ b/src/broadcasting.jl @@ -88,23 +88,20 @@ end # Adapted from Base.Broadcast.copyto!, base/broadcasting.jl, line 832 -# @inline function Base.copyto!(dest::Taylor1{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{AbstractArray} # 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) -# # I is the coefficients index -# @simd for I in eachindex(bc′) -# @inbounds dest[I-1] = getcoeff(bc′[I], I-1) -# end -# return dest -# end -# +@inline function Base.copyto!(dest::Taylor1{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{Taylor1{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 + # @inline function Base.copyto!(dest::HomogeneousPolynomial{T}, bc::Broadcasted) where {T<:NumberNotSeries} # 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 @@ -122,23 +119,24 @@ end # end # return dest # end -# -# @inline function Base.copyto!(dest::TaylorN{T}, bc::Broadcasted) where {T<:NumberNotSeries} -# 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{AbstractArray} # 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.coeffs, bc) -# # I is the coefficients index -# @simd for I in eachindex(bc′) -# @inbounds dest[I] = getcoeff(bc′[I], I) -# end -# return dest -# end + +@inline function Base.copyto!(dest::TaylorN{T}, bc::Broadcasted) where {T<:NumberNotSeries} + 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{TaylorN{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.coeffs, bc) + # I is the coefficients index + @simd for I in eachindex(bc′.args) + @show(I, bc′.args[I], dest[I]) + @inbounds dest[I] = getcoeff(bc′[I], I) + end + return dest +end # Broadcasted extensions diff --git a/test/broadcasting.jl b/test/broadcasting.jl index f769bfcb..22ec9eca 100644 --- a/test/broadcasting.jl +++ b/test/broadcasting.jl @@ -37,8 +37,13 @@ end @test xH .== xH @test yH .≈ yH + @test xH .== xH @test x[2] .== y[2] + @test 1 .* xH == xH + @test 1 .* [xH] == [xH] + @test [1] .* xH == [xH] + @test x .== x @test y .≈ y @test x .!= (1 + x) @@ -50,7 +55,7 @@ end @test typeof(1.0 .+ [y]) == Vector{TaylorN{Float64}} @test 1.0 .+ [x, 2y] == [1.0 + x, 1.0 + 2y] @test [1.0,2.0] .+ [x 2y] == [1.0+x 1.0+2y; 2.0+x 2.0+2y] - @test [1.0] .+ x == [1.0 + x] + @test [1.0] .+ x == x .+ [1.0] == [1.0 + x] @test 1.0 .* y == y @test typeof(1.0 .* x .* y) == TaylorN{Float64} # @@ -70,21 +75,27 @@ end @testset "Broadcasting for nested Taylor1's" begin t = Taylor1(Int, 3) - # @test typeof(similar(t)) == Taylor1{Int} - @test get_order(t) == 3 + ts = similar(t); + @test get_order(t) == get_order(ts) == 3 + @test typeof(ts) == Taylor1{Int} + @. ts = 3 * t^2 - 1 + @test ts == 3 * t^2 - 1 - tt = Taylor1([zero(t), one(t), zero(t)]) - # @test typeof(similar(tt)) == Taylor1{Taylor1{Int}} - # @test get_order(tt) == get_order(similar(tt)) == 2 + tt = Taylor1([zero(t), one(t)], 2) + tts = similar(tt); + @test typeof(tts) == Taylor1{Taylor1{Int}} + @test get_order(tt) == get_order(tts) == 2 @test tt .== tt + @. tts = 3 * tt^2 - 1 + @test tts == 3 * tt^2 - 1 + @test eltype(similar(1.0*tt)) == Taylor1{Float64} ttt = Taylor1([zero(tt), one(tt)]) - # @test typeof(similar(ttt)) == Taylor1{Taylor1{Taylor1{Int}}} - # @test get_order(ttt) == get_order(similar(ttt)) == 1 + ttts = similar(ttt); + @test typeof(ttts) == Taylor1{Taylor1{Taylor1{Int}}} + @test get_order(ttt) == get_order(ttts) == 1 @test ttt .≈ ttt + @. ttts = 3 * ttt^2 - 1 + @test ttts == -1 - tf = Taylor1([zero(1.0*t), one(t)]) - # @test typeof(similar(tf)) == Taylor1{Taylor1{Float64}} - # @test get_order(tf) == get_order(similar(tf)) == 1 - @test tf .== tt end diff --git a/test/manyvariables.jl b/test/manyvariables.jl index 5ed5870f..50f3a3f3 100644 --- a/test/manyvariables.jl +++ b/test/manyvariables.jl @@ -44,6 +44,8 @@ eeuler = Base.MathConstants.e @test length(get_variables()) == get_numvars() x, y = set_variables("x y", order=6) + @test axes(x) == axes(y) == () + @test axes(x[1]) == axes(y[2]) == () @test size(x) == (7,) @test firstindex(x) == 0 @test lastindex(y) == get_order() @@ -83,6 +85,9 @@ eeuler = Base.MathConstants.e yT = TaylorN(Int, 2, order=17) zeroT = zero( TaylorN([xH],1) ) @test zeroT.coeffs == zeros(HomogeneousPolynomial{Int}, 1) + @test size(xH) == (2,) + @test firstindex(xH) == 1 + @test lastindex(yH) == 2 @test length(zeros(HomogeneousPolynomial{Int}, 1)) == 2 @test one(HomogeneousPolynomial(1,1)) == HomogeneousPolynomial([1,1]) uT = one(convert(TaylorN{Float64},yT)) diff --git a/test/onevariable.jl b/test/onevariable.jl index 9b22fe0a..cc958c40 100644 --- a/test/onevariable.jl +++ b/test/onevariable.jl @@ -28,7 +28,8 @@ eeuler = Base.MathConstants.e @test iterate(t) == (0.0, 1) @test iterate(t, 1) == (1.0, 2) @test iterate(t, 16) == nothing - + @test axes(t) == () + @test axes([t]) == (Base.OneTo(1),) v = [1,2] @test typeof(TaylorSeries.resize_coeffs1!(v,3)) == Nothing From f6ce08ae219dfc646488c1f5127721666c196c69 Mon Sep 17 00:00:00 2001 From: Luis Benet Date: Thu, 25 Apr 2019 19:20:39 -0500 Subject: [PATCH 11/12] Further fixes and cleanup Inner mutating functions are not extended to be broadcasted --- src/TaylorSeries.jl | 2 +- src/auxiliary.jl | 106 +++++++++++++----------------- src/broadcasting.jl | 147 +++++++++++++++++++----------------------- test/broadcasting.jl | 86 ++++++++++++------------ test/manyvariables.jl | 34 +++++----- test/onevariable.jl | 20 +++--- 6 files changed, 185 insertions(+), 210 deletions(-) diff --git a/src/TaylorSeries.jl b/src/TaylorSeries.jl index c820d76c..188f85c5 100644 --- a/src/TaylorSeries.jl +++ b/src/TaylorSeries.jl @@ -60,7 +60,6 @@ include("hash_tables.jl") include("constructors.jl") include("conversion.jl") include("auxiliary.jl") -include("broadcasting.jl") include("arithmetic.jl") include("power.jl") include("functions.jl") @@ -68,6 +67,7 @@ include("other_functions.jl") include("evaluate.jl") include("calculus.jl") include("dictmutfunct.jl") +include("broadcasting.jl") include("printing.jl") function __init__() diff --git a/src/auxiliary.jl b/src/auxiliary.jl index cf8fe84f..dc42235a 100644 --- a/src/auxiliary.jl +++ b/src/auxiliary.jl @@ -179,33 +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 - @inline iterate(a::$T, state=0) = state > a.order ? nothing : (a.coeffs[state+1], state+1) + 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) - # Base.iterate(rS::Iterators.Reverse{$T}, state=rS.itr.order) = state < 0 ? nothing : (a.coeffs[state], state-1) @inline eltype(::$T{S}) where {S<:Number} = S - @inline length(a::$T) = length(a.coeffs) @inline size(a::$T) = size(a.coeffs) - @inline firstindex(a::$T) = 0 - @inline lastindex(a::$T) = a.order @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) -@inline eltype(::HomogeneousPolynomial{S}) where {S<:Number} = S -@inline length(a::HomogeneousPolynomial) = size_table[a.order+1]#length(a.coeffs) -@inline size(a::HomogeneousPolynomial) = size(a.coeffs) -@inline firstindex(a::HomogeneousPolynomial) = 1 -@inline lastindex(a::HomogeneousPolynomial) = length(a) -@inline get_order(a::HomogeneousPolynomial) = a.order -@inline axes(a::HomogeneousPolynomial) = () - ## fixorder ## for T in (:Taylor1, :TaylorN) @@ -240,62 +237,51 @@ end ## similar ## -similar(a::Taylor1) = Taylor1(similar(a.coeffs), a.order) -function similar(a::Array{Taylor1{T},1}) where {T} - ret = Vector{Taylor1{T}}(undef, size(a,1)) - fill!(ret, similar(a[1].coeffs)) - return ret -end # similar(a::Taylor1, R::Type) = Taylor1(similar(a.coeffs, R), a.order) # similar(a::Array{Taylor1{T},1}, R::Type) where {T} = # convert(Vector{promote_type(Taylor1{T},R)}, similar(a)) -similar(a::HomogeneousPolynomial) = HomogeneousPolynomial(similar(a.coeffs), a.order) -# function similar(a::Array{HomogeneousPolynomial{T},1}) where {T} -# ret = Vector{HomogeneousPolynomial{T}}(undef, size(a,1)) -# @simd for i in eachindex(a) -# @inbounds ret[i] = similar(a[i]) -# end -# return ret -# end # similar(a::Array{HomogeneousPolynomial{T},1}, R::Type{<:NumberNotSeriesN}) where {T} = # convert(Array{HomogeneousPolynomial{R},1}, similar(a)) -similar(a::TaylorN) = TaylorN(similar(a.coeffs), a.order) -# function similar(a::Array{TaylorN{T},1}) where {T} -# ret = Vector{TaylorN{T}}(undef, size(a,1)) -# @simd for i in eachindex(a) -# @inbounds ret[i] = similar(a[i]) -# end -# return ret -# end # similar(a::Array{taylor_expand{T},1}, R::Type) where {T} = # convert(Vector{promote_type(TaylorN{T},R)}, similar(a)) -## copyto! ## -# Inspired from base/abstractarray.jl, line 665 -function copyto!(dst::Taylor1{T}, src::Taylor1{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 +for T in (:Taylor1, :HomogeneousPolynomial, :TaylorN) + @eval begin + ## similar ## + similar(a::$T) = $T(similar(a.coeffs), a.order) + if $T == Taylor1 + function similar(a::Array{$T{T},1}) where {T} + ret = Vector{$T{T}}(undef, size(a,1)) + fill!(ret, similar(a[1].coeffs)) + return ret + end + else + function similar(a::Array{$T{T},1}) where {T} + ret = Vector{$T{T}}(undef, size(a,1)) + @simd for i in eachindex(a) + @inbounds ret[i] = similar(a[i]) + end + return ret + end + end -function copyto!(dst::TaylorN{T}, src::TaylorN{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]) + ## copyto! ## + # Inspired from base/abstractarray.jl, line 665 + 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 - return dst end diff --git a/src/broadcasting.jl b/src/broadcasting.jl index c0a723f6..27212120 100644 --- a/src/broadcasting.jl +++ b/src/broadcasting.jl @@ -54,97 +54,82 @@ find_taylor(a::HomogeneousPolynomial, rest) = a find_taylor(a::TaylorN, rest) = a find_taylor(::AbstractArray, rest) = find_taylor(rest) -# Extend Base.similar -function Base.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 +# 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 Base.similar(bc::Broadcasted{HomogeneousPolynomialStyle{S}}, ::Type{T}) where {S, T} - # Proper promotion - @show(Base.Broadcast.eltypes(bc.args)) - # combine_eltypes(f, args::Tuple) = Base._return_type(f, eltypes(args)) - R = Base.Broadcast.combine_eltypes(bc.f, bc.args) - @show(bc.f, bc.args) - # Scan the inputs for the HomogeneousPolynomial: - A = find_taylor(bc) - @show(A, R) - # Create the output - return HomogeneousPolynomial(similar(A.coeffs, R), A.order) -end +# function similar(bc::Broadcasted{HomogeneousPolynomialStyle{S}}, ::Type{T}) where {S, T} +# # Proper promotion +# @show(Base.Broadcast.eltypes(bc.args)) +# # combine_eltypes(f, args::Tuple) = Base._return_type(f, eltypes(args)) +# R = Base.Broadcast.combine_eltypes(bc.f, bc.args) +# @show(bc.f, bc.args) +# # Scan the inputs for the HomogeneousPolynomial: +# A = find_taylor(bc) +# @show(A, R) +# # Create the output +# return HomogeneousPolynomial(similar(A.coeffs, R), A.order) +# end -function Base.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 +# 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 -@inline function Base.copyto!(dest::Taylor1{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{Taylor1{T}} # only a single input argument to broadcast! - A = bc.args[1] - if axes(dest) == axes(A) - return copyto!(dest, A) +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 - bc′ = Base.Broadcast.preprocess(dest, bc) - copyto!(dest, bc′[1]) - return dest end -# @inline function Base.copyto!(dest::HomogeneousPolynomial{T}, bc::Broadcasted) where {T<:NumberNotSeries} -# 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{AbstractArray} # only a single input argument to broadcast! -# A = bc.args[1] -# if axes(dest) == axes(A) -# return copyto!(dest, A) + +# Broadcasted extensions +# 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 -# bc′= Base.Broadcast.preprocess(dest.coeffs, bc) -# # I is the coefficients index -# @simd for I in eachindex(bc′) -# # @inbounds dest[I] = getcoeff(bc′[I], I) -# @inbounds dest[I] = bc′[I].coeffs[I] +# 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 -# return dest # end - -@inline function Base.copyto!(dest::TaylorN{T}, bc::Broadcasted) where {T<:NumberNotSeries} - 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{TaylorN{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.coeffs, bc) - # I is the coefficients index - @simd for I in eachindex(bc′.args) - @show(I, bc′.args[I], dest[I]) - @inbounds dest[I] = getcoeff(bc′[I], I) - end - return dest -end - - -# Broadcasted extensions -# This prevents broadcasting being applied to `a` and `b` -# for the mutating functions, and to act only in `k` -function broadcasted(::Taylor1Style{T}, f!, a::Taylor1{T}, b::Taylor1{T}, k) where {T} - @inbounds for i in eachindex(k) - f!(a, b, k[i]) - end - nothing -end diff --git a/test/broadcasting.jl b/test/broadcasting.jl index 22ec9eca..5c420019 100644 --- a/test/broadcasting.jl +++ b/test/broadcasting.jl @@ -25,9 +25,37 @@ using Test st = sin(t) @test st .== st - @test st(pi/3) == evaluate(st, pi/3) + @test st == sin.(t) @test st.(pi/3) == evaluate(st, pi/3) + @test st(pi/3) == evaluate.(st, pi/3) @test st.([0.0, pi/3]) == evaluate(st, [0.0, pi/3]) + + # Nested Taylor1 tests + t = Taylor1(Int, 3) + ts = similar(t); + @test get_order(t) == get_order(ts) == 3 + @test typeof(ts) == Taylor1{Int} + ts .= t + @test ts == t + @. ts = 3 * t^2 - 1 + @test ts == 3 * t^2 - 1 + + tt = Taylor1([zero(t), one(t)], 2) + tts = similar(tt); + @test typeof(tts) == Taylor1{Taylor1{Int}} + @test get_order(tt) == get_order(tts) == 2 + @test tt .== tt + @. tts = 3 * tt^2 - 1 + @test tts == 3 * tt^2 - 1 + @test eltype(similar(1.0*tt)) == Taylor1{Float64} + + ttt = Taylor1([zero(tt), one(tt)]) + ttts = similar(ttt); + @test typeof(ttts) == Taylor1{Taylor1{Taylor1{Int}}} + @test get_order(ttt) == get_order(ttts) == 1 + @test ttt .≈ ttt + @. ttts = 3 * ttt^2 - 1 + @test ttts == -1 end @testset "Broadcasting for HomogeneousPolynomial and TaylorN" begin @@ -40,6 +68,14 @@ end @test xH .== xH @test x[2] .== y[2] + xHs = similar(xH) + @test typeof(xHs) == typeof(xH) + @test get_order(xHs) == get_order(xH) + xHs .= xH + @test xHs == xH + @. xHs = 2 * xH + yH + @test xHs == 2 * xH + yH + @test 1 .* xH == xH @test 1 .* [xH] == [xH] @test [1] .* xH == [xH] @@ -48,6 +84,14 @@ end @test y .≈ y @test x .!= (1 + x) + p = similar(x) + @test typeof(p) == typeof(x) + @test get_order(p) == get_order(x) + p .= x + @test p == x + @. p = 1 + 2*x + 3x^2 - x * y + @test p == 1 + 2*x + 3*x^2 - x * y + @test 1.0 .+ x == 1.0 + x @test y .+ x == y + x @test typeof(big"1.0" .+ x) == TaylorN{BigFloat} @@ -58,44 +102,4 @@ end @test [1.0] .+ x == x .+ [1.0] == [1.0 + x] @test 1.0 .* y == y @test typeof(1.0 .* x .* y) == TaylorN{Float64} - # - # @test 1.0 .+ t == 1.0 + t - # @test typeof(1.0 .+ t) == Taylor1{Float64} - # @test 1.0 .+ [t] == [1.0 + t] - # @test typeof(1.0 .+ [t]) == Vector{Taylor1{Float64}} - # @test 1.0 .+ [t, 2t] == [1.0 + t, 1.0 + 2t] - # @test [1.0,2.0] .+ [t 2t] == [1.0+t 1.0+2t; 2.0+t 2.0+2t] - # @test [1.0] .+ t == [1.0 + t] - # @test 1.0 .* t == t - # @test typeof(1.0 .* t) == Taylor1{Float64} - - # xH = HomogeneousPolynomial([one(t), zero(t)]) - # @assert typeof(similar(xHt1)) == HomogeneousPolynomial{Taylor1{Int64}} -end - -@testset "Broadcasting for nested Taylor1's" begin - t = Taylor1(Int, 3) - ts = similar(t); - @test get_order(t) == get_order(ts) == 3 - @test typeof(ts) == Taylor1{Int} - @. ts = 3 * t^2 - 1 - @test ts == 3 * t^2 - 1 - - tt = Taylor1([zero(t), one(t)], 2) - tts = similar(tt); - @test typeof(tts) == Taylor1{Taylor1{Int}} - @test get_order(tt) == get_order(tts) == 2 - @test tt .== tt - @. tts = 3 * tt^2 - 1 - @test tts == 3 * tt^2 - 1 - @test eltype(similar(1.0*tt)) == Taylor1{Float64} - - ttt = Taylor1([zero(tt), one(tt)]) - ttts = similar(ttt); - @test typeof(ttts) == Taylor1{Taylor1{Taylor1{Int}}} - @test get_order(ttt) == get_order(ttts) == 1 - @test ttt .≈ ttt - @. ttts = 3 * ttt^2 - 1 - @test ttts == -1 - end diff --git a/test/manyvariables.jl b/test/manyvariables.jl index 50f3a3f3..98f2477c 100644 --- a/test/manyvariables.jl +++ b/test/manyvariables.jl @@ -631,31 +631,31 @@ eeuler = Base.MathConstants.e @test Q() == evaluate(Q) @test Q[1:3]() == evaluate(Q[1:3]) - # dx = set_variables("x", numvars=4, order=10) - # for i in 1:4 - # @test deg2rad(180+dx[i]) == pi + deg2rad(1.0)dx[i] - # @test rad2deg(pi+dx[i]) == 180.0+rad2deg(1.0)dx[i] - # end - # p = sin(exp(dx[1]*dx[2])+dx[3]*dx[2])/(1.0+dx[4]^2) - # q = zero(p) - # TaylorSeries.deg2rad!(q, p, 0) - # @test q[0] == p[0]*(pi/180) + dx = set_variables("x", numvars=4, order=10) + for i in 1:4 + @test deg2rad(180+dx[i]) == pi + deg2rad(1.0)dx[i] + @test rad2deg(pi+dx[i]) == 180.0+rad2deg(1.0)dx[i] + end + p = sin(exp(dx[1]*dx[2])+dx[3]*dx[2])/(1.0+dx[4]^2) + q = zero(p) + TaylorSeries.deg2rad!(q, p, 0) + @test q[0] == p[0]*(pi/180) # TaylorSeries.deg2rad!.(q, p, [1,3,5]) # for i in [0,1,3,5] # @test q[i] == p[i]*(pi/180) # end - # TaylorSeries.rad2deg!(q, p, 0) - # @test q[0] == p[0]*(180/pi) + TaylorSeries.rad2deg!(q, p, 0) + @test q[0] == p[0]*(180/pi) # TaylorSeries.rad2deg!.(q, p, [1,3,5]) # for i in [0,1,3,5] # @test q[i] == p[i]*(180/pi) # end - # xT = 5+TaylorN(Int, 1, order=10) - # yT = TaylorN(2, order=10) - # TaylorSeries.deg2rad!(yT, xT, 0) - # @test yT[0] == xT[0]*(pi/180) - # TaylorSeries.rad2deg!(yT, xT, 0) - # @test yT[0] == xT[0]*(180/pi) + xT = 5+TaylorN(Int, 1, order=10) + yT = TaylorN(2, order=10) + TaylorSeries.deg2rad!(yT, xT, 0) + @test yT[0] == xT[0]*(pi/180) + TaylorSeries.rad2deg!(yT, xT, 0) + @test yT[0] == xT[0]*(180/pi) end @testset "Integrate for several variables" begin diff --git a/test/onevariable.jl b/test/onevariable.jl index cc958c40..e0bd97ba 100644 --- a/test/onevariable.jl +++ b/test/onevariable.jl @@ -469,22 +469,22 @@ eeuler = Base.MathConstants.e TaylorSeries.deg2rad!(b, a, 0) @test a == c @test a[0]*(pi/180) == b[0] - TaylorSeries.deg2rad!.(b, a, [0,1,2]) - @test a == c - for i in 0:2 - @test a[i]*(pi/180) == b[i] - end + # TaylorSeries.deg2rad!.(b, a, [0,1,2]) + # @test a == c + # for i in 0:2 + # @test a[i]*(pi/180) == b[i] + # end a = Taylor1(rand(10)) b = Taylor1(rand(10)) c = deepcopy(a) TaylorSeries.rad2deg!(b, a, 0) @test a == c @test a[0]*(180/pi) == b[0] - TaylorSeries.rad2deg!.(b, a, [0,1,2]) - @test a == c - for i in 0:2 - @test a[i]*(180/pi) == b[i] - end + # TaylorSeries.rad2deg!.(b, a, [0,1,2]) + # @test a == c + # for i in 0:2 + # @test a[i]*(180/pi) == b[i] + # end # Test additional Taylor1 constructors @test Taylor1{Float64}(true) == Taylor1([1.0]) From 3152be7278605535f9fe0180dde43cdc4f0de0e4 Mon Sep 17 00:00:00 2001 From: Luis Benet Date: Thu, 25 Apr 2019 21:14:13 -0500 Subject: [PATCH 12/12] Use zero in similar to avoid undef's (with mixtures) And fix the vector case of constant_term and linear_polynomial --- src/auxiliary.jl | 35 ++++++++++------------------------- src/broadcasting.jl | 3 --- test/onevariable.jl | 4 ++-- 3 files changed, 12 insertions(+), 30 deletions(-) diff --git a/src/auxiliary.jl b/src/auxiliary.jl index dc42235a..1112c67c 100644 --- a/src/auxiliary.jl +++ b/src/auxiliary.jl @@ -236,37 +236,22 @@ function Base.findlast(a::Taylor1{T}) where {T<:Number} end -## similar ## -# similar(a::Taylor1, R::Type) = Taylor1(similar(a.coeffs, R), a.order) -# similar(a::Array{Taylor1{T},1}, R::Type) where {T} = -# convert(Vector{promote_type(Taylor1{T},R)}, similar(a)) - -# similar(a::Array{HomogeneousPolynomial{T},1}, R::Type{<:NumberNotSeriesN}) where {T} = -# convert(Array{HomogeneousPolynomial{R},1}, similar(a)) - -# similar(a::Array{taylor_expand{T},1}, R::Type) where {T} = -# convert(Vector{promote_type(TaylorN{T},R)}, similar(a)) - - - +## similar and copyto! ## for T in (:Taylor1, :HomogeneousPolynomial, :TaylorN) @eval begin ## similar ## similar(a::$T) = $T(similar(a.coeffs), a.order) - if $T == Taylor1 - function similar(a::Array{$T{T},1}) where {T} - ret = Vector{$T{T}}(undef, size(a,1)) + # similar(a::$T, R::Type) = $T(similar(a.coeffs, R), a.order) + function similar(a::Array{$T{T},1}) where {T} + ret = Vector{$T{T}}(undef, size(a,1)) + if $T == Taylor1 fill!(ret, similar(a[1].coeffs)) - return ret - end - else - function similar(a::Array{$T{T},1}) where {T} - ret = Vector{$T{T}}(undef, size(a,1)) + else @simd for i in eachindex(a) - @inbounds ret[i] = similar(a[i]) + @inbounds ret[i] = zero(a[i]) # avoid `undef`s end - return ret end + return ret end ## copyto! ## @@ -296,7 +281,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 @@ -311,6 +296,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 diff --git a/src/broadcasting.jl b/src/broadcasting.jl index 27212120..89735df4 100644 --- a/src/broadcasting.jl +++ b/src/broadcasting.jl @@ -66,13 +66,10 @@ find_taylor(::AbstractArray, rest) = find_taylor(rest) # function similar(bc::Broadcasted{HomogeneousPolynomialStyle{S}}, ::Type{T}) where {S, T} # # Proper promotion -# @show(Base.Broadcast.eltypes(bc.args)) # # combine_eltypes(f, args::Tuple) = Base._return_type(f, eltypes(args)) # R = Base.Broadcast.combine_eltypes(bc.f, bc.args) -# @show(bc.f, bc.args) # # Scan the inputs for the HomogeneousPolynomial: # A = find_taylor(bc) -# @show(A, R) # # Create the output # return HomogeneousPolynomial(similar(A.coeffs, R), A.order) # end diff --git a/test/onevariable.jl b/test/onevariable.jl index e0bd97ba..ad3e57c7 100644 --- a/test/onevariable.jl +++ b/test/onevariable.jl @@ -114,11 +114,11 @@ eeuler = Base.MathConstants.e @test constant_term(2.0) == 2.0 @test constant_term(t) == 0 @test constant_term(tim) == complex(0, 0) - @test constant_term([zt, t]) == zt + @test constant_term([zt, t]) == [0, 0] @test linear_polynomial(2) == 2 @test linear_polynomial(t) == t @test linear_polynomial(tim^2) == zero(tim) - @test linear_polynomial([zero(tim), tim]) == zero(tim) + @test linear_polynomial([zero(tim), tim]) == [zero(tim), tim] @test ot == 1 @test 0.0 == zt