From 50f65ac82688ee5aa46174a72e649d88757d816d Mon Sep 17 00:00:00 2001 From: Ollie Farley Date: Mon, 21 Mar 2022 16:27:42 +0000 Subject: [PATCH 1/2] changed how random seeds work in phase screens, fix for #61 --- aotools/turbulence/infinitephasescreen.py | 18 +++++------------- aotools/turbulence/phasescreen.py | 16 +++++----------- 2 files changed, 10 insertions(+), 24 deletions(-) diff --git a/aotools/turbulence/infinitephasescreen.py b/aotools/turbulence/infinitephasescreen.py index 8c521e9..53ae6df 100644 --- a/aotools/turbulence/infinitephasescreen.py +++ b/aotools/turbulence/infinitephasescreen.py @@ -177,19 +177,16 @@ def make_initial_screen(self): # phase screen will make it *really* random if no seed at all given. # If a seed is here, screen must be repeatable wiht same seed - if self.random_seed is not None: - seed = numpy.random.randint(0, 2**32 -1) - else: - seed = None + self._R = numpy.random.default_rng(self.random_seed) self._scrn = phasescreen.ft_phase_screen( - self.r0, self.stencil_length, self.pixel_scale, self.L0, 1e-10, seed=seed + self.r0, self.stencil_length, self.pixel_scale, self.L0, 1e-10, seed=self._R ) self._scrn = self._scrn[:, :self.nx_size] def get_new_row(self): - random_data = numpy.random.normal(0, 1, size=self.nx_size) + random_data = self._R.normal(0, 1, size=self.nx_size) stencil_data = self._scrn[(self.stencil_coords[:, 0], self.stencil_coords[:, 1])] new_row = self.A_mat.dot(stencil_data) + self.B_mat.dot(random_data) @@ -270,6 +267,7 @@ class PhaseScreenVonKarman(PhaseScreen): pixel_scale(float): Size of each phase pixel in metres r0 (float): fried parameter (metres) L0 (float): Outer scale (metres) + random_seed (int, optional): seed for the random number generator n_columns (int, optional): Number of columns to use to continue screen, default is 2 """ def __init__(self, nx_size, pixel_scale, r0, L0, random_seed=None, n_columns=2): @@ -286,9 +284,6 @@ def __init__(self, nx_size, pixel_scale, r0, L0, random_seed=None, n_columns=2): self.random_seed = random_seed - if random_seed is not None: - numpy.random.seed(random_seed) - self.set_X_coords() self.set_stencil_coords() @@ -400,9 +395,6 @@ def __init__(self, nx_size, pixel_scale, r0, L0, random_seed=None, stencil_lengt self.stencil_length = stencil_length_factor * self.nx_size self.random_seed = random_seed - if random_seed is not None: - numpy.random.seed(random_seed) - # Coordinate of Fried's "reference point" that stops the screen diverging self.reference_coord = (1, 1) @@ -417,7 +409,7 @@ def __init__(self, nx_size, pixel_scale, r0, L0, random_seed=None, stencil_lengt self.make_initial_screen() def get_new_row(self): - random_data = numpy.random.normal(0, 1, size=self.nx_size) + random_data = self._R.normal(0, 1, size=self.nx_size) stencil_data = self._scrn[(self.stencil_coords[:, 0], self.stencil_coords[:, 1])] diff --git a/aotools/turbulence/phasescreen.py b/aotools/turbulence/phasescreen.py index b41f078..ab3ce2a 100644 --- a/aotools/turbulence/phasescreen.py +++ b/aotools/turbulence/phasescreen.py @@ -35,10 +35,7 @@ def ft_sh_phase_screen(r0, N, delta, L0, l0, FFT=None, seed=None): Returns: ndarray: numpy array representing phase screen in radians """ - R = random.SystemRandom(time.time()) - if seed is None: - seed = int(R.random()*100000) - numpy.random.seed(seed) + R = numpy.random.default_rng(seed) D = N * delta # high-frequency screen from FFT method @@ -71,8 +68,8 @@ def ft_sh_phase_screen(r0, N, delta, L0, l0, FFT=None, seed=None): PSD_phi[1,1] = 0 # random draws of Fourier coefficients - cn = ( (numpy.random.normal(size=(3,3)) - + 1j*numpy.random.normal(size=(3,3)) ) + cn = ( (R.normal(size=(3,3)) + + 1j*R.normal(size=(3,3)) ) * numpy.sqrt(PSD_phi)*del_f ) SH = numpy.zeros((N,N),dtype="complex") # loop over frequencies on this grid @@ -117,10 +114,7 @@ def ft_phase_screen(r0, N, delta, L0, l0, FFT=None, seed=None): L0 = float(L0) l0 = float(l0) - R = random.SystemRandom(time.time()) - if seed is None: - seed = int(R.random()*100000) - numpy.random.seed(seed) + R = numpy.random.default_rng(seed) del_f = 1./(N*delta) @@ -136,7 +130,7 @@ def ft_phase_screen(r0, N, delta, L0, l0, FFT=None, seed=None): PSD_phi[int(N/2), int(N/2)] = 0 - cn = ((numpy.random.normal(size=(N, N))+1j * numpy.random.normal(size=(N, N))) * numpy.sqrt(PSD_phi)*del_f) + cn = ((R.normal(size=(N, N))+1j * R.normal(size=(N, N))) * numpy.sqrt(PSD_phi)*del_f) phs = ift2(cn, 1, FFT).real From 4d609f2fb174c7014a86cec91e3d18a49f0f3952 Mon Sep 17 00:00:00 2001 From: Ollie Farley Date: Mon, 21 Mar 2022 17:03:51 +0000 Subject: [PATCH 2/2] added docstrings for seed --- aotools/turbulence/phasescreen.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/aotools/turbulence/phasescreen.py b/aotools/turbulence/phasescreen.py index ab3ce2a..b5d4f6c 100644 --- a/aotools/turbulence/phasescreen.py +++ b/aotools/turbulence/phasescreen.py @@ -31,6 +31,8 @@ def ft_sh_phase_screen(r0, N, delta, L0, l0, FFT=None, seed=None): delta (float): size in Metres of each pxl L0 (float): Size of outer-scale in metres l0 (float): inner scale in metres + seed (int, optional): seed for random number generator. If provided, + allows for deterministic screens Returns: ndarray: numpy array representing phase screen in radians @@ -99,6 +101,8 @@ def ft_phase_screen(r0, N, delta, L0, l0, FFT=None, seed=None): delta (float): size in Metres of each pxl L0 (float): Size of outer-scale in metres l0 (float): inner scale in metres + seed (int, optional): seed for random number generator. If provided, + allows for deterministic screens .. note:: The phase screen is returned as a 2d array, with each element representing the phase