I am using Measurements.jl for error propagation in my code and have come across an error when trying to use it with a Bessel function with a complex input.
This works
using Measurements
x = 1 ± 0.5
y = besselk(1. , x)
and this works
y = besselk(1., complex(1.,1.))
But this does not
using Measurements
x = 1 ± 0.5
y = besselk(1. , complex(1. , x))
and gives the following error
LoadError: MethodError: no method matching besselk(::Measurements.Measurement{Float64}, ::Complex{Measurements.Measurement{Float64}})
Closest candidates are:
besselk(::T<:AbstractFloat, ::Complex{T<:AbstractFloat}) where T<:AbstractFloat at C:\Users\xxxxxx.julia\v0.6\SpecialFunctions\src\bessel.jl:457
besselk(::Real, ::Complex) at C:\Users\xxxxxx.julia\v0.6\SpecialFunctions\src\bessel.jl:454
besselk(::Any…; kwargs…) at deprecated.jl:1294
…
while loading TestCase.jl, in expression starting on line 7
in besselk at SpecialFunctions\src\bessel.jl:455
I was under the impression that Measurements.jl was compatible with Julia functions that used AbstractFloat data types?
Is there something I am missing?
Measurements.jl
supports any function whose argument’s type is no more specific than AbstractFloat
and internally calls functions already supported, which is not the case for special functions with complex argument. I added support only for special functions with real argument, I should do the same for complex arguments as well. Patches are welcome
Many thanks. If I can figure out a patch I will be sure to submit it.
Look to definition of besselh in https://github.com/JuliaPhysics/Measurements.jl/blob/master/src/math.jl. It’s a function of real argument, but the output is complex. The result
function is extensively commented, at the beginning of the file. The first argument is the nominal value, the second is the derivative, the last one is the input argument
Uhm, looking better, I believe another result
method should be defined to support this case So far, I didn’t have to worry about functions of complex arguments because most of such functions works out of the box (because they internally call supported functions), the only problem are special functions of complex argument because they only wrap external C functions
@Jason_Kilpatrick I just added the missing method for result
function. The arguments are described in the comments before the new method.
In order to be able to add support to besselj0
you should compute the derivatives of the real and imaginary parts of the result with respect to the real and imaginary parts of the argument, so they’re 4 derivatives in total. If you find out those quantities, besselj0
with complex argument can be readily defined.
In the meantime, you can propagate the uncertainty of besselj0
with complex argument (or any complex-valued function of one complex argument) by using the following macro, similar to the @uncertain
macro which works only for functions of real arguments:
julia> using Measurements
julia> macro uncertain_complex(expr::Expr)
f = esc(expr.args[1]) # Function name
a = esc(expr.args[2]) # Argument, of Complex{Measurement} type
return :( Measurements.result($f(Measurements.value($a)),
vcat(Calculus.gradient(x -> real($f(complex(x[1], x[2]))),
[reim(Measurements.value($a))...]),
Calculus.gradient(x -> imag($f(complex(x[1], x[2]))),
[reim(Measurements.value($a))...])),
$a) )
end
@uncertain_complex (macro with 1 method)
julia> @uncertain_complex besselj0(complex(1, 1 ± 0.5))
(0.9376084768060292 ± 0.18251401441177362) - (0.4965299476091221 ± 0.30708016745937555)im
julia> @uncertain_complex besselj0(complex(1 ± 0.5))
(0.7651976865579666 ± 0.22002529286918884) + (0.0 ± 0.0)im
julia> besselj0(1 ± 0.5)
0.7651976865579666 ± 0.22002529287246675
The last two lines shows that the macro correctly works in the case of a complex measurement with null imaginary parts (within numerical accuracy of Calculus.gradient
), which is a good test case. This macro requires checking out master
branch of the package, but I’ll probably tag a new version next week.
Many thanks for the quick solution.