Trouble with Measurements.jl and Bessel Functions

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 :slight_smile:

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 :slight_smile: 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.