Skip to content

ADD numpy random state argument to synalm #1

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

Open
sukhdeep2 wants to merge 4 commits into
base: master
Choose a base branch
from
Open

ADD numpy random state argument to synalm #1

Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
c1b7400
ADD numpy random state argument to synalm
sukhdeep1989 Jan 28, 2019
799ef4e
fixed typo
sukhdeep1989 Jan 28, 2019
72d9d12
changed from caps to lower case
sukhdeep1989 Jan 28, 2019
b5f2add
removed stack overflow link
sukhdeep1989 Jan 28, 2019
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
30 changes: 24 additions & 6 deletions healpy/sphtfunc.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -347,7 +347,7 @@ def alm2map(
return np.array(output)


def synalm(cls, lmax=None, mmax=None, new=False, verbose=True):
def synalm(cls, lmax=None, mmax=None, new=False, rng=None, verbose=True):
"""Generate a set of alm given cl.
The cl are given as a float array. Corresponding alm are generated.
If lmax is None, it is assumed lmax=cl.size-1
Expand All @@ -369,6 +369,10 @@ def synalm(cls, lmax=None, mmax=None, new=False, verbose=True):
(e.g. TT, EE, BB, TE, EB, TB or TT, EE, BB, TE if 4 cl as input).
If False, use the old ordering, ie by row
(e.g. TT, TE, TB, EE, EB, BB or TT, TE, EE, BB if 4 cl as input).

rng: numpy random number generator state. Useful when generating multiple
maps in parallel. Typically, `rng=np.random.RandomState(seed)`


Returns
-------
Expand Down Expand Up @@ -408,8 +412,14 @@ def synalm(cls, lmax=None, mmax=None, new=False, verbose=True):
cls_list = [np.asarray(cls, dtype=np.float64)]
szalm = Alm.getsize(lmax, mmax)
alm = np.zeros(szalm, "D")
alm.real = np.random.standard_normal(szalm)
alm.imag = np.random.standard_normal(szalm)

if rng is None:
alm.real = np.random.standard_normal(szalm)
alm.imag = np.random.standard_normal(szalm)
else:
alm.real = rng.standard_normal(szalm)
alm.imag = rng.standard_normal(szalm)

alms_list = [alm]
sphtlib._synalm(cls_list, alms_list, lmax, mmax)
return alm
Expand Down Expand Up @@ -442,8 +452,14 @@ def synalm(cls, lmax=None, mmax=None, new=False, verbose=True):
alms_list = []
for i in six.moves.xrange(Nspec):
alm = np.zeros(szalm, "D")
alm.real = np.random.standard_normal(szalm)
alm.imag = np.random.standard_normal(szalm)

if rng is None:
alm.real = np.random.standard_normal(szalm)
alm.imag = np.random.standard_normal(szalm)
else:
alm.real = rng.standard_normal(szalm)
alm.imag = rng.standard_normal(szalm)

alms_list.append(alm)
if new: # new input order: input given by diagonal, should be given by row
cls_list = new_to_old_spectra_order(cls_list)
Expand All @@ -467,6 +483,7 @@ def synfast(
fwhm=0.0,
sigma=None,
new=False,
rng=None,
verbose=True,
):
"""Create a map(s) from cl(s).
Expand Down Expand Up @@ -497,6 +514,7 @@ def synfast(
sigma : float, scalar, optional
The sigma of the Gaussian used to smooth the map (applied on alm)
[in radians]
rng: Numpy random number generator state, passed to synalm.

Returns
-------
Expand Down Expand Up @@ -524,7 +542,7 @@ def synfast(
cls_lmax = cb.len_array_or_arrays(cls) - 1
if lmax is None or lmax < 0:
lmax = min(cls_lmax, 3 * nside - 1)
alms = synalm(cls, lmax=lmax, mmax=mmax, new=new, verbose=verbose)
alms = synalm(cls, lmax=lmax, mmax=mmax, new=new, rng=rng, verbose=verbose)
maps = alm2map(
alms,
nside,
Expand Down