Modify signature of evaluate for TaylorN's (#247)

* Add keyword parameter to some methods of evaluate

and add a keyword parameter to some methods of evaluate

* Add some tests

* Add documentation

* Remove broken tests

The removed  tests are  related to function-like evaluation, which
therefore implies using `sorting=true` in  `evaluate`

* Function-like evaluation  with/without sorting...

and related updates in the docs

* Minor fix related to display

* Remove a redundant method of `evaluate!`

docs/src/userguide.md

@@ -410,6 +410,35 @@ exy([.1,.02]) == exp(0.12)
 exy(:x, 0.0)
 ```
 
+Internally, `evaluate` for `TaylorN` considers separately
+the contributions of all `HomogeneousPolynomial`s by `order`,
+which are finally added up *after* sorting them in place (which is the default)
+in increasing order by `abs2`. This is done in order to
+use as many significant figures as possible of all terms
+in the final sum, which then should yield a more
+accurate result. This default can be changed to a non-sorting
+sum thought, which may be more performant or useful for
+certain subtypes of `Number` which, for instance, do not have `isless`
+defined. See
+[this issue](https://github.com/JuliaDiff/TaylorSeries.jl/issues/242)
+for a motivating example. This can be done using the keyword
+`sorting` in `evaluate`, which expects a `Bool`, or using a
+that boolean as the *first* argument in the function-like evaluation.
+
+```@repl userguide
+exy([.1,.02]) # default is `sorting=true`
+evaluate(exy, [.1,.02]; sorting=false)
+exy(false, [.1,.02])
+```
+
+In the examples shown above, the first entry corresponds to the
+default case (`sorting=true`), which yields the same result as
+`exp(0.12)`, and the remaining two illustrate
+turning off sorting the terms. Note that the results are not
+identical, since [floating point addition is not
+associative](https://en.wikipedia.org/wiki/Associative_property#Nonassociativity_of_floating_point_calculation),
+which may introduce rounding errors.
+
 The functions `taylor_expand` and `update!` work as well for `TaylorN`.
 
 ```@repl userguide

src/evaluate.jl

@@ -51,19 +51,9 @@ representing the Taylor expansion for the dependent variables
 of an ODE at *time* `δt`. It updates the vector `x0` with the
 computed values.
 """
-function evaluate!(x::AbstractArray{Taylor1{T}}, δt::T,
-        x0::AbstractArray{T}) where {T<:Number}
-
-    # @assert length(x) == length(x0)
-    @inbounds for i in eachindex(x, x0)
-        x0[i] = evaluate( x[i], δt )
-    end
-    nothing
-end
 function evaluate!(x::AbstractArray{Taylor1{T}}, δt::S,
         x0::AbstractArray{T}) where {T<:Number, S<:Number}
 
-    # @assert length(x) == length(x0)
     @inbounds for i in eachindex(x, x0)
         x0[i] = evaluate( x[i], δt )
     end
@@ -110,7 +100,6 @@ evaluate(p::Taylor1{T}, x::Array{S}) where {T<:Number, S<:Number} =
 
 #function-like behavior for Taylor1
 (p::Taylor1)(x) = evaluate(p, x)
-
 (p::Taylor1)() = evaluate(p)
 
 #function-like behavior for Vector{Taylor1}
@@ -146,11 +135,10 @@ function evaluate!(x::AbstractArray{TaylorN{T}}, δx::Array{Taylor1{T},1},
 end
 
 function evaluate!(x::AbstractArray{TaylorN{T}}, δx::Array{TaylorN{T},1},
-        x0::AbstractArray{TaylorN{T}}) where {T<:NumberNotSeriesN}
+        x0::AbstractArray{TaylorN{T}}; sorting::Bool=true) where {T<:NumberNotSeriesN}
 
-    # @assert length(x) == length(x0)
     @inbounds for i in eachindex(x, x0)
-        x0[i] = evaluate( x[i], δx )
+        x0[i] = _evaluate( x[i], δx, Val(sorting) )
     end
     nothing
 end
@@ -202,7 +190,7 @@ evaluate(a::HomogeneousPolynomial{T}, vals::Array{S,1} ) where
 
 evaluate(a::HomogeneousPolynomial, v, vals...) = evaluate(a, (v, vals...,))
 
-evaluate(a::HomogeneousPolynomial, v) = evaluate(a, v...)
+evaluate(a::HomogeneousPolynomial, v) = evaluate(a, [v...])
 
 function evaluate(a::HomogeneousPolynomial)
     a.order == 0 && return a[1]
@@ -217,31 +205,46 @@ end
 (p::HomogeneousPolynomial)() = evaluate(p)
 
 """
-    evaluate(a, [vals])
+    evaluate(a, [vals]; sorting::Bool=true)
 
 Evaluate the `TaylorN` polynomial `a` at `vals`.
-If `vals` is ommitted, it's evaluated at zero.
-Note that the syntax `a(vals)` is equivalent to `evaluate(a, vals)`; and `a()`
-is equivalent to `evaluate(a)`.
+If `vals` is ommitted, it's evaluated at zero. The
+keyword parameter `sorting` can be used to avoid
+sorting (in increasing order by `abs2`) the
+terms that are added.
+
+Note that the syntax `a(vals)` is equivalent to
+`evaluate(a, vals)`; and `a()` is equivalent to
+`evaluate(a)`. No extension exists that incorporates
+`sorting`.
 """
-function evaluate(a::TaylorN{T}, vals::NTuple{N,S}) where
-        {T<:Number,S<:NumberNotSeries, N}
+evaluate(a::TaylorN{T}, vals::NTuple; sorting::Bool=true) where {T<:Number} =
+    _evaluate(a, vals, Val(sorting))
 
-    @assert N == get_numvars()
+evaluate(a::TaylorN, vals; sorting::Bool=true) = _evaluate(a, (vals...,), Val(sorting))
 
-    R = promote_type(T,S)
+evaluate(a::TaylorN, v, vals...; sorting::Bool=true) =
+    _evaluate(a, (v, vals...,), Val(sorting))
+
+function _evaluate(a::TaylorN{T}, vals) where {T<:Number}
+    @assert get_numvars() == length(vals)
+    R = promote_type(T,typeof(vals[1]))
     a_length = length(a)
     suma = zeros(R, a_length)
     @inbounds for homPol in length(a):-1:1
         suma[homPol] = evaluate(a.coeffs[homPol], vals)
     end
+    return suma
+end
 
+function _evaluate(a::TaylorN{T}, vals::NTuple, ::Val{true}) where {T<:Number}
+    suma = _evaluate(a, vals)
     return sum( sort!(suma, by=abs2) )
 end
-
-evaluate(a::TaylorN, vals) = evaluate(a, (vals...,))
-
-evaluate(a::TaylorN, v, vals...) = evaluate(a, (v, vals...,))
+function _evaluate(a::TaylorN{T}, vals::NTuple, ::Val{false}) where {T<:Number}
+    suma = _evaluate(a, vals)
+    return sum( suma )
+end
 
 function evaluate(a::TaylorN{T}, vals::NTuple{N,Taylor1{S}}) where
         {T<:Number, S<:NumberNotSeries, N}
@@ -329,6 +332,8 @@ evaluate(A::AbstractArray{TaylorN{T}}) where {T<:Number} = evaluate.(A)
 (p::TaylorN)(s::Symbol, x) = evaluate(p, s, x)
 (p::TaylorN)(x::Pair) = evaluate(p, first(x), last(x))
 (p::TaylorN)(x, v...) = evaluate(p, (x, v...,))
+(p::TaylorN)(b::Bool, x) = evaluate(p, x, sorting=b)
+(p::TaylorN)(b::Bool, x, v...) = evaluate(p, (x, v...,), sorting=b)
 
 #function-like behavior for AbstractArray{TaylorN{T}}
 if VERSION > v"1.1"

src/printing.jl

@@ -181,7 +181,6 @@ function homogPol2str(a::HomogeneousPolynomial{Taylor1{T}}) where {T<:Number}
 end
 
 function numbr2str(zz, ifirst::Bool=false)
-    iszero(zz) && return string( zz )
     plusmin = ifelse( ifirst, string(""), string("+ ") )
     return string(plusmin, zz)
 end

test/manyvariables.jl

@@ -5,9 +5,10 @@ using TaylorSeries
 
 using Test
 using LinearAlgebra
-eeuler = Base.MathConstants.e
 
 @testset "Tests for HomogeneousPolynomial and TaylorN" begin
+    eeuler = Base.MathConstants.e
+
     @test HomogeneousPolynomial <: AbstractSeries
     @test HomogeneousPolynomial{Int} <: AbstractSeries{Int}
     @test TaylorN{Float64} <: AbstractSeries{Float64}
@@ -216,7 +217,7 @@ eeuler = Base.MathConstants.e
     @test evaluate(xH) == zero(eltype(xH))
     @test xH() == zero(eltype(xH))
     @test xH([1,1]) == evaluate(xH, [1,1])
-    @test xH((1,1)) == 1
+    @test xH((1,1)) == xH(1, 1.0) == evaluate(xH, (1, 1.0)) == 1
     hp = -5.4xH+6.89yH
     @test hp([1,1]) == evaluate(hp, [1,1])
     vr = rand(2)
@@ -376,6 +377,9 @@ eeuler = Base.MathConstants.e
     @test exy(0.1im, 0.01im) == exp(0.11im)
     @test evaluate(exy,(0.1im, 0.01im)) == exp(0.11im)
     @test exy((0.1im, 0.01im)) == exp(0.11im)
+    @test exy(true, (0.1im, 0.01im)) == exp(0.11im)
+    @test evaluate(exy, (0.1im, 0.01im), sorting=false) == exy(false, (0.1im, 0.01im))
+    @test evaluate(exy, (0.1im, 0.01im), sorting=false) == exy(false, 0.1im, 0.01im)
     @test evaluate(exy,[0.1im, 0.01im]) == exp(0.11im)
     @test exy([0.1im, 0.01im]) == exp(0.11im)
     @test isapprox(evaluate(exy, (1,1)), eeuler^2)
@@ -394,6 +398,8 @@ eeuler = Base.MathConstants.e
     v = zeros(Int, 2)
     @test evaluate!([xT, yT], ones(Int, 2), v) == nothing
     @test v == ones(2)
+    @test evaluate!([xT, yT][1:2], ones(Int, 2), v) == nothing
+    @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]
     @test evaluate(A_TN) == [0.0 0.0 0.0; 0.0 0.0 0.0]

test/onevariable.jl

@@ -5,7 +5,6 @@ using TaylorSeries
 
 using Test
 using LinearAlgebra, SparseArrays
-const eeuler = Base.MathConstants.e
 
 # This is used to check the fallack of pretty_print
 struct SymbNumber <: Number
@@ -14,6 +13,8 @@ end
 Base.iszero(::SymbNumber) = false
 
 @testset "Tests for Taylor1 expansions" begin
+    eeuler = Base.MathConstants.e
+
     ta(a) = Taylor1([a,one(a)],15)
     t = Taylor1(Int,15)
     tim = im*t