Skip to content

Commit

Permalink
Fixes in evaluate and pow (#369)
Browse files Browse the repository at this point in the history 
* Fixes in evaluate, and rename one method evaluate! -> _evaluate!

* Add function barrier _pow

* Fixes in pow!

* Fixes in _pow for Interval coeffs, and tests

* More fixes

* Metaprogramming and fixes in
pow!

for integer power, pow! uses power_by_squaring, is defined for TaylorN, or when coeffs are intervals

* Further minor fixes

* Another minor fix

* Bump patch version
  • Loading branch information
@lbenet
lbenet authored Jan 23, 2025
1 parent aa788e2 commit 84a007f
Show file tree
Hide file tree
Showing 6 changed files with 254 additions and 310 deletions.
2 changes: 1 addition & 1 deletion Project.toml
View file
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
name = "TaylorSeries"
uuid = "6aa5eb33-94cf-58f4-a9d0-e4b2c4fc25ea"
repo = "https://github.com/JuliaDiff/TaylorSeries.jl.git"
version = "0.18.2"
version = "0.18.3"

[deps]
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
Expand Down
123 changes: 66 additions & 57 deletions ext/TaylorSeriesIAExt.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,75 +8,84 @@ import TaylorSeries: evaluate, _evaluate, normalize_taylor, square

isdefined(Base, :get_extension) ? (using IntervalArithmetic) : (using ..IntervalArithmetic)

# Method used for Taylor1{Interval{T}}^n
for T in (:Taylor1, :TaylorN)
@eval begin
function ^(a::$T{Interval{S}}, n::Integer) where {S<:Real}
n == 0 && return one(a)
n == 1 && return copy(a)
n == 2 && return square(a)
n < 0 && return a^float(n)
return power_by_squaring(a, n)
@eval ^(a::$T{Interval{T}}, n::Integer) where {T<:Real} = TS.power_by_squaring(a, n)

@eval ^(a::$T{Interval{T}}, r::Rational) where {T<:Real} = a^float(r)

@eval function ^(a::$T{Interval{T}}, r::S) where {T<:Real, S<:Real}
isinteger(r) && r ≥ 0 && return TS.power_by_squaring(a, Integer(r))
a0 = constant_term(a) ∩ Interval(zero(T), T(Inf))
@assert !isempty(a0)
aux = one(a0^r)
a[0] = aux * a0
r == 0.5 && return sqrt(a)
if $T == TaylorN
if iszero(a0)
throw(DomainError(a,
"""The 0-th order TaylorN coefficient must be non-zero
in order to expand `^` around 0."""))
end
end
return TS._pow(a, r)
end
^(a::$T{Interval{S}}, r::Rational) where {S<:Real} = a^float(r)
end
end

function ^(a::Taylor1{Interval{T}}, r::S) where {T<:Real, S<:Real}
a0 = constant_term(a) ∩ Interval(zero(T), T(Inf))
aux = one(a0)^r

iszero(r) && return Taylor1(aux, a.order)
aa = one(aux) * a
aa[0] = one(aux) * a0
r == 1 && return aa
r == 2 && return square(aa)
r == 1/2 && return sqrt(aa)

l0 = findfirst(a)
lnull = trunc(Int, r*l0 )
if (a.order-lnull < 0) || (lnull > a.order)
return Taylor1( zero(aux), a.order)
end
c_order = l0 == 0 ? a.order : min(a.order, trunc(Int,r*a.order))
c = Taylor1(zero(aux), c_order)
for k = 0:c_order
TS.pow!(c, aa, c, r, k)
end
# _pow
function TS._pow(a::$T{Interval{S}}, n::Integer) where {S<:Real}
n < 0 && return TS._pow(a, float(n))
return TS.power_by_squaring(a, n)
end

return c
end
function ^(a::TaylorN{Interval{T}}, r::S) where {T<:Real, S<:Real}
a0 = constant_term(a) ∩ Interval(zero(T), T(Inf))
a0r = a0^r
aux = one(a0r)

iszero(r) && return TaylorN(aux, a.order)
aa = aux * a
aa[0] = aux * a0
r == 1 && return aa
r == 2 && return square(aa)
r == 1/2 && return sqrt(aa)
isinteger(r) && return aa^round(Int,r)

# @assert !iszero(a0)
iszero(a0) && throw(DomainError(a,
"""The 0-th order TaylorN coefficient must be non-zero
in order to expand `^` around 0."""))

c = TaylorN( a0r, a.order)
for ord in 1:a.order
TS.pow!(c, aa, c, r, ord)
end
function TS._pow(a::$T{Interval{T}}, r::S) where {T<:Real, S<:Real}
isinteger(r) && r ≥ 0 && return TS.power_by_squaring(a, Integer(r))
a0 = constant_term(a) ∩ Interval(zero(T), T(Inf))
@assert !isempty(a0)
aux = one(a0^r)
a[0] = aux * a0
r == 0.5 && return sqrt(a)
if $T == Taylor1
l0 = findfirst(a)
# Index of first non-zero coefficient of the result; must be integer
!isinteger(r*l0) && throw(DomainError(a,
"""The 0-th order Taylor1 coefficient must be non-zero
to raise the Taylor1 polynomial to a non-integer exponent."""))
lnull = trunc(Int, r*l0 )
(lnull > a.order) && return $T( zero(aux), a.order)
c_order = l0 == 0 ? a.order : min(a.order, trunc(Int, r*a.order))
else
if iszero(a0)
throw(DomainError(a,
"""The 0-th order TaylorN coefficient must be non-zero
in order to expand `^` around 0."""))
end
c_order = a.order
end
#
c = $T(zero(aux), c_order)
aux0 = zero(c)
for k in eachindex(c)
TS.pow!(c, a, aux0, r, k)
end
return c
end

return c
function TS.pow!(res::$T{Interval{T}}, a::$T{Interval{T}}, aux::$T{Interval{T}},
r::S, k::Int) where {T<:Real, S<:Integer}
(r == 0) && return TS.one!(res, a, k)
(r == 1) && return TS.identity!(res, a, k)
(r == 2) && return TS.sqr!(res, a, k)
TS.power_by_squaring!(res, a, aux, r)
return nothing
end
end
end


for T in (:Taylor1, :TaylorN)
@eval function log(a::$T{Interval{S}}) where {S<:Real}
iszero(constant_term(a)) && throw(DomainError(a,
"""The 0-th order coefficient must be non-zero in order to expand `log` around 0."""))
"""The 0-th order coefficient must be non-zero in order to expand `log` around 0."""))
a0 = constant_term(a) ∩ Interval(S(0.0), S(Inf))
order = a.order
aux = log(a0)
Expand Down
107 changes: 48 additions & 59 deletions src/evaluate.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -254,39 +254,32 @@ Note that the syntax `a(vals)` is equivalent to
`evaluate(a)`; use a(b::Bool, x) corresponds to
evaluate(a, x, sorting=b).
"""
function evaluate(a::TaylorN, vals::NTuple{N,<:Number};
sorting::Bool=true) where {N}
function evaluate(a::TaylorN, vals::NTuple{N,<:Number}; sorting::Bool=true) where {N}
@assert get_numvars() == N
return _evaluate(a, vals, Val(sorting))
end

function evaluate(a::TaylorN, vals::NTuple{N,<:AbstractSeries};
sorting::Bool=false) where {N}
function evaluate(a::TaylorN, vals::NTuple{N,<:AbstractSeries}; sorting::Bool=false) where {N}
@assert get_numvars() == N
return _evaluate(a, vals, Val(sorting))
end

evaluate(a::TaylorN{T}, vals::AbstractVector{<:Number};
sorting::Bool=true) where {T<:NumberNotSeries} =
evaluate(a, (vals...,); sorting=sorting)
sorting::Bool=true) where {T<:NumberNotSeries} = evaluate(a, (vals...,); sorting=sorting)

evaluate(a::TaylorN{T}, vals::AbstractVector{<:AbstractSeries};
sorting::Bool=false) where {T<:NumberNotSeries} =
evaluate(a, (vals...,); sorting=sorting)
sorting::Bool=false) where {T<:NumberNotSeries} = evaluate(a, (vals...,); sorting=sorting)

evaluate(a::TaylorN{Taylor1{T}}, vals::AbstractVector{S};
sorting::Bool=false) where {T, S} =
evaluate(a, (vals...,); sorting=sorting)
sorting::Bool=false) where {T, S} = evaluate(a, (vals...,); sorting=sorting)

function evaluate(a::TaylorN{T}, s::Symbol, val::S) where
{T<:Number, S<:NumberNotSeriesN}
function evaluate(a::TaylorN{T}, s::Symbol, val::S) where {T<:Number, S<:NumberNotSeriesN}
ind = lookupvar(s)
@assert (1 ≤ ind ≤ get_numvars()) "Symbol is not a TaylorN variable; see `get_variable_names()`"
return evaluate(a, ind, val)
end

function evaluate(a::TaylorN{T}, ind::Int, val::S) where
{T<:Number, S<:NumberNotSeriesN}
function evaluate(a::TaylorN{T}, ind::Int, val::S) where {T<:Number, S<:NumberNotSeriesN}
@assert (1 ≤ ind ≤ get_numvars()) "Invalid `ind`; it must be between 1 and `get_numvars()`"
R = promote_type(T,S)
return _evaluate(convert(TaylorN{R}, a), ind, convert(R, val))
Expand All @@ -305,8 +298,7 @@ function evaluate(a::TaylorN{T}, ind::Int, val::TaylorN) where {T<:Number}
return _evaluate(a, ind, val)
end

evaluate(a::TaylorN{T}, x::Pair{Symbol,S}) where {T, S} =
evaluate(a, first(x), last(x))
evaluate(a::TaylorN{T}, x::Pair{Symbol,S}) where {T, S} = evaluate(a, first(x), last(x))

evaluate(a::TaylorN{T}) where {T<:Number} = constant_term(a)

Expand Down Expand Up @@ -338,14 +330,6 @@ function _evaluate(a::TaylorN{T}, vals::NTuple{N,<:Number}) where {N,T<:Number}
return suma
end

function _evaluate!(res::Vector{TaylorN{T}}, a::TaylorN{T}, vals::NTuple{N,<:TaylorN},
valscache::Vector{TaylorN{T}}, aux::TaylorN{T}) where {N,T<:Number}
@inbounds for homPol in eachindex(a)
_evaluate!(res[homPol+1], a[homPol], vals, valscache, aux)
end
return nothing
end

function _evaluate(a::TaylorN{T}, vals::NTuple{N,<:TaylorN}) where {N,T<:Number}
R = promote_type(T, TS.numtype(vals[1]))
suma = [TaylorN(zero(R), vals[1].order) for _ in eachindex(a)]
Expand Down Expand Up @@ -375,27 +359,37 @@ function _evaluate(a::TaylorN{T}, ind::Int, val::TaylorN{T}) where {T<:NumberNot
return suma
end

function _evaluate!(res::Vector{TaylorN{T}}, a::TaylorN{T}, vals::NTuple{N,<:TaylorN},
valscache::Vector{TaylorN{T}}, aux::TaylorN{T}) where {N,T<:Number}
@inbounds for homPol in eachindex(a)
_evaluate!(res[homPol+1], a[homPol], vals, valscache, aux)
end
return nothing
end

function _evaluate!(suma::TaylorN{T}, a::HomogeneousPolynomial{T}, ind::Int, val::T) where
{T<:NumberNotSeriesN}
order = a.order
orderTN = get_order()
if order == 0
suma[0] = a[1]*one(val)
suma[0][1] = a[1]*one(val)
return nothing
end
vv = val .^ (0:order)
# ct = @isonethread coeff_table[order+1]
ct = deepcopy(coeff_table[order+1])
vct = zero(coeff_table[order+1][1])
zct = zero(coeff_table[order+1][1])
for (i, a_coeff) in enumerate(a.coeffs)
iszero(a_coeff) && continue
if ct[i][ind] == 0
vpow = coeff_table[order+1][i][ind]
if vpow == 0
suma[order][i] += a_coeff
continue
end
vpow = ct[i][ind]
vct .= coeff_table[order+1][i]
zct[ind] = vpow
red_order = order - vpow
ct[i][ind] -= vpow
kdic = in_base(get_order(), ct[i])
ct[i][ind] += vpow
kdic = in_base(orderTN, vct - zct)
zct[ind] = 0
pos = pos_table[red_order+1][kdic]
suma[red_order][pos] += a_coeff * vv[vpow+1]
end
Expand All @@ -406,32 +400,33 @@ function _evaluate!(suma::TaylorN{T}, a::HomogeneousPolynomial{T}, ind::Int,
val::TaylorN{T}, aux::TaylorN{T}) where {T<:NumberNotSeriesN}
order = a.order
if order == 0
suma[0] = a[1]
suma[0][1] = a[1]
return nothing
end
vv = zero(suma)
ct = coeff_table[order+1]
vvaux = zero(vv)
za = zero(a)
for (i, a_coeff) in enumerate(a.coeffs)
iszero(a_coeff) && continue
if ct[i][ind] == 0
vpow = coeff_table[order+1][i][ind]
if vpow == 0
suma[order][i] += a_coeff
continue
end
za[i] = a_coeff
zero!(aux)
_evaluate!(aux, za, ind, one(T))
za[i] = zero(T)
vpow = ct[i][ind]
# vv = val ^ vpow
if constant_term(val) == 0
vv = val ^ vpow
zero!(vvaux)
power_by_squaring!(vv, val, vvaux, vpow)
else
for ordQ in eachindex(val)
zero!(vv, ordQ)
pow!(vv, val, vv, vpow, ordQ)
pow!(vv, val, vvaux, vpow, ordQ)
end
end
za[i] = a_coeff
zero!(aux)
_evaluate!(aux, za, ind, one(T))
za[i] = zero(a_coeff)
for ordQ in eachindex(suma)
mul!(suma, vv, aux, ordQ)
end
Expand All @@ -456,7 +451,7 @@ evaluate(A::AbstractArray{TaylorN{T}}) where {T<:Number} = evaluate.(A)
#function-like behavior for TaylorN
(p::TaylorN)(x) = evaluate(p, x)
(p::TaylorN)() = evaluate(p)
(p::TaylorN)(s::S, x) where {S<:Union{Symbol, Int}}= evaluate(p, s, x)
(p::TaylorN)(s::S, x) where {S<:Union{Symbol, Int}} = evaluate(p, s, x)
(p::TaylorN)(x::Pair) = evaluate(p, first(x), last(x))
(p::TaylorN)(x, v::Vararg{T}) where {T} = evaluate(p, (x, v...,))
(p::TaylorN)(b::Bool, x) = evaluate(p, x, sorting=b)
Expand All @@ -470,39 +465,33 @@ evaluate(A::AbstractArray{TaylorN{T}}) where {T<:Number} = evaluate.(A)


"""
evaluate!(x, δt, x0)
evaluate!(x, δt, dest)
Evaluates each element of `x::AbstractArray{Taylor1{T}}`,
representing the Taylor expansion for the dependent variables
of an ODE at *time* `δt`. It updates the vector `x0` with the
of an ODE at *time* `δt`. It updates the vector `dest` with the
computed values.
"""
function evaluate!(x::AbstractArray{Taylor1{T}}, δt::S,
x0::AbstractArray{T}) where {T<:Number, S<:Number}
x0 .= evaluate.( x, δt )
dest::AbstractArray{T}) where {T<:Number, S<:Number}
dest .= evaluate.( x, δt )
return nothing
end
# function evaluate!(x::AbstractArray{Taylor1{Taylor1{T}}}, δt::Taylor1{T},
# x0::AbstractArray{Taylor1{T}}) where {T<:Number}
# x0 .= evaluate.( x, Ref(δt) )
# # x0 .= evaluate.( x, δt )
# return nothing
# end

## In place evaluation of multivariable arrays
function evaluate!(x::AbstractArray{TaylorN{T}}, δx::Array{T,1},
x0::AbstractArray{T}) where {T<:Number}
x0 .= evaluate.( x, Ref(δx) )
dest::AbstractArray{T}) where {T<:Number}
dest .= evaluate.( x, Ref(δx) )
return nothing
end

function evaluate!(x::AbstractArray{TaylorN{T}}, δx::Array{TaylorN{T},1},
x0::AbstractArray{TaylorN{T}}; sorting::Bool=true) where {T<:NumberNotSeriesN}
x0 .= evaluate.( x, Ref(δx), sorting = sorting)
dest::AbstractArray{TaylorN{T}}; sorting::Bool=true) where {T<:NumberNotSeriesN}
dest .= evaluate.( x, Ref(δx), sorting = sorting)
return nothing
end

function evaluate!(a::TaylorN{T}, vals::NTuple{N,TaylorN{T}}, dest::TaylorN{T},
function _evaluate!(a::TaylorN{T}, vals::NTuple{N,TaylorN{T}}, dest::TaylorN{T},
valscache::Vector{TaylorN{T}}, aux::TaylorN{T}) where {N,T<:Number}
@inbounds for homPol in eachindex(a)
_evaluate!(dest, a[homPol], vals, valscache, aux)
Expand All @@ -518,7 +507,7 @@ function evaluate!(a::AbstractArray{TaylorN{T}}, vals::NTuple{N,TaylorN{T}},
# loop over elements of `a`
for i in eachindex(a)
(!iszero(dest[i])) && zero!(dest[i])
evaluate!(a[i], vals, dest[i], valscache, aux)
_evaluate!(a[i], vals, dest[i], valscache, aux)
end
return nothing
end
Expand Down
Loading

2 comments on commit 84a007f

@lbenet
Copy link
Member Author

@lbenet lbenet commented on 84a007f Jan 23, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@JuliaRegistrator register

@JuliaRegistrator
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Registration pull request created: JuliaRegistries/General/123587

Tip: Release Notes

Did you know you can add release notes too? Just add markdown formatted text underneath the comment after the text
"Release notes:" and it will be added to the registry PR, and if TagBot is installed it will also be added to the
release that TagBot creates. i.e.

@JuliaRegistrator register

Release notes:

## Breaking changes

- blah

To add them here just re-invoke and the PR will be updated.

Tagging

After the above pull request is merged, it is recommended that a tag is created on this repository for the registered package version.

This will be done automatically if the Julia TagBot GitHub Action is installed, or can be done manually through the github interface, or via:

git tag -a v0.18.3 -m "<description of version>" 84a007faac0a62d610c5cce2be8e859f7d3e698e
git push origin v0.18.3

Please sign in to comment.