Giter Site home page Giter Site logo

baggepinnen / montecarlomeasurements.jl Goto Github PK

View Code? Open in Web Editor NEW
259.0 7.0 17.0 4.77 MB

Propagation of distributions by Monte-Carlo sampling: Real number types with uncertainty represented by samples.

Home Page: https://baggepinnen.github.io/MonteCarloMeasurements.jl/stable/

License: MIT License

Julia 100.00%
monte-carlo monte-carlo-sampling numeric-types uncertainty-propagation uncertainty-sampling error-propagation uncertainties physical-quantities error-analysis robust-optimization

montecarlomeasurements.jl's Introduction

logo Build Status codecov Documentation, stable Documentation, latest arXiv article

Imagine you had a type that behaved like your standard Float64 but it really represented a probability distribution like Gamma(0.5) or MvNormal(m, S). Then you could call y=f(x) and have y be the probability distribution y=p(f(x)). This package gives you such a type.

This package facilitates working with probability distributions by means of Monte-Carlo methods, in a way that allows for propagation of probability distributions through functions. This is useful for, e.g., nonlinear uncertainty propagation. A variable or parameter might be associated with uncertainty if it is measured or otherwise estimated from data. We provide two core types to represent probability distributions: Particles and StaticParticles, both <: Real. (The name "Particles" comes from the particle-filtering literature.) These types all form a Monte-Carlo approximation of the distribution of a floating point number, i.e., the distribution is represented by samples/particles. Correlated quantities are handled as well, see multivariate particles below.

Although several interesting use cases for doing calculations with probability distributions have popped up (see Examples), the original goal of the package is similar to that of Measurements.jl, to propagate the uncertainty from input of a function to the output. The difference compared to a Measurement is that Particles represent the distribution using a vector of unweighted particles, and can thus represent arbitrary distributions and handle nonlinear uncertainty propagation well. Functions like f(x) = x², f(x) = sign(x) at x=0 and long-time integration, are examples that are not handled well using linear uncertainty propagation ala Measurements.jl. MonteCarloMeasurements also support correlations between quantities.

A number of type Particles behaves just as any other Number while partaking in calculations. Particles also behave like a distribution, so after a calculation, an approximation to the complete distribution of the output is captured and represented by the output particles. mean, std etc. can be extracted from the particles using the corresponding functions pmean and pstd. Particles also interact with Distributions.jl, so that you can call, e.g., Normal(p) and get back a Normal type from distributions or fit(Gamma, p) to get a Gammadistribution. Particles can also be asked for maximum/minimum, quantile etc. using functions with a prefix p, i.e., pmaximum. If particles are plotted with plot(p), a histogram is displayed. This requires Plots.jl. A kernel-density estimate can be obtained by density(p) is StatsPlots.jl is loaded.

Below, we show an example where an input uncertainty is propagated through σ(x)

transformed densities

In the figure above, we see the probability-density function of the input p(x) depicted on the x-axis. The density of the output p(y) = f(x) is shown on the y-axis. Linear uncertainty propagation does this by linearizing f(x) and using the equations for an affine transformation of a Gaussian distribution, and hence produces a Gaussian approximation to the output density. The particles form a sampled approximation of the input density p(x). After propagating them through f(x), they form a sampled approximation to p(y) which correspond very well to the true output density, even though only 20 particles were used in this example. The figure can be reproduced by examples/transformed_densities.jl.

Quick start

using MonteCarloMeasurements, Plots
a = π ± 0.1 # Construct Gaussian uncertain parameters using ± (\\pm)
# Particles{Float64,2000}
#  3.14159 ± 0.1
b = 2  0.1 # ∓ (\\mp) creates StaticParticles (with StaticArrays)
# StaticParticles{Float64,100}
#  2.0 ± 0.0999
pstd(a)     # Ask about statistical properties
# 0.09999231528930486
sin(a)      # Use them like any real number
# Particles{Float64,2000}
#  1.2168e-16 ± 0.0995
plot(a)     # Plot them
b = sin.(1:0.1:5)  0.1; # Create multivariate uncertain numbers
plot(b)                   # Vectors of particles can be plotted
using Distributions
c = Particles(500, Poisson(3.)) # Create uncertain numbers distributed according to a given distribution
# Particles{Int64,500}
#  2.882 ± 1.7

For further help, see the documentation, the examples folder or the arXiv paper.

montecarlomeasurements.jl's People

Contributors

aplavin avatar baggepinnen avatar florisvankempen avatar gdalle avatar github-actions[bot] avatar joachimbrand avatar juliatagbot avatar masonprotter avatar sethaxen avatar simeonschaub avatar simontreillou avatar

Stargazers

 avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar

Watchers

 avatar  avatar  avatar  avatar  avatar  avatar  avatar

montecarlomeasurements.jl's Issues

Problem with custom Sampler/Distribution

I am trying to use MonteCarloMeasurements.jl with a custom Sampler/Distribution from Distributions.jl and I am having some problems. The first question I have is: can I use a simple Sampleable or is a fully fledged Distribution required to create Particles? The reason I would rather use a sampler is that it is very easy to sample from this circular distribution (Projected Normal) but as far as I know there are no closed form formulas for the CDF, quantiles, etc.

This is what I tried:

using Distributions

import Base: rand
import Random

struct ProjectedNormal <: Sampleable{Univariate,Continuous}
    μ::Float64
    ξ::Float64
end

function rand(rng::Random.AbstractRNG, pn::ProjectedNormal)
    r = rand(rng)
    return atan(pn.ξ*cos(pn.μ)+randn(), pn.ξ*sin(pn.μ)+randn())
end

pn1 = ProjectedNormal(1.0,2.0)

using MonteCarloMeasurements

y = MonteCarloMeasurements.Particles(500, pn1)

MethodError: no method matching Particles(::Int64, ::ProjectedNormal)
Closest candidates are:
  Particles(::Integer) at C:\Users\arrigo\.julia\packages\MonteCarloMeasurements\Jcg6F\src\types.jl:70
  Particles(::Integer, !Matched::Distribution{Multivariate,S} where S<:ValueSupport) at C:\Users\arrigo\.julia\packages\MonteCarloMeasurements\Jcg6F\src\types.jl:77
  Particles(::Integer, !Matched::Distribution{#s32,VS} where #s32; kwargs...) where VS at C:\Users\arrigo\.julia\packages\MonteCarloMeasurements\Jcg6F\src\types.jl:70
  ...

Stacktrace:
 [1] top-level scope at In[12]:2

Any ideas? A related question is: what is the minimum number of methods that I should implement for a Distribution? I used skellam.jl as an example but for some reason I can't even sample from it. This is the code:

module ProjectedNormals

using Random
using Distributions
using SpecialFunctions

struct ProjectedNormal{T<:Real} <: ContinuousUnivariateDistribution
	μ::T	# Mean angle
	ξ::T	# SNR (A/sigma)
	
    function ProjectedNormal{T}(μ::T, ξ::T) where {T <: Real}
        return new{T}(μ, ξ)
    end
end

ProjectedNormal(μ::Real, ξ::Real) = ProjectedNormal(promote(μ, ξ)...)
ProjectedNormal(μ::Integer, ξ::Integer) = Skellam(float(μ), float(ξ))
ProjectedNormal() = ProjectedNormal(0.0, 1.0, check_args=false)

#@distr_support ProjectedNormal -π π

#### Parameters
params(d::ProjectedNormal) = (d.μ, d.ξ)
partype(::ProjectedNormal{T}) where {T} = T

function pdf(d::ProjectedNormal, θ::Real)
	μ, ξ = params(d)

	1/2π*exp(-ξ^2/2*(1+sqrt(π/2)*ξ*cos(θ-μ)*exp(ξ^2*cos(θ-μ)^2/2)*(1+erf(ξ*cos(θ-μ)/sqrt(2)))))	
end

#### Sampling
#@rand_rdist(ProjectedNormal)
function rand(d::ProjectedNormal)
	μ, ξ = params(d)
	
	atan(ξ*cos(μ)+randn(), ξ*sin(μ)+randn())
end

function rand(rng::AbstractRNG, d::ProjectedNormal)
	μ, ξ = params(d)

	atan(cos(μ)+randn(rng),sin(μ)+randn(rng))
end

end

and this is the error that I get when I try to sample from a ProjectedNormal:

MethodError: no method matching iterate(::ProjectedNormals.ProjectedNormal{Float64})
Closest candidates are:
  iterate(!Matched::Core.SimpleVector) at essentials.jl:600
  iterate(!Matched::Core.SimpleVector, !Matched::Any) at essentials.jl:600
  iterate(!Matched::ExponentialBackOff) at error.jl:218
  ...

Stacktrace:
 [1] copyto!(::Array{Float64,1}, ::ProjectedNormals.ProjectedNormal{Float64}) at .\abstractarray.jl:722
 [2] _collect(::UnitRange{Int64}, ::ProjectedNormals.ProjectedNormal{Float64}, ::Base.HasEltype, ::Base.HasLength) at .\array.jl:566
 [3] collect(::ProjectedNormals.ProjectedNormal{Float64}) at .\array.jl:560
 [4] #quantile#53(::Bool, ::typeof(quantile), ::ProjectedNormals.ProjectedNormal{Float64}, ::Float64) at D:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.3\Statistics\src\Statistics.jl:966
 [5] quantile(::ProjectedNormals.ProjectedNormal{Float64}, ::Float64) at D:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.3\Statistics\src\Statistics.jl:966
 [6] rand(::Random._GLOBAL_RNG, ::ProjectedNormals.ProjectedNormal{Float64}) at C:\Users\arrigo\.julia\packages\Distributions\kPXE9\src\univariates.jl:175
 [7] rand(::ProjectedNormals.ProjectedNormal{Float64}) at C:\Users\arrigo\.julia\packages\Distributions\kPXE9\src\genericrand.jl:22
 [8] top-level scope at In[12]:1

InexactError

@show a = Particles(1000,Bernoulli(0.1))
@show a * 2
@show a * 2.05

returns:

a = Particles(1000, Bernoulli(0.1)) = 0.094 ± 0.29
a * 2 = 0.188 ± 0.58
InexactError: Int64(2.05)

Intuitively I feel like the third line should work similar to the second, but maybe I don't understand some of the limitations of the particle technique used here?

Place package in JuliaPhysics?

As has been suggested by @giordano on discourse already, maybe you want to have this package hosted under JuliaPhysics.

Although this package isn't strictly a physics package, Measurements.jl is part of JuliaPhysics so it might make sense to have them next to each other as they seem to serve the same purpose.

Feel free to close this issue if you prefer to leave it as a separate package.

Allow specification of random number generator everywhere

Right now there are a few rand calls embedded in functions where we can't pass a custom rng. This makes testing functionality that depends on MCM challenging because we can't set a seed.

Everywhere where rand is called, there should be a way to pass an rng::AbstractRNG parameter. Likwise, interface functions that might internally have a stochastic call should have a version that takes rng.

While this can be done with a keyword argument and slurping, I find it more useful to implement using the pattern

myfun(args...; kwargs...) = myfun(Random.GLOBAL_RNG, args...; kwargs...)
myfun(rng::AbstractRNG, args...; kwargs...) = ... # my functionality

How to make MonteCarloMeasurements work with Optim

using MonteCarloMeasurements,Optim
f(x)=(1.0 - x[1])^2 + 100.0 * (x[2] - x[1]^2)^2
x=[0.0±0.0,0.0±0.1]
optimize(f,x)

ERROR: MethodError: no method matching -(::Particles{Float64,2000}, ::Array{Float64,1})
Closest candidates are:
-(::Particles{T,N}) where {T, N} at C:\Users\Administrator.julia\packages\MonteCarloMeasurements\J2g8U\src\register_primitive.jl:95
-(::Real, ::Complex{Bool}) at complex.jl:304
-(::Real, ::Complex) at complex.jl:316
...
Stacktrace:
[1] (::Statistics.var"#8#9"{Array{Float64,1}})(::Particles{Float64,2000}) at D:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.4\Statistics\src\Statistics.jl:219
[2] _mapreduce(::Statistics.var"#8#9"{Array{Float64,1}}, ::typeof(+), ::IndexLinear, ::Array{Particles{Float64,2000},1}) at .\reduce.jl:400
[3] _mapreduce_dim(::Function, ::Function, ::NamedTuple{(),Tuple{}}, ::Array{Particles{Float64,2000},1}, ::Colon) at .\reducedim.jl:312
[4] #mapreduce#580 at .\reducedim.jl:307 [inlined]
[5] mapreduce(::Function, ::Function, ::Array{Particles{Float64,2000},1}) at .\reducedim.jl:307
[6] centralize_sumabs2(::Array{Particles{Float64,2000},1}, ::Array{Float64,1}) at D:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.4\Statistics\src\Statistics.jl:220
[7] _varm(::Array{Particles{Float64,2000},1}, ::Array{Float64,1}, ::Bool, ::Colon) at D:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.4\Statistics\src\Statistics.jl:307
[8] #varm#11 at D:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.4\Statistics\src\Statistics.jl:297 [inlined]
[9] _var at D:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.4\Statistics\src\Statistics.jl:340 [inlined]
[10] #var#15 at D:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.4\Statistics\src\Statistics.jl:335 [inlined]
[11] var at D:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.4\Statistics\src\Statistics.jl:335 [inlined]
[12] nmobjective(::Array{Particles{Float64,2000},1}, ::Int64, ::Int64) at C:\Users\Administrator.julia\packages\Optim\SFpsz\src\multivariate\solvers\zeroth_order\nelder_mead.jl:104
[13] initial_state(::NelderMead{Optim.AffineSimplexer,Optim.AdaptiveParameters}, ::Optim.Options{Float64,Nothing}, ::NonDifferentiable{Particles{Float64,2000},Array{Particles{Float64,2000},1}}, ::Array{Particles{Float64,2000},1}) at C:\Users\Administrator.julia\packages\Optim\SFpsz\src\multivariate\solvers\zeroth_order\nelder_mead.jl:171
[14] optimize(::NonDifferentiable{Particles{Float64,2000},Array{Particles{Float64,2000},1}}, ::Array{Particles{Float64,2000},1}, ::NelderMead{Optim.AffineSimplexer,Optim.AdaptiveParameters}, ::Optim.Options{Float64,Nothing}) at C:\Users\Administrator.julia\packages\Optim\SFpsz\src\multivariate\optimize\optimize.jl:33
[15] optimize(::Function, ::Array{Particles{Float64,2000},1}; inplace::Bool, autodiff::Symbol, kwargs::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}) at C:\Users\Administrator.julia\packages\Optim\SFpsz\src\multivariate\optimize\interface.jl:64
[16] optimize(::Function, ::Array{Particles{Float64,2000},1}) at C:\Users\Administrator.julia\packages\Optim\SFpsz\src\multivariate\optimize\interface.jl:58
[17] top-level scope at REPL[8]:1
[18] eval(::Module, ::Any) at .\boot.jl:331
[19] eval_user_input(::Any, ::REPL.REPLBackend) at D:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.4\REPL\src\REPL.jl:86
[20] run_backend(::REPL.REPLBackend) at C:\Users\Administrator.julia\packages\Revise\tV8FE\src\Revise.jl:1165
[21] top-level scope at none:0

Show of Vector{Particles}

The show method seems to to some rounding with particles in an array:

julia> [10000+i ± 0.5 for i=1:10:50]
5-element Array{Particles{Float64,2000},1}:
 10000.0 ± 0.5
 10000.0 ± 0.5
 10000.0 ± 0.5
 10000.0 ± 0.5
 10000.0 ± 0.5

julia> ans[1]
Particles{Float64,2000}
 10001.0 ± 0.5

Is this intentional? If so, then it should take into account the significant digits according to the uncertainty.

Using Particles{Bool,N} in a conditional

Thanks for the great package that Chad Scherrer introduced me to!

Currently,

julia> Bernoulli(Particles(10, Uniform()))
ERROR: TypeError: non-boolean (Particles{Bool,10}) used in boolean context

is an error due to the check_args call. You can get around this with Bernoulli(Particles(10, Uniform()), check_args=false), but it would be nice if that weren't needed. Would it be dangerous/too special-cased if a Particles{Bool} with all 1s behaved like a true in a conditional?

Particles as Distribution parameters

Really enjoying this package, thank you for your work on it.

For Soss.jl, here's a simple model:

m = @model begin
    x ~ Normal(0,1)
    y ~ Normal(x^2, 1) |> iid(10)
end

For simple examples like this, I can now use MonteCarloMeasurements to get

julia> particles(m)
(x = 0.00257 ± 1.0, y = Particles{Float64,1000}[1.0 ± 1.7, 1.0 ± 1.7, 1.0 ± 1.7, 0.998 ± 1.7, 1.0 ± 1.8, 1.0 ± 1.7, 0.998 ± 1.7, 0.997 ± 1.7, 0.999 ± 1.7, 0.996 ± 1.8])

The way I get there... it's not so pretty. Let's start with μ and σ:

julia> μ = Particles(1000,Normal(0,1))
Part1000(0.000235 ± 1.0)

julia> σ = Particles(1000,Normal(0,1))^2
Part1000(1.008 ± 1.5)

Surprisingly, this works just fine:

julia> Normal(μ,σ)
Normal{Particles{Float64,1000}}(
μ: 0.000235 ± 1.0
σ: 1.01 ± 1.5
)

Now trying the obvious thing,

julia> x = rand(Normal(μ,σ))
Part1000(-0.03248 ± 1.0)

It thinks it works, but it's really horribly broken:

julia> (x-μ)/σ
Part1000(-0.03246 ± 2.6e-8)

So I got it to work with a little helper function:

N = 1000
parts(x::Normal{P} where {P <: AbstractParticles}) = Particles(length(x.μ), Normal()) * x.σ + x.μ
parts(x::Sampleable{F,S}) where {F,S} = Particles(N,x)
parts(x::Integer) = parts(float(x))
parts(x::Real) = parts(repeat([x],N))
parts(x::AbstractArray) = Particles(x)
parts(p::Particles) = p 

So now I have

julia> x = parts(Normal(μ,σ))
Part1000(0.05171 ± 2.37)

julia> (x-μ)/σ
Part1000(0.0007729 ± 1.0)

Much better.

This is fine for one distribution, but most don't compose as nicely as a normal.

Many distributions break entirely, seemingly because of argument checking in Distributions.jl:

julia> Bernoulli(Particles(1000,Beta(2,3)))
ERROR: TypeError: non-boolean (Particles{Float64,1000}) used in boolean context
Stacktrace:
 [1] macro expansion at /home/chad/.julia/packages/Distributions/Iltex/src/utils.jl:5 [inlined]
 [2] Bernoulli{Particles{Float64,1000}}(::Particles{Float64,1000}) at /home/chad/.julia/packages/Distributions/Iltex/src/univariate/discrete/bernoulli.jl:31
 [3] Bernoulli(::Particles{Float64,1000}) at /home/chad/.julia/packages/Distributions/Iltex/src/univariate/discrete/bernoulli.jl:37
 [4] top-level scope at none:0

Do you see a path toward...

  1. Being able to build distributions with particle parameters in this way?
  2. Getting rand to work properly?

bymap does not support function with 2 output values

I have a function that takes two inputs, a numeric array and an array of Particles, that I call using bymap:

ret = @bymap decode(v,Y) X = ret[1] a = ret[2]
When I added the second (non-Particle) output value I got this error:

Output with dimension >2 is currently not supported by `bymap`.

Is there any workaround?

A few bugs with MonteCarloMeasurements + ControlSystems

Hey, sorry to dump multiple things on here, I can separate these into different issues if you'd like. I'm not even sure if this or ControlSystems.jl is the right place for this kind of thing. I'm trying to use MonteCarloMeasurements.jl + ControlSystems.jl to get something close to MATLAB's Robust Control Toolbox. Trying to replicate this example fails on a few parts.

  1. Here, trying to call bode on F, or G1 crashes my system sometimes.
using ControlSystems
using GenericLinearAlgebra
using MonteCarloMeasurements
using Plots

unsafe_comparisons(true, verbose=false)

s = tf("s")

C = 100 * ss((s + 1) / (0.001*s + 1))^3

k = 1.0 ± 0.2
m1 = 1.0 ± 0.2
m2 = 1.0 ± 0.2

G1 = 1 / s^2 / m1
G2 = 1 / s^2 / m2

F = [0; G1] * [1 -1] + [1 -1]' * [0 G2]

bode(G1)     # This crashes
bode(F)      # This crashes
bode(F[2,1]) # This also crashes

This gives me the error

Unreachable reached at 000000003CBF36DF

Please submit a bug report with steps to reproduce this fault, and any error messages that follow (in their entirety). Thanks.
Exception: EXCEPTION_ILLEGAL_INSTRUCTION at 0x3cbf36df -- getindex at .\array.jl:809
in expression starting at REPL[2]:1
getindex at .\array.jl:809
iterate at .\array.jl:785 [inlined]
iterate at .\array.jl:785 [inlined]
iterate at .\generator.jl:44 [inlined]
_collect at .\array.jl:699 [inlined]
collect_similar at .\array.jl:628 [inlined]
map at .\abstractarray.jl:2162 [inlined]
zpkdata at C:\Users\jdiegelm\.julia\packages\ControlSystems\H3Br2\src\analysis.jl:123
_bounds_and_features at C:\Users\jdiegelm\.julia\packages\ControlSystems\H3Br2\src\freqresp.jl:152
#146 at C:\Users\jdiegelm\.julia\packages\ControlSystems\H3Br2\src\freqresp.jl:138
iterate at .\generator.jl:47 [inlined]
_collect at .\array.jl:699 [inlined]
collect_similar at .\array.jl:628
unknown function (ip: 000000003CBF3D64)
map at .\abstractarray.jl:2162 [inlined]
_default_freq_vector at C:\Users\jdiegelm\.julia\packages\ControlSystems\H3Br2\src\freqresp.jl:138
_default_freq_vector at C:\Users\jdiegelm\.julia\packages\ControlSystems\H3Br2\src\freqresp.jl:145 [inlined]
bode at C:\Users\jdiegelm\.julia\packages\ControlSystems\H3Br2\src\freqresp.jl:108
unknown function (ip: 000000003CBF3640)
jl_apply at /cygdrive/d/buildbot/worker/package_win64/build/src\julia.h:1690 [inlined]
do_call at /cygdrive/d/buildbot/worker/package_win64/build/src\interpreter.c:117
eval_value at /cygdrive/d/buildbot/worker/package_win64/build/src\interpreter.c:206
eval_stmt_value at /cygdrive/d/buildbot/worker/package_win64/build/src\interpreter.c:157 [inlined]
eval_body at /cygdrive/d/buildbot/worker/package_win64/build/src\interpreter.c:548
jl_interpret_toplevel_thunk at /cygdrive/d/buildbot/worker/package_win64/build/src\interpreter.c:660
jl_toplevel_eval_flex at /cygdrive/d/buildbot/worker/package_win64/build/src\toplevel.c:840
jl_toplevel_eval_flex at /cygdrive/d/buildbot/worker/package_win64/build/src\toplevel.c:790
jl_toplevel_eval at /cygdrive/d/buildbot/worker/package_win64/build/src\toplevel.c:849 [inlined]
jl_toplevel_eval_in at /cygdrive/d/buildbot/worker/package_win64/build/src\toplevel.c:883
eval at .\boot.jl:331
eval_user_input at D:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.5\REPL\src\REPL.jl:134
repl_backend_loop at D:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.5\REPL\src\REPL.jl:195
start_repl_backend at D:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.5\REPL\src\REPL.jl:180
#run_repl#37 at D:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.5\REPL\src\REPL.jl:292
run_repl at D:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.5\REPL\src\REPL.jl:288
#807 at .\client.jl:399
jfptr_YY.807_38807 at C:\Users\jdiegelm\AppData\Local\Programs\Julia\Julia 1.5.2\lib\julia\sys.dll (unknown line)
jl_apply at /cygdrive/d/buildbot/worker/package_win64/build/src\julia.h:1690 [inlined]
do_apply at /cygdrive/d/buildbot/worker/package_win64/build/src\builtins.c:655
jl_f__apply at /cygdrive/d/buildbot/worker/package_win64/build/src\builtins.c:669 [inlined]
jl_f__apply_latest at /cygdrive/d/buildbot/worker/package_win64/build/src\builtins.c:705
#invokelatest#1 at .\essentials.jl:710 [inlined]
invokelatest at .\essentials.jl:709 [inlined]
run_main_repl at .\client.jl:383
exec_options at .\client.jl:313
_start at .\client.jl:506
jfptr__start_36391 at C:\Users\jdiegelm\AppData\Local\Programs\Julia\Julia 1.5.2\lib\julia\sys.dll (unknown line)
jl_apply at /cygdrive/d/buildbot/worker/package_win64/build/src\julia.h:1690 [inlined]
true_main at /cygdrive/d/buildbot/worker/package_win64/build/ui\repl.c:106
wmain at /cygdrive/d/buildbot/worker/package_win64/build/ui\repl.c:227
__tmainCRTStartup at /usr/src/debug/mingw64-x86_64-runtime-7.0.0-1/crt\crtexe.c:334
mainCRTStartup at /usr/src/debug/mingw64-x86_64-runtime-7.0.0-1/crt\crtexe.c:223
BaseThreadInitThunk at C:\Windows\System32\KERNEL32.DLL (unknown line)
RtlUserThreadStart at C:\Windows\SYSTEM32\ntdll.dll (unknown line)
Allocations: 114178914 (Pool: 114140923; Big: 37991); GC: 110

What's weird is if I put some other stuff before the bode calls, it works for G1.

using ControlSystems
using GenericLinearAlgebra
using MonteCarloMeasurements
using Plots

unsafe_comparisons(true, verbose=false)

s = tf("s")

k = 1.0 ± 0.2
m1 = 1.0 ± 0.2
m2 = 1.0 ± 0.2

G1 = 1 / (m1 * s^2)
G2 = 1 / (m2 * s^2)

F = [0; G1] * [1 -1] + [1 -1]' * [0 G2]

P = lft(F, ss(k))

C = 100 * ss((s + 1) / (0.001*s + 1))^3

L = P * C

T = feedback(L, 1)

bode(G1)
  1. Trying to convert T, P, or L to zpk form hits a stack overflow error. This one, I'm sure, is more on the ControlSystems side though. It looks like it's expecting the eltypes of the state space matrices to be AbstractFloats. I tried converting to tf first, but that gives me an error that the iteration limit for a schur decomposition is reached.
julia> tf(T)
ERROR: ArgumentError: iteration limit 330 reached
Stacktrace:
 [1] _schur!(::GenericLinearAlgebra.HessenbergFactorization{Particles{Float64,2000},Array{Particles{Float64,2000},2},GenericLinearAlgebra.Householder{Particles{Float64,2000},S} where S<:(StridedArray{T, 1} where 
T)}; tol::Float64, debug::Bool, shiftmethod::Symbol, maxiter::Int64, kwargs::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}) at C:\Users\jdiegelm\.julia\packages\GenericLinearAlgebra\JoKjw\src\eigenGeneral.jl:104
 [2] _schur! at C:\Users\jdiegelm\.julia\packages\GenericLinearAlgebra\JoKjw\src\eigenGeneral.jl:91 [inlined]
 [3] #_schur!#21 at C:\Users\jdiegelm\.julia\packages\GenericLinearAlgebra\JoKjw\src\eigenGeneral.jl:169 [inlined]
 [4] _schur! at C:\Users\jdiegelm\.julia\packages\GenericLinearAlgebra\JoKjw\src\eigenGeneral.jl:169 [inlined]
 [5] #_eigvals!#23 at C:\Users\jdiegelm\.julia\packages\GenericLinearAlgebra\JoKjw\src\eigenGeneral.jl:246 [inlined]
 [6] _eigvals! at C:\Users\jdiegelm\.julia\packages\GenericLinearAlgebra\JoKjw\src\eigenGeneral.jl:246 [inlined]
 [7] eigvals!(::Array{Particles{Float64,2000},2}; sortby::Nothing, kwargs::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}) at C:\Users\jdiegelm\.julia\packages\GenericLinearAlgebra\JoKjw\src\eigenGeneral.jl:260
 [8] #eigvals#73 at D:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.5\LinearAlgebra\src\eigen.jl:326 [inlined]
 [9] #eigvalsnosort#82 at C:\Users\jdiegelm\.julia\packages\ControlSystems\H3Br2\src\utilities.jl:35 [inlined]
 [10] eigvalsnosort at C:\Users\jdiegelm\.julia\packages\ControlSystems\H3Br2\src\utilities.jl:35 [inlined]
 [11] charpoly(::Array{Particles{Float64,2000},2}) at C:\Users\jdiegelm\.julia\packages\ControlSystems\H3Br2\src\types\conversion.jl:310
 [12] convert(::Type{TransferFunction{Continuous,ControlSystems.SisoRational{Particles{Float64,2000}}}}, ::StateSpace{Continuous,Particles{Float64,2000},Array{Particles{Float64,2000},2}}) at C:\Users\jdiegelm\.julia\packages\ControlSystems\H3Br2\src\types\conversion.jl:261
 [13] convert(::Type{TransferFunction{Continuous,ControlSystems.SisoRational}}, ::StateSpace{Continuous,Particles{Float64,2000},Array{Particles{Float64,2000},2}}) at C:\Users\jdiegelm\.julia\packages\ControlSystems\H3Br2\src\types\conversion.jl:270
 [14] convert at C:\Users\jdiegelm\.julia\packages\ControlSystems\H3Br2\src\types\conversion.jl:250 [inlined]
 [15] tf(::StateSpace{Continuous,Particles{Float64,2000},Array{Particles{Float64,2000},2}}) at C:\Users\jdiegelm\.julia\packages\ControlSystems\H3Br2\src\types\tf.jl:69
 [16] top-level scope at REPL[41]:1

What's weird here is that I can call eigvals(T.A) and it works fine.

0.3.5 release?

The register_primitive fix is a big help to me, so now I'm requiring >=0.3.5 of your package. But that seems to cause problems for at least some users. Any chance you could tag a v0.3.5 release?

Computing tail probabilities

I'm wondering if there is a way to compute tail probabilities from a Particles variable without converting them to a vector. Say that p(x) is the (empirical) probability density associated with a Particles P, I would like to compute for example \integral_{-\infty} ^b p(x) dx = Prob(P < b)
Maybe one could redefine

P < b 

to do that?

Error with Complex Particles

Hello everybody, it is still me... According to the MonteCarloMeasurements.jl paper it seems that Complex Particles are supported, however I am getting an error, see this MWE:

using MonteCarloMeasurements
using Distributions

Y = [Particles(100, Normal(10+3i,2)) for i=1:3]

a = sin(2π*Y[1])
b = cos(2π*Y[2])

@show typeof(a)
@show typeof(a)

z = a + im*b

gives

typeof(a) = Particles{Float64,100}
typeof(a) = Particles{Float64,100}
Comparison operators are not well defined for uncertain values and are currently turned off. Call `unsafe_comparisons(true)` to enable comparison operators for particles using the current reduction function mean. Change this function using `set_comparison_function(f)`.

Stacktrace:
 [1] error(::String) at .\error.jl:33

Are Complex Particles still not fully supported maybe?

register_primitive causes "eval into closed module" warnings

My library needs this line in order for things to work properly:

foreach(register_primitive, [<=, >=, <, >])

But this makes using Soss spit out lots of warnings like this:

WARNING: eval into closed module MonteCarloMeasurements:
Expr(:function, Expr(:call, Expr(:., Base, :(:>)), Expr(:::, :p, :StaticParticles)), Expr(:block, #= Symbol("/home/chad/.julia/packages/MonteCarloMeasurements/nd0Nh/src/register_primitive.jl"):37 =#, Expr(:call, :StaticParticles, Expr(:call, :map, Expr(:., Base, :(:>)), Expr(:., :p, :(:particles))))))
  ** incremental compilation may be fatally broken for this module **

Should I be calling register_primitive differently? Do you know a way around this?

Question about register_primitive()

I realized that I need to call

register_primitive(floor)

since floor() by default does not operate correctly on Particles. If I call it from the REPL then everything works fine. However when I call it from a function defined inside a module the code just hangs. I read in the help section:

  register_primitive(f, eval=eval)

  Register both single and multi-argument function so that it works with particles. If you want to
  register functions from within a module, you must pass the modules eval function.

however this is not 100% clear to me. An example would be great!

Wasserstein distance between particles

It would be neat to provide functionality for the Wasserstein-p distance between particles. Should be super simple for 1D, and doable for MvParticles of reasonable dimensions.

fft of Particles vector not working

I am getting an error when I try to compute the fft of a vector of Particles. Are there any workarounds? In my specific case this is a 3 point fft, so I could just use an explicit formula, however I'd like to do this in the more general case as well at some point.

set2expr only works for F->F functions

If the function returns an object of different type, set2expr will not work. This expression should be built based on the output type rather than the input type.

Strange behavior of approximate equality

MWE:

ulia> import MonteCarloMeasurements

julia> a = 12.345
12.345

julia> b = MonteCarloMeasurements.Particles([a])
Particles{Float64, 1}
 12.345

julia> a == b
true

julia> a  b
false

Rework Weighted Particles

I've been thinking about about the approach for WeightedParticles. When I've needed this I haven't been able to use the built-in type. As it stands, the weighting vector is per-Particles, but I've needed a shared weighting vector across several Particles.

Is the non-sharing an intentional design choice, or are you open to this sharing idea? I think it would look something like this:

struct LogWeights{T,N} <: AbstractParticles{T,N}
    particles::Vector{T}
    next 
end

struct WeightedParticles{T,N} <: AbstractParticles{T,N}
    particles::Vector{T}
    logweights::Ref{LogWeights{W,N} where W <: Real}
end

Functions like +, *, etc on WeightedParticles would require === on the logweights, to be sure that the entire system evolves as one. The next function is mostly for resampling. The idea is that when a resample is triggered, existing particles become "stale" and need to pass through this function to reindex the values.

I think there are a few ways of implementing this last bit, but I wanted to get your take on the shared weights idea first before digging into it much more.

StructArrays

Hi @baggepinnen , we've talked before about StructArrays...

I haven't used this package much before, but I did get this going and it seems to work well:

julia> using StructArrays, MonteCarloMeasurements, NamedTupleTools, BenchmarkTools

julia> function structarray(nt::NamedTuple)
           return StructArray(prototype(nt)(getproperty.(values(nt), :particles)))
       end
structarray (generic function with 1 method)

julia> nt = (x=Particles(5), y=Particles(5))
(x = -4.44e-17 ± 0.98, y = -4.44e-17 ± 0.98)

julia> @btime structarray($nt)
  2.063 ns (0 allocations: 0 bytes)
5-element StructArray(::Array{Float64,1}, ::Array{Float64,1}) with eltype NamedTuple{(:x, :y),Tuple{Float64,Float64}}:
 (x = -1.2815515655446004, y = -1.2815515655446004)
 (x = -0.5244005127080408, y = -0.5244005127080408)
 (x = 0.5244005127080407, y = 0.5244005127080407)
 (x = 1.2815515655446004, y = 0.0)
 (x = 0.0, y = 1.2815515655446004)

Not sure how this might fit in with your manipulations in the package, just passing it along in case it's helpful :)

Weights are not carried through map

If WeightedParticles are partaking in calculations, their weights are reset to 1/N. Break out function definition loops of type loop and define separate function definition loop for WeightedParticles.

register_primitive doesn't work for Uniful Quantity objects wrapping Particles

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

::Particles{Bool,N}

Can you help me understand what's needed to implement this? Here's the use case:

First, some imports, and some bivariate normal particles:

julia> using Distributions, MonteCarloMeasurements, LinearAlgebra, TransformVariables

julia> z = Particles(1000,MvNormal(zeros(2,2) + I))
2-element Array{Particles{Float64,1000},1}:
  0.052 ± 0.99
 -0.0192 ± 1.0

Now we define a transform. I'll show the typical use case as well:

julia> t = as((y=asℝ₊,p=as𝕀))
TransformVariables.TransformTuple{NamedTuple{(:y, :p),Tuple{TransformVariables.ShiftedExp{true,Float64},TransformVariables.ScaledShiftedLogistic{Float64}}}}((y = TransformVariables.ShiftedExp{true,Float64}(0.0), p = TransformVariables.ScaledShiftedLogistic{Float64}(1.0, 0.0)), 2)

julia> t([0.0,0.0])
(y = 1.0, p = 0.5)

julia> t([0.0,0.0]) |> inverse(t)
2-element Array{Float64,1}:
 0.0
 0.0

This works great, and it works fine on particles as well, going "forward". But the inverse transform breaks, because it expects a Boolean:

julia> foreach(register_primitive, [<=, >=, <, >])

julia> t(z)
(y = 1.71 ± 2.1, p = 0.496 ± 0.21)

julia> t(z) |> inverse(t)
ERROR: TypeError: non-boolean (Particles{Float64,1000}) used in boolean context
Stacktrace:
 [1] macro expansion at /home/chad/.julia/packages/ArgCheck/xX4DA/src/checks.jl:162 [inlined]
 [2] inverse at /home/chad/.julia/packages/TransformVariables/nMbbh/src/scalar.jl:75 [inlined]
 [3] inverse!(::SubArray{Float64,1,Array{Float64,1},Tuple{UnitRange{Int64}},true}, ::TransformVariables.ShiftedExp{true,Float64}, ::Particles{Float64,1000}) at /home/chad/.julia/packages/TransformVariables/nMbbh/src/scalar.jl:23
 [4] _inverse!_tuple(::Array{Float64,1}, ::Tuple{TransformVariables.ShiftedExp{true,Float64},TransformVariables.ScaledShiftedLogistic{Float64}}, ::Tuple{Particles{Float64,1000},Particles{Float64,1000}}) at /home/chad/.julia/packages/TransformVariables/nMbbh/src/aggregation.jl:196
 [5] inverse!(::Array{Float64,1}, ::TransformVariables.TransformTuple{NamedTuple{(:y, :p),Tuple{TransformVariables.ShiftedExp{true,Float64},TransformVariables.ScaledShiftedLogistic{Float64}}}}, ::NamedTuple{(:y, :p),Tuple{Particles{Float64,1000},Particles{Float64,1000}}}) at /home/chad/.julia/packages/TransformVariables/nMbbh/src/aggregation.jl:237
 [6] inverse(::TransformVariables.TransformTuple{NamedTuple{(:y, :p),Tuple{TransformVariables.ShiftedExp{true,Float64},TransformVariables.ScaledShiftedLogistic{Float64}}}}, ::NamedTuple{(:y, :p),Tuple{Particles{Float64,1000},Particles{Float64,1000}}}) at /home/chad/.julia/packages/TransformVariables/nMbbh/src/generic.jl:206
 [7] InverseTransform at /home/chad/.julia/packages/TransformVariables/nMbbh/src/generic.jl:109 [inlined]
 [8] |>(::NamedTuple{(:y, :p),Tuple{Particles{Float64,1000},Particles{Float64,1000}}}, ::TransformVariables.InverseTransform{TransformVariables.TransformTuple{NamedTuple{(:y, :p),Tuple{TransformVariables.ShiftedExp{true,Float64},TransformVariables.ScaledShiftedLogistic{Float64}}}}}) at ./operators.jl:813
 [9] top-level scope at none:0

Error with Arrayof Particles

I want to pass to a custom function an Array of Particles, so I redefined one argument from Array{Float64,1} to Array{Number,1}, however I am getting an error. With a non-Array argument everything works fine. Here is a MWE:

Y = [Particles(100, Normal(10+3i,2)) for i=1:3]

function g(x::Float64, Y::Array{Number,1})
    sum(Y) - x
end

g(1.0,Y)

MethodError: no method matching g(::Float64, ::Array{Particles{Float64,100},1})
Closest candidates are:
  g(::Float64, !Matched::Array{Number,1}) at In[21]:2

Stacktrace:
 [1] top-level scope at In[25]:1

What am I missing?

Why MonteCarlo and why Particles?

Hi,

first of all, great package!

I'm curious why you chose the name MonteCarloMeasurements.jl for this package. Could you elaborate the Monte Carlo context that you had in mind (apart from Monte Carlo generically representing randomness)?
Typically, say in a MCMC scenario, I do not know the distribution function of an observable in advance but would get individual, correlated "Particles" from a simulation one after another. AFAICS, I can't use this package in such a case as there

  • doesn't seem to be functionality to add particles one by one
  • and, more importantly, it doesn't respect correlations in the data, see for example here (it only considers correlations if I give the respective distribution function in the constructor in advance)

Is this assessment correct?

On a side note, I was confused by the name Particles at first (maybe that's because I'm a physicist). It seems to expose the technical detail that the uncertainty is represented by a "data cloud" to the user. Personally I find Measurements.jl's choice measurement a bit more natural. Note that the related python package seems to be called uncertainties, a general name as well.

EDIT:
FYI, my package BinningAnalysis.jl provides tools to estimate the standard error of the mean of correlated data.

ribbonplot goes beyond particle bounds

What makes MonteCarloMeasurements so apealling to me is the ability to approximate distributions that are very different from Gaussian. But some plotting methods still seem to approximate the approximation by putting things in terms of a Gaussian.

In particular, I'm thinking of

help?> ribbonplot
search: ribbonplot ribbonplot!

  ribbonplot(x,y,[q=2])

  Plots a vector of particles with a ribbon representing q*std(y). Default width is
  2σ

The 2σ approximation is especially problematic in cases with bounded support, where the plot often crosses into zero-probability land. The 2σ width is presumably to be a stand-in for a credible interval. But we already have something much more precise, so why not use it?

To me it would be much more natural to have two methods:

  • ribbonplot(x,y,[q::Real]) could give a ribbon spanning from the (1-q)/2 quantile to the (1+q)/2 quantile (so the coverage area is q)
  • ribbonplot(x,y,[q::(Real,Real)]) could give a ribbon spanning from the q[1] quantile to q[2].

Do you see any drawbacks to this approach for your use cases?

NaN Particles

This seems inconsistent:

julia> NaN ± NaN
ERROR: ArgumentError: Normal: the condition σ >= zero(σ) is not satisfied.

julia> NaN ± Inf
Particles{Float64,2000}
 NaN ± NaN

Multiplication with particles from Unitful.jl

Seeing that this worked:

julia> (1 ± 0.5)m * (1 ± 0)kg
Part10000(1.0 kg m ± 0.4999920218301007 kg m)

I was hoping that this would work, but it doesn't

julia> (1 ± 0.5)m * 1kg
ERROR: DimensionError: kg and 0.586026044547503 kg m are not dimensionally compatible.
Stacktrace:
 [1] convert(::Type{Quantity{Quantity{Float64,D,U} where U where D,�,Unitful.FreeUnits{(kg,),�,nothing}}}, ::Quantity{Float64,�*�,Unitful.FreeUnits{(kg, m),�*�,nothing}}) at C:\Users\jdiegelm\.julia\packages\Unitful\dZYmO\src\conversion.jl:108
 [2] setindex!(::Array{Quantity{Quantity{Float64,D,U} where U where D,�,Unitful.FreeUnits{(kg,),�,nothing}},1}, ::Quantity{Float64,�*�,Unitful.FreeUnits{(kg, m),�*�,nothing}}, ::Int64) at .\array.jl:825
 [3] copyto! at .\multidimensional.jl:962 [inlined]
 [4] Array{Quantity{Quantity{Float64,D,U} where U where D,�,Unitful.FreeUnits{(kg,),�,nothing}},1}(::Array{Quantity{Float64,�*�,Unitful.FreeUnits{(kg, m),�*�,nothing}},1}) at .\array.jl:541
 [5] convert at .\array.jl:533 [inlined]
 [6] Particles at C:\Users\jdiegelm\.julia\packages\MonteCarloMeasurements\FcMud\src\types.jl:18 [inlined]
 [7] *(::Particles{Quantity{Float64,�,Unitful.FreeUnits{(m,),�,nothing}},10000}, ::Quantity{Int64,�,Unitful.FreeUnits{(kg,),�,nothing}}) at C:\Users\jdiegelm\.julia\packages\MonteCarloMeasurements\FcMud\src\unitful.jl:28
 [8] top-level scope at REPL[23]:1

Not sure if this is a MonteCarloMeasurements.jl issue or a Unitful.jl issue, so I thought I'd put it here.

bymap on StaticParticles gives Particle

julia> using MonteCarloMeasurements

julia> p = StaticParticles(10)
SPart10(-4.441e-17 ± 0.989)

julia> h(x::Float64) = x + 1; p = StaticParticles(10)
SPart10(-8.882e-17 ± 0.989)

julia> @bymap h(p)
Part10(1.0 ± 0.989)

I would expect this to preserve in the input type.

Stop overloading `norm(::Particle)`?

The problem is that norm in Base works recursively, so I was quite surprised when norm(::Vector{Particle}, 3) returned just a float instead of a Particle.

Utilize Blas more

  • dot can be expressed like gemm or gemv for one or two vectors of particles.
  • In-place mul! and friends
  • axpy! and friends
  • more complicated expressions can utilize Einsum or Tullio etc.

Difficult functions

An idea using cassette or irtools: specify the function that is difficult, execute in a difficultContext and when the difficult function is called, change behavior of all particle methods to operate on a single particle only and repeat this for all particles.

Dynamo with context and dispatch seems to be able to do the job
https://mikeinnes.github.io/IRTools.jl/latest/dynamo/

Setfield mught be a useful component.
Metatheory.jl might also be useful.

Particles as indices

Hi @baggepinnen ,

For some Soss work, I have something like (simplified example)

julia> a = [Particles(), Particles()+1, Particles()+2]
3-element Array{Particles{Float64,2000},1}:
 -5.33e-18 ± 1.0
  1.0 ± 1.0
  2.0 ± 1.0

julia> i = Particles(DiscreteUniform(1,3))
Particles{Int64,2000}
 1.984 ± 0.821

and I'd like to be able to do

julia> a[i]
Particles{Float64,2000}
 0.995411 ± 1.29

My current approach is

function Base.getindex(a::AbstractArray{P}, i::Particles) where P <: AbstractParticles
    return Particles([a[i.particles[n]][n] for n in eachindex(i.particles)])
end

This works fine, but maybe I have 2, 3, or 4 indices instead of just one. And for n indices we might have some particles and some Ints, with 2^n possibilities.

Do you have a better way to set up this kind of dispatch? Ideally this would be in MCM, but I'm not sure how to even approach the general case

Recommend Projects

  • React photo React

    A declarative, efficient, and flexible JavaScript library for building user interfaces.

  • Vue.js photo Vue.js

    🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.

  • Typescript photo Typescript

    TypeScript is a superset of JavaScript that compiles to clean JavaScript output.

  • TensorFlow photo TensorFlow

    An Open Source Machine Learning Framework for Everyone

  • Django photo Django

    The Web framework for perfectionists with deadlines.

  • D3 photo D3

    Bring data to life with SVG, Canvas and HTML. 📊📈🎉

Recommend Topics

  • javascript

    JavaScript (JS) is a lightweight interpreted programming language with first-class functions.

  • web

    Some thing interesting about web. New door for the world.

  • server

    A server is a program made to process requests and deliver data to clients.

  • Machine learning

    Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.

  • Game

    Some thing interesting about game, make everyone happy.

Recommend Org

  • Facebook photo Facebook

    We are working to build community through open source technology. NB: members must have two-factor auth.

  • Microsoft photo Microsoft

    Open source projects and samples from Microsoft.

  • Google photo Google

    Google ❤️ Open Source for everyone.

  • D3 photo D3

    Data-Driven Documents codes.