Skip to content

Commit

Permalink
Mathematically consistent big O notation (JuliaDiff#226)
Browse files Browse the repository at this point in the history 
* Change fixorder towards the smallest order (except when it is <= 0).

* Fixes related to TaylorN and mixtures

* Skip tests related to intervals, for Julia 1.0

IntervalArithmetic.jl requires Julia v1.1

* Minor fix in TaylorN constructor

* Use broadcasting in evaluate methods

* Promote method involving nested Taylor1s

* Add some type constraints and small fixes to convert methods

* More tests

* Add documentation of bigO convention used
  • Loading branch information
@lbenet
lbenet authored Oct 5, 2019
1 parent ca7fbc5 commit 53a8652
Show file tree
Hide file tree
Showing 9 changed files with 149 additions and 107 deletions.
45 changes: 24 additions & 21 deletions docs/src/userguide.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -47,7 +47,9 @@ t = shift_taylor(0.0) # Independent variable `t`
```

Note that the information about the maximum order considered is displayed
using a big-𝒪 notation. In some cases, it is desirable to not display
using a big-𝒪 notation. The convention followed when different orders are
combined is consistent with the mathematics and the big-𝒪 notation, that is,
to propagate the lowest order. In some cases, it is desirable to not display
the big-𝒪 notation. The function [`displayBigO`](@ref) allows to
control whether it is displayed or not.
```@repl userguide
Expand All @@ -69,9 +71,9 @@ t
```

The definition of `shift_taylor(a)` uses the method
[`Taylor1([::Type{Float64}], [order::Int=1])`](@ref), which is a
[`Taylor1([::Type{Float64}], order::Int)`](@ref), which is a
shortcut to define the independent variable of a Taylor expansion,
of given type and order (defaults are `Float64` and `order=1`).
of given type and order (the default is `Float64`).
This is one of the easiest ways to work with the package.

The usual arithmetic operators (`+`, `-`, `*`, `/`, `^`, `==`) have been
Expand All @@ -87,10 +89,11 @@ t*(t^2-4)/(t+2)
tI = im*t
(1-t)^3.2
(1+t)^t
Taylor1(3) + Taylor1(5) == 2Taylor1(3) # big-𝒪 convention applies
t^6 # t is of order 5
```

If no valid Taylor expansion can be computed, an error is thrown, for instance
If no valid Taylor expansion can be computed an error is thrown, for instance,
when a derivative is not defined (or simply diverges):

```@repl userguide
Expand All @@ -104,11 +107,11 @@ are computed recursively. At the moment of this writing, these functions
are `exp`, `log`, `sqrt`, the trigonometric functions
`sin`, `cos` and `tan`, their inverses, as well as the hyperbolic functions
`sinh`, `cosh` and `tanh` and their inverses;
more will be added in the future. Note that this way of obtaining the
Taylor coefficients is not the *laziest* way, in particular for many independent
more functions will be added in the future. Note that this way of obtaining the
Taylor coefficients is not a *lazy* way, in particular for many independent
variables. Yet, it is quite efficient, especially for the integration of
ordinary differential equations, which is among the applications we have in
mind (see also
mind (see
[TaylorIntegration.jl](https://github.com/PerezHz/TaylorIntegration.jl)).

```@repl userguide
Expand Down Expand Up @@ -144,7 +147,7 @@ functions return the corresponding [`Taylor1`](@ref) expansions.
The last coefficient of a derivative is set to zero to keep the
same order as the original polynomial; for the integral, an
integration constant may be set by the user (the default is zero). The
order of the resulting polynomial is not changed. The value of the
order of the resulting polynomial is kept unchanged. The value of the
``n``-th (``n \ge 0``)
derivative is obtained using `derivative(n,a)`, where `a` is a Taylor series.

Expand Down Expand Up @@ -173,7 +176,7 @@ eBig = evaluate( exp(tBig), one(BigFloat) )
ℯ - eBig
```

Another way to obtain the value of a `Taylor1` polynomial `p` at a given value `x`, is to call `p` as if it were a function, i.e., `p(x)`:
Another way to obtain the value of a `Taylor1` polynomial `p` at a given value `x`, is to call `p` as if it was a function, i.e., `p(x)`:

```@repl userguide
t = Taylor1(15)
Expand All @@ -186,13 +189,12 @@ p() == p(0.0) # p() is a shortcut to obtain the 0-th order coefficient of `p`
```

Note that the syntax `p(x)` is equivalent to `evaluate(p, x)`, whereas `p()` is
equivalent to `evaluate(p)`. For more details about function-like behavior for a
given type in Julia, see the [Function-like objects](https://docs.julialang.org/en/stable/manual/methods/#Function-like-objects-1)
section of the Julia manual.
equivalent to `evaluate(p)`.

Useful shortcuts are [`taylor_expand`](@ref) are [`update!`](@ref).
Useful shortcuts are [`taylor_expand`](@ref) and [`update!`](@ref).
The former returns
the expansion of a function around a given value `t0`. In turn, `update!`
the expansion of a function around a given value `t0`, mimicking the use
of `shift_taylor` above. In turn, `update!`
provides an in-place update of a given Taylor polynomial, that is, it shifts
it further by the provided amount.

Expand Down Expand Up @@ -237,15 +239,15 @@ one can specify the variables using a vector of symbols.
set_variables([:x, :y], order=10)
```

Similarly, numbered variables are also available by specifying a single
Similarly, subindexed variables are also available by specifying a single
variable name and the optional keyword argument `numvars`:

```@repl userguide
set_variables("α", numvars=3)
```

Alternatively to `set_variables`, [`get_variables`](@ref) can be used if one
doesn't want to change internal dictionaries. `get_variables` returns a vector
does not want to change internal dictionaries. `get_variables` returns a vector
of `TaylorN` independent variables of a desired `order`
(lesser than `get_order` so the
internals doesn't have to change) with the length and variable names defined
Expand Down Expand Up @@ -300,13 +302,14 @@ the corresponding independent variable, i.e. ``x \to x+a``.

As before, the usual arithmetic operators (`+`, `-`, `*`, `/`, `^`, `==`)
have been extended to work with [`TaylorN`](@ref) objects, including the
appropriate
promotions to deal with numbers. (Some of the arithmetic operations have
been extended for
appropriate promotions to deal with numbers.
Note that some of the arithmetic operations have been extended for
[`HomogeneousPolynomial`](@ref), whenever the result is a
[`HomogeneousPolynomial`](@ref); division, for instance, is not extended.)
[`HomogeneousPolynomial`](@ref); division, for instance, is not extended.
The same convention used for `Taylor1` objects is used when combining
`TaylorN` polynomials of different order.

Also, the elementary functions have been
The elementary functions have also been
implemented, again by computing their coefficients recursively:

```@repl userguide
Expand Down
8 changes: 5 additions & 3 deletions src/auxiliary.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -243,9 +243,11 @@ for T in (:Taylor1, :TaylorN)
@eval begin
@inline function fixorder(a::$T, b::$T)
a.order == b.order && return a, b
a.order < b.order &&
return $T(copy(a.coeffs), b.order), $T(copy(b.coeffs), b.order)
return $T(copy(a.coeffs), a.order), $T(copy(b.coeffs), a.order)
minorder, maxorder = minmax(a.order, b.order)
if minorder ≤ 0
minorder = maxorder
end
return $T(copy(a.coeffs), minorder), $T(copy(b.coeffs), minorder)
end
end
end
Expand Down
38 changes: 22 additions & 16 deletions src/constructors.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -48,7 +48,6 @@ Taylor1(x::Taylor1{T}) where {T<:Number} = x
Taylor1(coeffs::Array{T,1}, order::Int) where {T<:Number} = Taylor1{T}(coeffs, order)
Taylor1(coeffs::Array{T,1}) where {T<:Number} = Taylor1(coeffs, length(coeffs)-1)
function Taylor1(x::T, order::Int) where {T<:Number}
# v = zeros(T, order+1)
v = fill(zero(x), order+1)
v[1] = x
return Taylor1(v, order)
Expand All @@ -69,8 +68,8 @@ julia> Taylor1(Rational{Int}, 4)
1//1 t + 𝒪(t⁵)
```
"""
Taylor1(::Type{T}, order::Int=1) where {T<:Number} = Taylor1( [zero(T), one(T)], order)
Taylor1(order::Int=1) = Taylor1(Float64, order)
Taylor1(::Type{T}, order::Int) where {T<:Number} = Taylor1( [zero(T), one(T)], order)
Taylor1(order::Int) = Taylor1(Float64, order)


######################### HomogeneousPolynomial
Expand Down Expand Up @@ -153,26 +152,31 @@ struct TaylorN{T<:Number} <: AbstractSeries{T}
order :: Int

function TaylorN{T}(v::Array{HomogeneousPolynomial{T},1}, order::Int) where T<:Number
m = maxorderH(v)
order = max( m, order )
coeffs = zeros(HomogeneousPolynomial{T}, order)
@simd for i in eachindex(v)
@inbounds ord = v[i].order
@inbounds coeffs[ord+1] += v[i]
@inbounds for i in eachindex(v)
ord = v[i].order
if ord ≤ order
coeffs[ord+1] += v[i]
end
end
new{T}(coeffs, order)
end
end

TaylorN(x::TaylorN{T}) where T<:Number = x
TaylorN(x::Array{HomogeneousPolynomial{T},1}, order::Int) where {T<:Number} =
TaylorN{T}(x, order)
TaylorN(x::Array{HomogeneousPolynomial{T},1}) where {T<:Number} = TaylorN(x,0)
TaylorN(x::HomogeneousPolynomial{T},order::Int) where {T<:Number} =
TaylorN([x], order)
TaylorN(x::HomogeneousPolynomial{T}) where {T<:Number} = TaylorN([x], x.order)
function TaylorN(x::Array{HomogeneousPolynomial{T},1}, order::Int) where {T<:Number}
if order == 0
order = maxorderH(x)
end
return TaylorN{T}(x, order)
end
TaylorN(x::Array{HomogeneousPolynomial{T},1}) where {T<:Number} =
TaylorN(x, maxorderH(x))
TaylorN(x::HomogeneousPolynomial{T}, order::Int) where {T<:Number} =
TaylorN( [x], order )
TaylorN(x::HomogeneousPolynomial{T}) where {T<:Number} = TaylorN(x, x.order)
TaylorN(x::T, order::Int) where {T<:Number} =
TaylorN([HomogeneousPolynomial([x], 0)], order)
TaylorN(HomogeneousPolynomial([x], 0), order)

# Shortcut to define TaylorN independent variables
"""
Expand All @@ -192,9 +196,11 @@ julia> TaylorN(Rational{Int},2)
```
"""
TaylorN(::Type{T}, nv::Int; order::Int=get_order()) where {T<:Number} =
return TaylorN( [HomogeneousPolynomial(T, nv)], order )
TaylorN( HomogeneousPolynomial(T, nv), order )
TaylorN(nv::Int; order::Int=get_order()) = TaylorN(Float64, nv, order=order)



# A `Number` which is not an `AbstractSeries`
const NumberNotSeries = Union{setdiff(subtypes(Number), [AbstractSeries])...}

Expand Down
63 changes: 43 additions & 20 deletions src/conversion.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -10,20 +10,16 @@
convert(::Type{Taylor1{T}}, a::Taylor1) where {T<:Number} =
Taylor1(convert(Array{T,1}, a.coeffs), a.order)

function convert(::Type{Taylor1{Rational{T}}}, a::Taylor1{S}) where
{T<:Integer, S<:AbstractFloat}
convert(::Type{Taylor1{T}}, a::Taylor1{T}) where {T<:Number} = a

la = length(a.coeffs)
v = Array{Rational{T}}(undef, la)
v .= rationalize.(a[0:la-1])
return Taylor1(v)
end
convert(::Type{Taylor1{Rational{T}}}, a::Taylor1{S}) where
{T<:Integer, S<:AbstractFloat} = Taylor1(rationalize.(a[:]), a.order)

convert(::Type{Taylor1{T}}, b::Array{T,1}) where {T<:Number} =
Taylor1(b, length(b)-1)

convert(::Type{Taylor1{T}}, b::Array{S,1}) where {T<:Number, S<:Number} =
Taylor1(convert(Array{T,1},b))
Taylor1(convert(Array{T,1},b), length(b)-1)

convert(::Type{Taylor1{T}}, b::S) where {T<:Number, S<:Number} =
Taylor1([convert(T,b)], 0)
Expand All @@ -36,6 +32,9 @@ convert(::Type{Taylor1}, a::T) where {T<:Number} = Taylor1(a, 0)
convert(::Type{HomogeneousPolynomial{T}}, a::HomogeneousPolynomial) where {T<:Number} =
HomogeneousPolynomial(convert(Array{T,1}, a.coeffs), a.order)

convert(::Type{HomogeneousPolynomial{T}}, a::HomogeneousPolynomial{T}) where {T<:Number} =
a

function convert(::Type{HomogeneousPolynomial{Rational{T}}},
a::HomogeneousPolynomial{S}) where {T<:Integer, S<:AbstractFloat}

Expand Down Expand Up @@ -63,11 +62,13 @@ convert(::Type{HomogeneousPolynomial}, a::Number) = HomogeneousPolynomial([a],0)
convert(::Type{TaylorN{T}}, a::TaylorN) where {T<:Number} =
TaylorN( convert(Array{HomogeneousPolynomial{T},1}, a.coeffs), a.order)

convert(::Type{TaylorN{T}}, a::TaylorN{T}) where {T<:Number} = a

convert(::Type{TaylorN{T}}, b::HomogeneousPolynomial{S}) where {T<:Number, S<:Number} =
TaylorN( [convert(HomogeneousPolynomial{T}, b)], b.order)

convert(::Type{TaylorN{T}}, b::Array{HomogeneousPolynomial{S},1}) where {T<:Number, S<:Number} =
TaylorN( convert(Array{HomogeneousPolynomial{T},1}, b), length(b)-1)
TaylorN( convert(Array{HomogeneousPolynomial{T},1}, b), maxorderH(b))

convert(::Type{TaylorN{T}}, b::S) where {T<:Number, S<:Number} =
TaylorN( [HomogeneousPolynomial([convert(T, b)], 0)], 0)
Expand All @@ -76,15 +77,15 @@ convert(::Type{TaylorN{T}}, b::HomogeneousPolynomial{T}) where {T<:Number} =
TaylorN( [b], b.order)

convert(::Type{TaylorN{T}}, b::Array{HomogeneousPolynomial{T},1}) where {T<:Number} =
TaylorN( b, length(b)-1)
TaylorN( b, maxorderH(b))

convert(::Type{TaylorN{T}}, b::T) where {T<:Number} =
TaylorN( [HomogeneousPolynomial([b], 0)], 0)

convert(::Type{TaylorN}, b::Number) = TaylorN( [HomogeneousPolynomial([b], 0)], 0)


function convert(::Type{TaylorN{Taylor1{T}}}, s::Taylor1{TaylorN{T}}) where {T<:Number}
function convert(::Type{TaylorN{Taylor1{T}}}, s::Taylor1{TaylorN{T}}) where {T<:NumberNotSeries}

orderN = get_order()
r = zeros(HomogeneousPolynomial{Taylor1{T}}, orderN)
Expand All @@ -103,7 +104,7 @@ function convert(::Type{TaylorN{Taylor1{T}}}, s::Taylor1{TaylorN{T}}) where {T<:
return TaylorN(r)
end

function convert(::Type{Taylor1{TaylorN{T}}}, s::TaylorN{Taylor1{T}}) where {T<:Number}
function convert(::Type{Taylor1{TaylorN{T}}}, s::TaylorN{Taylor1{T}}) where {T<:NumberNotSeries}

ordert = 0
for ordHP in eachindex(s)
Expand All @@ -115,12 +116,12 @@ function convert(::Type{Taylor1{TaylorN{T}}}, s::TaylorN{Taylor1{T}}) where {T<:
end

@inbounds for ordN in eachindex(s)
vHP = HomogeneousPolynomial(zeros(T,length(s[ordN])))
vHP = HomogeneousPolynomial(zeros(T, length(s[ordN])))
@inbounds for ihp in eachindex(s[ordN].coeffs)
@inbounds for ind in eachindex(s[ordN][ihp].coeffs)
c = s[ordN][ihp][ind-1]
vHP[ihp] = c
vT[ind] += TaylorN(vHP)
vT[ind] += TaylorN(vHP, s.order)
vHP[ihp] = zero(T)
end
end
Expand All @@ -129,21 +130,21 @@ function convert(::Type{Taylor1{TaylorN{T}}}, s::TaylorN{Taylor1{T}}) where {T<:
end

function convert(::Type{Array{TaylorN{Taylor1{T}},N}},
s::Array{Taylor1{TaylorN{T}},N}) where {T<:Number, N}
s::Array{Taylor1{TaylorN{T}},N}) where {T<:NumberNotSeries, N}

v = Array{TaylorN{Taylor1{T}}}(undef, size(s))
for ind in eachindex(s)
v[ind] = convert(TaylorN{Taylor1{T}}, s[ind])
@simd for ind in eachindex(s)
@inbounds v[ind] = convert(TaylorN{Taylor1{T}}, s[ind])
end
return v
end

function convert(::Type{Array{Taylor1{TaylorN{T}},N}},
s::Array{TaylorN{Taylor1{T}},N}) where {T<:Number, N}
s::Array{TaylorN{Taylor1{T}},N}) where {T<:NumberNotSeries, N}

v = Array{Taylor1{TaylorN{T}}}(undef, size(s))
for ind in eachindex(s)
v[ind] = convert(Taylor1{TaylorN{T}}, s[ind])
@simd for ind in eachindex(s)
@inbounds v[ind] = convert(Taylor1{TaylorN{T}}, s[ind])
end
return v
end
Expand All @@ -170,10 +171,14 @@ promote_rule(::Type{Taylor1{T}}, ::Type{S}) where {T<:Number, S<:Number} =
promote_rule(::Type{Taylor1{Taylor1{T}}}, ::Type{Taylor1{T}}) where {T<:Number} =
Taylor1{Taylor1{T}}


promote_rule(::Type{HomogeneousPolynomial{T}},
::Type{HomogeneousPolynomial{S}}) where {T<:Number, S<:Number} =
HomogeneousPolynomial{promote_type(T,S)}

promote_rule(::Type{HomogeneousPolynomial{T}},
::Type{HomogeneousPolynomial{T}}) where {T<:Number} = HomogeneousPolynomial{T}

promote_rule(::Type{HomogeneousPolynomial{T}},
::Type{Array{S,1}}) where {T<:Number, S<:Number} = HomogeneousPolynomial{promote_type(T,S)}

Expand All @@ -184,6 +189,8 @@ promote_rule(::Type{HomogeneousPolynomial{T}}, ::Type{S}) where
promote_rule(::Type{TaylorN{T}}, ::Type{TaylorN{S}}) where {T<:Number, S<:Number}=
TaylorN{promote_type(T,S)}

promote_rule(::Type{TaylorN{T}}, ::Type{TaylorN{T}}) where {T<:Number} = TaylorN{T}

promote_rule(::Type{TaylorN{T}}, ::Type{HomogeneousPolynomial{S}}) where
{T<:Number, S<:Number} = TaylorN{promote_type(T,S)}

Expand All @@ -203,3 +210,19 @@ promote_rule(::Type{S}, ::Type{T}) where

promote_rule(::Type{Taylor1{T}}, ::Type{TaylorN{S}}) where {T<:NumberNotSeries, S<:NumberNotSeries} =
throw(ArgumentError("There is no reasonable promotion among `Taylor1{$T}` and `TaylorN{$S}` types"))


# Nested Taylor1's
function promote(a::Taylor1{Taylor1{T}}, b::Taylor1{T}) where {T<:NumberNotSeriesN}
order_a = get_order(a)
order_b = get_order(b)
zb = zero(b)
new_bcoeffs = similar(a.coeffs)
new_bcoeffs[1] = b
@inbounds for ind in 2:order_a+1
new_bcoeffs[ind] = zb
end
return a, Taylor1(b, order_a)
end
promote(b::Taylor1{T}, a::Taylor1{Taylor1{T}}) where {T<:NumberNotSeriesN} =
reverse(promote(a, b))
Loading

0 comments on commit 53a8652

Please sign in to comment.