how to handle efficiently different probability distribution functions?

Question

I have a particle system that as usually creates new particles, updates them and destroys...

In the emitter module there is for loop that resets particles:

foreach p in particles
    p.position = rand()
    p.velocity = rand()

usually when using C's rand() function we get uniform distibution, but what when I would like to use some other distribution (gaussian for instance) ?

how to change that code so that it will handle several (or at least two) different ways of generating new particles' parameters?

Of course you can create some object: like RandomGenerator, and use some virtual function calls and handle those different behaviours. But this piece of code should be very fast (when updating thousands of particles), so using virtual functions is not good I think.

or maybe I should not care and simply write:

foreach p in particles
    p.position = useGaussian ? gausRand() : UniRand()
    p.velocity = useGaussian ? gausRand() : UniRand()

we can narrow number of different distribution and use only two or three of them...

please notice that my example is very simple, but in real code you have several particle parametrs configurations.

I would like to get some general advice on that subject.

Have you profiled this code to determine that the virtual call overhead is a problem? — Oliver Charlesworth, Jun 02 '12 at 10:22
I did not profile, but I wanted to get some general thoughts on that case. When there are thousands of particles the virtual method call does matter I think. — fen, Jun 02 '12 at 10:25
The general thoughts are that you should not optimise until profiling has determined something to be a problem. — Oliver Charlesworth, Jun 02 '12 at 10:31
@fen: Not really. Only if you redistribute them 3000 times/sec, say. — Puppy, Jun 02 '12 at 14:46

score 3 · Answer 1 · answered Jun 02 '12 at 11:06

Usually when using C's rand() function we get uniform distibution, but what when I would like to use some other distribution (gaussian for instance) ?

The Box Muller Transform is a very clever algorithm which uses trig functions to generate random numbers from a Gaussian distribution using a uniform distribution as input (ie. using rand()). You specify the mean and standard deviation, and invoke this function to generate new variates. The only drawback is that it is more expensive than simply calling rand, as it also calls sin() and cos() (though only every second invocation).

how to change that code so that it will handle several (or at least two) different ways of generating new particles' parameters?

I suggest you start with the RandomGenerator and virtual methods approach. It will be the easiest to maintain. Start with the simplest approach and profile before you try to optimise things.

Given the computational complexity of generating random numbers, the cost of producing a variate will far outweigh the overhead of a virtual method call versus a static function call.

If this really, really isn't fast enough you could always generate a pool of random numbers, generating more in the background as required.

score 3 · Accepted Answer · edited May 23 '17 at 12:04

While @gavinb's answer is a perfectly valid way to go, I'd suggest to avoid reinventing the wheel and go with the standard facilities: if you have c++11 support, use std::normal_distribution and its relatives (see, e.g. C++ TR1: how to use the normal_distribution?). Otherwise, use boost::random library.

Since these are header-only (the boost version at least), there are no polymorphic calls involved, so you don't have to worry about them. Which, of course, does not obviate the utmost relevance of @Oli Charlesworth's advice.

EDIT: if the overhead due to polymorphic calls is non-negligible, you can always template your functions on an enumerated type of distributions and specialize them as necessary.

In a nutshell, it's as simple as this:

#include<iostream>

// template on an int selector
template<int N> void foo(){ std::cout<<"42\n"; }
template<> void foo<1>() {std::cout<<"1\n";}

//now use an enum
enum  distr_types {UNIF, NORMAL, UNKNOWN}; 
template<distr_types T> void bar() {std::cout<<"fourty two\n";}
template<> void bar<UNIF>() {std::cout<<"UNIF\n";}
template<> void bar<NORMAL>(){std::cout<<"NORMAL\n";}

int main(){
  foo<3>();
  foo<1>();

  bar<UNIF>();
  bar<NORMAL>();
  bar<UNKNOWN>();
}

But if you find yourself doing things of this sort, it's worth having a look at one of good C++ books.

Very good point - rolling your own as I described is best done only when you can't add additional dependencies. Using `std::normal_distribution` or `boost::random` should probably be the first alternatives explored. +1 :) — gavinb, Jun 02 '12 at 12:38
both your answers are good... I simply need to write the code, profile it and then optimize (if needed). Probably I will have some performance issues only when using hundreds of thousands of particles. And maybe the easiest method will be good enough for me. — fen, Jun 02 '12 at 13:06
@Zhenya: can you provide me some more info (links) for that template functions? — fen, Jun 02 '12 at 14:04
@fen: you're welcome. By the way, you might consider editing the title of your post: is it really about some sort of particles? — ev-br, Jun 02 '12 at 14:56

how to handle efficiently different probability distribution functions?

2 Answers2