# Path: cora/util/hputil.py
# def ang_positions(nside):
# def nside_for_lmax(lmax, accuracy_boost=1):
# def unpack_alm(alm, lmax, fullm=False):
# def pack_alm(almarray, lmax=None):
# def _make_full_alm(alm_half, centered=False):
# def _make_half_alm(alm_full):
# def sphtrans_real(hpmap, lmax=None, lside=None):
# def sphtrans_complex(hpmap, lmax=None, centered=False, lside=None):
# def sphtrans_real_pol(hpmaps, lmax=None, lside=None):
# def sphtrans_complex_pol(hpmaps, lmax=None, centered=False, lside=None):
# def sphtrans_inv_real(alm, nside):
# def sphtrans_inv_real_pol(alm, nside):
# def sphtrans_inv_complex(alm, nside):
# def sphtrans_sky(skymap, lmax=None):
# def sphtrans_inv_sky(alm, nside):
# def coord_x2y(map, x, y):
# def coord_g2c(map_):
# def coord_c2g(map_):
# def sph_ps(map1, map2=None, lmax=None):
# 
# Path: cora/util/nputil.py
# def save_ndarray_list(fname, la):
# def load_ndarray_list(fname):
# def matrix_root_manynull(mat, threshold=1e-16, truncate=True):
# def complex_std_normal(shape, rng=None):


import numpy as np
import scipy.linalg as la
import scipy.integrate as si
import healpy
from cora.util import hputil, nputil

    Returns
    -------
    hpmaps : np.ndarray (numz, npix)
        The Healpix maps. hpmaps[i] is the i'th map.
    """

    numz = corr.shape[1]
    maxl = corr.shape[0] - 1

    if corr.shape[2] != numz:
        raise Exception("Correlation matrix is incorrect shape.")

    alm_array = np.zeros((numz, 1, maxl + 1, maxl + 1), dtype=np.complex128)

    # Generate gaussian deviates and transform to have correct correlation
    # structure
    for l in range(maxl + 1):
        # Add in a small diagonal to try and ensure positive definiteness
        cmax = corr[l].diagonal().max() * 1e-14
        corrm = corr[l] + np.identity(numz) * cmax

        trans = nputil.matrix_root_manynull(corrm, truncate=False)
        gaussvars = nputil.complex_std_normal((numz, l + 1), rng=rng)
        alm_array[:, 0, l, : (l + 1)] = np.dot(trans, gaussvars)

    if alms:
        return alm_array

    # Perform the spherical harmonic transform for each z

# Path: cora/util/hputil.py
# def ang_positions(nside):
# def nside_for_lmax(lmax, accuracy_boost=1):
# def unpack_alm(alm, lmax, fullm=False):
# def pack_alm(almarray, lmax=None):
# def _make_full_alm(alm_half, centered=False):
# def _make_half_alm(alm_full):
# def sphtrans_real(hpmap, lmax=None, lside=None):
# def sphtrans_complex(hpmap, lmax=None, centered=False, lside=None):
# def sphtrans_real_pol(hpmaps, lmax=None, lside=None):
# def sphtrans_complex_pol(hpmaps, lmax=None, centered=False, lside=None):
# def sphtrans_inv_real(alm, nside):
# def sphtrans_inv_real_pol(alm, nside):
# def sphtrans_inv_complex(alm, nside):
# def sphtrans_sky(skymap, lmax=None):
# def sphtrans_inv_sky(alm, nside):
# def coord_x2y(map, x, y):
# def coord_g2c(map_):
# def coord_c2g(map_):
# def sph_ps(map1, map2=None, lmax=None):
# 
# Path: cora/util/nputil.py
# def save_ndarray_list(fname, la):
# def load_ndarray_list(fname):
# def matrix_root_manynull(mat, threshold=1e-16, truncate=True):
# def complex_std_normal(shape, rng=None):


# Path: cora/core/skysim.py
import numpy as np
import scipy.linalg as la
import scipy.integrate as si
import healpy
from cora.util import hputil, nputil


    Returns
    -------
    hpmaps : np.ndarray (numz, npix)
        The Healpix maps. hpmaps[i] is the i'th map.
    """

    numz = corr.shape[1]
    maxl = corr.shape[0] - 1

    if corr.shape[2] != numz:
        raise Exception("Correlation matrix is incorrect shape.")

    alm_array = np.zeros((numz, 1, maxl + 1, maxl + 1), dtype=np.complex128)

    # Generate gaussian deviates and transform to have correct correlation
    # structure
    for l in range(maxl + 1):
        # Add in a small diagonal to try and ensure positive definiteness
        cmax = corr[l].diagonal().max() * 1e-14
        corrm = corr[l] + np.identity(numz) * cmax

        trans = nputil.matrix_root_manynull(corrm, truncate=False)
        gaussvars = nputil.complex_std_normal((numz, l + 1), rng=rng)
        alm_array[:, 0, l, : (l + 1)] = np.dot(trans, gaussvars)

    if alms:
        return alm_array

    # Perform the spherical harmonic transform for each z

    sky = hputil.sphtrans_inv_sky(alm_array, nside)