Add numtype to replace eltype (#289)

* Add numtype (not exported) to replace eltype

* Add numtype docstrings and add it to the docs

* Add more tests

* Export normalize_taylor from intervals.jl

* Add docstrings of nonlinear_polynomial to docs

* Bump patch version, and upgrade some compat entries

Project.toml

@@ -1,7 +1,7 @@
 name = "TaylorSeries"
 uuid = "6aa5eb33-94cf-58f4-a9d0-e4b2c4fc25ea"
 repo = "https://github.com/JuliaDiff/TaylorSeries.jl.git"
-version = "0.11.2"
+version = "0.11.3"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
@@ -10,7 +10,7 @@ Requires = "ae029012-a4dd-5104-9daa-d747884805df"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 
 [compat]
-IntervalArithmetic = "0.15, 0.16, 0.17, 0.18"
+IntervalArithmetic = "0.15, 0.16, 0.17, 0.18, 0.19, 0.20"
 Requires = "0.5.2, 1.0, 1.1"
 julia = "1"
 

docs/src/api.md

@@ -45,6 +45,7 @@ hessian
 hessian!
 constant_term
 linear_polynomial
+nonlinear_polynomial
 inverse
 abs
 norm
@@ -67,6 +68,7 @@ in_base
 make_inverse_dict
 resize_coeffs1!
 resize_coeffsHP!
+numtype
 mul!
 mul!(::HomogeneousPolynomial, ::HomogeneousPolynomial, ::HomogeneousPolynomial)
 mul!(::Vector{Taylor1{T}}, ::Union{Matrix{T},SparseMatrixCSC{T}},::Vector{Taylor1{T}}) where {T<:Number}

docs/src/index.md

@@ -11,6 +11,11 @@ A [Julia](http://julialang.org) package for Taylor expansions in one or more ind
 - [David P. Sanders](http://sistemas.fciencias.unam.mx/~dsanders/), Facultad
     de Ciencias, Universidad Nacional Autónoma de México (UNAM).
 
+### Citing
+
+If you find useful this package, please cite the paper:
+
+Benet, L., & Sanders, D. P. (2019). TaylorSeries.jl: Taylor expansions in one and several variables in Julia. Journal of Open Source Software, 4(36), 1–4. https://doi.org/10.5281/zenodo.2601941
 
 ### License
 

src/TaylorSeries.jl

@@ -45,7 +45,7 @@ import Base: zero, one, zeros, ones, isinf, isnan, iszero,
     power_by_squaring,
     rtoldefault, isfinite, isapprox, rad2deg, deg2rad
 
-export Taylor1, TaylorN, HomogeneousPolynomial, AbstractSeries
+export Taylor1, TaylorN, HomogeneousPolynomial, AbstractSeries, TS
 
 export getcoeff, derivative, integrate, differentiate,
     evaluate, evaluate!, inverse, set_taylor1_varname,
@@ -54,9 +54,10 @@ export getcoeff, derivative, integrate, differentiate,
     set_variables, get_variables,
     get_variable_names, get_variable_symbols,
     # jacobian, hessian, jacobian!, hessian!,
-    ∇, taylor_expand, update!, 
-    constant_term, linear_polynomial, nonlinear_polynomial,
-    normalize_taylor
+    ∇, taylor_expand, update!,
+    constant_term, linear_polynomial, nonlinear_polynomial
+
+const TS = TaylorSeries
 
 include("parameters.jl")
 include("hash_tables.jl")

src/arithmetic.jl

@@ -172,7 +172,7 @@ for (f, fc) in ((:+, :(add!)), (:-, :(subst!)))
         function ($f)(a::TaylorN{Taylor1{T}}, b::Taylor1{S}) where
                 {T<:NumberNotSeries,S<:NumberNotSeries}
             @inbounds aux = $f(a[0][1], b)
-            R = eltype(aux)
+            R = TS.numtype(aux)
             coeffs = Array{HomogeneousPolynomial{Taylor1{R}}}(undef, a.order+1)
             coeffs .= a.coeffs
             @inbounds coeffs[1] = aux
@@ -182,7 +182,7 @@ for (f, fc) in ((:+, :(add!)), (:-, :(subst!)))
         function ($f)(b::Taylor1{S}, a::TaylorN{Taylor1{T}}) where
                 {T<:NumberNotSeries,S<:NumberNotSeries}
             @inbounds aux = $f(b, a[0][1])
-            R = eltype(aux)
+            R = TS.numtype(aux)
             coeffs = Array{HomogeneousPolynomial{Taylor1{R}}}(undef, a.order+1)
             @__dot__ coeffs = $f(a.coeffs)
             @inbounds coeffs[1] = aux
@@ -192,7 +192,7 @@ for (f, fc) in ((:+, :(add!)), (:-, :(subst!)))
         function ($f)(a::Taylor1{TaylorN{T}}, b::TaylorN{S}) where
                 {T<:NumberNotSeries,S<:NumberNotSeries}
             @inbounds aux = $f(a[0], b)
-            R = eltype(aux)
+            R = TS.numtype(aux)
             coeffs = Array{TaylorN{R}}(undef, a.order+1)
             coeffs .= a.coeffs
             @inbounds coeffs[1] = aux
@@ -202,7 +202,7 @@ for (f, fc) in ((:+, :(add!)), (:-, :(subst!)))
         function ($f)(b::TaylorN{S}, a::Taylor1{TaylorN{T}}) where
                 {T<:NumberNotSeries,S<:NumberNotSeries}
             @inbounds aux = $f(b, a[0])
-            R = eltype(aux)
+            R = TS.numtype(aux)
             coeffs = Array{TaylorN{R}}(undef, a.order+1)
             @__dot__ coeffs = $f(a.coeffs)
             @inbounds coeffs[1] = aux

src/auxiliary.jl

@@ -229,12 +229,19 @@ for T in (:Taylor1, :HomogeneousPolynomial, :TaylorN)
             @inline lastindex(a::$T) = a.order
         end
         @inline eachindex(a::$T) = firstindex(a):lastindex(a)
-        @inline eltype(::$T{S}) where {S<:Number} = S
+        @inline numtype(::$T{S}) where {S<:Number} = S
         @inline size(a::$T) = size(a.coeffs)
         @inline get_order(a::$T) = a.order
         @inline axes(a::$T) = ()
     end
 end
+numtype(a) = eltype(a)
+
+@doc doc"""
+    numtype(a::AbstractSeries)
+
+Return the type of the elements of the coefficients of `a`.
+""" numtype
 
 
 ## fixorder ##
@@ -305,7 +312,7 @@ constant_term(a::Number) = a
 """
     linear_polynomial(a)
 
-Returns the linear part of `a` as a polynomial (`Taylor1` or `TaylorN`), 
+Returns the linear part of `a` as a polynomial (`Taylor1` or `TaylorN`),
 *without* the constant term. The fallback behavior is to return `a` itself.
 """
 linear_polynomial(a::Taylor1) = Taylor1([zero(a[1]), a[1]], a.order)

src/calculus.jl

@@ -34,7 +34,7 @@ const derivative = differentiate
     differentiate!(res, a) --> nothing
 
 In-place version of `differentiate`. Compute the `Taylor1` polynomial of the
-differential of `a::Taylor1` and return it as `res` (order of `res` remains 
+differential of `a::Taylor1` and return it as `res` (order of `res` remains
 unchanged).
 """
 function differentiate!(res::Taylor1, a::Taylor1)
@@ -127,8 +127,7 @@ to the `r`-th variable.
 """
 function differentiate(a::HomogeneousPolynomial, r::Int)
     @assert 1 ≤ r ≤ get_numvars()
-    # @show zero(a[1]) zero(eltype(a))
-    T = eltype(a)
+    T = TS.numtype(a)
     a.order == 0 && return HomogeneousPolynomial(zero(a[1]), 0)
     @inbounds num_coeffs = size_table[a.order]
     coeffs = zeros(T,num_coeffs)
@@ -158,7 +157,7 @@ to the `r`-th variable. The `r`-th variable may be also
 specified through its symbol.
 """
 function differentiate(a::TaylorN, r=1::Int)
-    T = eltype(a)
+    T = TS.numtype(a)
     coeffs = Array{HomogeneousPolynomial{T}}(undef, a.order)
 
     @inbounds for ord = 1:a.order
@@ -223,7 +222,7 @@ end
 Compute the gradient of the polynomial `f::TaylorN`.
 """
 function gradient(f::TaylorN)
-    T = eltype(f)
+    T = TS.numtype(f)
     numVars = get_numvars()
     grad = Array{TaylorN{T}}(undef, numVars)
     @inbounds for nv = 1:numVars
@@ -361,7 +360,7 @@ function integrate(a::HomogeneousPolynomial, r::Int)
     @inbounds posTb = pos_table[a.order+2]
     @inbounds num_coeffs = size_table[a.order+1]
 
-    T = promote_type(eltype(a), eltype(a[1]/1))
+    T = promote_type(TS.numtype(a), TS.numtype(a[1]/1))
     coeffs = zeros(T, size_table[a.order+2])
 
     @inbounds for i = 1:num_coeffs
@@ -388,7 +387,7 @@ where `x0` the integration constant and must be independent
 of the `r`-th variable; if `x0` is ommitted, it is taken as zero.
 """
 function integrate(a::TaylorN, r::Int)
-    T = promote_type(eltype(a), eltype(a[0]/1))
+    T = promote_type(TS.numtype(a), TS.numtype(a[0]/1))
     order_max = min(get_order(), a.order+1)
     coeffs = zeros(HomogeneousPolynomial{T}, order_max)
 
@@ -408,7 +407,7 @@ function integrate(a::TaylorN, r::Int, x0::TaylorN)
     return x0+res
 end
 integrate(a::TaylorN, r::Int, x0) =
-    integrate(a,r,TaylorN(HomogeneousPolynomial([convert(eltype(a),x0)], 0)))
+    integrate(a,r,TaylorN(HomogeneousPolynomial([convert(TS.numtype(a),x0)], 0)))
 
 integrate(a::TaylorN, s::Symbol) = integrate(a, lookupvar(s))
 integrate(a::TaylorN, s::Symbol, x0::TaylorN) = integrate(a, lookupvar(s), x0)

src/functions.jl

@@ -22,7 +22,7 @@ for T in (:Taylor1, :TaylorN)
         end
 
         function log(a::$T)
-            iszero(constant_term(a)) && throw(DomainError(a, 
+            iszero(constant_term(a)) && throw(DomainError(a,
                     """The 0-th order coefficient must be non-zero in order to expand `log` around 0."""))
 
             order = a.order
@@ -63,7 +63,7 @@ for T in (:Taylor1, :TaylorN)
 
         function asin(a::$T)
             a0 = constant_term(a)
-            a0^2 == one(a0) && throw(DomainError(a, 
+            a0^2 == one(a0) && throw(DomainError(a,
                 """Series expansion of asin(x) diverges at x = ±1."""))
 
             order = a.order
@@ -79,7 +79,7 @@ for T in (:Taylor1, :TaylorN)
 
         function acos(a::$T)
             a0 = constant_term(a)
-            a0^2 == one(a0) && throw(DomainError(a, 
+            a0^2 == one(a0) && throw(DomainError(a,
             """Series expansion of asin(x) diverges at x = ±1."""))
 
             order = a.order
@@ -100,7 +100,7 @@ for T in (:Taylor1, :TaylorN)
             aa = one(aux) * a
             c = $T( aux, order)
             r = $T(one(aux) + a0^2, order)
-            iszero(constant_term(r)) && throw(DomainError(a, 
+            iszero(constant_term(r)) && throw(DomainError(a,
                 """Series expansion of atan(x) diverges at x = ±im."""))
 
             for k in eachindex(a)
@@ -576,7 +576,7 @@ function inverse(f::Taylor1{T}) where {T<:Number}
     z = Taylor1(T,f.order)
     zdivf = z/f
     zdivfpown = zdivf
-    S = eltype(zdivf)
+    S = TS.numtype(zdivf)
     coeffs = zeros(S,f.order+1)
 
     @inbounds for n in 1:f.order

src/intervals.jl

@@ -105,6 +105,7 @@ are mapped by an affine transformation to the intervals `-1..1`
 normalize_taylor(a::TaylorN, I::IntervalBox{N,T}, symI::Bool=true) where {T<:Real,N} =
     _normalize(a, I, Val(symI))
 
+export normalize_taylor
 
 #  I -> -1..1
 function _normalize(a::Taylor1, I::Interval{T}, ::Val{true}) where {T<:Real}

test/intervals.jl

@@ -15,8 +15,12 @@ eeuler = Base.MathConstants.e
     ti = Taylor1(Interval{Float64}, 10)
     x, y = set_variables(Interval{Float64}, "x y")
 
-    @test eltype(ti) == Interval{Float64}
-    @test eltype(x) == Interval{Float64}
+    # @test eltype(ti) == Interval{Float64}
+    # @test eltype(x) == Interval{Float64}
+    @test eltype(ti) == Taylor1{Interval{Float64}}
+    @test eltype(x) == TaylorN{Interval{Float64}}
+    @test TS.numtype(ti) == Interval{Float64}
+    @test TS.numtype(x) == Interval{Float64}
 
     @test p4(ti,-a) == (ti-a)^4
     @test p5(ti,-a) == (ti-a)^5
@@ -85,12 +89,12 @@ eeuler = Base.MathConstants.e
     @test diam(g1(-1..1)) < diam(gt(ii))
 
     # Test display for Taylor1{Complex{Interval{T}}}
-    vc = [complex(1.5 .. 2, 0 ), complex(-2  .. -1, -1 .. 1 ), 
+    vc = [complex(1.5 .. 2, 0 ), complex(-2  .. -1, -1 .. 1 ),
         complex( -1 .. 1.5, -1 .. 1.5), complex( 0..0, -1 .. 1.5)]
     displayBigO(false)
-    @test string(Taylor1(vc, 5)) == 
+    @test string(Taylor1(vc, 5)) ==
         " ( [1.5, 2] + [0, 0]im ) - ( [1, 2] + [-1, 1]im ) t + ( [-1, 1.5] + [-1, 1.5]im ) t² + ( [0, 0] + [-1, 1.5]im ) t³ "
     displayBigO(true)
-    @test string(Taylor1(vc, 5)) == 
+    @test string(Taylor1(vc, 5)) ==
         " ( [1.5, 2] + [0, 0]im ) - ( [1, 2] + [-1, 1]im ) t + ( [-1, 1.5] + [-1, 1.5]im ) t² + ( [0, 0] + [-1, 1.5]im ) t³ + 𝒪(t⁶)"
 end

test/manyvariables.jl

@@ -175,7 +175,9 @@ using LinearAlgebra
     @test !iszero(uT)
     @test iszero(zeroT)
 
-    @test eltype(xH) == Int
+    @test convert(eltype(xH), xH) === xH
+    @test eltype(xH) == HomogeneousPolynomial{Int}
+    @test TS.numtype(xH) == Int
     @test length(xH) == 2
     @test zero(xH) == 0*xH
     @test one(yH) == xH+yH
@@ -192,7 +194,9 @@ using LinearAlgebra
     @test TaylorN(uT) == convert(TaylorN{Complex},1)
     @test get_numvars() == 2
     @test length(uT) == get_order()+1
-    @test eltype(convert(TaylorN{Complex{Float64}},1)) == Complex{Float64}
+    @test convert(eltype(xT), xT) === xT
+    @test eltype(convert(TaylorN{Complex{Float64}},1)) == TaylorN{Complex{Float64}}
+    @test TS.numtype(convert(TaylorN{Complex{Float64}},1)) == Complex{Float64}
 
     @test 1+xT+yT == TaylorN(1,1) + TaylorN([xH,yH],1)
     @test xT-yT-1 == TaylorN([-1,xH-yH])
@@ -222,7 +226,7 @@ using LinearAlgebra
     @test (xH/complex(0,BigInt(3)))' ==
         im*HomogeneousPolynomial([BigInt(1),0])/3
     @test evaluate(xH) == zero(eltype(xH))
-    @test xH() == zero(eltype(xH))
+    @test xH() == zero(TS.numtype(xH))
     @test xH([1,1]) == evaluate(xH, [1,1])
     @test xH((1,1)) == xH(1, 1.0) == evaluate(xH, (1, 1.0)) == 1
     hp = -5.4xH+6.89yH
@@ -653,7 +657,8 @@ using LinearAlgebra
     @test taylor_expand(f11, 1.0,2.0) == taylor_expand(f22, [1,2.0])
     @test evaluate(taylor_expand(x->x[1] + x[2], [1,2])) == 3.0
     f33(x,y) = 3x+y
-    @test eltype(taylor_expand(f33,1,1)) == eltype(1)
+    @test eltype(taylor_expand(f33,1,1)) == TaylorN{eltype(1)}
+    @test TS.numtype(taylor_expand(f33,1,1)) == eltype(1)
     x,y = get_variables()
     xysq = x^2 + y^2
     update!(xysq,[1.0,-2.0])