Hello all,

A request for feedback on the design of this code before I put it through some serious testing and submit a formal pull request to Phobos' std.random.

The goal here is to provide an interface that is close to the Boost/C++11 design while allowing for some extra interesting possibilities.  In particular the goal is to allow the possibility of different underlying generation algorithms or "engines" that can serve different requirements regarding speed, precision, memory consumption etc.  See: http://www.cse.cuhk.edu.hk/%7Ephwl/mt/public/archives/papers/grng_acmcs07.pdf

... for further details.

The present code contains both a function and struct interface for normal random number generation; the latter stores the parameters of the distribution (mean and standard deviation) and an instance of the underlying engine, but not the random number generator.

(An alternative design choice for the struct would have been to design it as a range with front/popFront, but for this to work I'd have had to also have it store the RNG, and this would run into the same problems as RandomSample has due to RNGs not being reference types.  Boost/C++11 also do not store the RNG internally but implement a variate_generator class which allows the coupling of a random-number distribution and an underlying generator; this is the approach I'd take for defining say a RandomVariate range.)

I've so far implemented one underlying engine, the Box-Muller method currently
used in Boost.  This is not the best engine possible -- an implementation of the
Ziggurat algorithm is desirable -- but it represents what is typically used in
other libraries.  The essence of Box-Muller is that it generates
(uniformly-distributed) numbers a pair at a time and then uses those to generate
a corresponding pair of normally-distributed numbers.  The method is described here:
https://en.wikipedia.org/wiki/Box-Muller_transform

I've stuck very closely to the implementation in Boost, to the point of reproducing a "trick" Boost uses to get round the fact that its uniform random number generator only supports the interval [a, b) whereas the canonical description of Box-Muller requires the random numbers to be distributed in the interval (0, 1].  Obviously in D this is avoidable, but I felt there was some advantage in reproducing identical results to Boost, at least within the limits of floating-point numerical rounding (and maybe also because D's PI variable goes to a greater number of decimal places).

 From a licensing point of view this close reproduction isn't an issue, as Boost
and Phobos both use the same licence, but there'll be a credit to the original
Boost authors in any pull request.  However, if there's a desire to
differentiate the code more strongly, that can be arranged :-)  I've attached a
small C++ code example to demonstrate the common output for the D version and
that of Boost.

Anyway, what I'd really appreciate feedback on is the interface design for these functions and structs.  A few concrete questions:

     (1) Does the order of template arguments seem optimal?  My inclination is
perhaps to tweak it so that the Uniform RNG type goes first, the underlying
engine second, and the input/output type T goes last, as these are most likely
the order in which people will want to change things.

     (2) Is there a desire to have a separate input and output type?
Alternatively, should I perhaps treat all input and internal processing as using
the real type and just use T for the output?

     (3) Given the default template arguments provided, is there any way to
ensure that an instance of the Normal struct can be implemented like this:

      auto normal01 = Normal(0, 1);

rather than having to do this:

      auto normal01 = Normal!()(0, 1);

Yes, I could create a function that returns an instance, but I'm not sure what to call it given that I want the normal() function to serve a similar purpose as uniform() ... :-)

Thanks in advance for any and all feedback.

Best wishes,

      -- Joe

Joseph Rushton Wakeling:

Thank you for the work on this. A normals generator is quite needed in Phobos.

> I've so far implemented one underlying engine, the Box-Muller method currently
> used in Boost.  This is not the best engine possible -- an implementation of the
> Ziggurat algorithm is desirable -- but it represents what is typically used in
> other libraries.

Using the Boost API is useful. And using D-specific features/syntax is equally good.
Using the Ziggurat algorithm is desired, otherwise people will have to implement it outside Phobos, because it's better if you have to generate many normal distributed values (and sometimes I need many of them). is it possible to have both algorithms? Maybe with a template argument?

Bye,
bearophile

On 10/23/2012 03:50 PM, bearophile wrote:
> Using the Boost API is useful. And using D-specific features/syntax is equally
> good.

Most of the existing std.random is based around the planned C++11 random API, which in turn derives from boost.  There are several ways in which D seems able to do things more elegantly, however ... :-)

> Using the Ziggurat algorithm is desired, otherwise people will have to implement
> it outside Phobos, because it's better if you have to generate many normal
> distributed values (and sometimes I need many of them). is it possible to have
> both algorithms? Maybe with a template argument?

Yes, that's exactly the plan.  You see how there is an "engine" template argument (NormalRandomNumberEngine) which is currently by default Box-Muller, but could be Ziggurat or others (Marsaglia etc.)  Other engines will be added later, Box-Muller is there as the first and easiest to implement.

The idea of a NormalRandomNumberEngine is that it is a struct which returns via opCall a number drawn from the distribution N(0, 1).  This is trivial to transform into N(mean, sigma).

> Using the Ziggurat algorithm is desired, otherwise people will have to implement it outside Phobos, because it's better if you have to generate many normal distributed values (and sometimes I need many of them). is it possible to have both algorithms? Maybe with a template argument?

I have an implementation of the Ziggurat algorithm at https://github.com/jerro/phobos/blob/master/std/random.d#L2035. It modified the Ziggurat algorithm a bit, so that it doesn't need as many layers to work well, which reduces the memory consumption and makes initialization faster. The cauchy, normal and exponential functions currently use global tables, but the zigguratAlgorithm function can also be used to implement a struct or class based API.

Joseph mentions having different engines for generating normally distributed random numbers, so maybe I could write an engine for his API based on what I already have. If we had multiple engines, I don't think Ziggurat algorithm should be the default, though, because it requires relatively expensive initialization.

On 10/23/2012 04:36 PM, jerro wrote:
> I have an implementation of the Ziggurat algorithm at
> https://github.com/jerro/phobos/blob/master/std/random.d#L2035. It modified the
> Ziggurat algorithm a bit, so that it doesn't need as many layers to work well,
> which reduces the memory consumption and makes initialization faster. The
> cauchy, normal and exponential functions currently use global tables, but the
> zigguratAlgorithm function can also be used to implement a struct or class based
> API.

Oh, nice!  I will take a look at this in the next days.

> Joseph mentions having different engines for generating normally distributed
> random numbers, so maybe I could write an engine for his API based on what I
> already have. If we had multiple engines, I don't think Ziggurat algorithm
> should be the default, though, because it requires relatively expensive
> initialization.

The caveat here is that the Ziggurat algorithm is both fast and AFAIK optimal from the statistical point of view -- yes, it's costly space-wise, but it's also safe.  The Polar method might be a nice alternative that's faster than Box-Muller, that has good statistical properties, and that requires much less space.  (Box-Muller actually isn't really that good statistically.)

The choice probably needs to depend on the typical number of normal variates we think people are going to want to be generating.  The way my code is set up you should only need to allocate one engine per thread, so it seems like a once-off hit that's trivial in the context of any moderately sized simulation.

On 10/23/2012 04:36 PM, jerro wrote:
> I have an implementation of the Ziggurat algorithm at
> https://github.com/jerro/phobos/blob/master/std/random.d#L2035.

OK, so I had a reasonable (not totally in-depth!) look through.  Looks good! And of course reminds of the other benefit of Ziggurat, that it can be used for multiple different random number distributions.

It looks like it should be readily possible to integrate the core Ziggurat functionality and to convert the normal() function in your code to be a NormalRandomNumberEngine struct -- so assuming you like my own architecture proposal, let's do this (I'm happy to hear suggestions for alternative designs, as I don't feel particularly confident in my API-designing abilities).

For the other distributions, my feeling is that in some cases there's a value in also having this "engine" approach, e.g. for exponentially-distributed numbers one could use Ziggurat or one could use the approach

    T u = uniform!("[)", T, T, UniformRandomNumberGenerator)(0.0, 1.0, urng);
    return -log(1 - u)/lambda;

... which is not as fast but has a much lower memory footprint.

> It modified the Ziggurat algorithm a bit, so that it doesn't need as many layers
> to work well, which reduces the memory consumption and makes initialization faster.

Can you expand on this, and maybe provide a reference?  I don't doubt your code's effectiveness but I think where RNGs are concerned we really need to be able to justify our algorithmic choices.  There's too much literature out there showing how commonly-used algorithms actually carry statistical flaws.

Bigger picture on my approach to non-uniform random number distributions.  The goal is to have the following:

    * Where useful, it should be possible to define and use multiple different
      internal "engines" for generating random numbers from the given
      distribution

    * For each distribution, there should be a function interface and a struct
      interface.

    * The struct implementation should store the distribution parameters and an
      instance of the internal engine (if any).

    * The function implementation should have 2 versions: one which allows the
      user to pass an engine of choice as input, one which contains a static
      instance of the specified engine (hence, thread-safe, distinguished
      according to both engine type and underlying uniform RNG type).

    * ... unless there's no call for distinct underlying engines, in which case
      the function version just takes parameters and uniform RNG :-)

    * The struct version should be useful to couple with an RNG instance to
      create an arbitrary random-number range à la Boost's variate_generator
      class.

So, a nice way to expand on our respective approaches might be to incorporate your Ziggurat, adapt it to my Normal engine API, and for me to write a similar setup for exponentially-distributed random numbers that uses the simple approach above and can also use your Ziggurat implementation.

What do you think?

Best wishes,

    -- Joe

> Can you expand on this, and maybe provide a reference?  I don't doubt your code's effectiveness but I think where RNGs are concerned we really need to be able to justify our algorithmic choices.  There's too much literature out there showing how commonly-used algorithms actually carry statistical flaws.

The first change I made to the algorithm is shifting layer boundaries. Basic Ziggurat algorithm uses n layers with area A. Because computing samples from the bottom layer and the top layer is typically more expensive than computing samples from other layers, I shifted the layer boundaries so that the top and the bottom layers have areas A / 2 and other n - 1 layers have area A. That way we only need to draw half as many samples from the top and the bottom layer.

To better describe the second change, I drew a picture and uploaded it to http://i.imgur.com/bDRpP.png . The basic Ziggurat algorithm first chooses the layer randomly and then chooses a random number x between 0 and xavg. If x is below xmin, it returns it. The algorithm I use does this too. When x is above xmin, the basic Ziggurat algorithm chooses a random point in the outer layer area (see the picture), checks if y < f(y), returns its x coordinate if it is, and chooses a new point otherwise. My algorithm uses precomputed low and high x offsets for each layer. It chooses a point below the high diagonal, and checks if it is below the low diagonal. If it is (and it is in most cases), it returns its x coordinate without computing f(x). If only checks if y < f(x) when y is above the low diagonal. That way it doesn't need to make as many calls to f(x).

> It looks like it should be readily possible to integrate the core Ziggurat functionality and to convert the normal() function in your code to be a NormalRandomNumberEngine struct

I already have a NormalDist struct at https://github.com/jerro/phobos/blob/new-api/std/random.d#L2363 - converting that should be trivial.

>
> For the other distributions, my feeling is that in some cases there's a value in also having this "engine" approach, e.g. for exponentially-distributed numbers one could use Ziggurat or one could use the approach
>
>     T u = uniform!("[)", T, T, UniformRandomNumberGenerator)(0.0, 1.0, urng);
>     return -log(1 - u)/lambda;
>
> ... which is not as fast but has a much lower memory footprint.

I agree that the it's a good thing to design the API so that we can use different engines

> Can you expand on this, and maybe provide a reference?  I don't doubt your code's effectiveness but I think where RNGs are concerned we really need to be able to justify our algorithmic choices.  There's too much literature out there showing how commonly-used algorithms actually carry statistical flaws.

I have only tested a distribution of the output samples, which seems to be correct. I agree that we should try to avoid statistical flaws as much as possible. Do you know of any good tests for nonuniform random number generators?

> Bigger picture on my approach to non-uniform random number distributions.  The goal is to have the following:
>
>     * Where useful, it should be possible to define and use multiple different
>       internal "engines" for generating random numbers from the given
>       distribution
>
>     * For each distribution, there should be a function interface and a struct
>       interface.
>
>     * The struct implementation should store the distribution parameters and an
>       instance of the internal engine (if any).
>
>     * The function implementation should have 2 versions: one which allows the
>       user to pass an engine of choice as input, one which contains a static
>       instance of the specified engine (hence, thread-safe, distinguished
>       according to both engine type and underlying uniform RNG type).

>     * ... unless there's no call for distinct underlying engines, in which case
>       the function version just takes parameters and uniform RNG :-)
>
>     * The struct version should be useful to couple with an RNG instance to
>       create an arbitrary random-number range à la Boost's variate_generator
>       class.

Seems OK to me. Maybe the function that uses a static engine should have one version that takes Rng as a parameter and one that uses rndGen, similar to how uniform() works now?

How would the static engines be initialized? They could be initialized on first use, but this would slow sample generation a bit. Another options is to initialize them in the static constructor. I think it would be best to make the engine a static field of a struct template, and initialize it in the structs static constructor (similar to what my ZigguratTable struct does). That way you don't need to check if the engine is initialized every time you generate a sample and the engine isn't initialized unless the function template that uses it is instantiated.

> So, a nice way to expand on our respective approaches might be to incorporate your Ziggurat, adapt it to my Normal engine API, and for me to write a similar setup for exponentially-distributed random numbers that uses the simple approach above and can also use your Ziggurat implementation.
>
> What do you think?

I think we should do that. Do you already have some code somewhere where I can see it? I need to know what API the NormalRandomNumberEngine struct should have.

On 10/24/2012 11:23 PM, jerro wrote:
> I already have a NormalDist struct at
> https://github.com/jerro/phobos/blob/new-api/std/random.d#L2363 - converting
> that should be trivial.

I'll have a go at this some time in the next days. :-)

> I have only tested a distribution of the output samples, which seems to be
> correct. I agree that we should try to avoid statistical flaws as much as
> possible. Do you know of any good tests for nonuniform random number generators?

I'm not expert on this, but there are various tests described in this paper which compares different normal random number generation techniques:
http://www.cse.cuhk.edu.hk/%7Ephwl/mt/public/archives/papers/grng_acmcs07.pdf

... and the topic is also discussed in this paper:
http://www.doornik.com/research/ziggurat.pdf

I don't have sufficient knowhow to really feel confident about how to test these things though, and part of me feels it might be worth at some point approaching some stats experts and asking them to help define a test suite for D's random number functionality.

> Seems OK to me. Maybe the function that uses a static engine should have one
> version that takes Rng as a parameter and one that uses rndGen, similar to how
> uniform() works now?

Well, as I defined the function, the UniformRNG input has a default value of rndGen, so if you call it just as normal(mean, sigma); it should work as intended.  But if there's some factor which means it would be better to define a separate function which doesn't receive the RNG as input, I can do that.

> How would the static engines be initialized? They could be initialized on first
> use, but this would slow sample generation a bit. Another options is to
> initialize them in the static constructor. I think it would be best to make the
> engine a static field of a struct template, and initialize it in the structs
> static constructor (similar to what my ZigguratTable struct does). That way you
> don't need to check if the engine is initialized every time you generate a
> sample and the engine isn't initialized unless the function template that uses
> it is instantiated.

I think it probably depends on the details of the technique.  E.g. with Box-Muller you cache 2 random values, and you need to re-generate them for every pair of normal random variates the user requests, so there's no problem with having it initialized in opCall on the first call to the engine.

In general I favour that approach of initializing in opCall because it means that you don't have to ever pass anything to the constructor and also if I have a simulation which uses a lot of random numbers, its output won't be changed simply by creating an instance of a normal random number generator -- only when I start actually _using_ that generator in place of some other source of randomness.

> I think we should do that. Do you already have some code somewhere where I can
> see it? I need to know what API the NormalRandomNumberEngine struct should have.

Did you get the code I attached to my original post on this topic?  Is there something extra that you need?

I'm going to turn that into patches for Phobos which I'll put up on GitHub in the next days, so we can pull and push and test/write together as needed.

Best wishes,

    -- Joe

> I'm not expert on this, but there are various tests described in this paper which compares different normal random number generation techniques:
> http://www.cse.cuhk.edu.hk/%7Ephwl/mt/public/archives/papers/grng_acmcs07.pdf

I'll implement the tests described in this paper.

> ... and the topic is also discussed in this paper:
> http://www.doornik.com/research/ziggurat.pdf

> I don't have sufficient knowhow to really feel confident about how to test these things though, and part of me feels it might be worth at some point approaching some stats experts and asking them to help define a test suite for D's random number functionality.

Maybe there are testsuites already written for this.

> Well, as I defined the function, the UniformRNG input has a default value of rndGen, so if you call it just as normal(mean, sigma); it should work as intended.  But if there's some factor which means it would be better to define a separate function which doesn't receive the RNG as input, I can do that.

I don't see any downside to this.

> I think it probably depends on the details of the technique.  E.g. with Box-Muller you cache 2 random values, and you need to re-generate them for every pair of normal random variates the user requests, so there's no problem with having it initialized in opCall on the first call to the engine.
>
> In general I favour that approach of initializing in opCall because it means that you don't have to ever pass anything to the constructor and also if I have a simulation which uses a lot of random numbers, its output won't be changed simply by creating an instance of a normal random number generator -- only when I start actually _using_ that generator in place of some other source of randomness.

I have only been thinking about the Ziggurat algorithm, but you are right, it does depend on the details of the technique. For Box-Muller (and other engines that cache samples) it only makes sense to compute the first samples in the opCall. But for the Ziggurat algorithm, tables that must be computed before you can start sampling aren't changed during sampling and computing the tables doesn't require any additional arguments. So it makes the most sense for those tables to be initialized in the struct's constructor in the struct based API.

In my previous post I was talking about initializing static instances of the engine used in the normal() function. The advantage of initializing in a static constructor is that you don't need an additional check every time the normal() function is called. But because we will also have a struct based API, that will not require such checks (at least not for all engines), this isn't really that important. So we can also initialize the global engine instance in a call to normal(), if this simplifies things.

> Did you get the code I attached to my original post on this topic?  Is there something extra that you need?

I didn't notice the attachments, sorry (but I see them now). It seems it will be easy enough to write NormalZigguratEngine, there's just one thing that I think should be added to the Normal struct. Engines should have a way to do their initialization in Normal's constructor (where possible), so I think something like this:

static if(is(typeof(_engine.initialize())))
    _engine.initialize();

should be added to Normal's constructor. Or maybe

_engine = NormalRandomNumberEngine();

so that a static opCall can be used.

> I'm going to turn that into patches for Phobos which I'll put up on GitHub in the next days, so we can pull and push and test/write together as needed.

Maybe we should also have a separate file in a separate branch for tests. There will probably be a lot of code needed to test this well and the tests could take a long time to run, so I don't think they should go into unit test blocks (we can put some simpler tests in unit test blocks later). It may also be useful to use external libraries such as dstats for tests.