Modernize random number generator setup

[hakonsbm]

Overview

This PR introduces the new random number generator setup based on C++11 RNG facilities. For background, please see the discussion in #1440.

Changes for users

• Changing the kernel random number generator type is very easy.
• Seeding the random number generators is very easy: a single 32-bit seed provides the basis for seeding all generators in the kernel.
• Random number sequences differ from previous NEST versions; it is not possible to obtain exactly the same random number sequences as in NEST 2.
• Random numbers generated may differ depending on which compiler (ie, which implementation of the C++ standard library) is used.
• double random numbers will have full 53-bit mantissa resolution, unsigned long random numbers full 64 bit range.
• User modules can add additional random number generator engines (must conform to C++11 interface).
• Random number generators are no longer available at the SLI level.

Example

In [12]: nest.GetKernelStatus('rng_types')                                                                                                                                              
Out[12]: ('mt19937', 'mt19937_64')

In [12]: nest.GetKernelStatus(['rng_type', 'rng_seed'])                                                                                                                                 
Out[12]: ('mt19937_64', 143202461)

In [13]: nest.SetKernelStatus({'rng_type': 'mt19937', 'rng_seed': 4453453})                                                                                                             

In [14]: nest.GetKernelStatus(['rng_type', 'rng_seed'])                                                                                                                                 
Out[14]: ('mt19937', 4453453)

Changes for developers

• Random number generators no longer depend on GSL
• Models using random numbers through distributions now use wrappers for C++11(-compatible) distribution classes. Wrappers are required to allow exchanging the underlying RNG engine at runtime.
• Three types of RNGs are now available from the KernelManager:
  • rank_synced: same as the global_rng in NEST 2
  • vp_specific: same as the thread_rng in NEST 2
  • vp_synced: one RNG instance per VP/thread, but all need to be used in sync corresponding to the rank-synced generators; this is useful eg for the parent process in mip_generator.

Issues fixed

Fixes #1440; and thus fixes #245, fixes #349, fixes #508, fixes #1296, fixes #1381, fixes #1405, fixes #1574.

Things left to do

• Revise documentation
• Code formatting
• Add Random123 generators
• Perform some benchmarks
• Membrane potential randomisation in SLI versions of hpc_benchmark.
• Detailed NEST 2-to-3 transition guide.
• Fix failing tests (mostly due to hard-coded expected spike times)
  • SLI tests
    • unittests/test_mip_corrdet.sli
    • unittests/test_ppd_sup_generator.sli
    • unittests/test_spike_dilutor.sli
    • regressiontests/ticket-156.sli
  • PyNEST tests
    • test_connect_all_patterns.py
      • testWithMPI
    • test_connect_fixed_outdegree.py
      • testStatistics
    • test_erfc_neuron.py
      • test_activation_function
    • test_glif_cond.py
      • test_lif_asc
      • test_lif_r_asc
      • test_lif_r_asc_a
    • test_random_parameter.py
      • test_normal

[jessica-mitchell]

Thanks, looks good on my end!

[stinebuu]

Looks mostly good to me now, there is just a missing link that needs fixing.

[stinebuu]

I am happy with this now, thanks!

[jougs]

As we meanwhile have three approvals, I assume we could technically go ahead and merge this.

@heplesser, @jasperalbers, @jhnnsnk, @pnbabu: speak now or forever hold your peace ;-)

[heplesser]

@jougs I would like to include the validations by @jhnnsnk and @jasperalbers properly here before merging provided all is well.

@jhnnsnk and @jasperalbers, could you do the following:

1. Run your tests on mesocircuit/MAM for 30 different seeds to get enough data for statistical tests.
2. Add a brief description of how you ran the tests here.
3. Attach the actual data here (or link to a publicly accessible resource with the data), in the format @jhnnsnk sent me earlier offline.
4. Run statistical tests for equality of distribution of the values observed. Relevant tests for non-parametric two-sample comparison of distributions available in SciPy are
   • Kolmogorov-Smirnov
   • Epps-Singleton
   • Anderson-Darling
   • Mann-Whitney, which is not a distribution tests but tests if data in one sample is significantly larger/smaller than in the other sample.

I can also take care of running the tests if you provide the data.

[pnbabu]

I have two minor comments. Otherwise, it looks good to me.

[jasperalbers]

@heplesser see the data and analysis of the MAM simulations here:
https://github.com/jasperalbers/rng-scientific_validation

The plotting class also allows you to look at individual distributions, but for this somehow the bins need to be set cleverly.
I am happy to rectify any possible mistakes I made, conduct further tests or clarify anything that is not clear!

[heplesser]

@jasperalbers Thank you very much! I forked your repo to https://github.com/heplesser/rng-scientific_validation and extended your script (in a rather messy way). It now plots the p-values and collects all figure for one statistical test in a single figure. I ran the Kolmogorov-Smirnov and the Epps-Singleton tests.

For each of the three quantities we look at (rates, cv, cc), we have 4 x 2 x 32 tests, i.e., 256 tests. If all is well, we would thus expect about 2.5 tests to show p-values below 0.01 and 1.25 tests to show p-values below 0.005 (and correspondingly for larger p-levels). For the KS-tests this seems to be the case, for Epps-Singleton it looks like we have some more small p-values than expected. On the other hand, the expected number of p-value depends on the assumption that all values are independent, which is not really the case, since rates will be somehow correlated between populations in a simulation.

A further test would be to run tests with 2.20.1 and master-rng using different random number generators to see what differences that makes.

Kolmogorov-Smirnov

Epps-Singleton

[heplesser]

@jasperalbers One more follow-up on my posting from last night: As a second-level test, one could collect all 256 p-values for rates (correspondingly for CV, CC) and look at their distribution. Under the null hypothesis (statistically identical results for NEST 2 and NEST 3 RNG schemes), the p-values should be uniformly distributed on [0, 1]. Knuth TAOCP vol 2 has something on this.

[jhnnsnk]

@heplesser @jasperalbers I have now created a PR that analyzes the 4x4 data (30 runs, earlier master vs. new RNG branch) in the same way as done above for the MAM. Plots are uploaded there. Note that there is only one "area" (column).

[jhnnsnk]

Based on the validation results with the 4x4 mm² model, I approve this PR.

[jhnnsnk]

Although the statistics looks good for all runs of the 4x4 mm² model, I noticed outliers in the measured wall-clock times.
The figure shows the times for the statistically analyzed master and RNG-branch simulations and in addition repeat simulations with the RNG-branch.
Presimulation was 0.5 s and simulation 10 s of biological time.
Simulations were run on JURECA DC using 4 compute nodes with 8 MPI processes per node and 16 threads per MPI process.
The spikes of all neurons were recorded in ASCII format. For the statistical analysis, only the data from the simulation (not the presimulation) were included.

Outliers may be to some extend due to the machine workload. It is not understood why the connect times of the repeated runs with the RNG-branch are rather on the level of the master branch.
The repeatedly shorter simulation times for some run ids (same seed) may require investigation beyond this PR.

[heplesser]

I have extended the statistical tests (see https://github.com/heplesser/rng-scientific_validation) to better be able to judge observed p-values for the KS and ES tests. This is in particular important since the p-values we collect from different populations within a simulation are not statistically independent from each other.

To get a reference, I split the the trials we had for Old and New RNGs, respectively, into two non-overlapping groups. I then compared the two groups each to each other, separately for Old and New RNGs ("Old vs Old" and "New vs New" in figures below). This provides a reference for the distribution of p-values that might be seen in situations where differences between simulations results do not differ because of differences in RNGs used but only due to different random seeds. I then compared simulations with Old and New RNGs, for all four possible combinations of the non-overlapping trial groups ("Old vs New" in figures below).

If all p-values collected had been from statistically independent data, and the null-hypothesis of each pair of samples compared were distributed according to the same underlying distribution, then p-values should be uniformly distributed on [0, 1]. Due to the correlations in our data, deviations are to be expected. The reference curves allow us to judge them. The main concern are too many small p-values (cumulative distribution functions for p-values above the diagonal), as these would indicate many cases for which the null-hypothesis of equal underlying distributions does not hold at reasonable significance levels such as 1% or 5%. CDFs are very step wise for KS and for ES for the 4x4 model due to limited number of samples, which affects KS more than ES.

Figures show that Old-vs-New CDFs (brown, four per plot) generally lie not systematically above the green (Old vs Old) or blue (New vs New) reference lines. Based on this observation, I conclude that there is no indication of problems introduced by the new RNGs compared to the old, and that we can now merge the new RNGs.

Thanks to @jasperalbers (MAM) and @jhnnsnk (4x4) for providing the test data!

[jasperalbers]

Based on the scientific validation of the MAM, I approve of this PR. Thanks @heplesser for conducting the extensive statistical tests!

[heplesser]

Just putting in a blocker so we can merge this in ceremonial style :).

[heplesser]

All ready for the merge party now!