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.

[carmignani]

Hello,
I did not follow this discussion, but today I did some tracking in the cluster with compiled matlab and with the quantum excitation on and I noticed that the random number generator seed was the same in each of the processes. So I wanted to change the seeds, but I noticed that the seed parameter in the quantum diffusion element has been removed.
Is there a way now to give a different seed to each process?
The field "seed" as we had a few years ago in the QuantDiffPass was not good?

[oscarxblanco]

Hi Nicola,
I have gone through the modifications on QuantDiff.c, and the removal of the seed dates back to Dec,5 2022.

8ef93db#diff-cb146f57a02faf5a944b035fc9e297c5f528e322d1ecd9a801b856f8a13602d2

This commit mentions an update of the pass function that does not longer takes a seed.

And there are several commits before and after referring to fixing the compilation with OpenMP.

The current version does not take Seed as argument and when called from matlab it uses a global generator described by @lfarv above.
static pcg32_random_t pcg32_global = PCG32_INITIALIZER

[carmignani]

Thanks @oscarxblanco. I think we need a way to pass a seed to the c code from matlab, for when we run the same identical process many times. In the past it was possible to specify a different seed in the quantum diffusion element for each process, but now it is not possible.
Can't we have a seed input in ringpass for example?

[oscarxblanco]

Maybe @lfarv or @swhite2401 can answer that ?
I don't know enough about those modification to give an answer right away.

[swhite2401]

The seed can be set in python through a dedicated function, we had a quick look this morning with Nicola but could no find the matlab equivalent. I don't whether it is easy to implement but it should be possible maybe it even exist... we have to wait for @lfarv.

A short description of random generators in the documentation would be helpful as well.

[lfarv]

Hi all. Here are some explanations, in addition to the previous ones:

The random generator has been taken out of elements because of parallel processing. In the case of python multiprocessing and MPI, the generators would be copied and would generate the same sequence. That's why there are now 2 different streams: a common one and a thread-specific one, handled in the C "main" function, atpass.c for Matlab and at.c for python.

Resetting the seed can be done in python because we can have several entry points into the same atpass.c C module: one entry atpass for tracking, and another one reset_rng for seeding the generators. That is not possible in Matlab because mex-functions have a single entry MexFunction and are independent of each other.

So the only solution in Matlab would be to add another input argument to the atpass call, to optionally give a seed when starting the tracking. This was also mentioned above for the python Lattice.track method, and could also be done for the Matlab atpass function.

The problem is tricky because of the multiple parallel processing options: OpenMP for both python and Matlab, python multiprocessing and MPI for python only, and Matlab parallel-processing toolbox. The additional argument on tracking startup should work for all except Matlab parallel toolbox, though the generated sequence will be different for different methods.
If yo agree, I'll open a pull request for adding a new argument to Lattice.track and atpass.

[swhite2401]

Sounds good to me! Would it be possible to write a bit of documentation as well? Thanks!

[carmignani]

An additional input to atpass in matlab would be perfect for me.
Can we use OpenMP from matlab?

[lfarv]

OK, let's go for that.

Can we use OpenMP from matlab?

Yes:
>> atmexall -openmp

[carmignani]

Is there some documentation for openMP in matlab AT? I was not aware about it.
Is the parallel toolbox needed?

[lfarv]

@carmignani
Nothing to do, no toolbox: it's pure C, so it's just a compilation flag. The behaviour described for python here is also valid for Matlab: OpneMP is active as soon as the number of particles is equal or more that 10. The preliminary installation step for MacOS also applies.
But I agree, this should be duplicated in the Matlab section…

[lfarv]

#905 allows giving a seed for random generators when starting a tracking in Lattice.track (python), ringpass and linepass (Matlab)