{'fraction_non_alphanumeric': 0.03037738053510924, 'fraction_numerical': 0.011858862016590723, 'mean_word_length': 4.989853044086774, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Vectorization of random number generation and reproducibility of concurrent particle transport simulation

Efficient random number generation with high quality statistical properties and exact reproducibility of Monte Carlo simulations are important requirements in many areas of computational science. VecRNG is a package providing pseudo-random number generation (pRNG) in the context of a new library VecMath. This library bundles up several general-purpose mathematical utilities, data structures, and algorithms having both SIMD and SIMT (GPUs) support based on VecCore. Several state-of-the-art RNG algorithms are implemented as kernels supporting parallel generation of random numbers in scalar, vector, and Cuda workflows. In this report, we will present design considerations, implementation details, and computing performance of parallel pRNG engines on both CPU and GPU. Reproducibility of propagating multiple particles in parallel for HEP event simulation is demonstrated, using GeantV-based examples, for both sequential and fine-grain track-level concurrent simulation workflows. Strategies for efficient uses of vectorized pRNG and non-overlapping streams of random number sequences in concurrent computing environments is discussed as well.


Introduction
The stochastic nature of many physical systems is often well modeled by Monte Carlo techniques, which require pseudorandom number generators (pRNGs) with good statistical properties and generally the ability to create independent sub-sequences. Ever-increasing computing capacity is being provided by systems with modern multi-core CPU processors and especially by HPC systems equipped with GPUs or other accelerators. Owing to this, the precision and speed of a large scale scientific simulation of complex systems have been improved significantly. To enable these, new types of pRNGs with high statistical quality and a long period have been recently developed [1,2,3]. The recent evolution of hardware architectures towards wider vector pipelines, GPUs/accelerators and many-threads opens new opportunities for concurrent simulation models taking advantage of both SIMD and SIMT (GPU). These require pRNGs suitable for massively parallel and scalable computing [4].
Despite the fact that there are many parallel pRNGs developed for specific hardware architectures, portable libraries of random number services which can be used across different architectures and in hybrid computing models are not commonly available. It is pertinent to use common pRNGs not only for the purpose of validation, but also for the full reproducibility of ACAT 2019 Journal of Physics: Conference Series 1525 (2020) 012054 IOP Publishing doi: 10.1088/1742-6596/1525/1/012054 2 simulation chains across different platforms, discussed in detail in Section 4. In this paper, we report on design considerations and an initial implementation of a small set of pRNG algorithms portable for SIMD and SIMT which rely on a common kernel. The reproducibility of propagating multiple particles in parallel for high energy physics (HEP) event simulation is also demonstrated.

Vectorization of pseudo random number generation
VecRNG is a part of VecMath which is a collection of vectorized math algorithms and utility functions for HEP applications, based on the VecCore library [5]. It provides parallel pRNGs implementations for both SIMD and SIMT workflows via architecture-independent common kernels. In this section, we describe the design and implementation of vectorized algorithms in the context of VecRNG.

Choice of Generators
We sought generators that meet strict quality requirements, belonging to families of generators which have been examined in depth [6], or have evidence from ergodic theory of exceptional decorrelation properties [3]. All must pass major crush-resistant tests such as DIEHARD [7] and BigCrush of TestU01 [8]. In addition we took into consideration constraints in size of the state and performance: 1) a very long period (≥ 2 200 ), obtained from a small state (in memory), 2) fast implementations and repeatability in sequence on the same hardware configuration, 3) efficient ways of splitting the sequence into long disjoint streams.
For the first phase of implementation, we selected the following representative generators from major classes of pRNG: MRG32k3a [1], Random123 [2], and MIXMAX [3]. Table 1 shows key properties of the generators selected for vectorization.

Design Consideration
VecRNG is designed to support concurrent computing models as well as SIMD and SIMT (GPU) with common kernels using backends provided by VecCore. The primary purpose of the use of backends is to provide access to the full performance of vector hardware, while hiding both hardware-specific details (i.e. width of vector registers) and the complexity of low-level instruction sets behind abstraction layers. In this way, generic kernels are highly portable across different hardware architectures or dependencies to external libraries. Today a vector backend can use either the Vc library [9] or UME::SIMD [10] for explicit SIMD vectorization. The CUDA and scalar backends share standard types; in addition there are cuda-specific extensions. Detailed descriptions of VecCore and its backends can be found elsewhere [11,12]. Another design choice was the exclusive use of static polymorphism, motivated by performance considerations. Every concrete generator inherits through the CRTP (curiously recurring template pattern) from the VecRNG base class, which defines mandatory methods and common interfaces. Last but not least, the implementation of VecRNG is header only and provides a minimal set of member methods. This approach allows more flexibility to higher level interfaces for specific computing applications, but minimizes the overhead compilation time.

Implementation Details
VecRNG is the base class for the static polymorphism and has only one protected data member (State t *fState) where State t is defined in each concrete generator and provided through RNG traits. An example for MRG32k3a is: The essential components of VecRNG interfaces are Uniform<Backend>() and Uniform<Backend>(State t& s) which generate the backend type of double precision u.i.i.d (uniformly independent and identically distributed) in [0,1), and update the internal fState and the given state s, respectively. The latter shown below is originally introduced for GPU, but also is useful for reproducibility which will be discussed in the next section. Each derived pRNG class is only responsible for implementing the generic Kernel method and its own auxiliary member functions, One of associated requirements for each generator in VecRNG is to provide an efficient skipahead algorithm, s n+p = f p (s n ) (i.e., advancing a state, s n , by p-sequences where p is the unit of the stream length or an arbitrary number) in order to assign disjointed multiple streams for different tasks. For an example, the mandatory method, Initialize(long n), moves the random state at the beginning of the given n th stream. Each generator supports both scalar and vector backends with a common kernel. Random123 has an extremely efficient stream assignment without any additional cost since the key serves as the stream index while MRG32k3a uses transition matrices (A) which recursively evaluate (A s mod m) using the binary decomposition of s. The vector backend uses N(=SIMD length) consecutive substreams and also supports the scalar return-type which corresponds to the first lane of the vector return-type.
For SIMT (GPU) applications, multiple independent streams are assigned to threads of all blocks with Initialize(State t *states, unsigned int n), where n is the multiplication of the number of blocks times the number of threads per block and state is a pre-allocated array of State t of which the size is n. Then the state is used for each thread of a block to generate a random number and update its state using Uniform<Scalar>(state[tid]) where tid is the thread index in CUDA kernel (i.e, tid = threadIdx.x + blockIdx.x * blockDim.x). This approach can be also used for applications on the host side which require statistically disjointed streams for n-parallel tasks.
Besides the Uniform method, some commonly used random variates are also provided. Implementation details of probability distribution functions are discussed separately [13].

Performance
Three important metrics to evaluate performance of pRNGs are speed (CPU), memory requirement, and the quality of random numbers. As pRNGs selected for this report are already satisfactory for two latter criteria, we concentrate on the relative performance between different backends. To have a minimal validation check for quality, the fluctuation of the sum within each sample is required to meet the condition | N/nv u − N/2| < N/12 where u is the vector of random numbers and n v is the size of vector lanes for double (e.g., 4 for AVX). Table 2 shows preliminary performance results of implemented pRNGs on Intel R Xeon(R) E5-2620 CPU (12 cores @ 2.00GHz). Performance and quality of the Intel MKL/VSL [14] random number generation depend on the size of output array (N) -we used N = 32 for this comparison.  Table 3 were measured using the kernel configuration with 26 blocks x 192 threads on NVidia Tesla K20M GPU (2496 CUDA cores @ 0.71GHz). The Word size (W) and round (R) used for Random123 were W4x32 R20 for Threefry and W4x32 R10 for Philox, respectively.
Philox uses the Advance Randomization System (ARS) which iterates bijection with rounds of the Feistel function and a couple of XOR operations. Poor vector performance of Philox is due to a lack of a native vector operation for conversion between 64 and 32 bit (SIMD) integers which forces it to use scalar operations. Performance of Random123 (Threefry and Philox) on GPU will be further optimized.

Reproducibility of simulation in parallel computing
HEP experiments have a goal for the reproducibility for detector simulation: simulations with the same initial configuration (primary particles and pRNG choice and seed) must give the same results. This must hold true even if different choices are made during a run, e.g. using vector kernels for a set of physics processes of selected tracks. A key practical reason is the need to reproduce and debug problems which occur during the simulation of a particular event or initial particle. In addition, the reliability of a simulation which cannot be repeated is more difficult to assess.
In the Geant Vector Prototype (GeantV) [16], baskets of tracks undergoing the same interaction are accumulated to enable computations on vectors of track properties. The aim is to use vectors for the bulk of the computational work. The remaining tail of tracks is treated with serial (non-vectorized) code, using the same algorithms as the vector code. Multi-threading is used to gather larger populations of tracks having similar properties and enable (wider) vectors, targeting better use of vector code and higher performance. Due to the out-of-order execution in multi-threading, different tracks are collected in baskets in each run. In addition, a different set of remaining 'unbasketised' tracks is run in scalar code in each run, in particular during phases of basket starvation and the ramping down of the simulation. So a particular algorithm must obtain the same pRNG output value (variate) for a track, whether it is processed as part of a vector in a basket of tracks (in 'vector' mode) or as a single track using the non-vectorized code (in 'scalar' mode).
To obtain the same results for a track's physics interactions (or other operations), the same sequence of output values of a pRNG is needed. This is accomplished by associating a single state of a pRNG with each track. When a new track is created either as a primary particle or in an interaction, a deterministically-defined new state of the pRNG must be generated and associated to it. This idea, called 'pseudo-random' trees, was first proposed in the 1980s in a particle transport application [17]. A first implementation was also created using linear congruential generators. Applications in other parallel and branching computations have been proposed since -see the recent review of Schaathun [18] for an overview and an evaluation of methods proposed. One such method, the pedigree method, for constructing seeds was developed by Leiserson et al. [19] exploiting deterministic parallel computations written in Cilk. This method was demonstrated in particle transport simulation [20] using Geant4 [21] as a testbed.
We demonstrate how such splittable/tree pRNGs can be used in practice within the constraints of a particle transport program which mixes vector and 'scalar' (non-vector) code which can be run in either a multi-threaded or serial implementation. We seek also to measure the overhead, compared to simulations which do not use these methods, and, as a result, do not reproduce the same results between runs. We discuss a number of considerations and key aspects of the implementation.
The approach adopted for GeantV depends on two parts: First, an initial seed for the scalar mode or a set of seeds for the vector mode are assigned by the run number or other units of parallel tasks. Second, a unique stream index is determined for each track. The stream index for the primary track consists of high precision bits set by the event number and low precision bits by the track index. For the secondary (or daughter) track, the stream index is generated in a collision resistant way using the current state of the pRNG carried by the (mother) track that undergoes an interaction.
To enable vector code for physics processes, a vector pRNG must be created in order to generate the output in each vector lane of the PRNG corresponding to that track. In our design, this is the role of an instance of a proxy class, which acts as a vector pRNG. The proxy both provides all the expected outputs in each vector lane (as through from the pRNG of its track) and advances the state of each track's pRNG accordingly. Our first implementation of a proxy class gathers the contents of the scalar pRNGs into an instance of the corresponding VecRng class (e.g. gathering MRG32k3a<double> into MRG32k3a<Double v>). The proxy instance is reusable, by explicitly attaching and detaching the set of track pRNG states. We have tested this capability using a limited set of GeantV physics processes, including in particular bremsstrahlung, ionisation, and Compton scattering, which undergo a self-contained e − − γ cascade process. The pRNG used in the tests is ThreeFry from the Random123 package. This counter based generator was chosen because the stream is easily split into, the size of the state is moderate (128 bits), and the method of initialization from a seed is trivial.
We compare the number of tracks and steps of a simulation of 1000 events, each comprised of ten 10 GeV electrons impinging on a 50 layer lead and liquid argon calorimeter. The number of tracks and steps of scalar ('non-vector') and vector configurations, with either 1 or 4 threads  each has been compared.
The simulation is run in two modes: the default mode in which a per-thread state of one serial pRNG and one vector pRNG are used in each thread, and the 'reproducible' mode in which our method is used. Using the values for the 'reproducible' mode run with 1 thread as baseline (Seq-1T), the ratios of the number of tracks and steps can be seen in Figure 1.
It is verified that the reproducible mode maintains the constant number of tracks and steps for all tested configurations, as required. In addition, those numbers are different from the single threaded mode for each of the 1-thread (Seq-1T) and vector 1-thread (Vec-1T) modes by 0.2 − 0.6%. In the multi-threaded mode, the number of tracks and steps fluctuates, as expected with averages (and variances) within one sigma of the values of the reproducible mode.
Reproducibility introduces an overhead in simulation time due to copying and assigning pRNG states during simulation workflows, gathering scalar states to a SIMD vector state or joining-splitting states for the proxy approach, and synchronizing the index of states in output (Random123 specific). Figure 2 shows an example of the CPU overhead as the fraction of Time(Reproducibility)/Time(Default Mode) using gathering scalar PRNG states to a vector state and splitting states visa versa. Another approach using the join-split method shows a similar (2-5%) performance degradation for the reproducibility mode.
Alternative proxy implementations are under development, including one that avoids the cost of copying the data. This is of most interest for the cases in which the average number of variates required is small and/or the PRNG state is large.

Conclusion
We implemented common kernels of pseudorandom number generation for SIMD and SIMT architecture using VecCore. We also demonstrated reproducibility of propagating multiple particles in parallel for HEP event simulation with concurrent workflows.