register_primitive doesn't work for Uniful Quantity objects wrapping Particles

[klaff]

Consider the functions f1 and f2:

using MonteCarloMeasurements
using Unitful

function f1(Vi)
    if Vi ≤ 0.0
        return 0.0
    elseif Vi ≥ 1.0
        return 1.0
    else
        Vi
    end
end

function f2(Vi)
    if Vi ≤ 0.0u"V"
        return 0.0u"V"
    elseif Vi ≥ 1.0u"V"
        return 1.0u"V"
    else
        return Vi
    end
end

then f1(0.0..1.0) succeeds, but f1(-0.5..1.5) fails. This is understandable as the second distribution spans the discontinuities. If I then register_primitive(f1) and subsequently execute f1(-0.5..1.5), it succeeds and inspection of the result shows it is as I expect.

However, if I try the same thing with f2 (which I intend to be the exact same function except that now I’m using Unitful and everything should be in units compatible with u"V"), f2 fails, regardless of register_primitive(f2). Specifically f2((0.0..1.0)u"V") works, but f2((-0.5..1.5)u"V") always fails.

@simeonschaub diagnosed this "The problem here is that (-0.5..1.5)u"V" is creating an object of type Quantity{Particles}, whereas register_primitive only overloads f2 for Particles, so that second method never gets called."

See https://discourse.julialang.org/t/problem-combining-unitful-and-montecarlomeasurements/35418

[baggepinnen]

Hey and thanks for the issue. You are not alone in realizing the difficulty of having Particles interacting with other spcial number types. Indeed, the same problems arises when dealing with Complex{Particles} / Particles{Complex}. Some discussion on the topic is held in #22

Simeon was right in asserting that the problem arises from the type dispatching in <= is not Particles. I thus tried the reverse representation, i.e., Particles{Quantity}, which did work, although the resulting Particles{Quantity} could not be displayed properly due to a missing method specifying how to round the result properly for display. Unfortunately, I have a very limited understanding of how Unitful works, but maybe you can spot the problem?

using Unitful, MonteCarloMeasurements
function f2(Vi)
    if Vi ≤ 0.0u"V"
        return 0.0u"V"
    elseif Vi ≥ 1.0u"V"
        return 1.0u"V"
    else
        return Vi
    end
end
p = Particles(rand(Uniform(-0.5,1.5), 500) .* 1u"V")
register_primitive(f2)
p2 = f2(p)

julia> mean(p2)
0.4808899851070955 V

julia> std(p2)
0.41605814592938534 V

[klaff]

It appears to be doing the right calculations based on extrema(p2), but plotting it suggests that the uniform distribution is no longer being performed systematically.

Compare histogram([x.val for x in p2.particles]) to plot(f1(-0.5..1.5)). (You probably need to do register_primitive(f1) first.

EDIT: Also, examining a few points of p1=f1(-.5..1.5) shows data points on .0002 intervals, while the p2 points look, well, random.

My current understanding of both Unitful and MonteCarloMeasurements is on the weak side, although I am motivated to get them working together. I'm using them for some electronics design work, and it's awesome to start my calculations with distributions rather than think nominally and then come through with tolerancing as a second step, and I want to have some help keeping out unit errors at the same time.

[baggepinnen]

The line Particles(rand(Uniform(-0.5,1.5), 500) .* 1u"V") creates particles by first sampling from the Uniform distribution. You can choose to call the systematic sampling instead like this

 Particles(systematic_sample(500,Uniform(-0.5,1.5)) .* 1u"V")

Standard sampling should work okay if the number of samples is large, but systematic sampling has lower variance, that's why it is the default internally in MonteCarloMeasurements.

[baggepinnen]

@klaff I've added some support for unitful in the PR linked above. If you are still interested in this it would be nice if you could try it out and see how it works for you. I'm not a user of Unitful myself so I'm not sure what to test, Your example should now work quite smoothly like this

using Unitful, MonteCarloMeasurements, Test
function f2(Vi)
    if Vi ≤ 0.0u"V"
        return 0.0u"V"
    elseif Vi ≥ 1.0u"V"
        return 1.0u"V"
    else
        return Vi
    end
end

p1 = Particles(500, Uniform(-0.5,1.5)) * 1u"V"
p2 = Particles(500, Uniform(-0.5,1.5)) * u"V"

@test typeof(p1) == typeof(p2)

register_primitive(f2)
p3 = f2(p1)
@test extrema(p3) == (0.0u"V", 1.0u"V")

[baggepinnen]

I'll close this issue now since the PR above has been merged. Feel free to comment if you run into any problems!