Merge branch 'master' into lb/IA_1.0-dev

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.7"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
@@ -11,26 +11,37 @@ SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 
 [weakdeps]
 IntervalArithmetic = "d1acc4aa-44c8-5952-acd4-ba5d80a2a253"
+JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
+StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
+RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd"
 
 [extensions]
 TaylorSeriesIAExt = "IntervalArithmetic"
+TaylorSeriesJLD2Ext = "JLD2"
+TaylorSeriesSAExt = "StaticArrays"
+TaylorSeriesRATExt = "RecursiveArrayTools"
 
 [compat]
 Aqua = "0.8"
 IntervalArithmetic = "0.22"
 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"
 Test = "<0.0.1, 1"
 julia = "1.9"
 
 [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"
+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", "LinearAlgebra", "SparseArrays", "Test", "Aqua"]
+test = ["IntervalArithmetic", "JLD2", "LinearAlgebra", "RecursiveArrayTools", "SparseArrays", "StaticArrays", "Test", "Aqua"]

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

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

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

src/TaylorSeries.jl

@@ -45,6 +45,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,
@@ -80,6 +82,12 @@ 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
+        @require JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819" begin
+            include("../ext/TaylorSeriesJLD2Ext.jl")
+        end
     end
 end
 

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}
@@ -258,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)
@@ -402,4 +403,3 @@ macro isonethread(expr)
         end
     end)
 end
-

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

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())

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)