Error in dblquad integration with Unitful and MonteCarloMeasurements

Specific Domains Modelling & Simulations
, , ,
1

Hello everyone,

This is my first post so apologies in advance for any mistakes in the question format.

I have been trying to calculate an integral of two variables using dblquad. The function to be integrated contains units from Unitful and errors from MonteCarloMeasurements.
The following function:

using MultiQuad
using Unitful
using MonteCarloMeasurements

f(y, x) = y*x * (0.5 Β± 0.0001)
res = dblquad((y,x) -> f(y,x), 2u"MeV", 3u"MeV", 30u"eV", 50u"eV")

is an example which however generates the same error as my actual function:

ERROR: LoadError: setindex!() with non-isbitstype eltype is not supported by StaticArrays. Consider using SizedArray.
Stacktrace:
  [1] error(s::String)
    @ Base ./error.jl:35
  [2] setindex!(v::MVector{2, Particles{Float64, 2000}}, val::Particles{Float64, 2000}, i::Int64)
    @ StaticArrays ~/.julia/packages/StaticArrays/oOCPP/src/MArray.jl:39
  [3] (::HCubature.GenzMalik{2, Float64})(f::MultiQuad.var"#integrand#39"{MultiQuad.var"#arg2#38"{Quantity{Int64, 𝐋^2 𝐌 𝐓^-2, Unitful.FreeUnits{(MeV,), 𝐋^2 𝐌 𝐓^-2, nothing}}, Quantity{Int64, 𝐋^2 𝐌 𝐓^-2, Unitful.FreeUnits{(MeV,), 𝐋^2 𝐌 𝐓^-2, nothing}}, MultiQuad.var"#arg1#37"{var"#50#51", Quantity{Int64, 𝐋^2 𝐌 𝐓^-2, Unitful.FreeUnits{(eV,), 𝐋^2 𝐌 𝐓^-2, nothing}}, Quantity{Int64, 𝐋^2 𝐌 𝐓^-2, Unitful.FreeUnits{(eV,), 𝐋^2 𝐌 𝐓^-2, nothing}}}, Unitful.FreeUnits{(eV^2, MeV^2), 𝐋^8 𝐌^4 𝐓^-8, nothing}}}, a::SVector{2, Float64}, b::SVector{2, Float64}, norm::typeof(norm))
    @ HCubature ~/.julia/packages/HCubature/ZIutt/src/genz-malik.jl:137
  [4] hcubature_(f::MultiQuad.var"#integrand#39"{MultiQuad.var"#arg2#38"{Quantity{Int64, 𝐋^2 𝐌 𝐓^-2, Unitful.FreeUnits{(MeV,), 𝐋^2 𝐌 𝐓^-2, nothing}}, Quantity{Int64, 𝐋^2 𝐌 𝐓^-2, Unitful.FreeUnits{(MeV,), 𝐋^2 𝐌 𝐓^-2, nothing}}, MultiQuad.var"#arg1#37"{var"#50#51", Quantity{Int64, 𝐋^2 𝐌 𝐓^-2, Unitful.FreeUnits{(eV,), 𝐋^2 𝐌 𝐓^-2, nothing}}, Quantity{Int64, 𝐋^2 𝐌 𝐓^-2, Unitful.FreeUnits{(eV,), 𝐋^2 𝐌 𝐓^-2, nothing}}}, Unitful.FreeUnits{(eV^2, MeV^2), 𝐋^8 𝐌^4 𝐓^-8, nothing}}}, a::SVector{2, Float64}, b::SVector{2, Float64}, norm::typeof(norm), rtol_::Int64, atol::Int64, maxevals::Int64, initdiv::Int64, buf::Nothing)
    @ HCubature ~/.julia/packages/HCubature/ZIutt/src/HCubature.jl:110
  [5] hcubature_(f::Function, a::Vector{Int64}, b::Vector{Int64}, norm::Function, rtol::Int64, atol::Int64, maxevals::Int64, initdiv::Int64, buf::Nothing)
    @ HCubature ~/.julia/packages/HCubature/ZIutt/src/HCubature.jl:179
  [6] #hcubature#4
    @ ~/.julia/packages/HCubature/ZIutt/src/HCubature.jl:235 [inlined]
  [7] hcubature
    @ ~/.julia/packages/HCubature/ZIutt/src/HCubature.jl:235 [inlined]
  [8] _hcubature_wrapper_2D(arg::var"#50#51", x1::Quantity{Int64, 𝐋^2 𝐌 𝐓^-2, Unitful.FreeUnits{(MeV,), 𝐋^2 𝐌 𝐓^-2, nothing}}, x2::Quantity{Int64, 𝐋^2 𝐌 𝐓^-2, Unitful.FreeUnits{(MeV,), 𝐋^2 𝐌 𝐓^-2, nothing}}, y1::Quantity{Int64, 𝐋^2 𝐌 𝐓^-2, Unitful.FreeUnits{(eV,), 𝐋^2 𝐌 𝐓^-2, nothing}}, y2::Quantity{Int64, 𝐋^2 𝐌 𝐓^-2, Unitful.FreeUnits{(eV,), 𝐋^2 𝐌 𝐓^-2, nothing}}; kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ MultiQuad ~/.julia/packages/MultiQuad/IdgOj/src/MultiQuad.jl:380
  [9] _hcubature_wrapper_2D(arg::Function, x1::Quantity{Int64, 𝐋^2 𝐌 𝐓^-2, Unitful.FreeUnits{(MeV,), 𝐋^2 𝐌 𝐓^-2, nothing}}, x2::Quantity{Int64, 𝐋^2 𝐌 𝐓^-2, Unitful.FreeUnits{(MeV,), 𝐋^2 𝐌 𝐓^-2, nothing}}, y1::Quantity{Int64, 𝐋^2 𝐌 𝐓^-2, Unitful.FreeUnits{(eV,), 𝐋^2 𝐌 𝐓^-2, nothing}}, y2::Quantity{Int64, 𝐋^2 𝐌 𝐓^-2, Unitful.FreeUnits{(eV,), 𝐋^2 𝐌 𝐓^-2, nothing}})
    @ MultiQuad ~/.julia/packages/MultiQuad/IdgOj/src/MultiQuad.jl:366
 [10] dblquad(arg::Function, x1::Quantity{Int64, 𝐋^2 𝐌 𝐓^-2, Unitful.FreeUnits{(MeV,), 𝐋^2 𝐌 𝐓^-2, nothing}}, x2::Quantity{Int64, 𝐋^2 𝐌 𝐓^-2, Unitful.FreeUnits{(MeV,), 𝐋^2 𝐌 𝐓^-2, nothing}}, y1::Quantity{Int64, 𝐋^2 𝐌 𝐓^-2, Unitful.FreeUnits{(eV,), 𝐋^2 𝐌 𝐓^-2, nothing}}, y2::Quantity{Int64, 𝐋^2 𝐌 𝐓^-2, Unitful.FreeUnits{(eV,), 𝐋^2 𝐌 𝐓^-2, nothing}}; method::Symbol, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ MultiQuad ~/.julia/packages/MultiQuad/IdgOj/src/MultiQuad.jl:453
 [11] dblquad(arg::Function, x1::Quantity{Int64, 𝐋^2 𝐌 𝐓^-2, Unitful.FreeUnits{(MeV,), 𝐋^2 𝐌 𝐓^-2, nothing}}, x2::Quantity{Int64, 𝐋^2 𝐌 𝐓^-2, Unitful.FreeUnits{(MeV,), 𝐋^2 𝐌 𝐓^-2, nothing}}, y1::Quantity{Int64, 𝐋^2 𝐌 𝐓^-2, Unitful.FreeUnits{(eV,), 𝐋^2 𝐌 𝐓^-2, nothing}}, y2::Quantity{Int64, 𝐋^2 𝐌 𝐓^-2, Unitful.FreeUnits{(eV,), 𝐋^2 𝐌 𝐓^-2, nothing}})
    @ MultiQuad ~/.julia/packages/MultiQuad/IdgOj/src/MultiQuad.jl:440
 [12] top-level scope
    @ ~/postdoc/G3/MC_CEvNS/scripts/reactor_cevns_events/reactor_cevns_events_v3.jl:139
 [13] include(fname::String)
    @ Base.MainInclude ./client.jl:478
 [14] top-level scope
    @ REPL[11]:1

Do you know if it’s possible to overcome this problem?

Thanks in advance

2

Hello and welcome to the community! :wave:

This appears to be a limitation of the integration package and StaticArrays. There is a workaround possible, which is to use particles with static arrays as storage instead of the standard-array particles.

If you switch the operator Β± to βˆ“ (\mp), an instance of StaticParticles{Float64, 100} is created

julia> (0.5 βˆ“ 0.0001)
0.5 Β± 9.99e-5 StaticParticles{Float64, 100}

notice that only 100 particles are used by default rather than the 2000 particles used by Β±. The reason for this is that large static arrays can incur a very large compilation time.

If you make this change, you will then hit the error

ERROR: Comparison of uncertain values using comparison mode
safe failed.
Comparison operators are not well defined for uncertain values.
Call `unsafe_comparisons(true)` to enable comparison operators
for particles using the current reduction function pmean.
Change this function using `set_comparison_function(f)`.
For safety reasons, the default safe comparison function
is maximally conservative and tests if the extreme values of
the distributions fulfil the comparison operator.

as the error says, you need to tell the system how to compare two uncertain values. If you deem that a naive comparison is okay (think this through), you can set

julia> unsafe_comparisons(true)
[ Info: Unsafe comparisons using the function `pmean`
has been enabled globally.
Use `@unsafe` to enable in a local expression only or
`unsafe_comparisons(false)` to turn off unsafe comparisons

julia> res = dblquad((y,x) -> f(y,x), 2u"MeV", 3u"MeV", 30u"eV", 50u"eV")

This makes the integration run, but it takes a very long time, making me question whether this integration method is suitable to combine with uncertain numbers.

Reducing the tolerance to

julia> res = dblquad((y,x) -> f(y,x), 2u"MeV", 3u"MeV", 30u"eV", 50u"eV", rtol=1e-3, atol=1e-5)
(1000.0000000000006 eVΒ² MeVΒ² Β± 0.0789374708719211 eVΒ² MeVΒ², 0.22170130367696628 eVΒ² MeVΒ² Β± 0.15892964329986264 eVΒ² MeVΒ²)

Gives a result very quickly, do please verify that it seems reasonable.

3

A somewhat safer approach is to register a function that operates with an uncertain value as input, for example like this

function uncertain_int(val)
    f(y, x) = y*x * val
    res = dblquad((y,x) -> f(y,x), 2u"MeV", 3u"MeV", 30u"eV", 50u"eV")#, rtol=1e-3, atol=1e-5)
    @assert pmaximum(res[2]) < 1e-10u"eV^2*MeV^2" # Verify small integration error
    res[1]
end
MonteCarloMeasurements.register_primitive(uncertain_int)

res = uncertain_int(0.5 Β± 0.0001)

julia> res = uncertain_int(0.5 Β± 0.0001)
1000.0000000000009 eVΒ² MeVΒ² Β± 0.19998463057861196 eVΒ² MeVΒ² Particles{Quantity{Float64, 𝐋⁸ 𝐌⁴ 𝐓⁻⁸, Unitful.FreeUnits{(eVΒ², MeVΒ²), 𝐋⁸ 𝐌⁴ 𝐓⁻⁸, nothing}}, 2000}

this results in a larger error estimate for the integral, which is an indication that the soundness of the initial attempt above was thwarted by the uncertain comparisons.

4

Thanks a lot for the quick reply, I appreciate it!
I managed to reproduce successfully your second solution! Although the given example works, when I tried to apply it to my function it gave me the same error. After playing around with the code, I found that if val is a higher dimensional object such as a Vector or Matrix as:

coef = [3.217 Β± 4.09e-2, -3.111 Β± 2.34e-2, 1.395 Β± 4.88e-3,
        -3.690e-1 Β± 6.08e-4, 4.445e-2 Β± 7.77e-5, -2.053e-3 Β± 6.79e-6]

Then the following:

function uncertain_int(val)
    f(y, x) = y*x * sum([value for value in val])
    res = dblquad((y,x) -> f(y,x), 2u"MeV", 3u"MeV", 30u"eV", 50u"eV")#, rtol=1e-3, atol=1e-5)
    @assert pmaximum(res[2]) < 1e-10u"eV^2*MeV^2" # Verify small integration error
    res[1]
end

MonteCarloMeasurements.register_primitive(uncertain_int)
res = uncertain_int(coef)

fails with the same error. And the function that I want to integrate contains vectors and matrices with uncertain values such as coeff which have to be calculated in conjunction with the x and y variables.

Please let me know if I need to provide more details! Thanks again!

5

Try

res = @bymap uncertain_int(coef)

In general, MonteCarloMeasurements will struggle if it cannot determine how to propagate particles through a particular function. If you have a function with problematic elements, such as uncertain control flow or too strict typing, you have to find a way to tell MCM how to handle it. register_primitive is one such way, but it relies on Particles appearing as input arguments to the function. In the latter example here, you instead have Array{Particles} which is not a recognized type. The macro @bymap does try a bit harder to look for β€œthings with particles in them”, and handle this correctly, but it can only do so much. I think it handles arrays, tuples and named tuples atm. Here are some more details on how to handle β€œdifficult functions”:
https://baggepinnen.github.io/MonteCarloMeasurements.jl/stable/overloading/

6

Ok, it took me a while but I managed to put everything together and make it work. The job is done using bymap as you mentioned and that really helped me, thank you! MCM is very β€œsensitive”, as you mentioned, when it comes to propagate particles in different functions/containers. In my case I have a vector of vectors such as the following:


# coefficients
coef_1 = [3.217 Β± 4.09e-2, -3.111 Β± 2.34e-2, 1.395 Β± 4.88e-3,
        -3.690e-1 Β± 6.08e-4, 4.445e-2 Β± 7.77e-5, -2.053e-3 Β± 6.79e-6]

coef_2 = [4.833e-1 Β± 1.24e-1, 1.927e-1 Β± 5.86e-2, -1.283e-1 Β± 1.11e-2,
        -6.762e-3 Β± 1.92e-3, 2.233e-3 Β± 2.84e-4, -1.536e-4 Β± 2.86e-5]

A_coef = [coef_1, coef_2]

# Q-values
Q_1 = (202.36 Β± 0.26)u"MeV"
Q_2 = (205.99 Β± 0.52)u"MeV"

Qvalues = [Q_1, Q_2]

# uncertain quantities
VV = [A_coef, [Qvalues]]

Then these quantities enter the double integral and that’s the messy part. I put here the solution I found to make it work because it may be useful for others:

        # function to be bymapped
        function uncertain_integral(V...)

            # function to be integrated
            function spectrum(EΞ½, Enr)
                
                # For my case, here I call my uncertain parameters from V
                # and define the calculations

            end

            res = dblquad((y, x) -> spectrum(y, x), Eβ‚™α΅£α΅’, Eβ‚™α΅£α΅’β‚Šβ‚, x->EΞ½min, x->EΞ½max,
                          rtol=1e-5, atol=1e-8)
            res[1]
        end

# splatting does the trick
res = bymap(uncertain_integral, (VV...)...)

What I learned is that bymap or Worspace are a must when it comes to propagation of particles with non-basic operations.

Thanks again