Transition to the new numpy random Generator infrastructure

[lucabaldini]

https://numpy.org/doc/stable/reference/random/new-or-different.html#new-or-different

[lucabaldini]

This is interesting, and the numpy documentation is a very good source of information about the matter.

I have made some experiments on a branch, and demonstrated a 10% speedup by using one of the new bit generators. The main thing, though, is the new principle that it is better to instantiate a new generator than re-seed a new one, and this would imply moving away from the global-state paradigm we are using. This, in turn, involves a significant refactoring.

I can see, though, how this would be beneficial (in terms, e.g., of test isolation) and this might be worth the effort.

[lucabaldini]

At a cursory look

./hexsample/digi.py:from hexsample.rng import rng
./hexsample/digi.py:            pha += rng.normal(0., self.enc, size=pha.shape)
grep: ./hexsample/__pycache__/mc.cpython-310.pyc: binary file matches
grep: ./hexsample/__pycache__/digi.cpython-310.pyc: binary file matches
grep: ./hexsample/__pycache__/source.cpython-310.pyc: binary file matches
grep: ./hexsample/__pycache__/sensor.cpython-310.pyc: binary file matches
grep: ./hexsample/__pycache__/rng.cpython-310.pyc: binary file matches
./hexsample/rng.py:    return np.random.default_rng(bit_generator)
./hexsample/rng.py:    rng = random_generator(seed=seed)
./hexsample/rng.py:rng = random_generator()
./hexsample/modeling.py:from hexsample.rng import rng
./hexsample/modeling.py:            return self.ppf(rng.random(size), *self.parameters)
./hexsample/source.py:from hexsample.rng import rng
./hexsample/source.py:        r = self.radius * np.sqrt(rng.uniform(size=size))
./hexsample/source.py:        theta = rng.uniform(0., 2. * np.pi, size=size)
./hexsample/source.py:        x = rng.normal(self.x0, self.sigma, size=size)
./hexsample/source.py:        y = rng.normal(self.y0, self.sigma, size=size)
./hexsample/source.py:        return rng.choice(self._energies, size, replace=True, p=self._probs)
./hexsample/source.py:        timestamp = rng.uniform(tmin, tmax, size)
./hexsample/bin/hxsim.py:from hexsample.rng import rng, set_random_seed
./hexsample/sensor.py:from hexsample.rng import rng
./hexsample/sensor.py:        return np.round(rng.normal(mean, sigma)).astype(int)
./hexsample/sensor.py:        return dist.ppf(rng.uniform(0., dist.cdf(self.thickness)))
./hexsample/mc.py:from hexsample.rng import rng
./hexsample/mc.py:        x = rng.normal(self.absx, sigma, size=self.num_pairs)
./hexsample/mc.py:        y = rng.normal(self.absy, sigma, size=self.num_pairs)

the places where the rng is used are

• digi.digitize()
• source.Beam***.rvs()
• source.LineForest.rvs()
• source.Source.rvs_timestamp()
• sensor.Material.rvs_num_pairs()
• sensor.Sensor.absorption_depth()
• mc.MonteCarloEvent.propagate()

It almost look like it would make sense to have a dedicated class in hexsample.mc, with a rng class member, deputed to do all these calculations, given the proper arguments? (Although it is not entirely clear to me how the class hierarchy below source.BeamBase would work out.)

[lucabaldini]

Another option would be to pass a rng object to the constructors of all the classes that need to throw random numbers---that is, all the classes that have some rvs*() method. This would amount to:

• digi.HexagonalReadout
• source.BeamBase and subclasses
• source.SpectrumBase and subclasses
• sensor.Material
• sensor.Sensor
• mc.MonteCarloEvent

One way or the other, this is a fairly intrusive change.

[lucabaldini]

Pull request: https://github.com/lucabaldini/hexsample/pull/35

[lucabaldini]

Shipped in 0.6.0