From a79b0221c200be65dfa1788d4f8457f176652c5d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jorge=20P=C3=A9rez?= Date: Sun, 31 Mar 2024 04:16:44 +0200 Subject: [PATCH 1/5] Add StaticArrays extension (#353) * Add StaticArrays extension * Bump patch version * Add promote_op for matrix-matrix product * Update runtests.jl * Add Base.float dispatches for TaylorSeries types and variables --- Project.toml | 8 ++++++-- ext/TaylorSeriesSAExt.jl | 14 ++++++++++++++ src/TaylorSeries.jl | 5 +++++ src/conversion.jl | 9 +++++++++ test/manyvariables.jl | 9 +++++++++ test/mixtures.jl | 11 +++++++++++ test/onevariable.jl | 11 +++++++++++ test/runtests.jl | 4 ++-- test/staticarrays.jl | 20 ++++++++++++++++++++ 9 files changed, 87 insertions(+), 4 deletions(-) create mode 100644 ext/TaylorSeriesSAExt.jl create mode 100644 test/staticarrays.jl diff --git a/Project.toml b/Project.toml index 1b5a38f9..5815c353 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "TaylorSeries" uuid = "6aa5eb33-94cf-58f4-a9d0-e4b2c4fc25ea" repo = "https://github.com/JuliaDiff/TaylorSeries.jl.git" -version = "0.17.2" +version = "0.17.3" [deps] LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" @@ -11,9 +11,11 @@ SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" [weakdeps] IntervalArithmetic = "d1acc4aa-44c8-5952-acd4-ba5d80a2a253" +StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" [extensions] TaylorSeriesIAExt = "IntervalArithmetic" +TaylorSeriesSAExt = "StaticArrays" [compat] Aqua = "0.8" @@ -22,6 +24,7 @@ LinearAlgebra = "<0.0.1, 1" Markdown = "<0.0.1, 1" Requires = "0.5.2, 1" SparseArrays = "<0.0.1, 1" +StaticArrays = "1" Test = "<0.0.1, 1" julia = "1" @@ -30,7 +33,8 @@ Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595" IntervalArithmetic = "d1acc4aa-44c8-5952-acd4-ba5d80a2a253" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" +StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [targets] -test = ["IntervalArithmetic", "LinearAlgebra", "SparseArrays", "Test", "Aqua"] +test = ["IntervalArithmetic", "LinearAlgebra", "SparseArrays", "StaticArrays", "Test", "Aqua"] diff --git a/ext/TaylorSeriesSAExt.jl b/ext/TaylorSeriesSAExt.jl new file mode 100644 index 00000000..0e5cadec --- /dev/null +++ b/ext/TaylorSeriesSAExt.jl @@ -0,0 +1,14 @@ +module TaylorSeriesSAExt + +using TaylorSeries + +import Base.promote_op +import LinearAlgebra: matprod + +isdefined(Base, :get_extension) ? (using StaticArrays) : (using ..StaticArrays) + +promote_op(::typeof(adjoint), ::Type{T}) where {T<:AbstractSeries} = T +promote_op(::typeof(matprod), ::Type{T}, ::Type{U}) where {T <: AbstractSeries, U <: AbstractFloat} = T +promote_op(::typeof(matprod), ::Type{T}, ::Type{T}) where {T <: AbstractSeries} = T + +end \ No newline at end of file diff --git a/src/TaylorSeries.jl b/src/TaylorSeries.jl index dca66013..23276b6b 100644 --- a/src/TaylorSeries.jl +++ b/src/TaylorSeries.jl @@ -49,6 +49,8 @@ import Base: zero, one, zeros, ones, isinf, isnan, iszero, isless, power_by_squaring, rtoldefault, isfinite, isapprox, rad2deg, deg2rad +import Base.float + export Taylor1, TaylorN, HomogeneousPolynomial, AbstractSeries, TS export getcoeff, derivative, integrate, differentiate, @@ -84,6 +86,9 @@ function __init__() @require IntervalArithmetic = "d1acc4aa-44c8-5952-acd4-ba5d80a2a253" begin include("../ext/TaylorSeriesIAExt.jl") end + @require StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" begin + include("../ext/TaylorSeriesSAExt.jl") + end end end diff --git a/src/conversion.jl b/src/conversion.jl index 154f5f17..53691240 100644 --- a/src/conversion.jl +++ b/src/conversion.jl @@ -201,3 +201,12 @@ function promote(a::Taylor1{Taylor1{T}}, b::Taylor1{T}) where {T<:NumberNotSerie end promote(b::Taylor1{T}, a::Taylor1{Taylor1{T}}) where {T<:NumberNotSeriesN} = reverse(promote(a, b)) + +# float +float(::Type{Taylor1{T}}) where T<:Number = Taylor1{float(T)} +float(::Type{HomogeneousPolynomial{T}}) where T<:Number = HomogeneousPolynomial{float(T)} +float(::Type{TaylorN{T}}) where T<:Number = TaylorN{float(T)} + +float(x::Taylor1{T}) where T<:Number = convert(Taylor1{float(T)}, x) +float(x::HomogeneousPolynomial{T}) where T<:Number = convert(HomogeneousPolynomial{float(T)}, x) +float(x::TaylorN{T}) where T<:Number = convert(TaylorN{float(T)}, x) diff --git a/test/manyvariables.jl b/test/manyvariables.jl index 6c91d637..5a519483 100644 --- a/test/manyvariables.jl +++ b/test/manyvariables.jl @@ -782,6 +782,15 @@ end # Lexicographic tests with 4 vars @test 1 > dx[1] > dx[2] > dx[3] > dx[4] @test dx[4]^2 < dx[3]*dx[4] < dx[3]^2 < dx[2]*dx[4] < dx[2]*dx[3] < dx[2]^2 < dx[1]*dx[4] < dx[1]*dx[3] < dx[1]*dx[2] < dx[1]^2 + + @testset "Test Base.float overloads for HomogeneousPolynomial and TaylorN" begin + @test float(HomogeneousPolynomial(-7, 2)) == HomogeneousPolynomial(-7.0, 2) + @test float(HomogeneousPolynomial(1+im, 2)) == HomogeneousPolynomial(float(1+im), 2) + @test float(TaylorN(Int, 2)) == TaylorN(2) + @test float(TaylorN(Complex{Rational}, 2)) == TaylorN(float(Complex{Rational}), 2) + @test float(HomogeneousPolynomial{Complex{Int}}) == float(HomogeneousPolynomial{Complex{Float64}}) + @test float(TaylorN{Complex{Rational}}) == float(TaylorN{Complex{Float64}}) + end end @testset "Integrate for several variables" begin diff --git a/test/mixtures.jl b/test/mixtures.jl index 2cb0078b..f5789bc4 100644 --- a/test/mixtures.jl +++ b/test/mixtures.jl @@ -495,6 +495,17 @@ using Test two = 2one(x[0]) @test two/x == 2/x == 2.0/x @test (2one(x))/x == 2/x + + @testset "Test Base.float overloads for Taylor1 and TaylorN mixtures" begin + q = get_variables(Int) + x1N = Taylor1(q) + @test float(x1N) == Taylor1(float.(q)) + xN1 = convert(TaylorN{Taylor1{Int}}, x1N) + @test float(xN1) == convert(TaylorN{Taylor1{Float64}}, Taylor1(float.(q))) + @test float(Taylor1{TaylorN{Int}}) == Taylor1{TaylorN{Float64}} + @test float(TaylorN{Taylor1{Int}}) == TaylorN{Taylor1{Float64}} + @test float(TaylorN{Taylor1{Complex{Int}}}) == TaylorN{Taylor1{Complex{Float64}}} + end end @testset "Tests with nested Taylor1s" begin diff --git a/test/onevariable.jl b/test/onevariable.jl index abd231d8..120a5a9f 100644 --- a/test/onevariable.jl +++ b/test/onevariable.jl @@ -625,6 +625,17 @@ Base.iszero(::SymbNumber) = false # Test fallback pretty_print st = Taylor1([SymbNumber(:x₀), SymbNumber(:x₁)]) @test string(st) == " SymbNumber(:x₀) + SymbNumber(:x₁) t + 𝒪(t²)" + + @testset "Test Base.float overloads for Taylor1" begin + @test float(Taylor1(-3, 2)) == Taylor1(-3.0, 2) + @test float(Taylor1(-1//2, 2)) == Taylor1(-0.5, 2) + @test float(Taylor1(3 - 0im, 2)) == Taylor1(3.0 - 0.0im, 2) + x = Taylor1(rand(5)) + @test float(x) == x + @test float(Taylor1{Int32}) == Taylor1{Float64} + @test float(Taylor1{Int}) == Taylor1{Float64} + @test float(Taylor1{Complex{Int}}) == Taylor1{ComplexF64} + end end @testset "Test inv for Matrix{Taylor1{Float64}}" begin diff --git a/test/runtests.jl b/test/runtests.jl index 707c02bd..afd7d3ab 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -3,7 +3,6 @@ # Tests for TaylorSeries testfiles = ( - "aqua.jl", "onevariable.jl", "manyvariables.jl", "mixtures.jl", @@ -11,7 +10,8 @@ testfiles = ( "intervals.jl", "broadcasting.jl", "identities_Euler.jl", - "fateman40.jl" + "fateman40.jl", + "staticarrays.jl" ) for file in testfiles diff --git a/test/staticarrays.jl b/test/staticarrays.jl new file mode 100644 index 00000000..ea84a8b0 --- /dev/null +++ b/test/staticarrays.jl @@ -0,0 +1,20 @@ +# This file is part of TaylorSeries.jl, MIT licensed +# + +using TaylorSeries, StaticArrays + +using Test +# eeuler = Base.MathConstants.e + +@testset "Tests TaylorSeries operations over StaticArrays types" begin + q = set_variables("q", order=2, numvars=2) + m = @SMatrix fill(Taylor1(rand(2).*q), 3, 3) + mt = m' + @test m isa SMatrix{3, 3, Taylor1{TaylorN{Float64}}, 9} + @test mt isa SMatrix{3, 3, Taylor1{TaylorN{Float64}}, 9} + v = @SVector [-1.1, 3.4, 7.62345e-1] + mtv = mt * v + @test mtv isa SVector{3, Taylor1{TaylorN{Float64}}} + mmt = m * mt + @test mmt isa SMatrix{3, 3, Taylor1{TaylorN{Float64}}, 9} +end From 87116db2213657ce53cf173d70321cf13631a881 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jorge=20P=C3=A9rez?= Date: Mon, 1 Apr 2024 20:12:05 +0200 Subject: [PATCH 2/5] Add JLD2 extension (#354) * Add JLD2 extension (move over from NEOs) * Add missing files * Bump patch version * Update Project.toml (order fields alphabetically) * Cleanup * Remove leftover comment --- Project.toml | 8 +++- ext/TaylorSeriesJLD2Ext.jl | 91 ++++++++++++++++++++++++++++++++++++++ src/TaylorSeries.jl | 3 ++ test/jld2.jl | 15 +++++++ test/runtests.jl | 3 +- test/staticarrays.jl | 1 - 6 files changed, 117 insertions(+), 4 deletions(-) create mode 100644 ext/TaylorSeriesJLD2Ext.jl create mode 100644 test/jld2.jl diff --git a/Project.toml b/Project.toml index 5815c353..8d9a2fb5 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "TaylorSeries" uuid = "6aa5eb33-94cf-58f4-a9d0-e4b2c4fc25ea" repo = "https://github.com/JuliaDiff/TaylorSeries.jl.git" -version = "0.17.3" +version = "0.17.4" [deps] LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" @@ -11,15 +11,18 @@ SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" [weakdeps] IntervalArithmetic = "d1acc4aa-44c8-5952-acd4-ba5d80a2a253" +JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" [extensions] TaylorSeriesIAExt = "IntervalArithmetic" +TaylorSeriesJLD2Ext = "JLD2" TaylorSeriesSAExt = "StaticArrays" [compat] Aqua = "0.8" IntervalArithmetic = "0.15 - 0.20" +JLD2 = "0.4" LinearAlgebra = "<0.0.1, 1" Markdown = "<0.0.1, 1" Requires = "0.5.2, 1" @@ -31,10 +34,11 @@ julia = "1" [extras] Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595" IntervalArithmetic = "d1acc4aa-44c8-5952-acd4-ba5d80a2a253" +JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [targets] -test = ["IntervalArithmetic", "LinearAlgebra", "SparseArrays", "StaticArrays", "Test", "Aqua"] +test = ["IntervalArithmetic", "JLD2", "LinearAlgebra", "SparseArrays", "StaticArrays", "Test", "Aqua"] diff --git a/ext/TaylorSeriesJLD2Ext.jl b/ext/TaylorSeriesJLD2Ext.jl new file mode 100644 index 00000000..6593060f --- /dev/null +++ b/ext/TaylorSeriesJLD2Ext.jl @@ -0,0 +1,91 @@ +module TaylorSeriesJLD2Ext + +import Base: convert +using TaylorSeries + +if isdefined(Base, :get_extension) + import JLD2: writeas +else + import ..JLD2: writeas +end + +@doc raw""" + TaylorNSerialization{T} + +Custom serialization struct to save a `TaylorN{T}` to a `.jld2` file. + +# Fields +- `vars::Vector{String}`: jet transport variables. +- `varorder::Int`: order of jet transport perturbations. +- `x::Vector{T}`: vector of coefficients. +""" +struct TaylorNSerialization{T} + vars::Vector{String} + varorder::Int + x::Vector{T} +end + +# Tell JLD2 to save TaylorN{T} as TaylorNSerialization{T} +writeas(::Type{TaylorN{T}}) where {T} = TaylorNSerialization{T} + +# Convert method to write .jld2 files +function convert(::Type{TaylorNSerialization{T}}, eph::TaylorN{T}) where {T} + # Variables + vars = TS.get_variable_names() + # Number of variables + n = length(vars) + # TaylorN order + varorder = eph.order + # Number of coefficients in each TaylorN + L = varorder + 1 + # Number of coefficients in each HomogeneousPolynomial + M = binomial(n + varorder, varorder) + + # Vector of coefficients + x = Vector{T}(undef, M) + + # Save coefficients + i = 1 + for i_1 in 0:varorder + # Iterate over i_1 order HomogeneousPolynomial + for i_2 in 1:binomial(n + i_1 - 1, i_1) + x[i] = eph.coeffs[i_1+1].coeffs[i_2] + i += 1 + end + end + + return TaylorNSerialization{T}(vars, varorder, x) +end + +# Convert method to read .jld2 files +function convert(::Type{TaylorN{T}}, eph::TaylorNSerialization{T}) where {T} + # Variables + vars = eph.vars + # Number of variables + n = length(vars) + # TaylorN order + varorder = eph.varorder + # Number of coefficients in each TaylorN + L = varorder + 1 + # Number of coefficients in each HomogeneousPolynomial + M = binomial(n + varorder, varorder) + + # Set variables + if TS.get_variable_names() != vars + TS.set_variables(T, vars, order = varorder) + end + + # Reconstruct TaylorN + i = 1 + TaylorN_coeffs = Vector{HomogeneousPolynomial{T}}(undef, L) + for i_1 in 0:varorder + # Reconstruct HomogeneousPolynomials + TaylorN_coeffs[i_1 + 1] = HomogeneousPolynomial(eph.x[i : i + binomial(n + i_1 - 1, i_1)-1], i_1) + i += binomial(n + i_1 - 1, i_1) + end + x = TaylorN{T}(TaylorN_coeffs, varorder) + + return x +end + +end \ No newline at end of file diff --git a/src/TaylorSeries.jl b/src/TaylorSeries.jl index 23276b6b..2e508e66 100644 --- a/src/TaylorSeries.jl +++ b/src/TaylorSeries.jl @@ -89,6 +89,9 @@ function __init__() @require StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" begin include("../ext/TaylorSeriesSAExt.jl") end + @require JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819" begin + include("../ext/TaylorSeriesJLD2Ext.jl") + end end end diff --git a/test/jld2.jl b/test/jld2.jl new file mode 100644 index 00000000..1cbf2e9d --- /dev/null +++ b/test/jld2.jl @@ -0,0 +1,15 @@ +# This file is part of TaylorSeries.jl, MIT licensed +# + +using TaylorSeries, JLD2 + +using Test + +@testset "Test TaylorSeries JLD2 extension" begin + dq = set_variables("q", order=4, numvars=6) + random_TaylorN = [cos(sum(dq .* rand(6))), sin(sum(dq .* rand(6))), tan(sum(dq .* rand(6)))] + jldsave("test.jld2"; random_TaylorN = random_TaylorN) + recovered_taylorN = JLD2.load("test.jld2", "random_TaylorN") + @test recovered_taylorN == random_TaylorN + rm("test.jld2") +end diff --git a/test/runtests.jl b/test/runtests.jl index afd7d3ab..9bdacdcb 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -11,7 +11,8 @@ testfiles = ( "broadcasting.jl", "identities_Euler.jl", "fateman40.jl", - "staticarrays.jl" + "staticarrays.jl", + "jld2.jl" ) for file in testfiles diff --git a/test/staticarrays.jl b/test/staticarrays.jl index ea84a8b0..188b6765 100644 --- a/test/staticarrays.jl +++ b/test/staticarrays.jl @@ -4,7 +4,6 @@ using TaylorSeries, StaticArrays using Test -# eeuler = Base.MathConstants.e @testset "Tests TaylorSeries operations over StaticArrays types" begin q = set_variables("q", order=2, numvars=2) From 3a7360f80fa83a8de723b87b7d437dd46f5484f0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jorge=20P=C3=A9rez?= Date: Wed, 17 Apr 2024 02:07:16 +0200 Subject: [PATCH 3/5] Add TaylorN constructor (#355) * Add TaylorN constructor * Bump patch version * Remove some tests added by mistake * Add test suggested by @lbenet --- Project.toml | 2 +- src/constructors.jl | 3 ++- test/manyvariables.jl | 7 +++++++ 3 files changed, 10 insertions(+), 2 deletions(-) diff --git a/Project.toml b/Project.toml index 8d9a2fb5..5289193f 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "TaylorSeries" uuid = "6aa5eb33-94cf-58f4-a9d0-e4b2c4fc25ea" repo = "https://github.com/JuliaDiff/TaylorSeries.jl.git" -version = "0.17.4" +version = "0.17.5" [deps] LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" diff --git a/src/constructors.jl b/src/constructors.jl index b6cd3740..8c8a5b3e 100644 --- a/src/constructors.jl +++ b/src/constructors.jl @@ -211,5 +211,6 @@ const NumberNotSeries = Union{Real,Complex} # A `Number` which is not `TaylorN` nor a `HomogeneousPolynomial` const NumberNotSeriesN = Union{Real,Complex,Taylor1} -## Additional Taylor1 outer constructor ## +## Additional Taylor1 and TaylorN outer constructor ## Taylor1{T}(x::S) where {T<:Number,S<:NumberNotSeries} = Taylor1([convert(T,x)], 0) +TaylorN{T}(x::S) where {T<:Number,S<:NumberNotSeries} = TaylorN(convert(T, x), TaylorSeries.get_order()) diff --git a/test/manyvariables.jl b/test/manyvariables.jl index 5a519483..40122b27 100644 --- a/test/manyvariables.jl +++ b/test/manyvariables.jl @@ -6,6 +6,12 @@ using Test # using LinearAlgebra @testset "Test hash tables" begin + a = TaylorN{Float64}(9//10) + a isa TaylorN{Float64} + @test constant_term(a) == Float64(9//10) + b = TaylorN{Complex{Float64}}(-4//7im) + @test b isa TaylorN{Complex{Float64}} + @test_throws MethodError TaylorN(-4//7im) # Issue #85 is solved! set_variables("x", numvars=66, order=1) @test TS._params_TaylorN_.order == get_order() == 1 @@ -791,6 +797,7 @@ end @test float(HomogeneousPolynomial{Complex{Int}}) == float(HomogeneousPolynomial{Complex{Float64}}) @test float(TaylorN{Complex{Rational}}) == float(TaylorN{Complex{Float64}}) end + end @testset "Integrate for several variables" begin From fdfa666f7e81ca2ff606d4a7336ba1e1b6056950 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jorge=20P=C3=A9rez?= Date: Tue, 23 Apr 2024 20:40:05 +0200 Subject: [PATCH 4/5] Add methods in support of PerezHz/TaylorIntegration.jl#192 (#356) * Add RecursiveArrayTools extension and setindex! dispatch for mixtures * Add RecursiveArrayTools extension and a Base.setindex! for mixtures * Fix in setindex! * Fix test * Add RAT compat entry (found by Aqua) --- Project.toml | 8 ++++++-- ext/TaylorSeriesRATExt.jl | 11 +++++++++++ src/auxiliary.jl | 2 +- test/mixtures.jl | 6 ++++++ test/rat.jl | 15 +++++++++++++++ test/runtests.jl | 3 ++- 6 files changed, 41 insertions(+), 4 deletions(-) create mode 100644 ext/TaylorSeriesRATExt.jl create mode 100644 test/rat.jl diff --git a/Project.toml b/Project.toml index 5289193f..1331fa67 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "TaylorSeries" uuid = "6aa5eb33-94cf-58f4-a9d0-e4b2c4fc25ea" repo = "https://github.com/JuliaDiff/TaylorSeries.jl.git" -version = "0.17.5" +version = "0.17.6" [deps] LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" @@ -13,11 +13,13 @@ SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" IntervalArithmetic = "d1acc4aa-44c8-5952-acd4-ba5d80a2a253" JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" +RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd" [extensions] TaylorSeriesIAExt = "IntervalArithmetic" TaylorSeriesJLD2Ext = "JLD2" TaylorSeriesSAExt = "StaticArrays" +TaylorSeriesRATExt = "RecursiveArrayTools" [compat] Aqua = "0.8" @@ -25,6 +27,7 @@ IntervalArithmetic = "0.15 - 0.20" JLD2 = "0.4" LinearAlgebra = "<0.0.1, 1" Markdown = "<0.0.1, 1" +RecursiveArrayTools = "2, 3" Requires = "0.5.2, 1" SparseArrays = "<0.0.1, 1" StaticArrays = "1" @@ -36,9 +39,10 @@ Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595" IntervalArithmetic = "d1acc4aa-44c8-5952-acd4-ba5d80a2a253" JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [targets] -test = ["IntervalArithmetic", "JLD2", "LinearAlgebra", "SparseArrays", "StaticArrays", "Test", "Aqua"] +test = ["IntervalArithmetic", "JLD2", "LinearAlgebra", "RecursiveArrayTools", "SparseArrays", "StaticArrays", "Test", "Aqua"] diff --git a/ext/TaylorSeriesRATExt.jl b/ext/TaylorSeriesRATExt.jl new file mode 100644 index 00000000..cc59f7f4 --- /dev/null +++ b/ext/TaylorSeriesRATExt.jl @@ -0,0 +1,11 @@ +module TaylorSeriesRATExt + +using TaylorSeries + +isdefined(Base, :get_extension) ? (import RecursiveArrayTools) : (import ..RecursiveArrayTools) + +function RecursiveArrayTools.recursivecopy(a::AbstractArray{<:AbstractSeries, N}) where N + deepcopy(a) +end + +end \ No newline at end of file diff --git a/src/auxiliary.jl b/src/auxiliary.jl index fd0406f5..7ab075dd 100644 --- a/src/auxiliary.jl +++ b/src/auxiliary.jl @@ -84,6 +84,7 @@ getindex(a::Taylor1{T}, u::StepRange{Int,Int}) where {T<:Number} = view(a.coeffs, u[:] .+ 1) setindex!(a::Taylor1{T}, x::T, n::Int) where {T<:Number} = a.coeffs[n+1] = x +setindex!(a::Taylor1{T}, x::T, n::Int) where {T<:AbstractSeries} = setindex!(a.coeffs, deepcopy(x), n+1) setindex!(a::Taylor1{T}, x::T, u::UnitRange{Int}) where {T<:Number} = a.coeffs[u .+ 1] .= x function setindex!(a::Taylor1{T}, x::Array{T,1}, u::UnitRange{Int}) where {T<:Number} @@ -405,4 +406,3 @@ macro isonethread(expr) end end) end - diff --git a/test/mixtures.jl b/test/mixtures.jl index f5789bc4..587889ab 100644 --- a/test/mixtures.jl +++ b/test/mixtures.jl @@ -496,6 +496,12 @@ using Test @test two/x == 2/x == 2.0/x @test (2one(x))/x == 2/x + dq = get_variables() + x = Taylor1(exp.(dq), 5) + x[1] = sin(dq[1]*dq[2]) + @test x[1] == sin(dq[1]*dq[2]) + @test x[1] !== sin(dq[1]*dq[2]) + @testset "Test Base.float overloads for Taylor1 and TaylorN mixtures" begin q = get_variables(Int) x1N = Taylor1(q) diff --git a/test/rat.jl b/test/rat.jl new file mode 100644 index 00000000..0fd42c63 --- /dev/null +++ b/test/rat.jl @@ -0,0 +1,15 @@ +# This file is part of TaylorSeries.jl, MIT licensed +# + +using TaylorSeries, RecursiveArrayTools + +using Test + +@testset "Tests TaylorSeries RecursiveArrayTools extension" begin + dq = get_variables() + x = Taylor1([0.9+2dq[1],-1.1dq[1], 0.7dq[2], 0.5dq[1]-0.45dq[2],0.9dq[1]]) + xx = [x,x] + yy = recursivecopy(xx) + @test yy == xx # yy and xx are equal... + @test yy !== xx # ...but they're not the same object in memory +end diff --git a/test/runtests.jl b/test/runtests.jl index 9bdacdcb..5991b446 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -12,7 +12,8 @@ testfiles = ( "identities_Euler.jl", "fateman40.jl", "staticarrays.jl", - "jld2.jl" + "jld2.jl", + "rat.jl" ) for file in testfiles From 466845e4b0b3901c29c3f952bfebac60819389ce Mon Sep 17 00:00:00 2001 From: Luis Benet Date: Wed, 22 May 2024 17:33:36 -0600 Subject: [PATCH 5/5] Fixes and additions in `evaluate`, including some in-line methods (#357) * Minor fixes in sqr_orderzero! and power_by_squaring * Fix methods of evaluate for Taylor1{TaylorN{T}} with tests. Some methods have been commented for now * Run Aqua tests only once at the end * Improvements when dealing with intervals * Add evaluate(::Taylor1{T}, ::TaylorN{T}), and other fixes * Add evaluate(::Taylor1{T}, ::Taylor1{TaylorN{T}) * Add _horner! as in-line evaluate method for certain mixtures * Add `fixorder` when needed for evaluate * Fix typo * Another fix (in fixorder) * Add methods for evaluate of a TaylorN by a value or TaylorN, including in-line methods * Fixes and tests * More tests a minor improvement * Avoid using `@isonethread` in evaluate, differentiate and integrate * Bump patch version --- Project.toml | 2 +- ext/TaylorSeriesIAExt.jl | 43 ++--- src/auxiliary.jl | 2 +- src/calculus.jl | 10 +- src/evaluate.jl | 343 +++++++++++++++++++++++++++++++++------ src/power.jl | 35 ++-- test/manyvariables.jl | 80 +++++---- test/mixtures.jl | 14 ++ test/runtests.jl | 6 +- 9 files changed, 406 insertions(+), 129 deletions(-) diff --git a/Project.toml b/Project.toml index 1331fa67..1105fae3 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "TaylorSeries" uuid = "6aa5eb33-94cf-58f4-a9d0-e4b2c4fc25ea" repo = "https://github.com/JuliaDiff/TaylorSeries.jl.git" -version = "0.17.6" +version = "0.17.7" [deps] LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" diff --git a/ext/TaylorSeriesIAExt.jl b/ext/TaylorSeriesIAExt.jl index 8fa11cdf..30ecc7c4 100644 --- a/ext/TaylorSeriesIAExt.jl +++ b/ext/TaylorSeriesIAExt.jl @@ -159,7 +159,7 @@ for T in (:Taylor1, :TaylorN) end -function evaluate(a::Taylor1{T}, dx::Interval{S}) where {T<:Real, S<:Real} +function evaluate(a::Taylor1, dx::Interval{S}) where {S<:Real} order = a.order uno = one(dx) dx2 = dx^2 @@ -183,34 +183,13 @@ function evaluate(a::TaylorN, dx::IntervalBox{N,T}) where {T<:Real,N} @assert N == get_numvars() a_length = length(a) suma = zero(constant_term(a)) + Interval{T}(0, 0) - @inbounds for homPol in length(a):-1:1 - suma += evaluate(a.coeffs[homPol], dx) + @inbounds for homPol in reverse(eachindex(a)) + suma += evaluate(a[homPol], dx) end return suma end -function evaluate(a::Taylor1{TaylorN{T}}, dx::Interval{S}) where {T<:Real, S<:Real} - order = a.order - uno = one(dx) - dx2 = dx^2 - if iseven(order) - kend = order-2 - @inbounds sum_even = a[end]*uno - @inbounds sum_odd = a[end-1]*zero(dx) - else - kend = order-3 - @inbounds sum_odd = a[end]*uno - @inbounds sum_even = a[end-1]*uno - end - @inbounds for k in kend:-2:0 - sum_odd = sum_odd*dx2 + a[k+1] - sum_even = sum_even*dx2 + a[k] - end - return sum_even + sum_odd*dx -end - - function evaluate(a::HomogeneousPolynomial, dx::IntervalBox{N,T}) where {T<:Real,N} @assert N == get_numvars() dx == IntervalBox(-1..1, Val(N)) && return _evaluate(a, dx, Val(true)) @@ -226,16 +205,16 @@ function _evaluate(a::HomogeneousPolynomial, dx::IntervalBox{N,T}, ::Val{true} ) @inbounds suma = a[1]*Interval{T}(0,0) Ieven = Interval{T}(0,1) - for (i,a_coeff) in enumerate(a.coeffs) + for (i, a_coeff) in enumerate(a.coeffs) iszero(a_coeff) && continue if isodd(sum(ct[i])) - suma += sum(a_coeff) * dx[1] - continue - end - @inbounds tmp = iseven(ct[i][1]) ? Ieven : dx[1] - for n in 2:N - @inbounds vv = iseven(ct[i][n]) ? Ieven : dx[1] - tmp *= vv + tmp = dx[1] + else + tmp = Ieven + for n in eachindex(ct[i]) + iseven(ct[i][n]) && continue + tmp *= dx[1] + end end suma += a_coeff * tmp end diff --git a/src/auxiliary.jl b/src/auxiliary.jl index 7ab075dd..22c1255c 100644 --- a/src/auxiliary.jl +++ b/src/auxiliary.jl @@ -259,7 +259,7 @@ end ## fixorder ## for T in (:Taylor1, :TaylorN) @eval begin - @inline function fixorder(a::$T{T}, b::$T{T}) where {T<:Number} + @inline function fixorder(a::$T, b::$T) a.order == b.order && return a, b minorder = _minorder(a, b) return $T(copy(a.coeffs), minorder), $T(copy(b.coeffs), minorder) diff --git a/src/calculus.jl b/src/calculus.jl index d3460230..466e81f8 100644 --- a/src/calculus.jl +++ b/src/calculus.jl @@ -141,9 +141,10 @@ function differentiate(a::HomogeneousPolynomial, r::Int) coeffs = zeros(T, num_coeffs) @inbounds posTb = pos_table[a.order] @inbounds num_coeffs = size_table[a.order+1] - + ct = deepcopy(coeff_table[a.order+1]) @inbounds for i = 1:num_coeffs - iind = @isonethread coeff_table[a.order+1][i] + # iind = @isonethread coeff_table[a.order+1][i] + iind = ct[i] n = iind[r] n == 0 && continue iind[r] -= 1 @@ -370,9 +371,10 @@ function integrate(a::HomogeneousPolynomial, r::Int) T = promote_type(TS.numtype(a), TS.numtype(a[1]/1)) coeffs = zeros(T, size_table[a.order+2]) - + ct = deepcopy(coeff_table[a.order+1]) @inbounds for i = 1:num_coeffs - iind = @isonethread coeff_table[a.order+1][i] + # iind = @isonethread coeff_table[a.order+1][i] + iind = ct[i] n = iind[r] n == order_max && continue iind[r] += 1 diff --git a/src/evaluate.jl b/src/evaluate.jl index 0cca2983..f65e4c9e 100644 --- a/src/evaluate.jl +++ b/src/evaluate.jl @@ -17,17 +17,19 @@ equivalent to `evaluate(a,dx)`, and `a()` is equivalent to `evaluate(a)`. function evaluate(a::Taylor1{T}, dx::T) where {T<:Number} @inbounds suma = zero(a[end]) @inbounds for k in reverse(eachindex(a)) - suma = suma*dx + a[k] + suma = suma * dx + a[k] end return suma end + function evaluate(a::Taylor1{T}, dx::S) where {T<:Number, S<:Number} suma = a[end]*zero(dx) @inbounds for k in reverse(eachindex(a)) - suma = suma*dx + a[k] + suma = suma * dx + a[k] end return suma end + evaluate(a::Taylor1{T}) where {T<:Number} = a[0] @@ -41,6 +43,7 @@ is equivalent to `evaluate(x)`. """ evaluate(x::AbstractArray{Taylor1{T}}, δt::S) where {T<:Number, S<:Number} = evaluate.(x, δt) + evaluate(a::AbstractArray{Taylor1{T}}) where {T<:Number} = getcoeff.(a, 0) @@ -58,62 +61,89 @@ function evaluate(a::Taylor1{T}, x::Taylor1{T}) where {T<:Number} a, x = fixorder(a, x) end @inbounds suma = a[end]*zero(x) - @inbounds for k in reverse(eachindex(a)) - suma = suma*x + a[k] - end + aux = zero(suma) + _horner!(suma, a, x, aux) return suma end function evaluate(a::Taylor1{Taylor1{T}}, x::Taylor1{T}) where {T<:NumberNotSeriesN} - @inbounds suma = a[end]*zero(x) - @inbounds for k in reverse(eachindex(a)) - suma = suma*x + a[k] + if a.order != x.order + a, x = fixorder(a, x) end + @inbounds suma = a[end]*zero(x) + aux = zero(suma) + _horner!(suma, a, x, aux) return suma end + function evaluate(a::Taylor1{T}, x::Taylor1{Taylor1{T}}) where {T<:NumberNotSeriesN} - @inbounds suma = a[end]*zero(x) - @inbounds for k in reverse(eachindex(a)) - suma = suma*x + a[k] + if a.order != x.order + a, x = fixorder(a, x) end + @inbounds suma = a[end]*zero(x) + aux = zero(suma) + _horner!(suma, a, x, aux) return suma end evaluate(p::Taylor1{T}, x::AbstractArray{S}) where {T<:Number, S<:Number} = evaluate.(Ref(p), x) -function evaluate(a::Taylor1{TaylorN{T}}, dx::S) where - {T<:NumberNotSeries, S<:NumberNotSeries} - @inbounds suma = TaylorN( zero(T)*constant_term(dx), a[0].order ) - @inbounds for k in reverse(eachindex(a)) - for ordQ in eachindex(a[k]) - mul!(suma, suma, dx, ordQ) - add!(suma, suma, a[k], ordQ) - end +# Substitute a TaylorN into a Taylor1 +function evaluate(a::Taylor1{T}, dx::TaylorN{T}) where {T<:NumberNotSeries} + suma = TaylorN( zero(T), dx.order) + aux = TaylorN( zero(T), dx.order) + _horner!(suma, a, dx, aux) + return suma +end + +function evaluate(a::Taylor1{T}, dx::Taylor1{TaylorN{T}}) where {T<:NumberNotSeries} + if a.order != dx.order + a, dx = fixorder(a, dx) end + suma = Taylor1( zero(dx[0]), a.order) + aux = Taylor1( zero(dx[0]), a.order) + _horner!(suma, a, dx, aux) return suma end -function evaluate(a::Taylor1{T}, dx::TaylorN{T}) where {T<:NumberNotSeries} - @inbounds suma = TaylorN( zero(T), get_order(dx) ) - @inbounds for k in reverse(eachindex(a)) - suma = suma*dx + a[k] + +# Evaluate a Taylor1{TaylorN{T}} on Vector{TaylorN} is interpreted +# as a substitution on the TaylorN vars +function evaluate(a::Taylor1{TaylorN{T}}, dx::Vector{TaylorN{T}}) where {T<:NumberNotSeries} + @assert length(dx) == get_numvars() + suma = Taylor1( zero(a[0]), a.order) + suma.coeffs .= evaluate.(a[:], Ref(dx)) + return suma +end + +function evaluate(a::Taylor1{TaylorN{T}}, ind::Int, dx::T) where {T<:NumberNotSeries} + @assert (1 ≤ ind ≤ get_numvars()) "Invalid `ind`; it must be between 1 and `get_numvars()`" + suma = Taylor1( zero(a[0]), a.order) + for ord in eachindex(suma) + for ordQ in eachindex(a[0]) + _evaluate!(suma[ord], a[ord][ordQ], ind, dx) + end end return suma end -function evaluate(a::Taylor1{TaylorN{T}}, dx::TaylorN{T}) where {T<:NumberNotSeries} - order = min(minimum(get_order.(a[:])), get_order(dx)) - @inbounds suma = TaylorN( zero(T), order ) - @inbounds for k in reverse(eachindex(a)) - suma = suma*dx + a[k] +function evaluate(a::Taylor1{TaylorN{T}}, ind::Int, dx::TaylorN{T}) where {T<:NumberNotSeries} + @assert (1 ≤ ind ≤ get_numvars()) "Invalid `ind`; it must be between 1 and `get_numvars()`" + suma = Taylor1( zero(a[0]), a.order) + aux = zero(dx) + for ord in eachindex(suma) + for ordQ in eachindex(a[0]) + _evaluate!(suma[ord], a[ord][ordQ], ind, dx, aux) + end end return suma end + #function-like behavior for Taylor1 (p::Taylor1)(x) = evaluate(p, x) -(p::Taylor1)() = evaluate(p) +(p::Taylor1)() = evaluate(p) #function-like behavior for Vector{Taylor1} (asumes Julia version >= 1.6) (p::Array{Taylor1{T}})(x) where {T<:Number} = evaluate.(p, x) @@ -153,14 +183,62 @@ function _evaluate(a::HomogeneousPolynomial{T}, vals::NTuple) where {T} a.order == 0 && return a[1]*one(vals[1]) ct = coeff_table[a.order+1] suma = zero(a[1])*vals[1] + vv = vals .^ ct[1] for (i, a_coeff) in enumerate(a.coeffs) iszero(a_coeff) && continue - @inbounds tmp = prod( vals .^ ct[i] ) + @inbounds vv .= vals .^ ct[i] + tmp = prod( vv ) suma += a_coeff * tmp end return suma end +function _evaluate(a::HomogeneousPolynomial{T}, vals::NTuple{N,<:TaylorN{T}}) where + {N,T<:NumberNotSeries} + # @assert length(vals) == get_numvars() + a.order == 0 && return a[1]*one(vals[1]) + ct = coeff_table[a.order+1] + suma = TaylorN(zero(T), vals[1].order) + # + # vv = power_by_squaring.(vals, ct[1]) + vv = vals .^ ct[1] + tmp = zero(suma) + aux = one(suma) + for (i, a_coeff) in enumerate(a.coeffs) + iszero(a_coeff) && continue + # @inbounds vv .= power_by_squaring.(vals, ct[i]) + vv .= vals .^ ct[i] + # tmp = prod( vv ) + for ord in eachindex(tmp) + @inbounds one!(aux, vv[1], ord) + end + for j in eachindex(vv) + for ord in eachindex(tmp) + zero!(tmp, ord) + @inbounds mul!(tmp, aux, vv[j], ord) + end + for ord in eachindex(tmp) + identity!(aux, tmp, ord) + end + end + # suma += a_coeff * tmp + for ord in eachindex(tmp) + for ordQ in eachindex(tmp[ord]) + zero!(aux[ord], ordQ) + aux[ord][ordQ] = a_coeff * tmp[ord][ordQ] + suma[ord][ordQ] += aux[ord][ordQ] + end + end + end + return suma +end + +function _evaluate(a::HomogeneousPolynomial{T}, ind::Int, val::T) where {T<:NumberNotSeries} + suma = TaylorN(zero(T), get_order()) + _evaluate!(suma, a, ind, val) + return suma +end + #function-like behavior for HomogeneousPolynomial (p::HomogeneousPolynomial)(x) = evaluate(p, x) @@ -209,23 +287,44 @@ evaluate(a::TaylorN{Taylor1{T}}, vals::AbstractVector{S}; function evaluate(a::TaylorN{T}, s::Symbol, val::S) where {T<:Number, S<:NumberNotSeriesN} - vars = get_variables(T) ind = lookupvar(s) - @inbounds vars[ind] = val + zero(vars[ind]) - return _evaluate(a, (vars...,), Val(false)) + @assert (1 ≤ ind ≤ get_numvars()) "Symbol is not a TaylorN variable; see `get_variable_names()`" + return evaluate(a, ind, val) end -evaluate(a::TaylorN{T}, x::Pair{Symbol,S}) where {T<:NumberNotSeries, S} = - evaluate(a, first(x), last(x)) +function evaluate(a::TaylorN{T}, ind::Int, val::S) where + {T<:Number, S<:NumberNotSeriesN} + @assert (1 ≤ ind ≤ get_numvars()) "Invalid `ind`; it must be between 1 and `get_numvars()`" + R = promote_type(T,S) + return _evaluate(convert(TaylorN{R}, a), ind, convert(R, val)) +end + +function evaluate(a::TaylorN{T}, s::Symbol, val::TaylorN) where {T<:Number} + ind = lookupvar(s) + @assert (1 ≤ ind ≤ get_numvars()) "Symbol is not a TaylorN variable; see `get_variable_names()`" + return evaluate(a, ind, val) +end + +function evaluate(a::TaylorN{T}, ind::Int, val::TaylorN) where {T<:Number} + @assert (1 ≤ ind ≤ get_numvars()) "Invalid `ind`; it must be between 1 and `get_numvars()`" + a, val = fixorder(a, val) + a, val = promote(a, val) + return _evaluate(a, ind, val) +end -evaluate(a::TaylorN{Taylor1{T}}, x::Pair{Symbol,S}) where {T<:NumberNotSeries, S} = +evaluate(a::TaylorN{T}, x::Pair{Symbol,S}) where {T, S} = evaluate(a, first(x), last(x)) evaluate(a::TaylorN{T}) where {T<:Number} = constant_term(a) # _evaluate -# Returns a vector with the evaluation of the HomogeneousPolynomials +_evaluate(a::TaylorN{T}, vals::NTuple, ::Val{true}) where {T<:NumberNotSeries} = + sum( sort!(_evaluate(a, vals), by=abs2) ) + +_evaluate(a::TaylorN{T}, vals::NTuple, ::Val{false}) where {T<:Number} = + sum( _evaluate(a, vals) ) + function _evaluate(a::TaylorN{T}, vals::NTuple{N,<:Number}) where {N,T<:Number} R = promote_type(T, typeof(vals[1])) a_length = length(a) @@ -246,11 +345,89 @@ function _evaluate(a::TaylorN{T}, vals::NTuple{N,<:TaylorN}) where {N,T<:Number} return suma end -_evaluate(a::TaylorN{T}, vals::NTuple, ::Val{true}) where {T<:NumberNotSeries} = - sum( sort!(_evaluate(a, vals), by=abs2) ) -_evaluate(a::TaylorN{T}, vals::NTuple, ::Val{false}) where {T<:Number} = - sum( _evaluate(a, vals) ) +function _evaluate(a::TaylorN{T}, ind::Int, val::T) where {T<:NumberNotSeriesN} + suma = TaylorN(zero(a[0]*val), a.order) + vval = convert(numtype(suma), val) + suma, a = promote(suma, a) + @inbounds for ordQ in eachindex(a) + _evaluate!(suma, a[ordQ], ind, vval) + end + return suma +end + +function _evaluate(a::TaylorN{T}, ind::Int, val::TaylorN{T}) where {T<:NumberNotSeriesN} + suma = TaylorN(zero(a[0]), a.order) + aux = zero(suma) + @inbounds for ordQ in eachindex(a) + _evaluate!(suma, a[ordQ], ind, val, aux) + end + return suma +end + +function _evaluate!(suma::TaylorN{T}, a::HomogeneousPolynomial{T}, ind::Int, val::T) where + {T<:NumberNotSeriesN} + order = a.order + if order == 0 + suma[0] = a[1]*one(val) + return nothing + end + vv = val .^ (0:order) + # ct = @isonethread coeff_table[order+1] + ct = deepcopy(coeff_table[order+1]) + for (i, a_coeff) in enumerate(a.coeffs) + iszero(a_coeff) && continue + if ct[i][ind] == 0 + suma[order][i] += a_coeff + continue + end + vpow = ct[i][ind] + red_order = order - vpow + ct[i][ind] -= vpow + kdic = in_base(get_order(), ct[i]) + ct[i][ind] += vpow + pos = pos_table[red_order+1][kdic] + suma[red_order][pos] += a_coeff * vv[vpow+1] + end + return nothing +end + +function _evaluate!(suma::TaylorN{T}, a::HomogeneousPolynomial{T}, ind::Int, + val::TaylorN{T}, aux::TaylorN{T}) where {T<:NumberNotSeriesN} + order = a.order + if order == 0 + suma[0] = a[1] + return nothing + end + vv = zero(suma) + ct = coeff_table[order+1] + za = zero(a) + for (i, a_coeff) in enumerate(a.coeffs) + iszero(a_coeff) && continue + if ct[i][ind] == 0 + suma[order][i] += a_coeff + continue + end + za[i] = a_coeff + zero!(aux) + _evaluate!(aux, za, ind, one(T)) + za[i] = zero(T) + vpow = ct[i][ind] + # vv = val ^ vpow + if constant_term(val) == 0 + vv = val ^ vpow + else + for ordQ in eachindex(val) + zero!(vv, ordQ) + pow!(vv, val, vpow, ordQ) + end + end + for ordQ in eachindex(suma) + mul!(suma, vv, aux, ordQ) + end + end + return nothing +end #High-dim array evaluation @@ -269,7 +446,7 @@ evaluate(A::AbstractArray{TaylorN{T}}) where {T<:Number} = evaluate.(A) #function-like behavior for TaylorN (p::TaylorN)(x) = evaluate(p, x) (p::TaylorN)() = evaluate(p) -(p::TaylorN)(s::Symbol, x) = evaluate(p, s, x) +(p::TaylorN)(s::S, x) where {S<:Union{Symbol, Int}}= evaluate(p, s, x) (p::TaylorN)(x::Pair) = evaluate(p, first(x), last(x)) (p::TaylorN)(x, v::Vararg{T}) where {T} = evaluate(p, (x, v...,)) (p::TaylorN)(b::Bool, x) = evaluate(p, x, sorting=b) @@ -300,16 +477,88 @@ end ## In place evaluation of multivariable arrays function evaluate!(x::AbstractArray{TaylorN{T}}, δx::Array{T,1}, x0::AbstractArray{T}) where {T<:Number} - @inbounds for i in eachindex(x, x0) - x0[i] = evaluate( x[i], δx ) - end + x0 .= evaluate.( x, Ref(δx) ) return nothing end function evaluate!(x::AbstractArray{TaylorN{T}}, δx::Array{TaylorN{T},1}, x0::AbstractArray{TaylorN{T}}; sorting::Bool=true) where {T<:NumberNotSeriesN} - @inbounds for i in eachindex(x, x0) - x0[i] = _evaluate( x[i], δx, Val(sorting) ) + x0 .= _evaluate.( x, Ref(δx), Ref(Val(sorting)) ) + return nothing +end + + +# In-place Horner methods, used when the result of an evaluation (substitution) +# is Taylor1{} +function _horner!(suma::Taylor1{T}, a::Taylor1{T}, x::Taylor1{T}, + aux::Taylor1{T}) where {T<:Number} + @inbounds for k in reverse(eachindex(a)) + for ord in eachindex(aux) + mul!(aux, suma, x, ord) + end + for ord in eachindex(aux) + identity!(suma, aux, ord) + end + add!(suma, suma, a[k], 0) + end + return nothing +end + +function _horner!(suma::Taylor1{T}, a::Taylor1{Taylor1{T}}, x::Taylor1{T}, + aux::Taylor1{T}) where {T<:Number} + @inbounds for k in reverse(eachindex(a)) + for ord in eachindex(aux) + mul!(aux, suma, x, ord) + end + for ord in eachindex(aux) + identity!(suma, aux, ord) + add!(suma, suma, a[k], ord) + end end return nothing end + +function _horner!(suma::Taylor1{Taylor1{T}}, a::Taylor1{T}, x::Taylor1{Taylor1{T}}, + aux::Taylor1{Taylor1{T}}) where {T<:Number} + @inbounds for k in reverse(eachindex(a)) + for ord in eachindex(aux) + mul!(aux, suma, x, ord) + end + for ord in eachindex(aux) + identity!(suma, aux, ord) + add!(suma, suma, a[k], ord) + end + end + return nothing +end + +function _horner!(suma::TaylorN{T}, a::Taylor1{T}, dx::TaylorN{T}, + aux::TaylorN{T}) where {T<:NumberNotSeries} + @inbounds for k in reverse(eachindex(a)) + for ordQ in eachindex(suma) + zero!(aux, ordQ) + mul!(aux, suma, dx, ordQ) + end + for ordQ in eachindex(suma) + identity!(suma, aux, ordQ) + end + add!(suma, suma, a[k], 0) + end + return nothing +end + +function _horner!(suma::Taylor1{TaylorN{T}}, a::Taylor1{T}, dx::Taylor1{TaylorN{T}}, + aux::Taylor1{TaylorN{T}}) where {T<:NumberNotSeries} + @inbounds for k in reverse(eachindex(a)) + # aux = suma * dx + for ord in eachindex(aux) + zero!(aux, ord) + mul!(aux, suma, dx, ord) + end + for ord in eachindex(aux) + identity!(suma, aux, ord) + end + add!(suma, suma, a[k], 0) + end + return suma +end \ No newline at end of file diff --git a/src/power.jl b/src/power.jl index 9dc98acc..7e9788d0 100644 --- a/src/power.jl +++ b/src/power.jl @@ -64,16 +64,18 @@ end # Licensed under MIT "Expat" for T in (:Taylor1, :HomogeneousPolynomial, :TaylorN) @eval function power_by_squaring(x::$T, p::Integer) - p == 1 && return copy(x) - p == 0 && return one(x) - p == 2 && return square(x) + if p == 1 + return copy(x) + elseif p == 0 + return one(x) + elseif p == 2 + return square(x) + end t = trailing_zeros(p) + 1 p >>= t - while (t -= 1) > 0 x = square(x) end - y = x while p > 0 t = trailing_zeros(p) + 1 @@ -83,7 +85,6 @@ for T in (:Taylor1, :HomogeneousPolynomial, :TaylorN) end y *= x end - return y end end @@ -338,8 +339,8 @@ Return `a^2`; see [`TaylorSeries.sqr!`](@ref). for T in (:Taylor1, :TaylorN) @eval function square(a::$T) - c = $T( constant_term(a)^2, a.order) - for k in 1:a.order + c = $T( zero(constant_term(a)), a.order) + for k in eachindex(a) sqr!(c, a, k) end return c @@ -365,19 +366,29 @@ end #auxiliary function to avoid allocations @inline function sqr_orderzero!(c::Taylor1{T}, a::Taylor1{T}) where {T<:NumberNotSeries} - @inbounds c[0] = constant_term(a)^2 + @inbounds c[0] = a[0]^2 return nothing end @inline function sqr_orderzero!(c::TaylorN{T}, a::TaylorN{T}) where {T<:NumberNotSeries} - @inbounds c[0][1] = constant_term(a)^2 + @inbounds c[0][1] = a[0][1]^2 return nothing end @inline function sqr_orderzero!(c::Taylor1{TaylorN{T}}, a::Taylor1{TaylorN{T}}) where {T<:NumberNotSeries} - @inbounds c[0][0][1] = constant_term(a)^2 + @inbounds for ord in eachindex(c[0]) + sqr!(c[0], a[0], ord) + end return nothing end @inline function sqr_orderzero!(c::TaylorN{Taylor1{T}}, a::TaylorN{Taylor1{T}}) where {T<:NumberNotSeries} - @inbounds c[0][1][0] = a[0][1][0]^2 + @inbounds for ord in eachindex(c[0][1]) + sqr!(c[0][1], a[0][1], ord) + end + return nothing +end +@inline function sqr_orderzero!(c::Taylor1{Taylor1{T}}, a::Taylor1{Taylor1{T}}) where {T<:NumberNotSeriesN} + @inbounds for ord in eachindex(c[0]) + sqr!(c[0], a[0], ord) + end return nothing end diff --git a/test/manyvariables.jl b/test/manyvariables.jl index 40122b27..1b7f243b 100644 --- a/test/manyvariables.jl +++ b/test/manyvariables.jl @@ -77,7 +77,7 @@ end @test eachindex(x) == 0:6 @test iterate(x) == (HomogeneousPolynomial([0.0], 0), 1) @test iterate(y, 1) == (HomogeneousPolynomial([0.0, 1.0], 1), 2) - @test iterate(x, 7) == nothing + @test isnothing(iterate(x, 7)) @test x.order == 6 @test TS.name_taylorNvar(1) == " x" @@ -458,10 +458,13 @@ end @test ptxy(:x₁ => -1.0) == -1 + yT + (-1.0+yT^3)/3 + yT - yT^2 @test evaluate(ptxy, :x₁ => -1.0) == -1 + yT + (-1.0+yT^3)/3 + yT - yT^2 @test evaluate(ptxy, :x₁, -1.0) == -1 + yT + (-1.0+yT^3)/3 + yT - yT^2 + @test isa(evaluate(ptxy, :x₁, 1), TaylorN{Float64}) + @test evaluate(ptxy, :x₁, xT) == ptxy + @test evaluate(ptxy, 1, 1-yT) ≈ 4/3 + zero(yT) v = zeros(Int, 2) - @test evaluate!([xT, yT], ones(Int, 2), v) == nothing + @test isnothing(evaluate!([xT, yT], ones(Int, 2), v)) @test v == ones(2) - @test evaluate!([xT, yT][1:2], ones(Int, 2), v) == nothing + @test isnothing(evaluate!([xT, yT][1:2], ones(Int, 2), v)) @test v == ones(2) A_TN = [xT 2xT 3xT; yT 2yT 3yT] @test evaluate(A_TN, ones(2)) == [1.0 2.0 3.0; 1.0 2.0 3.0] @@ -813,7 +816,7 @@ end end -@testset "Consistency of coeff_table (differentiate)" begin +@testset "Consistency of coeff_table" begin order = 20 x, y, z, w = set_variables(Int128, "x y z w", numvars=4, order=2order) ctab = deepcopy(TS.coeff_table); @@ -840,8 +843,7 @@ end end f = fun(order); - df_exact = 4*20*fun(order-1); - fint_exact = 4*fun(order+1)/(order+1); + df_exact = 4*order*fun(order-1); df1 = diffs1(f); @test ctab == TS.coeff_table @@ -851,28 +853,48 @@ end @test ctab == TS.coeff_table @test df1 == df_exact - # function integ1(f) - # f1 = integrate(f, 1) - # f2 = integrate(f, 2) - # f3 = integrate(f, 3) - # f4 = integrate(f, 4) - # return f1 + f2 + f3 + f4 - # end - # - # function integ2(f) - # vder = zeros(typeof(f), get_numvars()) - # Threads.@threads for i = 1:get_numvars() - # vder[i] = integrate(f, i) - # end - # return sum(vder) - # end - # - # ii1 = integ1(f); - # @test ctab == TS.coeff_table - # @test ii1 == fint_exact - # - # ii2 = integ2(f); - # @test ctab == TS.coeff_table - # @test ii2 == fint_exact + function integ1(f) + f1 = integrate(f, 1) + f2 = integrate(f, 2) + f3 = integrate(f, 3) + f4 = integrate(f, 4) + return f1 + f2 + f3 + f4 + end + + function integ2(f) + T = typeof(integrate(f, 1)) + vinteg = zeros(T, get_numvars()) + Threads.@threads for i = 1:get_numvars() + vinteg[i] = integrate(f, i) + end + return sum(vinteg) + end + + ii1 = integ1(f); + ii2 = integ2(f); + @test ctab == TS.coeff_table + @test ii1 == ii2 + + function ev1(f) + f1 = f(1, 1.0) + f2 = f(2, 1.0) + f3 = f(3, 1.0) + f4 = f(4, 1.0) + return f1 + f2 + f3 + f4 + end + + function ev2(f) + T = typeof(evaluate(f, 1, 1.0)) + veval = zeros(T, get_numvars()) + Threads.@threads for i = 1:get_numvars() + veval[i] = evaluate(f, i, 1.0) + end + return sum(veval) + end + + ee1 = ev1(f); + ee2 = ev2(f); + @test ctab == TS.coeff_table + @test ee1 == ee2 end diff --git a/test/mixtures.jl b/test/mixtures.jl index 587889ab..1bb24412 100644 --- a/test/mixtures.jl +++ b/test/mixtures.jl @@ -162,6 +162,12 @@ using Test @test v == [TaylorN(1), TaylorN(1)^2] @test tN1() == t @test evaluate(tN1, :x₁ => 1.0) == TaylorN([HomogeneousPolynomial([1.0+t]), zero(xHt), yHt^2]) + @test evaluate(tN1, 1, 1.0) == TaylorN([HomogeneousPolynomial([1.0+t]), zero(xHt), yHt^2]) + @test evaluate(t, t1N) == t1N + @test evaluate(t1N, 0.5) == t1N[0] + t1N[1]/2 + t1N[2]/4 + @test evaluate(t1N, [t1N[0], zero(t1N[0])]) == Taylor1([t1N[0], t1N[1]], t.order) + @test evaluate(t1N, 2, 0.0) == Taylor1([t1N[0], t1N[1]], t.order) + @test evaluate(t1N, 1, 0.0) == Taylor1([zero(t1N[0]), t1N[1], t1N[2]], t.order) # Tests for functions of mixtures t1N = Taylor1([zero(TaylorN(Float64,1)), one(TaylorN(Float64,1))], 6) @@ -326,6 +332,14 @@ using Test @test tanh(atanh(xx*tt)) ≈ xx * tt @test_throws DomainError atanh(xx+tt) + # pp = xx*yy*(1+tt)^4 + @test evaluate(pp, 1, 0.0) == yy*(1+tt)^4 + @test evaluate(pp, 2, 0.0) == xx*(1+tt)^4 + @test evaluate(t, tt) == tt + @test evaluate(tt, t) == tt + @test evaluate(xx, 2, δy) == xx + @test evaluate(xx, 1, δy) == yy + #testing evaluate and function-like behavior of Taylor1, TaylorN for mixtures: t = Taylor1(25) p = cos(t) diff --git a/test/runtests.jl b/test/runtests.jl index 5991b446..ef447b14 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -13,11 +13,11 @@ testfiles = ( "fateman40.jl", "staticarrays.jl", "jld2.jl", - "rat.jl" + "rat.jl", + # Run Aqua tests at the very end + "aqua.jl", ) for file in testfiles include(file) end -# After `using intervalArithmetic` new ambiguities may arise -include("aqua.jl")