GitLab

Broadcastability

If a and b are instances of Series, one would expect [a, a] .+ b to return [a+b, a+b]. This can be achieved by defining:

import Base: broadcastable
using FormalSeries

Base.broadcastable(s::Series) = Ref(s)

a = Series((1.0, 1.0))
b = Series((1.0, 0.0))

[a, a] .+ b

This appears to be the recommended approach if the type should behave like a scalar (Julia Documentation on Broadcasting). I have not come across a situation where this should not be the case, similar to uwreal objects in ADerrors.

Derivatives of numerical intergals: compatibility with QuadGK

I needed to take derivatives of numerical integrals with respect to a variable. For the case of unidimensional integrals, I have used QuadGK, and I had to define these expansions:

import Base: iterate, isnan, isinf, isless
import LinearAlgebra: norm
using FormalSeries

iterate(s::FormalSeries.Series) = (s, nothing)

norm(s::FormalSeries.Series{T,N}) where {T <: Real, N} = s
norm(s::FormalSeries.Series{T,N}) where {T <: Complex, N} = sqrt(real(s)^2 + imag(s)^2)

isnan(s::FormalSeries.Series) = any(isnan.(s.c))
isinf(s::FormalSeries.Series) = any(isinf.(s.c))

isless(s1::FormalSeries.Series, s2::FormalSeries.Series) = isless(s1.c[1], s2.c[1])

I would say all but isless are reasonable. I defined iterate by looking at the definition for a Number. @edit iterate(2) leads to:

iterate(x::Number) = (x, nothing)

Then, for norm, I defined the operation that would be done at order 0. And isnan and isinf simply check each element.

The most problematic definition is isless, but I believe it is necessary for routines where a target precision is sought, like in numerical integration. This definition would allow stochastic things like the accept-reject step of Metropolis to "work" without error, and perhaps more cases where the derivative is not defined... but maybe it's more useful that FormalSeries works out of the box for relatively simple things like numerical integrals, and it should be up to the end user to decide if they have used FormalSeries correctly?

Anyway, with these definitions, expansions of numerical integrals with QuadGK can already be calculated:

using FormalSeries
using QuadGK

function ana_num(I, theta)
    return phi -> (I * (1 - cos(phi)) + theta * sin(phi) / (2pi)) * exp(-I * (1 - cos(phi)) - theta * sin(phi) / (2pi))
end

function ana_den(I, theta)
    return phi -> exp(-I * (1 - cos(phi)) - theta * sin(phi) / (2pi))
end

function integral_ana(I, theta)
    return quadgk(ana_num(I, theta), -pi, pi)[1] / quadgk(ana_den(I, theta), -pi, pi)[1]
end

theta = Series((0.0, 1.0, 0.0))

integral_ana(5.0, theta)