Support < and > through isless (#323)

* Add findfirst to TaylorN and HomogeneousPolynomial, and tests

* Add < and > for Taylor1, TaylorN and HomPolyn and `Real`s

* Add tests for < and >

* Simplify code and  improvements

* Overload isless (instead of <, >) and increase test coverage

* Add comments to documentation

* Support isless for mixtures, but not nested Taylor1s

* Small changes in tests (considers promotions)

* Small changes in tests

* More info on the lexicographic order (docs)

* Include in comments iss #326

* Remove repeated Taylor1 tests

* Small fix in docs

* Include ordering tests in a `@testset` sub-block

* Bump minor version

Project.toml

@@ -1,7 +1,7 @@
 name = "TaylorSeries"
 uuid = "6aa5eb33-94cf-58f4-a9d0-e4b2c4fc25ea"
 repo = "https://github.com/JuliaDiff/TaylorSeries.jl.git"
-version = "0.14.0"
+version = "0.15.0"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

docs/src/api.md

@@ -50,6 +50,7 @@ inverse
 abs
 norm
 isapprox
+isless
 isfinite
 displayBigO
 use_show_default

docs/src/userguide.md

@@ -105,6 +105,16 @@ of order 5; this is be consistent with the number of known coefficients of the
 returned series, since the result corresponds to factorize `t` in the numerator
 to cancel the same factor in the denominator.
 
+`Taylor1` is also equiped with a total order, provided by overloading [`isless`](@ref).
+The ordering is consistent with the interpretation that there are infinitessimal
+elements in the algebra; for details see M. Berz, "Automatic Differentiation as
+Nonarchimedean Analysis", Computer Arithmetic and Enclosure Methods, (1992), Elsevier,
+439-450. This is illustrated by:
+
+```@repl userguide
+1 > t > 2*t^2 > 100*t^3 >= 0
+```
+
 If no valid Taylor expansion can be computed an error is thrown, for instance,
 when a derivative is not defined at the expansion point, or it simply diverges.
 
@@ -336,12 +346,29 @@ Note that some of the arithmetic operations have been extended for
 [`HomogeneousPolynomial`](@ref); division, for instance, is not extended.
 The same convention used for `Taylor1` objects is used when combining
 `TaylorN` polynomials of different order.
+Both `HomogeneousPolynomial` and `TaylorN` are equiped with a total *lexicographical*
+order, provided by overloading [`isless`](@ref).
+The ordering is consistent with the interpretation that there are infinitessimal
+elements in the algebra; for details see M. Berz, "Automatic Differentiation as
+Nonarchimedean Analysis", Computer Arithmetic and Enclosure Methods, (1992), Elsevier,
+439-450.
+Essentially, the lexicographic order works as follows: smaller order monomials
+are *larger* than higher order monomials; when the order is the same, *larger* monomials
+appear before in the hash-tables; the function [`show_monomials`](@ref).
+```@repl userguide
+x, y = set_variables("x y", order=10);
+
+show_monomials(2)
+```
+Then, the following then holds:
+```@repl userguide
+0 < 1e8 * y^2 < x*y < x^2  < y < x/1e8 < 1.0
+```
 
 The elementary functions have also been
 implemented, again by computing their coefficients recursively:
 
 ```@repl userguide
-x, y = set_variables("x y", order=10);
 exy = exp(x+y)
 ```
 

src/TaylorSeries.jl

@@ -39,7 +39,7 @@ import Base: ==, +, -, *, /, ^
 import Base: iterate, size, eachindex, firstindex, lastindex,
     eltype, length, getindex, setindex!, axes, copyto!
 
-import Base: zero, one, zeros, ones, isinf, isnan, iszero,
+import Base: zero, one, zeros, ones, isinf, isnan, iszero, isless,
     convert, promote_rule, promote, show,
     real, imag, conj, adjoint,
     rem, mod, mod2pi, abs, abs2,

src/arithmetic.jl

@@ -10,7 +10,6 @@
 
 ## Equality ##
 for T in (:Taylor1, :TaylorN)
-
     @eval begin
         ==(a::$T{T}, b::$T{S}) where {T<:Number,S<:Number} = ==(promote(a,b)...)
 
@@ -34,11 +33,127 @@ function ==(a::HomogeneousPolynomial, b::HomogeneousPolynomial)
     return iszero(a.coeffs) && iszero(b.coeffs)
 end
 
+## Total ordering ##
+for T in (:Taylor1, :TaylorN)
+    @eval begin
+        @inline function isless(a::$T{<:Number}, b::Real)
+            a0 = constant_term(a)
+            a0 != b && return isless(a0, b)
+            nz = findfirst(a-b)
+            if nz == -1
+                return isless(zero(a0), zero(b))
+            else
+                return isless(a[nz], zero(b))
+            end
+        end
+        @inline function isless(b::Real, a::$T{<:Number})
+            a0 = constant_term(a)
+            a0 != b && return isless(b, a0)
+            nz = findfirst(b-a)
+            if nz == -1
+                return isless(zero(b), zero(a0))
+            else
+                return isless(zero(b), a[nz])
+            end
+        end
+        #
+        @inline isless(a::$T{T}, b::$T{S}) where {T<:Number, S<:Number} = isless(promote(a,b)...)
+        @inline isless(a::$T{T}, b::$T{T}) where {T<:Number} = isless(a - b, zero(constant_term(a)))
+    end
+end
+
+@inline function isless(a::HomogeneousPolynomial{<:Number}, b::Real)
+    orda = get_order(a)
+    if orda == 0
+        return isless(a[1], b)
+    else
+        !iszero(b) && return isless(zero(a[1]), b)
+        nz = max(findfirst(a), 1)
+        return isless(a[nz], b)
+    end
+end
+@inline function isless(b::Real, a::HomogeneousPolynomial{<:Number})
+    orda = get_order(a)
+    if orda == 0
+        return isless(b, a[1])
+    else
+        !iszero(b) && return isless(b, zero(a[1]))
+        nz = max(findfirst(a),1)
+        return isless(b, a[nz])
+    end
+end
+#
+@inline isless(a::HomogeneousPolynomial{T}, b::HomogeneousPolynomial{S}) where {T<:Number, S<:Number} = isless(promote(a,b)...)
+@inline function isless(a::HomogeneousPolynomial{T}, b::HomogeneousPolynomial{T}) where {T<:Number}
+    orda = get_order(a)
+    ordb = get_order(b)
+    if orda == ordb
+        return isless(a-b, zero(a[1]))
+    elseif orda < ordb
+        return isless(a, zero(a[1]))
+    else
+        return isless(-b, zero(a[1]))
+    end
+end
+
+# Mixtures
+@inline isless(a::Taylor1{TaylorN{T}}, b::Taylor1{TaylorN{T}}) where {T<:NumberNotSeries} =
+    isless(a - b, zero(T))
+@inline function isless(a::HomogeneousPolynomial{Taylor1{T}}, b::HomogeneousPolynomial{Taylor1{T}}) where {T<:NumberNotSeries}
+    orda = get_order(a)
+    ordb = get_order(b)
+    if orda == ordb
+        return isless(a-b, zero(T))
+    elseif orda < ordb
+        return isless(a, zero(T))
+    else
+        return isless(-b, zero(T))
+    end
+end
+@inline isless(a::TaylorN{Taylor1{T}}, b::TaylorN{Taylor1{T}}) where {T<:NumberNotSeries} = isless(a - b, zero(T))
+
+#= TODO: Nested Taylor1s; needs careful thinking; iss #326. The following works:
+@inline isless(a::Taylor1{Taylor1{T}}, b::Taylor1{Taylor1{T}}) where {T<:Number} = isless(a - b, zero(T))
+# Is the following correct?
+# ti = Taylor1(3)
+# to = Taylor1([zero(ti), one(ti)], 9)
+# tito = ti * to
+# ti > to > 0 # ok
+# to^2 < toti < ti^2 # ok
+# ti > ti^2 > to # is this ok?
+=#
+
+@doc doc"""
+    isless(a::Taylor1{<:Real}, b::Real)
+    isless(a::TaylorN{<:Real}, b::Real)
+
+Compute `isless` by comparing the `constant_term(a)` and `b`. If they are equal,
+returns `a[nz] < 0`, with `nz` the first
+non-zero coefficient after the constant term. This defines a total order.
+
+For many variables, the ordering includes a lexicographical convention in order to be
+total. We have opted for the simplest one, where the *larger* variable appears *before*
+when the `TaylorN` variables are defined (e.g., through [`set_variables`](@ref)).
+
+Refs:
+- M. Berz, AIP Conference Proceedings 177, 275 (1988); https://doi.org/10.1063/1.37800
+- M. Berz, "Automatic Differentiation as Nonarchimedean Analysis", Computer Arithmetic and
+    Enclosure Methods, (1992), Elsevier, 439-450.
+
+---
+
+    isless(a::Taylor1{<:Real}, b::Taylor1{<:Real})
+    isless(a::TaylorN{<:Real}, b::Taylor1{<:Real})
+
+Return `isless(a - b, zero(b))`.
+""" isless
+
+
+## zero and one ##
 for T in (:Taylor1, :HomogeneousPolynomial, :TaylorN)
     @eval iszero(a::$T) = iszero(a.coeffs)
 end
 
-## zero and one ##
 for T in (:Taylor1, :TaylorN)
     @eval zero(a::$T) = $T(zero.(a.coeffs))
     @eval function one(a::$T)

src/auxiliary.jl

@@ -268,11 +268,35 @@ end
 # Finds the first non zero entry; extended to Taylor1
 function Base.findfirst(a::Taylor1{T}) where {T<:Number}
     first = findfirst(x->!iszero(x), a.coeffs)
-    isa(first, Nothing) && return -1
+    isnothing(first) && return -1
     return first-1
 end
 # Finds the last non-zero entry; extended to Taylor1
 function Base.findlast(a::Taylor1{T}) where {T<:Number}
+    last = findlast(x->!iszero(x), a.coeffs)
+    isnothing(last) && return -1
+    return last-1
+end
+
+# Finds the first non zero entry; extended to HomogeneousPolynomial
+function Base.findfirst(a::HomogeneousPolynomial{T}) where {T<:Number}
+    first = findfirst(x->!iszero(x), a.coeffs)
+    isa(first, Nothing) && return -1
+    return first
+end
+function Base.findfirst(a::TaylorN{T}) where {T<:Number}
+    first = findfirst(x->!iszero(x), a.coeffs)
+    isa(first, Nothing) && return -1
+    return first-1
+end
+# Finds the last non-zero entry; extended to HomogeneousPolynomial
+function Base.findlast(a::HomogeneousPolynomial{T}) where {T<:Number}
+    last = findlast(x->!iszero(x), a.coeffs)
+    isa(last, Nothing) && return -1
+    return last
+end
+# Finds the last non-zero entry; extended to TaylorN
+function Base.findlast(a::TaylorN{T}) where {T<:Number}
     last = findlast(x->!iszero(x), a.coeffs)
     isa(last, Nothing) && return -1
     return last-1

test/manyvariables.jl

@@ -93,17 +93,23 @@ end
     @test v == [1,2,0]
     @test_throws AssertionError TS.resize_coeffsHP!(v,1)
     hpol_v = HomogeneousPolynomial(v)
+    @test findfirst(hpol_v) == 1
+    @test findlast(hpol_v) == 2
     hpol_v[3] = 3
     @test v == [1,2,3]
     hpol_v[1:3] = 3
     @test v == [3,3,3]
     hpol_v[1:2:2] = 0
     @test v == [0,3,3]
+    @test findfirst(hpol_v) == 2
+    @test findlast(hpol_v) == 3
     hpol_v[1:1:2] = [1,2]
     @test all(hpol_v[1:1:2] .== [1,2])
     @test v == [1,2,3]
     hpol_v[:] = zeros(Int, 3)
     @test hpol_v == 0
+    @test findfirst(hpol_v) == -1
+    @test findlast(hpol_v) == -1
 
     tn_v = TaylorN(HomogeneousPolynomial(zeros(Int, 3)))
     tn_v[0] = 1
@@ -132,9 +138,12 @@ end
     @test (@inferred real(xH)) == xH
     xT = TaylorN(xH, 17)
     yT = TaylorN(Int, 2, order=17)
+    @test findfirst(xT) == 1
+    @test findlast(yT) == 1
     @test (@inferred conj(xT)) == (@inferred adjoint(xT))
     @test (@inferred real(xT)) == (xT)
     zeroT = zero( TaylorN([xH],1) )
+    @test findfirst(zeroT) == findlast(zeroT) == -1
     @test (@inferred imag(xT)) == (zeroT)
     @test zeroT.coeffs == zeros(HomogeneousPolynomial{Int}, 1)
     @test size(xH) == (2,)
@@ -156,6 +165,16 @@ end
     @test TS.fixorder(xH,yH) == (xH,yH)
     @test_throws AssertionError TS.fixorder(zeros(xH,0)[1],yH)
 
+    @testset "Lexicographic order" begin
+        @test HomogeneousPolynomial([2.1]) > 0.5 * xH > yH > xH^2
+        @test -1.0*xH < -yH^2 < 0 < HomogeneousPolynomial([1.0])
+        @test -3 < HomogeneousPolynomial([-2]) < 0.0
+        @test !(zero(yH^2) > 0)
+        @test 1 > 0.5 * xT > yT > xT^2 > 0 > TaylorN(-1.0, 3)
+        @test -1 < -0.25 * yT < -xT^2 < 0.0 < TaylorN(1, 3)
+        @test !(zero(xT) > 0)
+        @test !(zero(yT^2) < 0)
+    end
 
     @test constant_term(xT) == 0
     @test constant_term(uT) == 1.0
@@ -761,6 +780,9 @@ end
     @test yT[0] == xT[0]*(pi/180)
     TS.rad2deg!(yT, xT, 0)
     @test yT[0] == xT[0]*(180/pi)
+    # 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
 end
 
 @testset "Integrate for several variables" begin