[ANN] MonteCarloMeasurements.jl

Thinking a bit more along those lines… This type could be used for black-box optimization with particle methods, e.g., particle-swarm optimization. The weighting would then be provided by the (inverse) cost function, and resampling would focus on promising areas of the parameter space. Maybe it would be a good idea to extract the particle type into a separate package to facilitate all these crazy ideas :stuck_out_tongue:

2 Likes

Not the most urgent thing to decide but Swarm.jl or ParticleSwarm.jl sounds cool.

ParticleSwarms.jl is clearer and reflects use in the literature (edited to pick up next comment). There are some repos named ParticleSwarmOptimization.jl

At the risk of exacerbating the way too early naming…if I understand the goal of the proposal correctly, I think ParticleSwarms.jl would align with the plural convention for things like Distributions.jl and StaticArrays.jl (if each subtype is a type of ParticleSwarm)

2 Likes

The package is now registrated and can be added with

] add MonteCarloMeasurements

(make sure to ] update first)

The readme is extended with further documentation, more plotting functions are implemented and I also added an example file where I visualize the Bode and Nyquist diagrams of an uncertain transfer-function model, as well as perform a time-domain simulation of the system and look at the uncertainties in the response.

I verified that the package works with DifferentialEquations.jl. The tutorial for solving differential equations using Measurement works for Particles as well. In the second pendulum example, I had to call solve(..., dt=0.01, adaptive=false) , otherwise the solver tried to estimate a suitable dt using the uncertain state, which created some problems.

On a side note, I implemented a very primitive CMA optimization algorithm as a proof of concept. It’s available here and seems to work okay. I do not intend to take this further though, as there is already the well developed BlackBoxOptim.jl

1 Like

Updates in v0.1.4

I implemented the weighted particles with resampling functionality (using numerically stable logsumexp).
I also added two macros, @bymap, @bypmap, detailed below.

Weighted particles

The type WeightedParticles contains an additional field logweights. You may modify this field as you see fit, e.g.

reweight(p,y) = (p.logweights .+= logpdf.(Normal(0,1), y .- p.particles))

where y would be some measurement. After this you can resample the particles using resample!(p). This performs a systematic resample with replacement, where each particle is sampled proportionally to exp.(logweights).

Monte-Carlo simulation by map/pmap

Some functions will not work when the input arguments are of type Particles. For this kind of function, we provide a fallback onto a traditional map(f,p.particles). The only thing you need to do is to decorate the function call with the macro @bymap like so:

f(x) = 3x^2
p = 1 ± 0.1
r = @bymap f(p)

We further provide the macro @bypmap which does exactly the same thing, but with a pmap (parallel map) instead, allowing you to run several invocations of f in a distributed fashion.

These macros will map the function f over each element of p::Particles{T,N}, such that f is only called with arguments of type T, e.g., Float64. This handles arguments that are multivaiate particles <: Vector{<:AbstractParticles} as well.

These macros will typically be slower than calling f(p), but allow for Monte-Carlo simulation of things like calls to Optim.optimize etc., which fail if called like optimize(f,p). If f is very expensive, @bypmap might prove prove faster than calling f with p, it’s worth a try. The usual caveats for distributed computing applies, all code must be loaded on all workers etc.

4 Likes

Add an ODE_DEFAULT_NORM overload in the init file of DiffEqBase so that the error norm returns a float if time is a float. That will make it work more nicely with adaptive time stepping. Otherwise time itself would need to be a swarm

Thanks for the tip! I created a WIP PR, WIP since it’s not quite working yet… :confused:

Some updates:

We now support

Further, the documentation in the readme is extended quite a lot

New examples

Robust probabilistic optimization

Here, we use MonteCarloMeasurements to perform robust optimization. With robust and probabilistic, we mean that we place some kind of bound on a quantile of an uncertain value, or otherwise make use of the probability distribution of some value that depend on the optimized parameters.

The application we consider is optimization of a PID controller. Normally, we are interested in controller performance and robustness against uncertainty. The robustness is often introduced by placing an upper bound on the, so called, sensitivity function. When the system to be controlled is parameterized by Particles, we can penalize both variance in the performance measure as well as the 90:th quantile of the maximum of the sensitivity function. This example illustrates how easy it is to incorporate probabilistic constrains or cost functions in an optimization problem using Particles.

Autodiff and Robust optimization

Another example using MonteCarloMeasurements to perform robust optimization, this time with automatic differentiation. We use Optim.jl to solve a linear program with probabilistic constraints using 4 different methods, two gradient free, one first-order and one second-order method.

Control systems

This example shows how to simulate control systems (using ControlSystems.jl
) with uncertain parameters. We calculate and display Bode diagrams, Nyquist diagrams and time-domain responses. We also illustrate how the package ControlSystemIdentification.jl interacts with MonteCarloMeasurements to facilitate the creation and analysis of uncertain systems.

We also perform some limited benchmarks.

Lantin Hypercube Sampling

We show how to initialize particles with LHS and how to make sure the sample gets the desired moments. We also visualize the statistics of the sample.

How MC uncertainty propagation works

We produce the first figure in this readme and explain in visual detail how different forms of uncertainty propagation propagates a probability distribution through a nonlinear function.

New benchmark

The benchmark results below come from examples/controlsystems.jl The benchmark consists of calculating the Bode curves for a linear system with uncertain parameters

w  = exp10.(LinRange(-1,1,200)) # Frequency vector
p  = 1 ± 0.1
ζ  = 0.3 ± 0.1
ω  = 1 ± 0.1
G  = tf([p*ω], [1, 2ζ*ω, ω^2])
t1 = @belapsed bode($G,$w)
   ⋮
Benchmark Results
Time with 500 particles 11.5937ms
Time with regular floating point 0.3456ms
Time with Measurements 0.5794ms
Time with 100 static particles 1.8212ms
Time with static sigmapoints 0.7252ms
500×floating point time 172.7780ms
Speedup particles vs. Manual 14.9x
Slowdown particles vs. Measurements 20.0x
Slowdown static vs. Measurements 3.1x
Slowdown sigma vs. Measurements 1.3x
15 Likes

Very cool package! Are there any plans to implement multi-threading as an alternative to pmap?

Thanks!

I have been considering it, but I think most functions would not benefit very much from naive threading. The speedup you get by using SIMD on almost every internal operation together with not redoing computation on deterministic input arguments is quite hard to beat. Possibly a blocking strategy where the computations are divided into exactly nthreads blocks could make use of both threads and SIMD. I’ll think about it :slight_smile:

Instead of just doing multithreading, can we try and get this finished up for automated GPU goodies? :slight_smile: https://github.com/baggepinnen/MonteCarloMeasurements.jl/pull/6 . Is there another place in the code which is calling maps?

Yeah, I can’t really help out with that, can never become friends with the gpu drivers… :confused:
I would be nice to have though.

map is used in the macro @bymap, which acts as a fallback for a top-level function f that fails when called with particles f(p1, p2, ...).

@bymap f(p1, p2, ...) just transforms the call to a map over f and assembles the result into a Particles in the end.

I don’t think @bymap/@bypmap needs GPU support though.

New version of Julia, new version of MonteCarloMeasurements :slight_smile:

Since I last updated this post, there has been quite a lot of usability changes, some highlights

  • Documentation has been improved and moved from the README to a proper website. https://baggepinnen.github.io/MonteCarloMeasurements.jl/dev/
  • We have a number of new contributers, thanks to all of you :pray:
  • Plotting is now more robust
  • Discrete distributions (i.e., distributions where a sample is an Int or a Bool etc.) work much better, thanks to @simeonschaub
  • Particles and missing should now play nicely together, thanks @simeonschaub.
  • Performance has increased, in the updated benchmark below by about 30% (not counting improvements performance improvements to julia), details here.
  • Conversion from particles and arrays of particles to regular arrays is much smoother, thanks to @sethaxen
  • MonteCarloMeasurements have been given support in the visualization library ArviZ.jl, thanks to @sethaxen
  • @cscherrer has written a nice blog post on how to make use of Particles for variational inference. MCM has also been incorporated into Soss.jl
  • The slack channel #monte-carlo-measurements have been created for discussions related to the package.
  • Both Zygote and ForwardDiff can differentiate stuff that contains Particles, robust optimization example.
  • Better handling of NamedTuple with particle fields.

Benchmark

The new results below are obtained with a faster computer and newer version of Julia. The relative timings in the bottom are however indicative of some nice performance improvements:

Benchmark Old Results New Results
Time with 500 particles 11.5937ms 1.3632ms
Time with regular floating point 0.3456ms 0.0821ms
Time with Measurements.jl 0.5794ms 0.1132ms
Time with 100 static particles 1.8212ms 0.2375ms
Time with static sigmapoints 0.7252ms 0.0991ms
500×floating point time 172.7780ms 41.0530ms
Speedup particles vs. Manual 14.9x 30.1x
Slowdown particles vs. Measurements.jl 20.0x 12.0x
Slowdown static vs. Measurements.jl 3.1x 2.1x
Slowdown sigma vs. Measurements.jl 1.3x 0.9x
13 Likes

I’d like to add support for something like this to BAT.jl, and I was wondering: what’s the relationship (if any) between MonteCarloMeasurements.jl and @carstenbauer’s MonteCarloObservable.jl](https://github.com/crstnbr/MonteCarloObservable.jl)? If they’re similar approaches to the problem, I may try to add support for both.

Based on the Readme of MonteCarloObservable it seems we are solving different problems.

I’m not sure how you envision support from MCM, it seems BAT is somewhat similar to Soss with some functionality overlapping with Arviz as well, so the support present in those packages may perhaos serve as inspiration for your package?

Ah, yes - I took a closer look at MonteCarloObservable, I think it goes in a different direction.

BAT is somewhat similar to Soss with some functionality overlapping with Arviz as well, so the support present in those packages may perhaos serve as inspiration for your package

Thanks, I’ll definitely have a look! BAT different from Soss in that it’s about use cases where the user has a custom likelihood distribution, it doesn’t provide a probabilistic programming DSL, but it also doesn’t require a forward model. MonteCarloMeasurements should fit in very well, though.

I also want to see if I can mesh MonteCarloMeasurements with ValueShapes.jl.

1 Like

I’m looking forward to see what you have in mind :slightly_smiling_face:

Ah, what I mean is, I want to make BAT support MCM - as in, being able to convert the BAT MCMC output to MCM Particles.