TLDR: I think Dual should be <:Real and that Complex{Dual} leads to more generic code than Dual{Complex}.
Hi all,
I am toying with mixed complex and dual numbers.
My use case is numerical integration of a nonlinear PDE (2D Navier-Stokes) and the variational equations, i.e. equations that describe the evolution of infinitesimal around the reference trajectory. The two equations can be marched simultaneously time by letting the solution be described using dual numbers, where the \epsilon component the perturbation. With julia type generic code, there is not much to be done to get the equations solved when the code for the main PDE is available. Quite neat!
I solve the PDE using a Fourier spectral technique, and thus the state variables are complex coefficients. FFTs are used to transform the solution to the physical domain, evaluate nonlinearity there, and then backtransform to the spectral domain.
When the solution is represented using dual numbers, what you really want to do is to perform transforms of the real and perturbation parts separately. This can be achieved using guru FFTW interface, which allows specifying the stride pattern and the number of to perform (it requires hacking into FFTW.jl, at the moment, but it can be done).
The main issue (an annoyance, really) I have is that, currently [1], dual numbers are defined like this
const ReComp = Union{Real,Complex}
immutable Dual{T<:ReComp} <: Number
value::T
epsilon::T
end
This definition does not play nice with the signature of methods that you would define in type-generic code, and requires (AFAICS), specialised signatures for dual types. For instance, what would be written naturally as
# this method signature would work for `T` begin a Float64 or a Dual{Float64} if Dual <: Real
function rfft!(U::AbstractVector{Complex{T}}, u::AbstractVector{T}) where {T<:Real}
# transform u to U
# if T < Dual, transform both value and epsilon parts separately, using
# FFTW guru interface with custom strides pattern
end
needs custom Dual number methods such as
function rfft!(U::AbstractVector{Dual{Complex{T}}}, u::AbstractVector{Dual{T}}) where {T<:Real}
# do stuff
end
This, and many other use cases where you have method signature such as
foo(a::Complex{T}, b::T) where {T} = # ...
as well as in type definitions, e.g.
struct Foo{T<:Real, A::AbstraArray{T}}
data::A
end
struct CFoo{T<:Complex, A::AbstraArray{T}}
data::A
end
would be solved gracefully if dual numbers where defined as
struct Dual{T<:Real} <: Real
value::T
epsilon::T
end
i.e. they where subtypes of real. One would then use Complex{Dual{T}} for automatic differentiation of functions with complex numbers, or other purposes.
I would be interested in knowing if there are cases where the current Dual{Complex} is better than Complex{Dual}.
Best,
Davide
.