Random Generator and AT elements

[oscarxblanco]

Dear all,
I have several questions about random generators in AT, you might have already discussed.
I see that AT uses a C implementation of the random generator https://www.pcg-random.org/using-pcg-c.html saved in the parameters struct. Two streams are created common_rng and thread_rng.

at/atintegrators/attypes.h

Lines 24 to 25 in 8620504

	struct pcg_state_setseq_64 *common_rng;
	struct pcg_state_setseq_64 *thread_rng;

I will separate the questions in three parts : PYAT, AT MATLAB, and Elements using random numbers.
Please, let me know what I got wrong about how it works.

---

PYAT
Both streams are initialized in pyat using the C function pcg32_srandom_r, where AT_RNG_ING, seed and rank are set before tracking using the static constants and the MPI settings. So, there is no way the user could set the two values used for initialization of each stream, right ?

at/pyat/at.c

Lines 799 to 800 in 8620504

	pcg32_srandom_r(&common_state, seed, AT_RNG_INC);
	pcg32_srandom_r(&thread_state, seed, rank);

at/atintegrators/atrandom.c

Line 28 in 8620504

#define AT_RNG_INC 0xda3e39cb94b95bdbULL

at/pyat/at/tracking/track.py

Line 43 in 8620504

seed = random.common.integers(0, high=_imax, dtype=int)

at/pyat/at.c

Lines 787 to 793 in 8620504

	#ifdef MPI
	int irank;
	MPI_Comm_rank(MPI_COMM_WORLD, &irank);
	uint64_t rank = irank;
	#else
	uint64_t rank = 0;
	#endif /* MPI */

---

AT MATLAB
In AT MATLAB it seems that they are only initialized once in the MATLAB session (not on every track as in pyat) and it uses static values, right ?

at/atmat/attrack/atpass.c

Lines 73 to 74 in 8620504

	static pcg32_random_t common_state = COMMON_PCG32_INITIALIZER;
	static pcg32_random_t thread_state = THREAD_PCG32_INITIALIZER;

at/atintegrators/atrandom.c

Lines 30 to 31 in 8620504

	#define COMMON_PCG32_INITIALIZER { AT_RNG_STATE, AT_RNG_INC, 0.0, false }
	#define THREAD_PCG32_INITIALIZER { AT_RNG_STATE, 1ULL, 0.0, false }

---

Elements using random numbers
There are a few elements that use random numbers: atQuantDiff and atVariableMultipole.
Both state that a seed could be used. However, the QuantDiffPass method uses the Param->thread_rng stream and the VariableThinMultipole uses a static initializer PCG32_INITIALIZER inside atrandn, but, we just saw that the user has no control over the seed, so, the user seed is actually never used, right ?

(QuantDiff)

at/atmat/lattice/element_creation/atQuantDiff.m

Lines 17 to 19 in 8620504

	% The optional field Seed can be added. In that case, the seed of the
	% random number generator is set at the first turn.
	% ELEM=ATQUANTDIFF(FAMNANE,RING,'Seed',4)

at/atintegrators/QuantDiffPass.c

Line 59 in 8620504

QuantDiffPass(r_in, Elem->Lmatp, nturn, Param->thread_rng, num_particles);

(VariableThinMultipole)

at/pyat/at/lattice/variable_elements.py

Lines 49 to 50 in 8620504

	Seed(int): Seed of the random number generator for white
	noise excitation. Default datetime.now()

at/atmat/lattice/element_creation/atvariablemultipole.m

Line 22 in 8620504

% SEED Input seed for the random number generator

at/atintegrators/VariableThinMPolePass.c

Line 69 in 8620504

val = atrandn(0.0, 1.0);

at/atintegrators/atrandom.c

Lines 153 to 156 in 8620504

	static inline double atrandn(double mean, double stdDev)
	{
	return atrandn_r(&pcg32_global, mean, stdDev);
	}

at/atintegrators/atrandom.c

Line 46 in 8620504

static pcg32_random_t pcg32_global = PCG32_INITIALIZER;

at/atintegrators/atrandom.c

Lines 42 to 44 in 8620504

	// If you *must* statically initialize it, here's one.
	#define PCG32_INITIALIZER { AT_RNG_STATE, AT_RNG_INC, 0.0, false }

at/atintegrators/atrandom.c

Lines 27 to 28 in 8620504

	#define AT_RNG_STATE 0x853c49e6748fea9bULL
	#define AT_RNG_INC 0xda3e39cb94b95bdbULL

---

Moreover, and assuming all I asked before is correct. The C implementation of the generator seems to be compatible with pyat through numpy https://numpy.org/doc/2.0/reference/random/bit_generators/pcg64.html , but, I couldn't find the same for Matlab. I guess it means that the generator can not belong to the element because it needs to work in both python and matlab, therefore, it needs to be kept at the C level.

If one desires to use the element seed in tracking, a new stream should be initialized before starting tracking, which means checking the ring for random elements and create every stream.

Alternatively, one could forget about the element seed and create an input for the user seed (initstate,initseq) to initialize the streams in the parameter struct before the tracking, and make sure the elements use only one of the two streams (common_rng or thread_rng). The element might be modified to use one or the other.

A third solution would be to remove any seed input on the element. It would only change the documentation and help messages.

One other option would be getting PCG to work inside MATLAB, so, the element could have its own stream when defined, and the tracking would only the element own stream.

At last, maybe a more complicated way, one could change to a generator that exist in both python and matlab, and again implement the stream when creating the element.

[lfarv]

I'll put things in another way: random generators used in python and C are distinct, and as far as I can see, there is no specific random generator in Matlab, only the Matlab standard ones. The motivation is to avoid having the C engine depending on any external library. The whole implementation is in atrandom.c, and the behaviour is identical when running under python or Matlab. the C generators can be accessed from python but not from Matlab. For element passmethods, the generators which matter are the C generators, and the will work for both Matlab and python. There are 2 streams:

• the "thread" stream Param->thread_rng will generate different streams in different parallel threads, either for python parallel execution, MPI or OpenMP. This is the one used in QuantumDiff, since particles tracked in different streams must be independent.
• the "common" stream Param->common_rng will generate the same stream for different parallel threads. This one should be used in variable elements for instance, so that particles tracked in different parallel threads experience the same element.

The buffers for these streams are defined in at.c:

at/pyat/at.c

Lines 73 to 75 in 8620504

	/* state buffers for RNGs */
	static pcg32_random_t common_state = COMMON_PCG32_INITIALIZER;
	static pcg32_random_t thread_state = THREAD_PCG32_INITIALIZER;

There is a third stream using a static private buffer, with simpler function calls without the buffer argument. Its buffer is defined here:

at/atintegrators/atrandom.c

Line 46 in 8620504

static pcg32_random_t pcg32_global = PCG32_INITIALIZER;

It should be used only in regions without parallel processing.

There is no provision to have a specific stream for each element: element random values should be taken from one of the 2 tracking streams. So you are right, the seed specified in the variable multipole is not used. But the present implementation in VariableThinMPolePass is wrong: on should use atrandn_r(Param->common_rng, 0.0, 1.0) instead of atrandn(0.0, 1.0). As it is, particles tracked in parallel threads will not see the same multipole.

Now comes the question of explicitly seeding the generators. It is possible in python, but not in Matlab. For the moment, a python random generator is used to provide the seed of C generators at the beginning of each tracking. It would be easy to add an optional seed argument to Lattice.track to allow a user-defined seed. If I understand well, that's what you propose here:

Alternatively, one could forget about the element seed and create an input for the user seed (initstate,initseq) to initialize the streams in the parameter struct before the tracking, and make sure the elements use only one of the two streams (common_rng or thread_rng).

That's more difficult in Matlab: the tracking streams must be declared in atpass.c to be accessed by element passmethods. Then they cannot be seen by other mex-functions allowing explicit seeding. Defining them in a shared library accessed by several Mex-function is a solution, but how to compile and link a shared library with the Mex command?

Finally about the "python" generators defined in lattice.options.py: they are completely independent form the C ones, and I don't think they are used anywhere except seeding the C generators.

[oscarxblanco]

Dear @lfarv , thank you for the answer it is clearer now.

From what I see, there is no need to add another stream right now.
I will remove the seed parameter from the VariableThinMultipole element that I have been editing.
Should I also remove the seed parameter info in the help messages of QuantDiff ?

And one more question about atrandn_r(Param->common_rng, 0.0, 1.0). Would it be acceptable to pass a mean and std defined by the user ? i.e. atrandn_r(Param->common_rng, user_mean, user_std) ?

[lfarv]

Would it be acceptable to pass a mean and std defined by the user ? i.e. atrandn_r(Param->common_rng, user_mean, user_std) ?

Of course you can, the function is built for that. It depends on what you want to do.