Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Change how random seeds work in phase screens, fix for #61 #81

Merged
merged 2 commits into from
Mar 22, 2022
Merged

Change how random seeds work in phase screens, fix for #61 #81

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
50f65ac
changed how random seeds work in phase screens, fix for #61
ojdf Mar 21, 2022
4d609f2
added docstrings for seed
ojdf Mar 21, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
18 changes: 5 additions & 13 deletions aotools/turbulence/infinitephasescreen.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -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):
Expand All @@ -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()

Expand Down Expand Up @@ -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)

Expand All @@ -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])]

Expand Down
20 changes: 9 additions & 11 deletions aotools/turbulence/phasescreen.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -31,14 +31,13 @@ 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
"""
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
Expand Down Expand Up @@ -71,8 +70,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
Expand Down Expand Up @@ -102,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
Expand All @@ -117,10 +118,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)

Expand All @@ -136,7 +134,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

Expand Down