) -> Generator[Iternorm, None, None]:

'''return the vector a and variance sigma^2 for iterative normal sampling'''

n: Tuple[int, ...]

if size is None:

n = ()

elif isinstance(size, int):

n = (size,)

else:

n = size

m = np.zeros((*n, k, k))

a = np.zeros((*n, k))

s = np.zeros((*n,))

q = (*n, k+1)

j = 0 if k > 0 else None

for i, x in enumerate(cov):

x = np.asanyarray(x)

if x.shape != q:

try:

x = np.broadcast_to(x, q)

except ValueError:

raise TypeError(f'covariance row {i}: shape {x.shape} cannot be broadcast to {q}') from None

# only need to update matrix A if there are correlations

if j is not None:

# compute new entries of matrix A

m[..., :, j] = 0

m[..., j:j+1, :] = np.matmul(a[..., np.newaxis, :], m)

m[..., j, j] = np.where(s != 0, -1, 0)

np.divide(m[..., j, :], -s[..., np.newaxis], where=(m[..., j, :] != 0), out=m[..., j, :])

# compute new vector a

c = x[..., 1:, np.newaxis]

a = np.matmul(m[..., :j], c[..., k-j:, :])

a += np.matmul(m[..., j:], c[..., :k-j, :])

a = a.reshape(*n, k)

# next rolling index

j = (j - 1) % k

# compute new standard deviation

s = x[..., 0] - np.einsum('...i,...i', a, a)

if np.any(s < 0):

raise ValueError('covariance matrix is not positive definite')

s = np.sqrt(s)

# yield the next index, vector a, and standard deviation s

) -> Generator[Array, None, None]:

'''Return array of cls as a covariance matrix for iterative sampling.'''

cov = np.zeros((nl, nc+1))

end = 0

for j in range(nf):

begin, end = end, end + j + 1

for i, cl in enumerate(cls[begin:end][:nc+1]):

if cl is None:

cov[:, i] = 0

else:

if i == 0 and np.any(np.less(cl, 0)):

raise ValueError('negative values in cl')

n = len(cl)

cov[:n, i] = cl

cov[n:, i] = 0

cov /= 2

'''multiply alm by bl'''

n = len(bl)

if inplace:

out = np.asanyarray(alm)

else:

out = np.copy(alm)

for m in range(n):

out[m*n-m*(m-1)//2:(m+1)*n-m*(m+1)//2] *= bl[m:]

) -> Cls:

'''Transform Cls to Gaussian Cls.'''

gls = []

for cl in cls:

if cl is not None and len(cl) > 0:

if cl[0] == 0:

monopole = 0.

else:

monopole = None

gl, info, err, niter = gaussiancl(cl, tfm, pars, monopole=monopole)

if info == 0:

warnings.warn('Gaussian cl did not converge, inexact transform')

else:

gl = []

gls.append(gl)

ncorr: Optional[int] = None, nside: Optional[int] = None

) -> Cls:

'''Compute Gaussian Cls for a Gaussian random field.

if ncorr is not None:

n = int((2*len(cls))**0.5)

if n*(n+1)//2 != len(cls):

raise ValueError('length of cls array is not a triangle number')

cls = [cls[i*(i+1)//2+j] if j <= ncorr else [] for i in range(n) for j in range(i+1)]

if nside is not None:

pw = hp.pixwin(nside, lmax=lmax)

gls = []

for cl in cls:

if cl is not None and len(cl) > 0:

if lmax is not None:

cl = cl[:lmax+1]

if nside is not None:

n = min(len(cl), len(pw))

cl = cl[:n] * pw[:n]**2

gls.append(cl)

ncorr: Optional[int] = None, nside: Optional[int] = None

) -> Cls:

'''Compute Gaussian Cls for a lognormal random field.'''

gls = gaussian_gls(cls, lmax=lmax, ncorr=ncorr, nside=nside)

rng: Optional[np.random.Generator] = None

) -> Generator[Array, None, None]:

'''Iteratively sample Gaussian random fields from Cls.

# get the default RNG if not given

if rng is None:

rng = np.random.default_rng()

# number of gls and number of fields

ngls = len(gls)

ngrf = int((2*ngls)**0.5)

# number of correlated fields if not specified

if ncorr is None:

ncorr = ngrf - 1

# number of modes

n = max((len(gl) for gl in gls if gl is not None), default=0)

if n == 0:

raise ValueError('all gls are empty')

# generates the covariance matrix for the iterative sampler

cov = cls2cov(gls, n, ngrf, ncorr)

# working arrays for the iterative sampling

z = np.zeros(n*(n+1)//2, dtype=np.complex128)

y = np.zeros((n*(n+1)//2, ncorr), dtype=np.complex128)

# generate the conditional normal distribution for iterative sampling

conditional_dist = iternorm(ncorr, cov, size=n)

# sample the fields from the conditional distribution

for j, a, s in conditional_dist:

# standard normal random variates for alm

# sample real and imaginary parts, then view as complex number

rng.standard_normal(n*(n+1), np.float64, z.view(np.float64))

# scale by standard deviation of the conditional distribution

# variance is distributed over real and imaginary part

alm = multalm(z, s)

# add the mean of the conditional distribution

for i in range(ncorr):

alm += multalm(y[:, i], a[:, i])

# store the standard normal in y array at the indicated index

if j is not None:

y[:, j] = z

# modes with m = 0 are real-valued and come first in array

alm[:n].real += alm[:n].imag

alm[:n].imag[:] = 0

# transform alm to maps

# can be performed in place on the temporary alm array

ncorr: Optional[int] = None,

rng: Optional[np.random.Generator] = None

) -> Generator[Array, None, None]:

'''Iterative sample lognormal random fields from Gaussian Cls.'''

for i, m in enumerate(generate_gaussian(gls, nside, ncorr=ncorr, rng=rng)):

# compute the variance of the auto-correlation

gl = gls[i*(i+1)//2]

ell = np.arange(len(gl))

var = np.sum((2*ell + 1)*gl)/(4*np.pi)

# fix mean of the Gaussian random field for lognormal transformation

m -= var/2

# exponentiate values in place and subtract 1 in one operation

np.expm1(m, out=m)

# lognormal shift, unless unity

if shift != 1:

m *= shift

# yield the lognormal map

"""Return a specific angular power spectrum from an array.

if j > i:

i, j = j, i

cl = cls[i*(i+1)//2+i-j]

if lmax is not None:

if len(cl) > lmax + 1:

cl = cl[:lmax+1]

else:

cl = np.pad(cl, (0, lmax + 1 - len(cl)))

"""Compute effective angular power spectra from weights.

from itertools import combinations_with_replacement, product

# this is the number of fields

n = int((2*len(cls))**0.5)

if n*(n+1)//2 != len(cls):

raise ValueError("length of cls is not a triangle number")

# find lmax if not given

if lmax is None:

lmax = max(map(len, cls), default=0) - 1

# broadcast weights1 such that its shape ends in n

weights1 = np.asanyarray(weights1)

weights2 = np.asanyarray(weights2) if weights2 is not None else weights1

shape1, shape2 = weights1.shape, weights2.shape

for i, shape in enumerate((shape1, shape2)):

if not shape or shape[0] != n:

raise ValueError(f"shape mismatch between fields and weights{i+1}")

# get the iterator over leading weight axes

# auto-spectra do not repeat identical computations

if weights2 is weights1:

pairs = combinations_with_replacement(np.ndindex(shape1[1:]), 2)

else:

pairs = product(np.ndindex(shape1[1:]), np.ndindex(shape2[1:]))

# create the output array: axes for all input axes plus lmax+1

out = np.empty(shape1[1:] + shape2[1:] + (lmax+1,))

# helper that will grab the entire first column (i.e. shells)

c = (slice(None),)

# compute all combined cls from pairs

# if weights2 is weights1, set the transpose elements in one pass

for j1, j2 in pairs:

w1, w2 = weights1[c + j1], weights2[c + j2]

cl = sum(

w1[i1] * w2[i2] * getcl(cls, i1, i2, lmax=lmax)

for i1, i2 in np.ndindex(n, n)

)

out[j1 + j2] = cl

if weights2 is weights1 and j1 != j2:

out[j2 + j1] = cl