# Error in dblquad integration with Unitful and MonteCarloMeasurements

**URL:** <https://discourse.julialang.org/t/error-in-dblquad-integration-with-unitful-and-montecarlomeasurements/110420>\
**Category:** Modelling & Simulations\
**Tags:** question, monte-carlo, quadgk, unitful\
**Created:** [February 19, 2024, 8:38pm UTC](https://discourse.julialang.org/t/error-in-dblquad-integration-with-unitful-and-montecarlomeasurements/110420 "2024-02-19T20:38:24Z")\
**Posts on this page:** 6\
**Page:** 1

<div class="post-metadata">

**Author:** ![GSavvidis](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gsavvidis/32/207145_2.png) [@GSavvidis](https://discourse.julialang.org/u/GSavvidis)\
**Post date:** [February 19, 2024, 8:38pm UTC](https://discourse.julialang.org/t/error-in-dblquad-integration-with-unitful-and-montecarlomeasurements/110420/1 "2024-02-19T20:38:24Z")

</div>

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:

```julia
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:

```julia
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

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [February 20, 2024, 6:23am UTC](https://discourse.julialang.org/t/error-in-dblquad-integration-with-unitful-and-montecarlomeasurements/110420/2 "2024-02-20T06:23:01Z")

</div>

Hello and welcome to the community! 👋

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-repl
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

```julia
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-repl
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
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.

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [February 20, 2024, 6:31am UTC](https://discourse.julialang.org/t/error-in-dblquad-integration-with-unitful-and-montecarlomeasurements/110420/3 "2024-02-20T06:31:15Z")

</div>

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

```julia
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-repl
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.

---

<div class="post-metadata">

**Author:** ![GSavvidis](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gsavvidis/32/207145_2.png) [@GSavvidis](https://discourse.julialang.org/u/GSavvidis)\
**Post date:** [February 21, 2024, 1:09am UTC](https://discourse.julialang.org/t/error-in-dblquad-integration-with-unitful-and-montecarlomeasurements/110420/4 "2024-02-21T01:09:20Z")

</div>

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:

```julia
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:

```julia
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!

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [February 21, 2024, 1:46am UTC](https://discourse.julialang.org/t/error-in-dblquad-integration-with-unitful-and-montecarlomeasurements/110420/5 "2024-02-21T01:46:03Z")

</div>

Try

```julia
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/](https://baggepinnen.github.io/MonteCarloMeasurements.jl/stable/overloading/)

---

<div class="post-metadata">

**Author:** ![GSavvidis](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gsavvidis/32/207145_2.png) [@GSavvidis](https://discourse.julialang.org/u/GSavvidis)\
**Post date:** [February 22, 2024, 7:00pm UTC](https://discourse.julialang.org/t/error-in-dblquad-integration-with-unitful-and-montecarlomeasurements/110420/6 "2024-02-22T19:00:05Z")

</div>

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:

```julia

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

```julia
        # 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