simulation.py 9.42 KB
Newer Older
1
"""
Syl's avatar
Syl committed
2
Set of routines to generate basic simulations
3 4 5 6 7 8 9 10 11 12

Author: Vanneste
"""
from __future__ import division

import numpy as np
import healpy as hp

from scipy import sparse

Syl's avatar
Syl committed
13

14 15 16 17 18 19
def muKarcmin2var(muKarcmin, nside):
    """
    Return pixel variance for a given nside and noise level [1e-6 K . arcmin]

    Parameters
    ----------
Syl's avatar
Syl committed
20 21 22 23
    muKarcmin : float
        Pixel noise [muK . arcmin]
    nside : int
        Healpix map resolution (power of 2)
24 25 26

    Returns
    ----------
Syl's avatar
Syl committed
27 28
    varperpix : float
        Variance per pixel
Julien's avatar
Julien committed
29 30 31 32 33 34

    Example
    ----------
    >>> var = muKarcmin2var(10.0, 16)
    >>> print(round(var * 1e16))
    21.0
35
    """
Julien's avatar
Julien committed
36 37
    pixarea = hp.nside2pixarea(nside, degrees=True)
    varperpix = (muKarcmin * 1e-6 / 60.)**2 / pixarea
38 39
    return varperpix

Syl's avatar
Syl committed
40

41 42
def pixvar2nl(pixvar, nside):
    """
Julien's avatar
Julien committed
43
    Return noise spectrum level for a given nside and pixel variance
44 45 46

    Parameters
    ----------
Syl's avatar
Syl committed
47 48 49 50
    pixvar : float
        Variance per pixel
    nside : int
        Healpix map resolution (power of 2)
51 52 53

    Returns
    ----------
Syl's avatar
Syl committed
54 55
    nl : float
        Noise spectrum level
Julien's avatar
Julien committed
56 57 58 59 60 61 62 63

    Example
    ----------
    >>> nside = 16
    >>> pixvar = muKarcmin2var(10.0, nside)
    >>> nl = pixvar2nl(pixvar, nside)
    >>> print(round(nl * 1e18))
    8.0
64
    """
Syl's avatar
Syl committed
65 66 67
    nl = pixvar*4.*np.pi/(12*nside**2.)
    return nl

68 69 70 71 72 73 74

def getNl(pixvar, nside, nbins):
    """
    Return noise spectrum for a given nside and pixel variance

    Parameters
    ----------
Syl's avatar
Syl committed
75 76 77 78 79 80
    pixvar : float
        Variance per pixel
    nside : int
        Healpix map resolution (power of 2)
    nbins : int
        Number of bins
81 82 83

    Returns
    ----------
Syl's avatar
Syl committed
84 85
    Nl : 1D array of float
        Noise spectrum
Julien's avatar
Julien committed
86 87 88 89 90 91 92 93

    Example
    ----------
    >>> nside = 16
    >>> pixvar = muKarcmin2var(10.0, nside)
    >>> nl = getNl(pixvar, nside, 2)
    >>> print(round(nl[0] * 1e18))
    8.0
94
    """
Syl's avatar
Syl committed
95 96
    Nl = pixvar*4.*np.pi/(12*nside**2.)*np.ones((nbins))
    return Nl
97 98


Syl's avatar
Syl committed
99
def getstokes(spec=None, temp=False, polar=False, corr=False):
100
    """
Syl's avatar
Syl committed
101
    Get the Stokes parameters number and name(s)
102 103 104

    Parameters
    ----------
Syl's avatar
Syl committed
105 106
    spec : bool
        If True, get Stokes parameters for polar (default: True)
Syl's avatar
Syl committed
107 108 109 110
    polar : bool
        If True, get Stokes parameters for polar (default: True)
    temp : bool
        If True, get Stokes parameters for temperature (default: False)
Syl's avatar
Syl committed
111
    corr : bool
Syl's avatar
Syl committed
112
        If True, get Stokes parameters for EB and TB (default: False)
113 114 115

    Returns
    ----------
Syl's avatar
Syl committed
116
    stokes : list of string
Syl's avatar
Syl committed
117 118 119
        Stokes variables names
    spec : int
        Spectra names
Syl's avatar
Syl committed
120
    istokes : list
Syl's avatar
Syl committed
121
        Indexes of power spectra
Julien's avatar
Julien committed
122 123 124

    Example
    ----------
Syl's avatar
Syl committed
125
    >>> getstokes(polar=True, temp=False, corr=False)
Julien's avatar
Julien committed
126
    (['Q', 'U'], ['EE', 'BB'], [1, 2])
Syl's avatar
Syl committed
127
    >>> getstokes(polar=True, temp=True, corr=False)
Julien's avatar
Julien committed
128
    (['I', 'Q', 'U'], ['TT', 'EE', 'BB', 'TE'], [0, 1, 2, 3])
Syl's avatar
Syl committed
129
    >>> getstokes(polar=True, temp=True, corr=True)
Julien's avatar
Julien committed
130
    (['I', 'Q', 'U'], ['TT', 'EE', 'BB', 'TE', 'EB', 'TB'], [0, 1, 2, 3, 4, 5])
131
    """
Syl's avatar
Syl committed
132 133 134 135 136 137 138
    if spec is not None:
        _temp = "TT" in spec or "TE" in spec or "TB" in spec or temp
        _polar = "EE" in spec or "BB" in spec or "TE" in spec or "TB" in \
            spec or "EB" in spec or polar
        _corr = "TE" in spec or "TB" in spec or "EB" in spec or corr
        if not _temp and not _polar and not _corr:
            print("invalid spectra list")
139
    else:
Syl's avatar
Syl committed
140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168
        _temp = temp
        _polar = polar
        _corr = corr

    speclist = []
    if _temp or (spec is None and corr):
        speclist.extend(["TT"])
    if _polar:
        speclist.extend(["EE", "BB"])
    if spec is not None and not corr:
        if 'TE' in spec:
            speclist.extend(["TE"])
        if 'EB' in spec:
            speclist.extend(["EB"])
        if 'TB' in spec:
            speclist.extend(["TB"])

    elif _corr:
        speclist.extend(["TE", "EB", "TB"])

    stokes = []
    if _temp:
        stokes.extend(["I"])
    if _polar:
        stokes.extend(["Q", "U"])

    ispecs = [['TT', 'EE', 'BB', 'TE', 'EB', 'TB'].index(s) for s in speclist]
    istokes = [['I', 'Q', 'U'].index(s) for s in stokes]
    return stokes, speclist, istokes, ispecs
Syl's avatar
Syl committed
169

170

Julien's avatar
Julien committed
171 172
def GetBinningMatrix(
        ellbins, lmax, norm=False, polar=True,
Syl's avatar
Syl committed
173
        temp=False, corr=False):
174
    """
Syl's avatar
Syl committed
175 176 177 178 179
    Return P (m,n) and Q (n,m) binning matrices such that
    Cb = P.Cl and Vbb = P.Vll.Q with m the number of bins and
    n the number of multipoles.
    In addition, returns ell (total non-binned multipole range)
    and ellval (binned multipole range)
180 181 182

    Parameters
    ----------
Syl's avatar
Syl committed
183 184 185 186 187 188 189 190 191 192
    ellbins : list of integers
        Bins lower bound
    lmax : int
        Maximum multipole
    norm : bool (default: False)
        If True, weight the binning scheme such that P = l*(l+1)/(2*pi)
    polar : bool
        If True, get Stokes parameters for polar (default: True)
    temp : bool
        If True, get Stokes parameters for temperature (default: False)
Syl's avatar
Syl committed
193
    corr : bool
Syl's avatar
Syl committed
194
        If True, get Stokes parameters for EB and TB (default: False)
195 196 197

    Returns
    ----------
Syl's avatar
Syl committed
198 199 200 201 202 203 204 205 206
    P : array of float (m,n)
        Binning matrix such that Cb = P.Cl
    Q : array of int (n,m)
        Binning matrix such that P.Q = I
        or P.Q = I * l(l+1)/(2pi) if norm=True
    ell : array of int (n)
        Multipoles range
    ellvall : array of float (m)
        Bins pivot range
Julien's avatar
Julien committed
207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248

    Example
    ----------
    >>> bins = np.array([2.0, 5.0, 10.0])
    >>> P, Q, ell, ellval = GetBinningMatrix(
    ...     bins, 10.0)
    >>> print(P) # doctest: +NORMALIZE_WHITESPACE
    [[ 0.33333333  0.33333333  0.33333333  0.          0.          0.       0.
       0.          0.          0.          0.          0.          0.       0.
       0.          0.          0.          0.        ]
     [ 0.          0.          0.          0.2         0.2         0.2      0.2
       0.2         0.          0.          0.          0.          0.       0.
       0.          0.          0.          0.        ]
     [ 0.          0.          0.          0.          0.          0.       0.
       0.          0.          0.33333333  0.33333333  0.33333333  0.       0.
       0.          0.          0.          0.        ]
     [ 0.          0.          0.          0.          0.          0.       0.
       0.          0.          0.          0.          0.          0.2      0.2
       0.2         0.2         0.2         0.        ]]
    >>> print(Q) # doctest: +NORMALIZE_WHITESPACE
    [[ 1.  0.  0.  0.]
     [ 1.  0.  0.  0.]
     [ 1.  0.  0.  0.]
     [ 0.  1.  0.  0.]
     [ 0.  1.  0.  0.]
     [ 0.  1.  0.  0.]
     [ 0.  1.  0.  0.]
     [ 0.  1.  0.  0.]
     [ 0.  0.  0.  0.]
     [ 0.  0.  1.  0.]
     [ 0.  0.  1.  0.]
     [ 0.  0.  1.  0.]
     [ 0.  0.  0.  1.]
     [ 0.  0.  0.  1.]
     [ 0.  0.  0.  1.]
     [ 0.  0.  0.  1.]
     [ 0.  0.  0.  1.]
     [ 0.  0.  0.  0.]]
    >>> print(ell)
    [  2.   3.   4.   5.   6.   7.   8.   9.  10.  11.]
    >>> print(ellval)
    [ 3.  7.]
249
    """
Syl's avatar
Syl committed
250
    # ### define Stokes
Syl's avatar
Syl committed
251
    nspec = 1
252

Julien's avatar
Julien committed
253 254 255 256 257
    nbins = len(ellbins) - 1
    ellmin = np.array(ellbins[0: nbins])
    ellmax = np.array(ellbins[1: nbins + 1]) - 1
    ell = np.arange(np.min(ellbins), lmax + 2)
    maskl = (ell[:-1] < (lmax + 2)) & (ell[:-1] > 1)
258

Julien's avatar
Julien committed
259 260 261 262 263
    # define min
    minell = np.array(ellbins[0: nbins])
    # and max of a bin
    maxell = np.array(ellbins[1: nbins + 1]) - 1
    ellval = (minell + maxell) * 0.5
264

Julien's avatar
Julien committed
265
    masklm = []
266
    for i in np.arange(nbins):
Julien's avatar
Julien committed
267
        masklm.append(((ell[:-1] >= minell[i]) & (ell[:-1] <= maxell[i])))
268

Syl's avatar
Syl committed
269
    allmasklm = nspec*[list(masklm)]
270
    masklM = np.array(sparse.block_diag(allmasklm[:]).toarray())
Julien's avatar
Julien committed
271
    binsnorm = np.array(
Syl's avatar
Syl committed
272
        nspec * [list(np.arange(minell[0], np.max(ellbins)))]).flatten()
273 274 275 276

    binsnorm = binsnorm*(binsnorm+1)/2./np.pi
    P = np.array(masklM)*1.
    Q = P.T
Julien's avatar
Julien committed
277
    P = P / np.sum(P, 1)[:, None]
278
    if norm:
Julien's avatar
Julien committed
279
        P *= binsnorm
280 281 282 283

    return P, Q, ell, ellval


Syl's avatar
Syl committed
284
def extrapolpixwin(nside, Slmax, pixwin=True):
Syl's avatar
Syl committed
285
    '''
286 287
    Parameters
    ----------
Syl's avatar
Syl committed
288 289 290 291
    nside : int
        Healpix map resolution
    Slmax : int
        Maximum multipole value computed for the pixel covariance pixel matrix
Syl's avatar
Syl committed
292
    pixwin : bool
Syl's avatar
Syl committed
293 294
        If True, multiplies the beam window function by the pixel
        window function. Default: True
295 296 297

    Returns
    ----------
Syl's avatar
Syl committed
298
    fpixwin : array of floats
299

Syl's avatar
Syl committed
300
    Example :
301
    ----------
Syl's avatar
Syl committed
302 303 304 305 306 307 308 309 310 311 312 313
    >>> print(hp.pixwin(2))
    [ 1.          0.977303    0.93310702  0.86971852  0.79038278  0.69905215
      0.60011811  0.49813949  0.39760902]
    >>> print(extrapolpixwin(2, 20, True))
    [ 1.          0.977303    0.93310702  0.86971852  0.79038278  0.69905215
      0.60011811  0.49813949  0.39760902  0.30702636  0.22743277  0.16147253
      0.10961864  0.07098755  0.04374858  0.02559774  0.01418623  0.00742903
      0.00366749  0.00170274]
    >>> print(extrapolpixwin(2, 20, False))
    [ 1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.
      1.  1.]
    '''
Syl's avatar
Syl committed
314
    if pixwin:
Syl's avatar
Syl committed
315 316 317 318 319 320 321 322
        prepixwin = np.array(hp.pixwin(nside))
        poly = np.polyfit(np.arange(len(prepixwin)), np.log(prepixwin),
                          deg=3, w=np.sqrt(prepixwin))
        y_int = np.polyval(poly, np.arange(Slmax))
        fpixwin = np.exp(y_int)
        fpixwin = np.append(prepixwin, fpixwin[len(prepixwin):])[: Slmax]
    else:
        fpixwin = np.ones((Slmax))
323

Syl's avatar
Syl committed
324
    return fpixwin
Julien's avatar
Julien committed
325 326 327


if __name__ == "__main__":
Julien's avatar
Julien committed
328 329 330 331 332 333 334 335
    """
    Run the doctest using

    python simulation.py

    If the tests are OK, the script should exit gracefuly, otherwise the
    failure(s) will be printed out.
    """
Julien's avatar
Julien committed
336 337 338
    import doctest
    if np.__version__ >= "1.14.0":
        np.set_printoptions(legacy="1.13")
Syl's avatar
Syl committed
339
    doctest.testmod()