In [17]:
%matplotlib inline

In [33]:
import matplotlib.pyplot as plt
import numpy as np
import healpy as hp

In [1]:
import numpy as np
import seaborn as sns
import healpy as hp
from astropy.io import fits
from matplotlib import gridspec
from matplotlib import ticker
from scipy.integrate import quad
import lib.cmb
import lib.noise
import calculations
import util.plot_funcs
import time
import pickle

%pylab inline

pylab.rcParams['figure.figsize'] = (10.0, 8.0)


Populating the interactive namespace from numpy and matplotlib
/usr/local/lib/python2.7/site-packages/pytz/__init__.py:35: UserWarning: Module argparse was already imported from /usr/local/Cellar/python/2.7.3/Frameworks/Python.framework/Versions/2.7/lib/python2.7/argparse.pyc, but /usr/local/lib/python2.7/site-packages is being added to sys.path
  from pkg_resources import resource_stream

NOTES/TODO:

Make hits map generator

- Do for single pixels in 4 corners of array
- Also do for full sky
- See how much overlap we have of full arrays

Inject gain noise and measure B-mode contamination

- Start with r = 0 cmb realization
- Make gain transfer function (in C_l space, constant function with a tick up at the desired ell)
- Make M realizations of gain transfer function
- Multiply cmb realization with gain transfer function realization
- Generate ILC with each cmb*gain map
- Estimate C_l's from each ILC map
- Repeat for N cmb realizations (N*M total runs)
- Repeat for L ell bins (how many???)

In [56]:
from IPython import parallel as p
rc = p.Client()
cluster = rc[:]
cluster.block = True

In [13]:
reload(calculations.util)


================================================================================
Importing calibration_gain.
cldict_normal_test already exists at /Users/jlazear/github/cmb/data/cldict_normal_test.hd5. Loading it.
cldict_noise already exists at /Users/jlazear/github/cmb/data/cldict_noise.hd5. Loading it.
Finished importing calibration_gain.
================================================================================
Out[13]:
<module 'calculations.util' from 'calculations/util.py'>

In [25]:
import calculations.util, calculations.calibration_gain
read_dict_from_hd5 = calculations.calibration_gain.read_dict_from_hd5

In [26]:
a4ln = read_dict_from_hd5('data/all_4_lownoise.hd5')['all_4_lownoise']

In [27]:
testdata = a4ln['ell_5']

In [30]:
fig = calculations.util.myplot(testdata, 'testdata', 'testdatalabel')
ax = fig.axes[0]
ax.set_ylim(ymax=.0001)
ax.set_xscale('linear')
ax.set_xlim(xmax=10)


Out[30]:
(2.0, 10)

In [36]:
npix = hp.nside2npix(512)
calmap = np.array([1.]*npix)
np.unique(hp.anafast(calmap))


Out[36]:
array([ 0.])

In [82]:
cls = np.zeros(400)
cls[2] = 1
map1 = hp.synfast(cls, 512)
map1 = 1. + map1*0.05/map1.std()
hp.mollview(map1)


Sigma is 0.000000 arcmin (0.000000 rad) 
-> fwhm is 0.000000 arcmin
The order of the input cl's will change in a future release.
Use new=True keyword to start using the new order.
See documentation of healpy.synalm.

In [114]:
cls_gain = np.zeros(1500)
cls_gain[500] = 1
map_gain = hp.synfast(cls_gain, 512)
gain_error = np.sqrt(2)
map_gain = 1. + map_gain*gain_error/map_gain.std()

cls_cmb = np.arange(1500, dtype='float')
map_cmb = hp.synfast(cls_cmb, 512)

map_combined = map_gain*map_cmb
hp.mollview(map_combined)
cls_recon = hp.anafast(map_combined, lmax=len(cls_gain)-1)
plt.figure()
plt.plot(cls_recon - cls_cmb)
print (cls_recon - cls_cmb).mean()
#plt.plot(cls_cmb)
#plt.yscale('log')


The order of the input cl's will change in a future release.
Use new=True keyword to start using the new order.
See documentation of healpy.synalm.
The order of the input cl's will change in a future release.
Use new=True keyword to start using the new order.
See documentation of healpy.synalm.
Sigma is 0.000000 arcmin (0.000000 rad) 
-> fwhm is 0.000000 arcmin
Sigma is 0.000000 arcmin (0.000000 rad) 
-> fwhm is 0.000000 arcmin
1423.33038761

In [46]:
retvalues = cluster.apply_sync(testfunc, range(5), a=a)
retvalues[:]


Out[46]:
[{'a': array([2])}, {'a': array([2])}]

In [795]:
tau = 3.
T = .1
N = 20

def filter(f, tau, T, N):
    ns = np.arange(0, N) - N//2
    return (1./tau)*np.sum(T*np.sinc(T*(f - ns/tau)))
filter = np.vectorize(filter)

fs = np.linspace(-50, 50, 1000)
# fs = np.linspace(-2, 2, 1000)

plot(fs, filter(fs, tau=tau, T=T, N=N), label='spread out')
# plot(fs, N*T*np.sinc(N*T*fs), label='all together')
legend()


Out[795]:
<matplotlib.legend.Legend at 0x11a0c2090>

In [725]:
print A_single
print Omega


1.288225e-06
0.242406840555

In [664]:
reload(util.plot_funcs)


Out[664]:
<module 'util.plot_funcs' from 'util/plot_funcs.py'>

In [57]:
myfname = 'data/lcdm_planck_2013_cl_lensed.dat'
first, _ = myfname.rsplit('.dat')
second, _ = first.rsplit('_lensed', 1)
second.rsplit('_cl', 1)


Out[57]:
['data/lcdm_planck_2013', '']

In [66]:
print lib.cmb.generate_cls(input='README.md', output='junk/explanatory_')
cldict = lib.cmb.get_dls(fname='junk/explanatory_cl.dat', lensed=False)

plot(cldict['ell'], cldict['BB'])



Out[66]:
[<matplotlib.lines.Line2D at 0x10fae6910>]

In [53]:
print lib.cmb.generate_cls(input='explanatory.ini', output='junk/explanatory_')
pdict = lib.cmb.get_lcdm_dls(lensed=False, fname='junk/explanatory_cl.dat')
plot(pdict['ell'], pdict['BB'])
xlim([2, 500])


Running CLASS version v2.3.1
Computing background
 -> age = 13.461689 Gyr
 -> conformal age = 13893.996602 Mpc
Computing thermodynamics with Y_He=0.2487
 -> recombination at z = 1086.807512
    corresponding to conformal time = 278.732613 Mpc
    with comoving sound horizon = 142.256131 Mpc
    angular diameter distance = 12.516244 Mpc
    and sound horizon angle 100*theta_s = 1.044828
 -> baryon drag stops at z = 1064.650406
    corresponding to conformal time = 283.022684 Mpc
    with comoving sound horizon rs = 144.161626 Mpc
 -> reionization with optical depth = 0.084456
    corresponding to conformal time = 4458.025551 Mpc
Computing sources
Computing primordial spectra (analytic spectrum)
No non-linear spectra requested. Nonlinear module skipped.
Computing transfers
Computing unlensed linear spectra
Computing lensed spectra (fast mode)
Writing output files in junk/explanatory_... 

Out[53]:
(2, 500)

In [41]:
reload(calculations.ilc_foreground_cleaning)
import calculations
# retdict2 = calculations.ilc_foreground_cleaning.polarized_ilc_reconstruction(freqs2, _debug=True,
#                  fname='junk/explanatory_cl.dat', lensed=False, regnoise=700.)

In [1151]:
Omega_pixel


Out[1151]:
3.994741635118857e-06

In [1114]:
Omega_single = 7.45446e-6  # sr, from bjohnson report
Omega_array = 9.54171e-3  # sr, from bjohnson report
m = 81  # uK_CMB sqrt(sr)
print m/np.sqrt(Omega_single), "uK per detector beamspot"
print m/np.sqrt(Omega_array), 'uK per array beamspot'

Omega_pixel = hp.nside2pixarea(512)

print m/np.sqrt(Omega_pixel), 'uK per Nside=512 pixel'

regnoise = 3.*np.sqrt(Omega_array/Omega_pixel)  # uK/pixel
print regnoise
# regnoise = 3.*np.sqrt(Omega_single/Omega_array)  # uK/pixel
# regnoise


29667.2249427 uK per detector beamspot
829.224093802 uK per array beamspot
40526.6467473 uK per Nside=512 pixel
146.618919

In [1154]:
print 'array/pixel', 3.*np.sqrt(Omega_array/Omega_pixel)  # uK/pixel
print 'pixel/array', 3.*np.sqrt(Omega_pixel/Omega_array)  # uK/pixel
print 'single/pixel', 3.*np.sqrt(Omega_single/Omega_pixel)  # uK/pixel
print 'pixel/single', 3.*np.sqrt(Omega_pixel/Omega_single)  # uK/pixel


array/pixel 146.618919
pixel/array 0.0613836199407
single/pixel 4.09812311319
pixel/single 2.19612728838

In [20]:
cl0p1dict = lib.cmb.get_cls(lensed=False, fname='junk/r0p1_cl.dat')

In [1225]:
# Load data
from lib.wmap import load_wmap_maps, load_wmap_masks
import lib.ilc

bands = ['K', 'Ka', 'Q', 'V', 'W']
fnames = ['data/wmap_band_smth_imap_r9_9yr_{0}_v5.fits'.format(band)
          for band in bands]
maskdata = 'data/wmap_ilc_rgn_defn_9yr_v5.fits'

print "Loading WMAP data..."
maps = load_wmap_maps(fnames)
masks, region = load_wmap_masks(maskdata)

# Compute weights
print "Computing ILC weights..."
nregions = masks.shape[0]
import time
t0 = time.time()
ws = lib.ilc.compute_ilc_weights(maps, masks)
ws2 = lib.ilc.compute_ilc_weights(maps, masks, noisefactor=True)
tf = time.time()
print "Computed ILC weights for {0} regions in {1} seconds.".format(
    nregions, tf - t0)

# Print weights
print "ILC weights"
linestr = "{0:^10} {1:^10} {2:^10} {3:^10} {4:^10} {5:^10}"
print linestr.format(*(['GROUP'] + bands))
for i in range(len(ws)):
    linelist = [i] + [str(w).ljust(8, '0')[:8] for w in ws[i]]
    print linestr.format(*linelist)
print
print linestr.format(*(['GROUP'] + bands))
for i in range(len(ws)):
    linelist = [i] + [str(w).ljust(8, '0')[:8] for w in ws2[i]]
    print linestr.format(*linelist)


Loading WMAP data...
Computing ILC weights...
Computed ILC weights for 12 regions in 34.8607280254 seconds.
ILC weights
  GROUP        K          Ka         Q          V          W     
    0       0.155986   -0.76841   -0.24397   2.257494   -0.40108 
    1       0.055803   -0.64892   0.208869   1.945169   -0.56091 
    2       0.048720   -0.43987   -0.23540   1.851762   -0.22519 
    3       -0.09964   0.205965   -0.64152   1.560084   -0.02487 
    4       -0.07787   0.094598   -0.42622   0.967506   0.441998 
    5       0.209118   -0.80956   -0.41398   2.482037   -0.46760 
    6       -0.08915   0.169148   -0.51964   1.000231   0.439416 
    7       0.159301   -0.82604   -0.05642   2.157844   -0.43467 
    8       0.201858   -0.23001   -1.69478   3.628966   -0.90602 
    9       -0.09946   -0.10505   -0.05691   1.156837   0.104587 
    10      0.176710   -0.92922   -0.06826   2.746804   -0.92602 
    11      0.237338   -0.85279   -0.59759   2.827484   -0.61443 

  GROUP        K          Ka         Q          V          W     
    0       -0.10318   -0.18066   -0.20556   0.130548   1.358861 
    1       -0.38896   -0.65386   -0.51836   1.406527   1.154672 
    2       -0.14007   -0.25004   -0.27500   0.313937   1.351184 
    3       -0.13936   -0.24997   -0.27583   0.386437   1.278731 
    4       -0.11944   -0.21280   -0.24024   0.259758   1.312741 
    5       -0.18967   -0.33002   -0.34178   0.582569   1.278923 
    6       -0.11695   -0.20821   -0.23559   0.243201   1.317563 
    7       -0.14933   -0.26612   -0.28859   0.414595   1.289470 
    8       -0.15318   -0.27656   -0.29779   0.495749   1.231793 
    9       -0.08794   -0.15301   -0.17776   0.055234   1.363483 
    10      -0.07211   -0.14106   -0.16644   0.209397   1.170225 
    11      -0.11033   -0.19839   -0.22449   0.186716   1.346503 

In [1227]:
# Construct ILC map
print "Constructing ILC map from weights..."
T_hat = lib.ilc.ilc_map_from_weights(maps, ws, region, sigma=None)
T_hat2 = lib.ilc.ilc_map_from_weights(maps, ws2, region, sigma=None)

# Display ILC map
print "Plotting ILC map. Close figure to exit."
hp.mollview(T_hat, title='ILC foreground-cleaned map using WMAP9 data',
            unit='mK', min=-.452123, max=.426628)
hp.mollview(T_hat2, title='ILC foreground-cleaned map using WMAP9 data',
            unit='mK', min=-.452123, max=.426628)
hp.mollview(T_hat2 - T_hat, title='Diff map', unit='mK',
            min=-.452123, max=.426628)


Constructing ILC map from weights...
Plotting ILC map. Close figure to exit.

In [26]:
plt.subplots?

In [204]:
reload(lib.cmb)
reload(lib.ilc)
reload(calculations.ilc_foreground_cleaning)
import calculations

In [74]:
from __future__ import print_function

def many_realizations(freqs, N=100, xxs=['BB'], fname=None, regnoise=0., lensed=False,
                      modcov=False, verbose=False, **kwargs):
    print("Finished: {0} of {1}".format(0, N), end='\r')
    reconfunc = calculations.ilc_foreground_cleaning.polarized_ilc_reconstruction
    cldict = {key: [] for key in xxs}
    cldict['weights_Q'] = []
    cldict['weights_U'] = []
    for i in range(N):
        regendust = True if (i == 0) else False
        retdict = reconfunc(freqs, _debug=False, fname=fname, regnoise=regnoise, 
                            lensed=lensed, verbose=verbose, modcov=modcov,
                            regeneratedust=regendust)
        cl_out = retdict['cl_out']
        for key in xxs:
            cldict[key].append(cl_out[key])
        cldict['weights_Q'].append(retdict['weights_Q'])
        cldict['weights_U'].append(retdict['weights_U'])
        
        print("Finished: {0} of {1}".format(i+1, N), end='\r')
    print("\n")

    for key, value in cldict.items():
        cldict[key] = np.array(value)
        cldict[key + '_mean'] = np.mean(cldict[key], axis=0)
        cldict[key + '_std'] = np.std(cldict[key], axis=0)

    cldict['ell'] = cl_out['ell']
    cldict['regnoise'] = retdict['regnoise']

    return cldict

cl0p1dict = lib.cmb.get_cls(lensed=False, fname='junk/r0p1_cl.dat')


Omega_single = 7.45446e-6  # sr, from bjohnson report
Omega_pixel = hp.nside2pixarea(512)  # pixel solid angle in sr
Omega_array = 9.54171e-3  # sr, from bjohnson report
regnoise = 3./1.3*np.array([1.3, 1.9, 6.7, 110.])*np.sqrt(Omega_single/Omega_pixel)  # uK/pixel
fname = 'junk/explanatory_cl.dat'
lensed = False
modcov = False

cldict_normal_test = calculations.ilc_foreground_cleaning.polarized_ilc_reconstruction(
                        frequencies=np.array([200., 270., 350., 600.])*1.e9, fname=fname, 
                        regnoise=regnoise, lensed=lensed, modcov=modcov, verbose=True)

noisemap = cldict_normal_test['noisemaps'][0][0]
cls_noise = hp.anafast([noisemap, noisemap, noisemap])
labels = ['TT', 'EE', 'BB', 'TE', 'EB', 'TB']
cldict_noise = {labels[i]: cls_noise[i] for i in range(len(labels))}
cldict_noise['ell'] = np.arange(len(cldict_noise['TT']), dtype='float')

prefix_noise = cldict_noise['ell']*(cldict_noise['ell'] + 1)/(2*np.pi)

def myplot(cld, name, label, nf=1.):
    figure()
    ell = cld['ell']
    prefix = ell*(ell+1)/(2*np.pi)
    
    for bb in cld['BB']:
        plot(ell, prefix*bb, 'k-', alpha=.015)

    bbmean = cld['BB_mean']
    bbstd = cld['BB_std']
    errorbar(ell, prefix*bbmean, yerr=prefix*bbstd, label='Output')
    plot(ell, 0*ell, 'g--', label='Input')

    nmax = where(ell<20.)[0][-1]
    # bbmax = prefix[nmax]*(bbmean[nmax] + 5*bbstd[nmax])*2
    bbmax = 0.1 # 0.007
    # ylim(ymax=bbmax)
    ylim([0, bbmax])
    # xlim(xmin=2, xmax=20)
    xlim(xmin=2, xmax=150)
    xscale('log')

    xlabel(r'$\ell$', fontsize=24)
    ylabel(r'$D_\ell^{BB}$ ($\mu K^2$)', fontsize=24)

    ell0p1 = cl0p1dict['ell']
    bb0p1 = cl0p1dict['BB']
    prefix = ell0p1*(ell0p1 + 1)/(2*np.pi)

    # Plot base noise level
    prefix_noise = cldict_noise['ell']*(cldict_noise['ell'] + 1)/(2*np.pi)
    plot(cldict_noise['ell'], prefix_noise*cldict_noise['BB']*nf**2, 'r-', label='Noise')

    # Plot weight-scaled noise level
    wQs = cld['weights_Q']
    wUs = cld['weights_U']
    wavg = (0.5*(wQs + wUs)).mean(axis=0)
    rnfact = cld['regnoise']/cld['regnoise'].min()
#     wfact = np.sum(wavg**2)
#     wfact = np.sum(np.abs(wavg)*rnfact)
#     wfact = np.sqrt(np.sum(wavg**2 * rnfact))
    wfact = np.sqrt(np.sum((wavg*rnfact)**2))*nf**2
    
    plot(cldict_noise['ell'], prefix_noise*cldict_noise['BB']*wfact, 'y-', label='Weighted Noise')
    
    plot(ell0p1, prefix*bb0p1, 'k--', label='r = 0.1')
    legend(loc='upper left', fontsize=20)
    title(label, fontsize=24)
    savefig(name + '.png', dpi=150) 
    f = gcf()
    return f

def save_data(results, fname='results.pickle'):
    copydict = results.copy()
    for key in copydict.keys():
        del copydict[key]['figure']
    with open(fname, 'w') as f:
        pickle.dump(copydict, f)


--------------------------------------------------------------------------------
Performing ILC simulation with frequencies: [ 200.  270.  350.  600.] GHz
Constructing CMB temperature and polarization maps.
Loaded dust from file: /Users/jlazear/github/cmb/data/dust.npy
Combining maps.
Instrument/regularization noise with std dev [   4.09812311    5.98956455   21.12109604  346.76426342] uK
Performing ILC on Q/U maps separately.
Reconstructing ILC Q/U maps.
Computing C_l's from reconstructed maps.
Done with ILC with frequencies [ 200.  270.  350.  600.] GHz!
--------------------------------------------------------------------------------

In [75]:
Omega_single = 7.45446e-6  # sr, from bjohnson report
Omega_pixel = hp.nside2pixarea(512)  # pixel solid angle in sr
noiselist = 3./1.3*np.array([1.3, 1.9, 6.7, 110.])*np.sqrt(Omega_single/Omega_pixel)  # uK/pixel
noisedict = dict(zip([200., 270., 350., 600.], noiselist))
# regnoise = 3.*np.sqrt(Omega_single/Omega_pixel)  # uK/pixel
fname = 'junk/explanatory_cl.dat'
lensed = False    # Include lensing?
noisefactor = 1.  # Noise level relative to PIPER target
modcov = False    # Efstathiou 2009 modified covariance?
verbose = False

baseargdict = {'lensed': lensed,
               'noisefactor': noisefactor,
               'modcov': modcov,
               'fname': fname,
               'verbose': verbose,
               'N': 100}

argdicts = {}
arglist = [('low_2_lownoise', [200., 270.], 0.1, '(200, 270) GHz, Low Noise'),
           #('all_4', [200., 270., 350., 600.], 1., 'All 4 Frequencies, Normal Noise'),
           #('low_high_2', [200., 600.], 1., '(200, 600) GHz, Normal Noise'),
           #('all_4_lownoise', [200., 270., 350., 600.], 0.1, 'All 4 Frequencies, Low Noise'),
           #('200_350', [200., 350.], 1., '(200, 350) GHz, Normal Noise'),
           #('350_600', [350., 600.], 1., '(350, 600) GHz, Normal Noise')
           ]

for i in arglist:
    name, freqs, noisefact, label = i
    regnoise = np.array([noisedict[float(freq)] for freq in freqs])
    freqs = np.array(freqs)*1.e9
    argdict = baseargdict.copy()
    argdict.update({'freqs': freqs, 'regnoise': regnoise*noisefact, 'noisefactor': noisefact})
    toadd = {'label': label, 'args': argdict}
    argdicts[name] = toadd

# ad = argdicts['low_2']  #DELME?
# results = {}

for name, ad in argdicts.items():
    t0 = time.time()
    print("Computing case: {0}".format(name))
    cld = many_realizations(**(ad['args']))
    results[name] = {'cldict': cld}
    f = myplot(cld, name, ad['label'])
    results[name]['figure'] = f
    tf = time.time()
    print("Finished computing ({1} s): {0}".format(name, tf - t0))
    print('-'*30)


Computing case: low_2_lownoise
freq = 200.0 GHz, factor = 0.40291038777
freq = 270.0 GHz, factor = 0.651242138951
Finished: 100 of 100

Finished computing (1238.85599995 s): low_2_lownoise
------------------------------

In [199]:
argdicts = {}
arglist = [('low_2', [200., 270.], 1., '(200, 270) GHz, Normal Noise'),
           ('all_4', [200., 270., 350., 600.], 1., 'All 4 Frequencies, Normal Noise'),
           ('low_high_2', [200., 600.], 1., '(200, 600) GHz, Normal Noise'),
           ('all_4_lownoise', [200., 270., 350., 600.], 0.1, 'All 4 Frequencies, Low Noise')]

for i in arglist:
    name, freqs, noisefact, label = i
    regnoise = np.array([noisedict[float(freq)] for freq in freqs])
    freqs = np.array(freqs)*1.e9
    argdict = baseargdict.copy()
    argdict.update({'freqs': freqs, 'regnoise': regnoise*noisefact, 'noisefactor': noisefact})
    toadd = {'label': label, 'args': argdict}
    argdicts[name] = toadd

for name, ad in argdicts.items():
    t0 = time.time()
#    print("Computing case: {0}".format(name))
#    cld = many_realizations(**(ad['args']))
#    results[name] = {'cldict': cld}
    cld = results[name]['cldict']
    nf = ad['args']['noisefactor']
    f = myplot(cld, name, ad['label'], nf)
    results[name]['figure'] = f
    tf = time.time()
#    print("Finished computing ({1} s): {0}".format(name, tf - t0))
#    print('-'*30)



In [25]:
results.keys()


Out[25]:
['low_2', 'low_high_2', 'all_4']

In [134]:
def save_data(results, fname='results.pickle'):
    copydict = results.copy()
    for key in copydict.keys():
        del copydict[key]['figure']
    with open(fname, 'w') as f:
        pickle.dump(copydict, f)

save_data(results)
        
with open('results.pickle', 'r') as f:
    newresults = pickle.load(f)
    
newresults.k
#results


Out[134]:
{'all_4': {'cldict': {'BB': array([[  0.00000000e+00,   0.00000000e+00,   1.10927848e-04, ...,
             1.28002212e-04,   1.26190556e-04,   1.30112182e-04],
          [  0.00000000e+00,   0.00000000e+00,   6.98910184e-05, ...,
             1.19857821e-04,   1.28158988e-04,   1.25698127e-04],
          [  0.00000000e+00,   0.00000000e+00,   2.14791197e-04, ...,
             1.18761402e-04,   1.20339676e-04,   1.17116768e-04],
          ..., 
          [  0.00000000e+00,   0.00000000e+00,   3.30133397e-04, ...,
             1.20306429e-04,   1.24181432e-04,   1.24175947e-04],
          [  0.00000000e+00,   0.00000000e+00,   3.93900239e-04, ...,
             1.17574703e-04,   1.19235089e-04,   1.22494882e-04],
          [  0.00000000e+00,   0.00000000e+00,   1.75447753e-04, ...,
             1.18423398e-04,   1.21951409e-04,   1.28853495e-04]]),
   'BB_mean': array([ 0.        ,  0.        ,  0.0002477 , ...,  0.00012041,
           0.00012335,  0.00012354]),
   'BB_std': array([  0.00000000e+00,   0.00000000e+00,   1.42030093e-04, ...,
            2.64588973e-06,   2.89649324e-06,   3.23606707e-06]),
   'ell': array([  0.00000000e+00,   1.00000000e+00,   2.00000000e+00, ...,
            1.53300000e+03,   1.53400000e+03,   1.53500000e+03]),
   'regnoise': array([   4.09812311,    5.98956455,   21.12109604,  346.76426342]),
   'weights_Q': array([[ 0.94369616,  0.1284236 , -0.06144784, -0.01067192],
          [ 0.94299502,  0.12937583, -0.06169792, -0.01067293],
          [ 0.94295233,  0.12955392, -0.06183658, -0.01066966],
          [ 0.94231651,  0.13000893, -0.06164558, -0.01067987],
          [ 0.9432238 ,  0.12953728, -0.0621007 , -0.01066038],
          [ 0.94387688,  0.12842127, -0.06163133, -0.01066683],
          [ 0.94253617,  0.12981912, -0.06167766, -0.01067763],
          [ 0.94289872,  0.12944271, -0.06166657, -0.01067487],
          [ 0.94352108,  0.12889681, -0.06175094, -0.01066695],
          [ 0.94240558,  0.13000893, -0.06173668, -0.01067783],
          [ 0.94283699,  0.12947873, -0.06164026, -0.01067546],
          [ 0.94302487,  0.12914996, -0.06149749, -0.01067734],
          [ 0.94361641,  0.12874367, -0.06169264, -0.01066745],
          [ 0.94369842,  0.12859274, -0.06162208, -0.01066908],
          [ 0.94282205,  0.12967328, -0.06182372, -0.01067161],
          [ 0.94386941,  0.12847682, -0.06168133, -0.0106649 ],
          [ 0.94373546,  0.12865356, -0.06172242, -0.0106666 ],
          [ 0.94310735,  0.1293484 , -0.06178601, -0.01066973],
          [ 0.94286007,  0.12950399, -0.06168922, -0.01067484],
          [ 0.94333172,  0.12868317, -0.06133522, -0.01067967],
          [ 0.94315122,  0.12913002, -0.06160716, -0.01067407],
          [ 0.9424938 ,  0.12991902, -0.06173647, -0.01067635],
          [ 0.9428682 ,  0.1295348 , -0.06173033, -0.01067267],
          [ 0.94263917,  0.1297468 , -0.06171051, -0.01067547],
          [ 0.94370886,  0.1286308 , -0.0616717 , -0.01066796],
          [ 0.94359711,  0.12912159, -0.06205966, -0.01065904],
          [ 0.94309755,  0.12924014, -0.0616648 , -0.01067288],
          [ 0.94351075,  0.12868105, -0.06151923, -0.01067257],
          [ 0.94300548,  0.12944118, -0.06177605, -0.01067061],
          [ 0.9424412 ,  0.12978193, -0.06154156, -0.01068157],
          [ 0.94348438,  0.12900674, -0.06182635, -0.01066477],
          [ 0.94311678,  0.12918754, -0.06163089, -0.01067342],
          [ 0.94316127,  0.12959673, -0.0620969 , -0.01066111],
          [ 0.94251748,  0.12970482, -0.06154244, -0.01067987],
          [ 0.94185508,  0.1305407 , -0.0617131 , -0.01068268],
          [ 0.94329   ,  0.12913559, -0.06175742, -0.01066817],
          [ 0.94323496,  0.12898672, -0.06154678, -0.01067489],
          [ 0.94294624,  0.12929754, -0.06156723, -0.01067655],
          [ 0.9425288 ,  0.13001989, -0.06187622, -0.01067247],
          [ 0.94379618,  0.12848498, -0.06161315, -0.01066802],
          [ 0.94284739,  0.12924244, -0.06140824, -0.01068159],
          [ 0.94291453,  0.12922535, -0.06145993, -0.01067995],
          [ 0.94281165,  0.12971071, -0.06185231, -0.01067005],
          [ 0.94295005,  0.12956072, -0.06184152, -0.01066925],
          [ 0.94271564,  0.1299565 , -0.06200489, -0.01066724],
          [ 0.94214646,  0.1304495 , -0.06192221, -0.01067375],
          [ 0.94322565,  0.12928955, -0.06184821, -0.01066699],
          [ 0.94348982,  0.12879177, -0.06161084, -0.01067074],
          [ 0.94334983,  0.12904649, -0.06172769, -0.01066864],
          [ 0.94227492,  0.13035689, -0.06195873, -0.01067307],
          [ 0.9435686 ,  0.1287916 , -0.06169161, -0.01066859],
          [ 0.94249067,  0.1298918 , -0.06170545, -0.01067703],
          [ 0.94259414,  0.12977728, -0.0616951 , -0.01067632],
          [ 0.94271425,  0.12966637, -0.06170556, -0.01067506],
          [ 0.94261881,  0.12978529, -0.0617293 , -0.0106748 ],
          [ 0.94348243,  0.12883593, -0.06164811, -0.01067025],
          [ 0.94330711,  0.12927335, -0.06191615, -0.01066431],
          [ 0.94265092,  0.12977092, -0.06174642, -0.01067541],
          [ 0.94298736,  0.129157  , -0.06146501, -0.01067935],
          [ 0.94352599,  0.12890106, -0.06176005, -0.01066699],
          [ 0.942815  ,  0.12940655, -0.06154302, -0.01067853],
          [ 0.94384354,  0.12871546, -0.0618989 , -0.0106601 ],
          [ 0.94348225,  0.12890734, -0.06172173, -0.01066786],
          [ 0.94296333,  0.12924275, -0.06152902, -0.01067706],
          [ 0.94353367,  0.12894707, -0.06181549, -0.01066525],
          [ 0.94371049,  0.12845512, -0.06149356, -0.01067204],
          [ 0.94249596,  0.12964339, -0.06145636, -0.01068299],
          [ 0.94338979,  0.12888759, -0.06160534, -0.01067204],
          [ 0.9429799 ,  0.12946631, -0.06177504, -0.01067118],
          [ 0.9431876 ,  0.12919124, -0.06170833, -0.01067052],
          [ 0.94273583,  0.12981342, -0.06187785, -0.01067139],
          [ 0.94365461,  0.12901363, -0.06200916, -0.01065907],
          [ 0.94302745,  0.12942993, -0.06178741, -0.01066998],
          [ 0.94427114,  0.12816546, -0.06177688, -0.01065972],
          [ 0.9423755 ,  0.12992832, -0.06162371, -0.01068011],
          [ 0.94402978,  0.12831987, -0.06168579, -0.01066386],
          [ 0.94272859,  0.12914959, -0.06118945, -0.01068873],
          [ 0.94368869,  0.12847064, -0.06148682, -0.01067251],
          [ 0.94315095,  0.12930097, -0.06178307, -0.01066885],
          [ 0.94222644,  0.1302109 , -0.06175996, -0.01067737],
          [ 0.94240692,  0.12977026, -0.06149438, -0.01068279],
          [ 0.94369566,  0.1287186 , -0.0617485 , -0.01066577],
          [ 0.9434407 ,  0.12893552, -0.06170782, -0.01066839],
          [ 0.94339641,  0.12861478, -0.06133275, -0.01067845],
          [ 0.94375215,  0.12838697, -0.0614664 , -0.01067272],
          [ 0.94369988,  0.12877384, -0.06180996, -0.01066376],
          [ 0.94270592,  0.12974966, -0.06178208, -0.0106735 ],
          [ 0.94254676,  0.12968974, -0.06155647, -0.01068003],
          [ 0.94237975,  0.12992359, -0.06162446, -0.01067889],
          [ 0.94310447,  0.12911711, -0.06154644, -0.01067515],
          [ 0.9429183 ,  0.12969882, -0.0619494 , -0.01066772],
          [ 0.94307406,  0.12920302, -0.06160243, -0.01067465],
          [ 0.94419482,  0.1280909 , -0.06162128, -0.01066444],
          [ 0.94173712,  0.13088173, -0.0619407 , -0.01067814],
          [ 0.94411333,  0.12840105, -0.06185487, -0.01065951],
          [ 0.94261015,  0.1296785 , -0.06161035, -0.0106783 ],
          [ 0.94332711,  0.1290201 , -0.06167731, -0.01066991],
          [ 0.94288677,  0.129121  , -0.06132491, -0.01068286],
          [ 0.94348772,  0.12864689, -0.06146047, -0.01067414],
          [ 0.94357091,  0.12878835, -0.06169142, -0.01066785]]),
   'weights_Q_mean': array([ 0.94309775,  0.12926338, -0.06168897, -0.01067216]),
   'weights_Q_std': array([  5.21718440e-04,   5.55970288e-04,   1.68449693e-04,
            6.13305330e-06]),
   'weights_U': array([[ 0.94334432,  0.12903235, -0.0617101 , -0.01066658],
          [ 0.94385932,  0.12835985, -0.06155486, -0.0106643 ],
          [ 0.94342926,  0.12893313, -0.0616963 , -0.01066608],
          [ 0.94287596,  0.12955435, -0.06176126, -0.01066906],
          [ 0.94364979,  0.12887269, -0.06186278, -0.01065971],
          [ 0.94369887,  0.1284091 , -0.06143942, -0.01066855],
          [ 0.94268089,  0.12972598, -0.06173633, -0.01067054],
          [ 0.9430142 ,  0.12932176, -0.06166358, -0.01067238],
          [ 0.94422225,  0.12802757, -0.06158852, -0.01066131],
          [ 0.94179448,  0.13033928, -0.06144898, -0.01068478],
          [ 0.94343062,  0.128907  , -0.06167187, -0.01066575],
          [ 0.94242886,  0.12975473, -0.06150592, -0.01067767],
          [ 0.94383884,  0.12834002, -0.06151277, -0.01066608],
          [ 0.94332639,  0.12907715, -0.06173838, -0.01066515],
          [ 0.94324686,  0.12927672, -0.06186096, -0.01066262],
          [ 0.94279187,  0.12984101, -0.06196864, -0.01066424],
          [ 0.94331296,  0.12878656, -0.06142743, -0.0106721 ],
          [ 0.94276578,  0.12972626, -0.06182544, -0.0106666 ],
          [ 0.9427984 ,  0.12959143, -0.06171993, -0.0106699 ],
          [ 0.94302605,  0.12936542, -0.06172492, -0.01066655],
          [ 0.94302171,  0.12929873, -0.06165017, -0.01067027],
          [ 0.9430724 ,  0.129043  , -0.06143993, -0.01067547],
          [ 0.94282131,  0.12933864, -0.06148366, -0.01067628],
          [ 0.94314152,  0.12907173, -0.06154019, -0.01067307],
          [ 0.94321218,  0.12924228, -0.06179046, -0.01066401],
          [ 0.94254673,  0.12979984, -0.06167402, -0.01067255],
          [ 0.94313899,  0.1296278 , -0.06210958, -0.01065721],
          [ 0.94374657,  0.12872431, -0.06181181, -0.01065907],
          [ 0.94255072,  0.12968524, -0.06156037, -0.01067559],
          [ 0.94252019,  0.12991459, -0.06176536, -0.01066941],
          [ 0.9429026 ,  0.12947628, -0.06170996, -0.01066892],
          [ 0.94386752,  0.12832404, -0.06152488, -0.01066668],
          [ 0.94347381,  0.12890018, -0.06171081, -0.01066318],
          [ 0.942738  ,  0.12952838, -0.06159421, -0.01067217],
          [ 0.94403411,  0.1282466 , -0.06162181, -0.0106589 ],
          [ 0.94307957,  0.12925945, -0.06167043, -0.01066858],
          [ 0.94342346,  0.12889705, -0.06165482, -0.01066568],
          [ 0.94373468,  0.12838075, -0.0614474 , -0.01066803],
          [ 0.94268985,  0.12973948, -0.06175839, -0.01067094],
          [ 0.94303188,  0.12900123, -0.06135624, -0.01067687],
          [ 0.94248596,  0.12982592, -0.06163752, -0.01067436],
          [ 0.94254437,  0.12981286, -0.0616851 , -0.01067213],
          [ 0.94297639,  0.12933018, -0.06163565, -0.01067092],
          [ 0.94243212,  0.12985671, -0.06161344, -0.0106754 ],
          [ 0.94335407,  0.12879905, -0.06148316, -0.01066996],
          [ 0.94277666,  0.12982799, -0.06194041, -0.01066424],
          [ 0.94246199,  0.1298285 , -0.0616153 , -0.01067519],
          [ 0.94310904,  0.12946382, -0.06190897, -0.0106639 ],
          [ 0.94299323,  0.12949499, -0.06182358, -0.01066465],
          [ 0.94319304,  0.12930964, -0.06183838, -0.0106643 ],
          [ 0.94352389,  0.12883394, -0.06169313, -0.0106647 ],
          [ 0.94316464,  0.12919088, -0.06168838, -0.01066714],
          [ 0.94209819,  0.13050458, -0.06193215, -0.01067062],
          [ 0.94371242,  0.12827415, -0.06131483, -0.01067173],
          [ 0.94307401,  0.12928762, -0.06169411, -0.01066753],
          [ 0.94298539,  0.12956282, -0.06188471, -0.01066349],
          [ 0.94300999,  0.12959073, -0.06193815, -0.01066256],
          [ 0.94238988,  0.12981025, -0.06152172, -0.01067842],
          [ 0.94290176,  0.12938983, -0.06162151, -0.01067008],
          [ 0.94282548,  0.12933512, -0.06148499, -0.01067561],
          [ 0.94282578,  0.12970932, -0.0618692 , -0.0106659 ],
          [ 0.9427736 ,  0.12960299, -0.06170794, -0.01066865],
          [ 0.94282814,  0.12937116, -0.06152596, -0.01067334],
          [ 0.94263854,  0.13002972, -0.0620037 , -0.01066455],
          [ 0.94221677,  0.13025566, -0.06180085, -0.01067158],
          [ 0.94363997,  0.12868094, -0.06165684, -0.01066407],
          [ 0.94285902,  0.12942424, -0.06161062, -0.01067265],
          [ 0.94289809,  0.12964572, -0.06187836, -0.01066545],
          [ 0.94256999,  0.12985874, -0.0617568 , -0.01067193],
          [ 0.94249585,  0.12985705, -0.06167887, -0.01067404],
          [ 0.94278787,  0.12974511, -0.06186656, -0.01066642],
          [ 0.94294236,  0.12942188, -0.06169573, -0.01066852],
          [ 0.94236208,  0.12986235, -0.06154793, -0.0106765 ],
          [ 0.94270951,  0.12951894, -0.06155305, -0.0106754 ],
          [ 0.94383198,  0.12853094, -0.0617036 , -0.01065932],
          [ 0.94323064,  0.12907047, -0.06163177, -0.01066934],
          [ 0.94239426,  0.13029721, -0.06202444, -0.01066703],
          [ 0.94281834,  0.12921968, -0.06135776, -0.01068027],
          [ 0.94287384,  0.12946427, -0.06166761, -0.0106705 ],
          [ 0.94325568,  0.12893644, -0.06152121, -0.0106709 ],
          [ 0.94267139,  0.12972375, -0.06172373, -0.01067142],
          [ 0.94256955,  0.1300684 , -0.06197318, -0.01066477],
          [ 0.9431137 ,  0.12906342, -0.0615055 , -0.01067162],
          [ 0.942461  ,  0.12975674, -0.06154166, -0.01067607],
          [ 0.94312282,  0.12918372, -0.06163678, -0.01066976],
          [ 0.94241058,  0.13012249, -0.06186322, -0.01066985],
          [ 0.94259721,  0.12977196, -0.06169917, -0.01067   ],
          [ 0.94264211,  0.12976504, -0.06173562, -0.01067152],
          [ 0.94305406,  0.12919467, -0.061578  , -0.01067073],
          [ 0.94238663,  0.13006091, -0.06177519, -0.01067235],
          [ 0.94311585,  0.12923864, -0.06168794, -0.01066656],
          [ 0.94272599,  0.12963122, -0.06168625, -0.01067097],
          [ 0.94176682,  0.13037028, -0.06145094, -0.01068616],
          [ 0.94334083,  0.12898298, -0.06165737, -0.01066644],
          [ 0.9430238 ,  0.12971767, -0.06208326, -0.01065822],
          [ 0.9430733 ,  0.1292248 , -0.06162946, -0.01066863],
          [ 0.94195056,  0.13043232, -0.06170442, -0.01067846],
          [ 0.94418391,  0.12837097, -0.06190286, -0.01065202],
          [ 0.94405633,  0.12823821, -0.061633  , -0.01066154],
          [ 0.94307054,  0.12934115, -0.06174566, -0.01066603]]),
   'weights_U_mean': array([ 0.94298586,  0.12936834, -0.0616852 , -0.01066899]),
   'weights_U_std': array([  4.97134444e-04,   5.41480643e-04,   1.61504763e-04,
            5.65606825e-06])}},
 'low_2': {'cldict': {'BB': array([[ 0.        ,  0.        ,  0.01238223, ...,  0.00059811,
            0.00063392,  0.00060923],
          [ 0.        ,  0.        ,  0.0043329 , ...,  0.00062302,
            0.00060765,  0.00065675],
          [ 0.        ,  0.        ,  0.00958019, ...,  0.00062638,
            0.00061814,  0.00060912],
          ..., 
          [ 0.        ,  0.        ,  0.00613781, ...,  0.00060467,
            0.00062157,  0.00063638],
          [ 0.        ,  0.        ,  0.0088833 , ...,  0.00060083,
            0.00064474,  0.00063356],
          [ 0.        ,  0.        ,  0.011994  , ...,  0.00061703,
            0.00062863,  0.00062524]]),
   'BB_mean': array([ 0.        ,  0.        ,  0.00867021, ...,  0.00061485,
           0.00063212,  0.00063016]),
   'BB_std': array([  0.00000000e+00,   0.00000000e+00,   2.16748380e-03, ...,
            1.51237205e-05,   1.69248680e-05,   1.56445324e-05]),
   'ell': array([  0.00000000e+00,   1.00000000e+00,   2.00000000e+00, ...,
            1.53300000e+03,   1.53400000e+03,   1.53500000e+03]),
   'regnoise': array([ 4.09812311,  5.98956455]),
   'weights_Q': array([[ 2.34698757, -1.34698757],
          [ 2.34693675, -1.34693675],
          [ 2.34689688, -1.34689688],
          [ 2.34698667, -1.34698667],
          [ 2.34699983, -1.34699983],
          [ 2.34680768, -1.34680768],
          [ 2.3470197 , -1.3470197 ],
          [ 2.34692822, -1.34692822],
          [ 2.34692292, -1.34692292],
          [ 2.34686445, -1.34686445],
          [ 2.34696179, -1.34696179],
          [ 2.34688741, -1.34688741],
          [ 2.34697348, -1.34697348],
          [ 2.34697376, -1.34697376],
          [ 2.34697171, -1.34697171],
          [ 2.34692167, -1.34692167],
          [ 2.3469758 , -1.3469758 ],
          [ 2.3469158 , -1.3469158 ],
          [ 2.34680869, -1.34680869],
          [ 2.34690168, -1.34690168],
          [ 2.34694906, -1.34694906],
          [ 2.34700865, -1.34700865],
          [ 2.34692865, -1.34692865],
          [ 2.34695844, -1.34695844],
          [ 2.34695755, -1.34695755],
          [ 2.34691871, -1.34691871],
          [ 2.34694363, -1.34694363],
          [ 2.34696298, -1.34696298],
          [ 2.34702692, -1.34702692],
          [ 2.34685686, -1.34685686],
          [ 2.34688632, -1.34688632],
          [ 2.34696636, -1.34696636],
          [ 2.34698517, -1.34698517],
          [ 2.34698388, -1.34698388],
          [ 2.34702904, -1.34702904],
          [ 2.34695427, -1.34695427],
          [ 2.34688389, -1.34688389],
          [ 2.34695924, -1.34695924],
          [ 2.34694039, -1.34694039],
          [ 2.34684067, -1.34684067],
          [ 2.34707517, -1.34707517],
          [ 2.34689034, -1.34689034],
          [ 2.34688945, -1.34688945],
          [ 2.34684408, -1.34684408],
          [ 2.34694717, -1.34694717],
          [ 2.3468558 , -1.3468558 ],
          [ 2.34694496, -1.34694496],
          [ 2.34691363, -1.34691363],
          [ 2.3469017 , -1.3469017 ],
          [ 2.34695044, -1.34695044],
          [ 2.34687456, -1.34687456],
          [ 2.34691865, -1.34691865],
          [ 2.3469287 , -1.3469287 ],
          [ 2.3469155 , -1.3469155 ],
          [ 2.34684457, -1.34684457],
          [ 2.34690975, -1.34690975],
          [ 2.34685509, -1.34685509],
          [ 2.34690611, -1.34690611],
          [ 2.34693989, -1.34693989],
          [ 2.34691364, -1.34691364],
          [ 2.34696981, -1.34696981],
          [ 2.34695942, -1.34695942],
          [ 2.34703367, -1.34703367],
          [ 2.34685939, -1.34685939],
          [ 2.34686221, -1.34686221],
          [ 2.34690287, -1.34690287],
          [ 2.34694532, -1.34694532],
          [ 2.34686389, -1.34686389],
          [ 2.34683016, -1.34683016],
          [ 2.3469235 , -1.3469235 ],
          [ 2.34693362, -1.34693362],
          [ 2.34689598, -1.34689598],
          [ 2.34695249, -1.34695249],
          [ 2.3469642 , -1.3469642 ],
          [ 2.34680722, -1.34680722],
          [ 2.34696869, -1.34696869],
          [ 2.3469361 , -1.3469361 ],
          [ 2.3469341 , -1.3469341 ],
          [ 2.34694253, -1.34694253],
          [ 2.34690549, -1.34690549],
          [ 2.34689136, -1.34689136],
          [ 2.34692378, -1.34692378],
          [ 2.34694352, -1.34694352],
          [ 2.34685968, -1.34685968],
          [ 2.34695315, -1.34695315],
          [ 2.34690483, -1.34690483],
          [ 2.34697503, -1.34697503],
          [ 2.34685538, -1.34685538],
          [ 2.34698928, -1.34698928],
          [ 2.34694504, -1.34694504],
          [ 2.3468014 , -1.3468014 ],
          [ 2.34702666, -1.34702666],
          [ 2.34696827, -1.34696827],
          [ 2.34694543, -1.34694543],
          [ 2.34690281, -1.34690281],
          [ 2.34697157, -1.34697157],
          [ 2.34689627, -1.34689627],
          [ 2.34702029, -1.34702029],
          [ 2.34693297, -1.34693297],
          [ 2.34694568, -1.34694568]]),
   'weights_Q_mean': array([ 2.34692845, -1.34692845]),
   'weights_Q_std': array([  5.42960748e-05,   5.42960748e-05]),
   'weights_U': array([[ 2.3430214 , -1.3430214 ],
          [ 2.34316641, -1.34316641],
          [ 2.3429486 , -1.3429486 ],
          [ 2.34327939, -1.34327939],
          [ 2.34324946, -1.34324946],
          [ 2.3432683 , -1.3432683 ],
          [ 2.34314484, -1.34314484],
          [ 2.34325397, -1.34325397],
          [ 2.34297659, -1.34297659],
          [ 2.3431606 , -1.3431606 ],
          [ 2.34319846, -1.34319846],
          [ 2.34306596, -1.34306596],
          [ 2.34293181, -1.34293181],
          [ 2.34295199, -1.34295199],
          [ 2.34316348, -1.34316348],
          [ 2.34318176, -1.34318176],
          [ 2.34304264, -1.34304264],
          [ 2.34314739, -1.34314739],
          [ 2.34295455, -1.34295455],
          [ 2.34334379, -1.34334379],
          [ 2.34324696, -1.34324696],
          [ 2.34294045, -1.34294045],
          [ 2.34307137, -1.34307137],
          [ 2.34293981, -1.34293981],
          [ 2.34330345, -1.34330345],
          [ 2.34326837, -1.34326837],
          [ 2.3429323 , -1.3429323 ],
          [ 2.34311542, -1.34311542],
          [ 2.34325552, -1.34325552],
          [ 2.34289921, -1.34289921],
          [ 2.34309324, -1.34309324],
          [ 2.34296269, -1.34296269],
          [ 2.34315996, -1.34315996],
          [ 2.34306122, -1.34306122],
          [ 2.3429786 , -1.3429786 ],
          [ 2.34313106, -1.34313106],
          [ 2.34305198, -1.34305198],
          [ 2.34320704, -1.34320704],
          [ 2.34306278, -1.34306278],
          [ 2.34297594, -1.34297594],
          [ 2.34316853, -1.34316853],
          [ 2.34324895, -1.34324895],
          [ 2.34303978, -1.34303978],
          [ 2.34289738, -1.34289738],
          [ 2.34324258, -1.34324258],
          [ 2.34315576, -1.34315576],
          [ 2.34331913, -1.34331913],
          [ 2.3431167 , -1.3431167 ],
          [ 2.34295293, -1.34295293],
          [ 2.34312295, -1.34312295],
          [ 2.3430468 , -1.3430468 ],
          [ 2.34303489, -1.34303489],
          [ 2.34302272, -1.34302272],
          [ 2.34318278, -1.34318278],
          [ 2.34316843, -1.34316843],
          [ 2.34328212, -1.34328212],
          [ 2.34306281, -1.34306281],
          [ 2.34315034, -1.34315034],
          [ 2.34308087, -1.34308087],
          [ 2.34316486, -1.34316486],
          [ 2.34325885, -1.34325885],
          [ 2.34324488, -1.34324488],
          [ 2.3431462 , -1.3431462 ],
          [ 2.34310915, -1.34310915],
          [ 2.34311353, -1.34311353],
          [ 2.34308004, -1.34308004],
          [ 2.34321874, -1.34321874],
          [ 2.34304742, -1.34304742],
          [ 2.34311126, -1.34311126],
          [ 2.34300953, -1.34300953],
          [ 2.34301009, -1.34301009],
          [ 2.34319976, -1.34319976],
          [ 2.34307215, -1.34307215],
          [ 2.34311881, -1.34311881],
          [ 2.3431589 , -1.3431589 ],
          [ 2.34312956, -1.34312956],
          [ 2.34311746, -1.34311746],
          [ 2.34307319, -1.34307319],
          [ 2.34315228, -1.34315228],
          [ 2.34303394, -1.34303394],
          [ 2.34309561, -1.34309561],
          [ 2.34312965, -1.34312965],
          [ 2.34315072, -1.34315072],
          [ 2.34295273, -1.34295273],
          [ 2.3433356 , -1.3433356 ],
          [ 2.34296175, -1.34296175],
          [ 2.34326563, -1.34326563],
          [ 2.34308535, -1.34308535],
          [ 2.34320403, -1.34320403],
          [ 2.34323621, -1.34323621],
          [ 2.34319075, -1.34319075],
          [ 2.34308608, -1.34308608],
          [ 2.34314402, -1.34314402],
          [ 2.34307337, -1.34307337],
          [ 2.34306195, -1.34306195],
          [ 2.34320964, -1.34320964],
          [ 2.34315747, -1.34315747],
          [ 2.34307887, -1.34307887],
          [ 2.34314339, -1.34314339],
          [ 2.34307671, -1.34307671]]),
   'weights_U_mean': array([ 2.34311615, -1.34311615]),
   'weights_U_std': array([ 0.00010603,  0.00010603])}},
 'low_high_2': {'cldict': {'BB': array([[ 0.        ,  0.        ,  0.00046466, ...,  0.00013153,
            0.00012928,  0.00013465],
          [ 0.        ,  0.        ,  0.00025477, ...,  0.000126  ,
            0.00013542,  0.00013684],
          [ 0.        ,  0.        ,  0.00024039, ...,  0.00012919,
            0.00013981,  0.00012934],
          ..., 
          [ 0.        ,  0.        ,  0.00032159, ...,  0.00012989,
            0.0001335 ,  0.00013197],
          [ 0.        ,  0.        ,  0.00014186, ...,  0.00013046,
            0.00012708,  0.00013407],
          [ 0.        ,  0.        ,  0.00027527, ...,  0.00013379,
            0.00012917,  0.00012842]]),
   'BB_mean': array([ 0.        ,  0.        ,  0.0003022 , ...,  0.00013063,
           0.00013388,  0.00013275]),
   'BB_std': array([  0.00000000e+00,   0.00000000e+00,   1.44203832e-04, ...,
            3.16416383e-06,   3.49995776e-06,   2.96253581e-06]),
   'ell': array([  0.00000000e+00,   1.00000000e+00,   2.00000000e+00, ...,
            1.53300000e+03,   1.53400000e+03,   1.53500000e+03]),
   'regnoise': array([   4.09812311,  346.76426342]),
   'weights_Q': array([[ 1.01161664, -0.01161664],
          [ 1.01161653, -0.01161653],
          [ 1.01161658, -0.01161658],
          [ 1.01161653, -0.01161653],
          [ 1.01161663, -0.01161663],
          [ 1.01161696, -0.01161696],
          [ 1.01161693, -0.01161693],
          [ 1.01161637, -0.01161637],
          [ 1.01161623, -0.01161623],
          [ 1.01161709, -0.01161709],
          [ 1.01161689, -0.01161689],
          [ 1.01161668, -0.01161668],
          [ 1.01161678, -0.01161678],
          [ 1.01161688, -0.01161688],
          [ 1.01161709, -0.01161709],
          [ 1.01161679, -0.01161679],
          [ 1.01161651, -0.01161651],
          [ 1.0116173 , -0.0116173 ],
          [ 1.01161667, -0.01161667],
          [ 1.01161672, -0.01161672],
          [ 1.01161673, -0.01161673],
          [ 1.01161628, -0.01161628],
          [ 1.01161684, -0.01161684],
          [ 1.01161681, -0.01161681],
          [ 1.01161697, -0.01161697],
          [ 1.01161681, -0.01161681],
          [ 1.01161617, -0.01161617],
          [ 1.01161751, -0.01161751],
          [ 1.01161646, -0.01161646],
          [ 1.01161699, -0.01161699],
          [ 1.01161739, -0.01161739],
          [ 1.01161768, -0.01161768],
          [ 1.01161658, -0.01161658],
          [ 1.01161705, -0.01161705],
          [ 1.01161697, -0.01161697],
          [ 1.01161665, -0.01161665],
          [ 1.01161685, -0.01161685],
          [ 1.01161631, -0.01161631],
          [ 1.01161679, -0.01161679],
          [ 1.01161725, -0.01161725],
          [ 1.01161649, -0.01161649],
          [ 1.0116168 , -0.0116168 ],
          [ 1.0116161 , -0.0116161 ],
          [ 1.01161641, -0.01161641],
          [ 1.01161697, -0.01161697],
          [ 1.01161727, -0.01161727],
          [ 1.0116167 , -0.0116167 ],
          [ 1.01161674, -0.01161674],
          [ 1.0116159 , -0.0116159 ],
          [ 1.01161602, -0.01161602],
          [ 1.01161628, -0.01161628],
          [ 1.01161635, -0.01161635],
          [ 1.01161704, -0.01161704],
          [ 1.01161569, -0.01161569],
          [ 1.01161699, -0.01161699],
          [ 1.011616  , -0.011616  ],
          [ 1.01161688, -0.01161688],
          [ 1.01161714, -0.01161714],
          [ 1.01161575, -0.01161575],
          [ 1.01161613, -0.01161613],
          [ 1.01161694, -0.01161694],
          [ 1.01161659, -0.01161659],
          [ 1.01161668, -0.01161668],
          [ 1.01161642, -0.01161642],
          [ 1.01161699, -0.01161699],
          [ 1.01161667, -0.01161667],
          [ 1.01161648, -0.01161648],
          [ 1.01161706, -0.01161706],
          [ 1.01161767, -0.01161767],
          [ 1.01161657, -0.01161657],
          [ 1.01161643, -0.01161643],
          [ 1.01161671, -0.01161671],
          [ 1.01161678, -0.01161678],
          [ 1.01161699, -0.01161699],
          [ 1.01161691, -0.01161691],
          [ 1.01161621, -0.01161621],
          [ 1.0116168 , -0.0116168 ],
          [ 1.01161645, -0.01161645],
          [ 1.0116162 , -0.0116162 ],
          [ 1.01161687, -0.01161687],
          [ 1.01161689, -0.01161689],
          [ 1.01161709, -0.01161709],
          [ 1.01161687, -0.01161687],
          [ 1.01161704, -0.01161704],
          [ 1.01161701, -0.01161701],
          [ 1.01161651, -0.01161651],
          [ 1.01161682, -0.01161682],
          [ 1.0116161 , -0.0116161 ],
          [ 1.01161688, -0.01161688],
          [ 1.01161725, -0.01161725],
          [ 1.01161635, -0.01161635],
          [ 1.01161756, -0.01161756],
          [ 1.01161681, -0.01161681],
          [ 1.01161719, -0.01161719],
          [ 1.01161665, -0.01161665],
          [ 1.0116168 , -0.0116168 ],
          [ 1.01161633, -0.01161633],
          [ 1.01161567, -0.01161567],
          [ 1.01161715, -0.01161715],
          [ 1.01161618, -0.01161618]]),
   'weights_Q_mean': array([ 1.01161671, -0.01161671]),
   'weights_Q_std': array([  4.05013725e-07,   4.05013725e-07]),
   'weights_U': array([[ 1.01161397, -0.01161397],
          [ 1.01161222, -0.01161222],
          [ 1.0116121 , -0.0116121 ],
          [ 1.01161281, -0.01161281],
          [ 1.01161161, -0.01161161],
          [ 1.01161032, -0.01161032],
          [ 1.01161193, -0.01161193],
          [ 1.01161177, -0.01161177],
          [ 1.01161297, -0.01161297],
          [ 1.01161312, -0.01161312],
          [ 1.01161243, -0.01161243],
          [ 1.0116115 , -0.0116115 ],
          [ 1.01161115, -0.01161115],
          [ 1.01161264, -0.01161264],
          [ 1.01161323, -0.01161323],
          [ 1.01161307, -0.01161307],
          [ 1.01161154, -0.01161154],
          [ 1.01161263, -0.01161263],
          [ 1.01161215, -0.01161215],
          [ 1.01161246, -0.01161246],
          [ 1.01161126, -0.01161126],
          [ 1.01161198, -0.01161198],
          [ 1.01161232, -0.01161232],
          [ 1.01161263, -0.01161263],
          [ 1.01161296, -0.01161296],
          [ 1.01161114, -0.01161114],
          [ 1.01161231, -0.01161231],
          [ 1.01161054, -0.01161054],
          [ 1.01161102, -0.01161102],
          [ 1.01161227, -0.01161227],
          [ 1.01161051, -0.01161051],
          [ 1.01161207, -0.01161207],
          [ 1.01161383, -0.01161383],
          [ 1.01161233, -0.01161233],
          [ 1.01161218, -0.01161218],
          [ 1.01161162, -0.01161162],
          [ 1.0116098 , -0.0116098 ],
          [ 1.01161215, -0.01161215],
          [ 1.01161273, -0.01161273],
          [ 1.01161176, -0.01161176],
          [ 1.01161181, -0.01161181],
          [ 1.01161302, -0.01161302],
          [ 1.01161279, -0.01161279],
          [ 1.01161219, -0.01161219],
          [ 1.01161187, -0.01161187],
          [ 1.01161421, -0.01161421],
          [ 1.01161244, -0.01161244],
          [ 1.01161256, -0.01161256],
          [ 1.01161369, -0.01161369],
          [ 1.01161311, -0.01161311],
          [ 1.01161297, -0.01161297],
          [ 1.01161272, -0.01161272],
          [ 1.01161337, -0.01161337],
          [ 1.01161152, -0.01161152],
          [ 1.01161176, -0.01161176],
          [ 1.01161291, -0.01161291],
          [ 1.01161245, -0.01161245],
          [ 1.01161263, -0.01161263],
          [ 1.01161294, -0.01161294],
          [ 1.01161289, -0.01161289],
          [ 1.01161321, -0.01161321],
          [ 1.01161194, -0.01161194],
          [ 1.01161252, -0.01161252],
          [ 1.01161254, -0.01161254],
          [ 1.01161148, -0.01161148],
          [ 1.01161179, -0.01161179],
          [ 1.01161299, -0.01161299],
          [ 1.01161155, -0.01161155],
          [ 1.01161258, -0.01161258],
          [ 1.01161195, -0.01161195],
          [ 1.01161201, -0.01161201],
          [ 1.01161232, -0.01161232],
          [ 1.01161174, -0.01161174],
          [ 1.01161162, -0.01161162],
          [ 1.01161281, -0.01161281],
          [ 1.01161382, -0.01161382],
          [ 1.01161083, -0.01161083],
          [ 1.01161236, -0.01161236],
          [ 1.01161167, -0.01161167],
          [ 1.01161189, -0.01161189],
          [ 1.01161168, -0.01161168],
          [ 1.01161125, -0.01161125],
          [ 1.01161109, -0.01161109],
          [ 1.01161263, -0.01161263],
          [ 1.01161196, -0.01161196],
          [ 1.01161054, -0.01161054],
          [ 1.0116109 , -0.0116109 ],
          [ 1.01161286, -0.01161286],
          [ 1.0116121 , -0.0116121 ],
          [ 1.01161204, -0.01161204],
          [ 1.01161198, -0.01161198],
          [ 1.0116094 , -0.0116094 ],
          [ 1.01161347, -0.01161347],
          [ 1.01161212, -0.01161212],
          [ 1.01161144, -0.01161144],
          [ 1.01161125, -0.01161125],
          [ 1.01161278, -0.01161278],
          [ 1.01161353, -0.01161353],
          [ 1.01161173, -0.01161173],
          [ 1.01161113, -0.01161113]]),
   'weights_U_mean': array([ 1.01161216, -0.01161216]),
   'weights_U_std': array([  8.76344367e-07,   8.76344367e-07])}}}

In [184]:
#nms = cldict_normal_test['noisemaps'][:, 0, :]
#noisecls = np.array([hp.anafast(nm) for nm in nms])
print(noisecls.mean(axis=1)/noisecls.mean(axis=1)[0])
print((regnoise/regnoise[0])**2)


[  1.00000000e+00   2.13424501e+00   2.65486713e+01   7.17723346e+03]
[  1.00000000e+00   2.13609467e+00   2.65621302e+01   7.15976331e+03]

In [160]:
noisecls = hp.anafast(cldict_normal_test['noisemap'][0])
noiseells = np.arange(len(noisecls))
noiseprefix = noiseells*(noiseells + 1)/(2*np.pi)
plot(noiseells, noiseprefix*noisecls)
xlim(xmax=20)
ylim(ymax=10)


Out[160]:
(0.0, 10)

In [146]:
f = results['low_2']['figure']
ax = f.get_axes()[0]
ax.set_ylim(ymax=10)
f


Out[146]:

In [25]:
N=100
freqs = np.array([200., 270.])*1.e9

Omega_single = 7.45446e-6  # sr, from bjohnson report
Omega_pixel = hp.nside2pixarea(512)  # pixel solid angle in sr
Omega_array = 9.54171e-3  # sr, from bjohnson report
regnoise = 3.*np.sqrt(Omega_single/Omega_pixel)  # uK/pixel
fname = 'junk/explanatory_cl.dat'
lensed = False
noisefactor = False

cldict_normal_test = calculations.ilc_foreground_cleaning.polarized_ilc_reconstruction(freqs, fname=fname, regnoise=regnoise,
                                      lensed=lensed, noisefactor=noisefactor)


--------------------------------------------------------------------------------
Performing ILC simulation with frequencies: [ 200.  270.] GHz
Constructing CMB temperature and polarization maps.
Loaded dust from file: /Users/jlazear/github/cmb/data/dust.npy
Combining maps.
Instrument/regularization noise with std dev 4.09812311319 uK
Performing ILC on Q/U maps separately.
Reconstructing ILC Q/U maps.
Computing C_l's from reconstructed maps.
Done with ILC with frequencies [ 200.  270.] GHz!
--------------------------------------------------------------------------------

In [131]:
noisemap = cldict_normal_test['noisemap']
cls_noise_low = hp.anafast([noisemap, noisemap, noisemap])
labels = ['TT', 'EE', 'BB', 'TE', 'EB', 'TB']
cldict_noise_low = {labels[i]: cls_noise_low[i] for i in range(len(labels))}
cldict_noise_low['ell'] = np.arange(len(cldict_noise_low['TT']), dtype='float')

prefix_noise_low = cldict_noise_low['ell']*(cldict_noise_low['ell'] + 1)/(2*np.pi)
plot(cldict_noise_low['ell'], prefix_noise_low*cldict_noise_low['BB'])


Out[131]:
[<matplotlib.lines.Line2D at 0x1116aea90>]

In [29]:
noisemap = cldict_normal_test['noisemap'][0]
cls_noise = hp.anafast([noisemap, noisemap, noisemap])
labels = ['TT', 'EE', 'BB', 'TE', 'EB', 'TB']
cldict_noise = {labels[i]: cls_noise[i] for i in range(len(labels))}
cldict_noise['ell'] = np.arange(len(cldict_noise['TT']), dtype='float')

prefix_noise = cldict_noise['ell']*(cldict_noise['ell'] + 1)/(2*np.pi)
plot(cldict_noise['ell'], prefix_noise*cldict_noise['BB'])


Out[29]:
[<matplotlib.lines.Line2D at 0x11495d5d0>]

In [59]:
cld = cldict_low
ell = cld['ell']
prefix = ell*(ell+1)/(2*np.pi)
# for bb in cld['BB']:
#     plot(ell, prefix*bb, 'k-', alpha=.015)

# plot(ell, prefix*bbmean, 'b-')
# errorbar(ell, prefix*bbmean, fmt=None, yerr=prefix*bbstd, color='b')
bbmean = cld['BB_mean']
bbstd = cld['BB_std']
errorbar(ell, prefix*bbmean, yerr=prefix*bbstd, label='Output')
plot(ell, 0*ell, 'g--', label='Input')

nmax = where(ell<20.)[0][-1]
# bbmax = prefix[nmax]*(bbmean[nmax] + 5*bbstd[nmax])*2
bbmax = 0.007
# ylim(ymax=bbmax)
ylim([0, bbmax])
# xlim(xmin=2, xmax=20)
xlim(xmin=2, xmax=150)
xscale('log')

xlabel(r'$\ell$', fontsize=24)
ylabel(r'$D_\ell^{BB}$ ($\mu K^2$)', fontsize=24)

ell0p1 = cl0p1dict['ell']
bb0p1 = cl0p1dict['BB']
prefix = ell0p1*(ell0p1 + 1)/(2*np.pi)

prefix_noise_low = cldict_noise_low['ell']*(cldict_noise_low['ell'] + 1)/(2*np.pi)
plot(cldict_noise_low['ell'], prefix_noise_low*cldict_noise_low['BB'], 'r-', label='Noise')
plot(ell0p1, prefix*bb0p1, 'k--', label='r = 0.1')
legend(loc='upper left', fontsize=20)
title('Pixel Noise = PIPER/10 = 0.4 uK/pixel')
savefig('low_noise.png', dpi=150)


---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-59-913462cb726a> in <module>()
     28 prefix = ell0p1*(ell0p1 + 1)/(2*np.pi)
     29 
---> 30 prefix_noise_low = cldict_noise_low['ell']*(cldict_noise_low['ell'] + 1)/(2*np.pi)
     31 plot(cldict_noise_low['ell'], prefix_noise_low*cldict_noise_low['BB'], 'r-', label='Noise')
     32 plot(ell0p1, prefix*bb0p1, 'k--', label='r = 0.1')

NameError: name 'cldict_noise_low' is not defined

In [89]:
nf = np.sqrt(np.sum(np.abs(cldict_normal_test['weights_Q'])**2))

In [90]:
cld = cldict_normal_test['cl_out']
ell = cld['ell']
prefix = ell*(ell+1)/(2*np.pi)
# for bb in cld['BB']:
plot(ell, prefix*bb, 'k-', alpha=1)

# plot(ell, prefix*bbmean, 'b-')
# errorbar(ell, prefix*bbmean, fmt=None, yerr=prefix*bbstd, color='b')
# bbmean = cld['BB_mean']
# bbstd = cld['BB_std']
# errorbar(ell, prefix*bbmean, yerr=prefix*bbstd, label='Output')
plot(ell, 0*ell, 'g--', label='Input')

prefix_noise = cldict_noise['ell']*(cldict_noise['ell'] + 1)/(2*np.pi)
plot(cldict_noise['ell'], nf*prefix_noise*cldict_noise['BB'], 'r-', label='Noise')

nmax = where(ell<20.)[0][-1]
# bbmax = prefix[nmax]*(bbmean[nmax] + 5*bbstd[nmax])*2
bbmax = 0.007
# ylim(ymax=bbmax)
ylim([0, bbmax])
# xlim(xmin=2, xmax=20)
xlim(xmin=2, xmax=150)
xscale('log')

xlabel(r'$\ell$', fontsize=24)
ylabel(r'$D_\ell^{BB}$ ($\mu K^2$)', fontsize=24)

ell0p1 = cl0p1dict['ell']
bb0p1 = cl0p1dict['BB']
prefix = ell0p1*(ell0p1 + 1)/(2*np.pi)

plot(ell0p1, prefix*bb0p1, 'k--', label='r = 0.1')
legend(loc='upper left', fontsize=20)
title('Pixel Noise = PIPER = 4 uK/pixel')
# savefig('piper_noise.png', dpi=150)


Out[90]:
<matplotlib.text.Text at 0x11f06f650>

In [104]:
cld = cldict_normal_mod
ell = cld['ell']
prefix = ell*(ell+1)/(2*np.pi)
for bb in cld['BB']:
    plot(ell, prefix*bb, 'k-', alpha=.015)

# plot(ell, prefix*bbmean, 'b-')
# errorbar(ell, prefix*bbmean, fmt=None, yerr=prefix*bbstd, color='b')
bbmean = cld['BB_mean']
bbstd = cld['BB_std']
errorbar(ell, prefix*bbmean, yerr=prefix*bbstd, label='Output')
plot(ell, 0*ell, 'g--', label='Input')

nmax = where(ell<20.)[0][-1]
# bbmax = prefix[nmax]*(bbmean[nmax] + 5*bbstd[nmax])*2
bbmax = 0.007
# ylim(ymax=bbmax)
# ylim([0, bbmax])
# xlim(xmin=2, xmax=20)
xlim(xmin=2, xmax=150)
xscale('log')

xlabel(r'$\ell$', fontsize=24)
ylabel(r'$D_\ell^{BB}$ ($\mu K^2$)', fontsize=24)

ell0p1 = cl0p1dict['ell']
bb0p1 = cl0p1dict['BB']
prefix = ell0p1*(ell0p1 + 1)/(2*np.pi)

plot(ell0p1, prefix*bb0p1, 'k--', label='r = 0.1')
legend(loc='upper left', fontsize=20)
title('Pixel Noise = PIPER = 4 uK/pixel; Modified covariance')
savefig('piper_noise_mod.png', dpi=150)



In [1259]:
for bb in cldict['BB']:
    ell = cldict['ell']
    plot(ell, ell*(ell+1)/(2*np.pi)*bb)

xscale('log')
xlim(xmax=20)
ylim(ymax=.0002)


Out[1259]:
(0.0, 0.0002)

In [1243]:
reload(calculations.ilc_foreground_cleaning)
import calculations

freqs2 = np.array([200., 270.])*1.e9
freqs4 = np.array([200., 270., 350., 600.])*1.e9
freqs14 = np.array([200., 600.])*1.e9

Omega_single = 7.45446e-6  # sr, from bjohnson report
Omega_pixel = hp.nside2pixarea(512)  # pixel solid angle in sr
Omega_array = 9.54171e-3  # sr, from bjohnson report
regnoise = 3.*np.sqrt(Omega_single/Omega_pixel)/10.  # uK/pixel
print "pixel noise = {0} uK".format(regnoise)

nf = False
retdict2 = calculations.ilc_foreground_cleaning.polarized_ilc_reconstruction(freqs2, _debug=True,
                 fname='junk/explanatory_cl.dat', lensed=False, regnoise=regnoise, noisefactor=nf)
retdict4 = calculations.ilc_foreground_cleaning.polarized_ilc_reconstruction(freqs4, _debug=True,
                 fname='junk/explanatory_cl.dat', lensed=False, regnoise=regnoise, noisefactor=nf)
retdict14 = calculations.ilc_foreground_cleaning.polarized_ilc_reconstruction(freqs14, _debug=True,
                 fname='junk/explanatory_cl.dat', lensed=False, regnoise=regnoise, noisefactor=nf)
print '-'*30, 'ALL DONE', '-'*30


pixel noise = 0.409812311319 uK
--------------------------------------------------------------------------------
Performing ILC simulation with frequencies: [ 200.  270.] GHz
Constructing CMB temperature and polarization maps.
Sigma is 0.000000 arcmin (0.000000 rad) 
-> fwhm is 0.000000 arcmin
Loaded dust from file: /Users/jlazear/github/cmb/data/dust.npy
Combining maps.
Instrument/regularization noise with std dev 0.409812311319 uK
Performing ILC on Q/U maps separately.
Reconstructing ILC Q/U maps.
Computing C_l's from reconstructed maps.
Performing analytical ILC simulation in foreground/background space
Done with ILC with frequencies [ 200.  270.] GHz!
--------------------------------------------------------------------------------
--------------------------------------------------------------------------------
Performing ILC simulation with frequencies: [ 200.  270.  350.  600.] GHz
Constructing CMB temperature and polarization maps.
Sigma is 0.000000 arcmin (0.000000 rad) 
-> fwhm is 0.000000 arcmin
Loaded dust from file: /Users/jlazear/github/cmb/data/dust.npy
Combining maps.
Instrument/regularization noise with std dev 0.409812311319 uK
Performing ILC on Q/U maps separately.
Reconstructing ILC Q/U maps.
Computing C_l's from reconstructed maps.
Performing analytical ILC simulation in foreground/background space
Done with ILC with frequencies [ 200.  270.  350.  600.] GHz!
--------------------------------------------------------------------------------
--------------------------------------------------------------------------------
Performing ILC simulation with frequencies: [ 200.  600.] GHz
Constructing CMB temperature and polarization maps.
Sigma is 0.000000 arcmin (0.000000 rad) 
-> fwhm is 0.000000 arcmin
Loaded dust from file: /Users/jlazear/github/cmb/data/dust.npy
Combining maps.
Instrument/regularization noise with std dev 0.409812311319 uK
Performing ILC on Q/U maps separately.
Reconstructing ILC Q/U maps.
Computing C_l's from reconstructed maps.
Performing analytical ILC simulation in foreground/background space
Done with ILC with frequencies [ 200.  600.] GHz!
--------------------------------------------------------------------------------
------------------------------ ALL DONE ------------------------------

In [1244]:
rdict = retdict2
cmbQUmaps = rdict['cmbQUmaps']
dustmaps = rdict['dustmaps']
noisemap = rdict['noisemap']
ilcQmap = rdict['ilcQmap']

hp.mollview(cmbQUmaps[0])
hp.mollview(dustmaps[0][0])
hp.mollview(noisemap)
hp.mollview(ilcQmap - cmbQUmaps[0])
print cmbQUmaps.shape
print dustmaps.shape
print noisemap.shape

print np.mean(cmbQUmaps[0]), np.std(cmbQUmaps[0]), np.std(cmbQUmaps[0])/np.sqrt(len(noisemap) - 1)
print np.mean(noisemap), np.std(noisemap), np.std(noisemap)/np.sqrt(len(noisemap) - 1)


(2, 3145728)
(2, 2, 3145728)
(3145728,)
0.0677824099768 4.25970807387 0.00240170311916
-0.000690470076247 0.409825854064 0.00023106748513

In [1242]:
rdict = retdict2
cl_in = rdict['cl_in']
cl_out = rdict['cl_out']

xxs = ['EE', 'BB']  # Available: ['TT', 'EE', 'TE', 'BB']
figs = {}
axes = {}
for xx in xxs:
    f, a = util.plot_funcs.Dl_resid_plot(cl_out['ell'], cl_out[xx], cl_in['ell'], cl_in[xx],
                                         label1='Output {0}'.format(xx),
                                         label2='Input {0}'.format(xx))
    figs[xx] = f
    axes[xx] = a

try:
    axesBB = axes['BB']
    axesBB[0].set_xlim([2, 50])
    axesBB[1].set_xlim([2, 50])
    axesBB[0].set_ylim([-0.001,.008])
#     axesBB[1].set_ylim([-.1,.1])
    axesBB[1].set_ylim([-.002,.002])
    axesBB[0].set_xscale('log')
    axesBB[1].set_xscale('log')
except KeyError:
    pass

util.plot_funcs.print_summary(rdict)


frequencies = [ 200.  270.] GHz
Q weights          = [ 2.34771765 -1.34771765]
Q weights expected = [ 2.84933324 -1.84933324]
log10(Q condition number) = 9.31228278679
var(ILC Q)          = 34.9654529043
var(ILC Q) expected = 30307.0632065
U weights          = [ 2.3475355 -1.3475355]
U weights expected = [ 2.04436231 -1.04436231]
log10(U condition number) = 11.266231145
var(ILC U)          = 34.8920765989
var(ILC U) expected = 1797.5338851

In [1245]:
rdict = retdict2
cl_in = rdict['cl_in']
cl_out = rdict['cl_out']

xxs = ['EE', 'BB']  # Available: ['TT', 'EE', 'TE', 'BB']
figs = {}
axes = {}
for xx in xxs:
    f, a = util.plot_funcs.Dl_resid_plot(cl_out['ell'], cl_out[xx], cl_in['ell'], cl_in[xx],
                                         label1='Output {0}'.format(xx),
                                         label2='Input {0}'.format(xx))
    figs[xx] = f
    axes[xx] = a

try:
    axesBB = axes['BB']
    axesBB[0].set_xlim([2, 50])
    axesBB[1].set_xlim([2, 50])
    axesBB[0].set_ylim([-0.001,.008])
#     axesBB[1].set_ylim([-.1,.1])
    axesBB[1].set_ylim([-.002,.002])
    axesBB[0].set_xscale('log')
    axesBB[1].set_xscale('log')
except KeyError:
    pass

util.plot_funcs.print_summary(rdict)


frequencies = [ 200.  270.] GHz
Q weights          = [ 2.34769091 -1.34769091]
Q weights expected = [ 2.63847333 -1.63847333]
log10(Q condition number) = 9.12629454212
var(ILC Q)          = 18.313651587
var(ILC Q) expected = 10191.0258483
U weights          = [ 2.34762116 -1.34762116]
U weights expected = [ 2.26804996 -1.26804996]
log10(U condition number) = 10.6217093567
var(ILC U)          = 18.2982360384
var(ILC U) expected = 139.719364522

In [1156]:
rdict = retdict4
cl_in = rdict['cl_in']
cl_out = rdict['cl_out']

xxs = ['BB']  # Available: ['TT', 'EE', 'TE', 'BB']
figs = {}
axes = {}
for xx in xxs:
    f, a = util.plot_funcs.Dl_resid_plot(cl_out['ell'], cl_out[xx], cl_in['ell'], cl_in[xx],
                                         label1='Output {0}'.format(xx),
                                         label2='Input {0}'.format(xx))
    figs[xx] = f
    axes[xx] = a

try:
    axesBB = axes['BB']
    axesBB[0].set_xlim([2, 500])
    axesBB[1].set_xlim([2, 500])
    axesBB[0].set_ylim([0,.008])
    axesBB[1].set_ylim([-.002,.002])
    axesBB[0].set_xscale('log')
    axesBB[1].set_xscale('log')
except KeyError:
    pass

util.plot_funcs.print_summary(rdict)


frequencies = [ 200.  270.  350.  600.] GHz
Q weights          = [ 0.5  0.5]
Q weights expected = [ 2.913994 -1.913994]
log10(Q condition number) = 9.20307105801
var(ILC Q)          = 34.9322660896
var(ILC Q) expected = 34.9322660245
U weights          = [ 0.5  0.5]
U weights expected = [ 2.31382355 -1.31382355]
log10(U condition number) = 10.1098003801
var(ILC U)          = 35.0413982038
var(ILC U) expected = 35.0413982038

In [667]:
rdict = retdict14
cl_in = rdict['cl_in']
cl_out = rdict['cl_out']

xxs = ['BB']  # Available: ['TT', 'EE', 'TE', 'BB']
figs = {}
axes = {}
for xx in xxs:
    f, a = util.plot_funcs.Dl_resid_plot(cl_out['ell'], cl_out[xx], cl_in['ell'], cl_in[xx],
                                         label1='Output {0}'.format(xx),
                                         label2='Input {0}'.format(xx))
    figs[xx] = f
    axes[xx] = a

try:
    axesBB = axes['BB']
    axesBB[0].set_xlim([2, 500])
    axesBB[1].set_xlim([2, 500])
    axesBB[0].set_ylim([0,.008])
    axesBB[1].set_ylim([-.002,.002])
    axesBB[0].set_xscale('log')
    axesBB[1].set_xscale('log')
except KeyError:
    pass

util.plot_funcs.print_summary(rdict)


frequencies = [ 200.  600.] GHz
Q weights          = [ 1.01161714 -0.01161714]
Q weights expected = [ 1.00638745 -0.00638745]
log10(Q condition number) = 11.4642386639
var(ILC Q)          = 18.1563461245
var(ILC Q) expected = 44295.6164144
U weights          = [ 1.01161754 -0.01161754]
U weights expected = [ 1.01149118 -0.01149118]
log10(U condition number) = 12.4831330348
var(ILC U)          = 18.1584997828
var(ILC U) expected = 22.2784486219

Checking the function for errors...


In [90]:
#     frequencies = np.array(frequencies)

#     # Construct CMB maps from spectra
#     (Tmap, Qmap, Umap), cldict = lib.cmb.generate_maps(Nside=Nside,
#                         n2r=True, return_cls=True, fname=fname, lensed=lensed)
#     cmbQUmaps = np.vstack([Qmap, Umap])  # uK

#     # 3 deg smoothed QU polarized synchrotron-tracking dust maps
#     MJypsr = lib.conversions.MJypsr
#     convs = [1.e6/(lib.conversions.dBdT(f)/MJypsr)  # uK_CMB/(MJy/sr)
#              for f in frequencies]

#     dustmaps = [lib.foregrounds.generate_synchro_traced_dust_QU_map(nu=nu,
#                         Nside=Nside,  n2r=True, sigma=3.*np.pi/180)
#                 for nu in frequencies]
#     dustmaps = [dustmaps[i]*convs[i] for i in range(len(dustmaps))]   # uK_CMB

#     # Construct CMB + foreground maps
#     totalQUmaps = [cmbQUmaps + dustmap for dustmap in dustmaps]

#     # Perform ILC on each set of Q and U maps
#     #
#     # # Get WMAP masks and regions
#     # maskdata = 'data/wmap_ilc_rgn_defn_9yr_v5.fits'
#     # masks, region = load_wmap_masks(maskdata)
#     # ws = compute_ilc_weights(maps, masks)

#     totalQmaps = np.vstack([totalQUmaps[i][0] for i in range(len(totalQUmaps))])
#     totalUmaps = np.vstack([totalQUmaps[i][1] for i in range(len(totalQUmaps))])

#     Qweights = lib.ilc.compute_ilc_weights(totalQmaps)
#     Uweights = lib.ilc.compute_ilc_weights(totalUmaps)

#     ilcQmap = np.dot(Qweights, totalQmaps)
#     ilcUmap = np.dot(Uweights, totalUmaps)

#     recon_cls = hp.anafast([Tmap, ilcQmap, ilcUmap])  # TT, EE,  BB, TE, EB, TB
#     labels = ['TT', 'EE', 'BB', 'TE', 'EB', 'TB']
#     recon_cldict = {labels[i]: recon_cls[i] for i in range(len(labels))}
#     recon_cldict['ell'] = np.arange(len(recon_cldict['TT']), dtype='float')
    
#     retdict = {'cl_in': cldict, 'cl_out': recon_cldict,
#                'cmbTmap': Tmap, 'cmbQUmaps': cmbQUmaps, 'dustmaps': dustmaps,
#                'totalQUmaps': totalQUmaps,
#                'weights_Q': Qweights, 'weights_U': Uweights,
#                'ilcQmap': ilcQmap, 'ilcUmap': ilcUmap}

#     return retdict

In [414]:
#     # Construct CMB maps from spectra
#     (Tmap, Qmap, Umap), cldict = lib.cmb.generate_maps(Nside=Nside,
#                         n2r=True, return_cls=True, fname=fname, lensed=lensed)
#     cmbQUmaps = np.vstack([Qmap, Umap])  # uK


# Testing this line
(Tmap, Qmap, Umap), cldict = lib.cmb.generate_maps(fname='junk/explanatory_cl.dat', 
                                                   lensed=False, return_cls=True, n2r=True)
cmbQUmaps = np.vstack([Qmap, Umap])  # uK


# Verify that the generated maps are sensible and match the input spectra
results = hp.anafast([Tmap, Qmap, Umap])
ells = np.arange(len(results[0]))

xxs = ['TT', 'EE', 'BB', 'TE', 'EB', 'TB']

for i, result in enumerate(results):
    xx = xxs[i]
    figure()
    plot(ells, ells*(ells+1)*result/(2*np.pi))
    eo = cldict['ell']
    f = eo*eo/(2*np.pi)
    try:
        plot(eo, f*cldict[xx])
    except KeyError:
        pass
    
print np.shape(cmbQUmaps)


Sigma is 0.000000 arcmin (0.000000 rad) 
-> fwhm is 0.000000 arcmin
(2, 3145728)

In [524]:
#     # 3 deg smoothed QU polarized synchrotron-tracking dust maps
#     MJypsr = lib.conversions.MJypsr
#     convs = [1.e6/(lib.conversions.dBdT(f)/MJypsr)  # uK_CMB/(MJy/sr)
#              for f in frequencies]

#     dustmaps = [lib.foregrounds.generate_synchro_traced_dust_QU_map(nu=nu,
#                         Nside=Nside,  n2r=True, sigma=3.*np.pi/180)
#                 for nu in frequencies]
#     dustmaps = [dustmaps[i]*convs[i] for i in range(len(dustmaps))]   # uK_CMB

    # 3 deg smoothed QU polarized synchrotron-tracking dust maps
frequencies = np.array([200.,270.,])*1.e9  # Hz
MJypsr = lib.conversions.MJypsr
convs = [1.e6/(lib.conversions.dBdT(f)/MJypsr)  # uK_CMB/(MJy/sr)
         for f in frequencies]

dustmaps = [lib.foregrounds.generate_synchro_traced_dust_QU_map(nu=nu, n2r=True, sigma=3.*np.pi/180,
                                                                Nside=512)
            for nu in frequencies]  # in MJy/sr

lefthalf = where(hp.pix2ang(512, arange(len(dm)))[1] >= np.pi)

for i in range(len(dustmaps)):
    for j in range(len(dustmaps[0])):
        dustmaps[i][j][lefthalf] *= 1#.5#*(i+1)

dustmaps = [dustmaps[i]*convs[i] for i in range(len(dustmaps))]   # MJy/sr -> uK_CMB

for i in range(len(dustmaps)):
    dustmap = dustmaps[i]
    for j in range(len(dustmap)):
        d = dustmaps[i][j]# /dustmaps[2][j]
        hp.mollview(d, norm='hist', unit='uK')
        print np.shape(where(np.isnan(d)))
        d2 = d[np.where(~np.isnan(d))[0]]
        print np.mean(d2), np.std(d2)
        
print np.shape(dustmaps)


freq = 200.0 GHz, factor = 0.40291038777
Sigma is 0.000000 arcmin (0.000000 rad) 
-> fwhm is 0.000000 arcmin
Sigma is 0.000000 arcmin (0.000000 rad) 
-> fwhm is 0.000000 arcmin
freq = 270.0 GHz, factor = 0.651242138951
Sigma is 0.000000 arcmin (0.000000 rad) 
-> fwhm is 0.000000 arcmin
Sigma is 0.000000 arcmin (0.000000 rad) 
-> fwhm is 0.000000 arcmin
(1, 0)
80.5932044722 467.446556851
(1, 0)
4.87127159843 186.626355317
(1, 0)
140.395339567 814.303371066
(1, 0)
8.48587461593 325.107689722
(2, 2, 3145728)

In [525]:
def naninfmean(d):
    return np.mean(notnaninf(d))

def notnaninf(d):
    d2 = d[~np.isnan(d)]
    d3 = d2[~np.isinf(d2)]
    return d3

# beta = 1.6
# (frequencies/353.e9)**beta

meanfracs = [[], []]
for i in range(len(dustmaps)):
    for j in range(len(cmbQUmaps)):
        dm = np.abs(dustmaps[i][j])
        cm = np.abs(cmbQUmaps[j])
        mf = dm.mean()/cm.mean()
        meanfracs[j].append(mf)
        print "mean fraction = ", mf

meanfracs = np.array(meanfracs)        
        
for i, meanfrac in enumerate(meanfracs):
    print "meanfrac[{0}] = ".format(i), meanfrac

for meanfrac in meanfracs:
    print 1/(meanfrac/meanfrac.min())
        
# dm = np.abs(dustmaps[0][0])
# cm = np.abs(cmbQUmaps[0])
# print dm.mean(), dm.std()
# print cm.mean(), cm.std()
# hist(dm, log=True)
# figure()
# hist(cm, log=True)


mean fraction =  37.6169648036
mean fraction =  18.8981144476
mean fraction =  65.529675625
mean fraction =  32.9209790253
meanfrac[0] =  [ 37.6169648   65.52967563]
meanfrac[1] =  [ 18.89811445  32.92097903]
[ 1.          0.57404473]
[ 1.          0.57404473]

In [526]:
#     # Construct CMB + foreground maps
#     totalQUmaps = [cmbQUmaps + dustmap for dustmap in dustmaps]

totalQUmaps = [cmbQUmaps + dustmap for dustmap in dustmaps]
np.shape(totalQUmaps)

# Try for various combinations of indices... Matching indices should give ~0 maps,
# others should give maps that look like the odd component out. 
# They all seem fine...
# hp.mollview((totalQUmaps[2][0] - dustmaps[2][0]) - cmbQUmaps[0], norm='hist')


Out[526]:
(2, 2, 3145728)

In [527]:
#     # Perform ILC on each set of Q and U maps
#     #
#     # # Get WMAP masks and regions
#     # maskdata = 'data/wmap_ilc_rgn_defn_9yr_v5.fits'
#     # masks, region = load_wmap_masks(maskdata)
#     # ws = compute_ilc_weights(maps, masks)

#     totalQmaps = np.vstack([totalQUmaps[i][0] for i in range(len(totalQUmaps))])
#     totalUmaps = np.vstack([totalQUmaps[i][1] for i in range(len(totalQUmaps))])

totalQmaps = np.vstack([totalQUmaps[i][0] for i in range(len(totalQUmaps))])
totalUmaps = np.vstack([totalQUmaps[i][1] for i in range(len(totalQUmaps))])

print np.shape(totalQmaps), "should be (4, 3145728)"
print np.shape(totalUmaps), "should be (4, 3145728)"

# Just checking that we re-grouped the maps correctly. Should all be True.
print [np.std(totalQmaps[i] - totalQUmaps[i][0]) == 0. for i in range(len(totalQmaps))]
print [np.std(totalUmaps[i] - totalQUmaps[i][1]) == 0. for i in range(len(totalUmaps))]

# Cross-checking the maps. Should all be non-zero.
print [np.std(totalQmaps[i] - totalQUmaps[i][1]) for i in range(len(totalQmaps))]
print [np.std(totalUmaps[i] - totalQUmaps[i][0]) for i in range(len(totalUmaps))]


(2, 3145728) should be (4, 3145728)
(2, 3145728) should be (4, 3145728)
[True, True]
[True, True]
[454.45839079957824, 791.64635836094124]
[454.45839079957824, 791.64635836094124]

In [528]:
#     Qweights = lib.ilc.compute_ilc_weights(totalQmaps)
#     Uweights = lib.ilc.compute_ilc_weights(totalUmaps)

#     ilcQmap = np.dot(Qweights, totalQmaps)
#     ilcUmap = np.dot(Uweights, totalUmaps)

Qweights = lib.ilc.compute_ilc_weights(totalQmaps)
Uweights = lib.ilc.compute_ilc_weights(totalUmaps)

ilcQmap = np.dot(Qweights, totalQmaps)
ilcUmap = np.dot(Uweights, totalUmaps)

print frequencies/1.e9
print Qweights
print Uweights

dQ = np.abs((ilcQmap - cmbQUmaps[0])/cmbQUmaps[0]) + .01
print dQ.min(), dQ.max()
dU = np.abs((ilcUmap - cmbQUmaps[1])/cmbQUmaps[1]) + .01

print ilcQmap.mean(), ilcQmap.std()
print cmbQUmaps[0].mean(), cmbQUmaps[0].std()
print dQ.mean(), dQ.std()

hp.mollview(ilcQmap - cmbQUmaps[0], norm='hist')
hp.mollview(dustmaps[0][0], norm='hist')

# hist(dQ, log=True)
# figure()
# hp.mollview(dQ, norm='log', min=.01, max=.1)


[ 200.  270.]
[ 2.3476266 -1.3476266]
[ 2.34769775 -1.34769775]
0.01 1510.96521336
0.0252528613484 4.25720372916
0.0229961589751 4.25722386856
0.021203997594 1.79885374545

In [529]:
S = cmbQUmaps[0]

for i in range(len(dustmaps)):
    F = dustmaps[i][0]
    print "bias factor @ {0} GHz, <SF>/<FF> =".format(frequencies[i]/1.e9), np.mean(S*F)/np.mean(F*F)
    Savg = np.mean(S)
    Favg = np.mean(F)
    Thatavg = np.mean(ilcQmap)
    print "bias factor2 @ {0} GHz, (<S> - <That>)/<F> =".format(frequencies[i]/1.e9), (Savg - Thatavg)/Favg


bias factor @ 200.0 GHz, <SF>/<FF> = -1.89555601586e-05
bias factor2 @ 200.0 GHz, (<S> - <That>)/<F> = -2.80011495771e-05
bias factor @ 270.0 GHz, <SF>/<FF> = -1.08813392826e-05
bias factor2 @ 270.0 GHz, (<S> - <That>)/<F> = -1.60739122843e-05

In [573]:
sigmaB2 = np.var(Qmap)
ones = np.ones([len(dustmaps), 1])

X = [covcoeff(Qmap, dustmaps[i][0]) for i in range(len(dustmaps))]
X = np.array([X]).T  # Make sure N_freqs x 1

F = np.cov(np.array(dustmaps)[..., 0])  # cov(R, R)

C = sigmaB2*ones*ones.T + X*ones.T + ones*X.T + F
Cinv = np.linalg.pinv(C)

print "det C =", np.linalg.det(C)
print "cond C =", np.linalg.cond(C) #np.linalg.norm(C)*np.linalg.norm(Cinv)
print "cond C =", np.linalg.norm(C)*np.linalg.norm(Cinv)
print "C =\n", C
print "C^-1 =\n", Cinv


det C = 172.13051595
cond C = 35.7338073468
cond C = 35.7617920502
C =
[[ 25.20186903  34.99383031]
 [ 34.99383031  55.42044022]]
C^-1 =
[[ 0.32196755 -0.20329824]
 [-0.20329824  0.1464114 ]]

In [574]:
print ones*X.T - np.outer(ones, X)
print X*ones.T - np.outer(X, ones)


[[ 0.  0.]
 [ 0.  0.]]
[[ 0.  0.]
 [ 0.  0.]]

In [597]:
print np.einsum("ij,jk,kl", ones.T, G_Q, ones)

def magnitude(A):
    wons = np.ones(len(A))
    return np.einsum("i,ij,j", wons, G_Q, wons)

print magnitude(G_Q)


[[ 111.66756077]]
111.667560772

In [578]:
def covcoeff(x, y):
    return np.mean(x*y) - np.mean(x)*np.mean(y)

X_Q = [covcoeff(Qmap, dustmaps[i][0]) for i in range(len(dustmaps))]
X_Q = np.array([X_Q]).T  # Make sure 1 x N_freqs

F_Q = np.cov(np.array(dustmaps)[..., 0])  # cov(R, R)

ones = np.ones([len(dustmaps), 1])

G_Q = F_Q + ones*X_Q.T
G_Qinv = np.linalg.pinv(G_Q)
F_Qinv = np.linalg.pinv(F_Q)

print "det =\n", np.linalg.det(G_Q)
print "G.G^-1 =\n", np.dot(G_Q, G_Qinv)
print "G =\n", G_Q
print "G^-1 =\n", G_Qinv
print "condition number =\n", np.linalg.cond(G_Q)

print "det =\n", np.linalg.det(F_Q)
print "F.F^-1 =\n", np.dot(F_Q, F_Qinv)
print "F =\n", F_Q
print "F^-1 =\n", F_Qinv
print "condition number =\n", np.linalg.cond(F_Q)


det =
4.5387727773e-06
G.G^-1 =
[[  1.00000009e+00  -2.98023224e-08]
 [  5.96046448e-08   9.99999970e-01]]
G =
[[ 13.19627848  22.98823977]
 [ 27.5282163   47.95482621]]
G^-1 =
[[ 10565593.60757547  -5064858.29142792]
 [ -6065123.55854286   2907455.34092085]]
condition number =
828433005.686
det =
0.0
F.F^-1 =
[[ 0.24785302  0.43176602]
 [ 0.43176602  0.75214698]]
F =
[[ 19.31464301  33.64658083]
 [ 33.64658083  58.61316727]]
F^-1 =
[[ 0.00318055  0.00554059]
 [ 0.00554059  0.00965184]]
condition number =
1.10147632603e+16

In [420]:
#     recon_cls = hp.anafast([Tmap, ilcQmap, ilcUmap])  # TT, EE,  BB, TE, EB, TB
#     labels = ['TT', 'EE', 'BB', 'TE', 'EB', 'TB']
#     recon_cldict = {labels[i]: recon_cls[i] for i in range(len(labels))}
#     recon_cldict['ell'] = np.arange(len(recon_cldict['TT']), dtype='float')
    
#     retdict = {'cl_in': cldict, 'cl_out': recon_cldict,
#                'cmbTmap': Tmap, 'cmbQUmaps': cmbQUmaps, 'dustmaps': dustmaps,
#                'totalQUmaps': totalQUmaps,
#                'weights_Q': Qweights, 'weights_U': Uweights,
#                'ilcQmap': ilcQmap, 'ilcUmap': ilcUmap}

#     return retdict

recon_cls = hp.anafast([Tmap, ilcQmap, ilcUmap])  # TT, EE,  BB, TE, EB, TB
labels = ['TT', 'EE', 'BB', 'TE', 'EB', 'TB']
recon_cldict = {labels[i]: recon_cls[i] for i in range(len(labels))}
recon_cldict['ell'] = np.arange(len(recon_cldict['TT']), dtype='float')

retdict = {'cl_in': cldict, 'cl_out': recon_cldict,
           'cmbTmap': Tmap, 'cmbQUmaps': cmbQUmaps, 'dustmaps': dustmaps,
           'totalQUmaps': totalQUmaps,
           'weights_Q': Qweights, 'weights_U': Uweights,
           'ilcQmap': ilcQmap, 'ilcUmap': ilcUmap}

In [421]:
clin = retdict['cl_in']
infact = clin['ell']*(clin['ell'] + 1)/(2*np.pi)
clout = retdict['cl_out']
outfact = clout['ell']*(clout['ell'] + 1)/(2*np.pi)


figtest, axestest = resid_plot(clout['ell'], outfact*clout['BB'], clin['ell'], infact*clin['BB'])
axestest[0].set_xscale('log')
axestest[0].set_xlim([1, 500])


Out[421]:
(1, 500)

Messing around with CLASS


In [3]:
reload(lib.cmb)
import lib.cmb

In [161]:
ells_lensed = lib.cmb.get_lcdm_cls(True)['ell']
TTs_lensed = lib.cmb.get_lcdm_cls(True)['TT']
ells_unlensed = lib.cmb.get_lcdm_cls(False)['ell']
TTs_unlensed = lib.cmb.get_lcdm_cls(False)['TT']
# plot(ells_unlensed, TTs_unlensed, label='unlensed')
# plot(ells_lensed, TTs_lensed, label='lensed')

ells_diff = ells_lensed
TTs_diff = TTs_lensed - TTs_unlensed[:len(ells_diff)]

plot(ells_diff, TTs_diff, label='diff')


Out[161]:
[<matplotlib.lines.Line2D at 0x1154f3d90>]

In [179]:
ilcmap = lib.wmap.load_wmap_maps(['data/wmap_ilc_9yr_v5.fits'])[0]

# cldict = lib.cmb.get_lcdm_cls(True)
# ells = cldict['ell']   # 1/rad
# TTs = cldict['TT']  # uK^2
# EEs = cldict['EE']  # uK^2
# TEs = cldict['TE']  # uK^2
# BBs = cldict['BB']  # uK^2

# TTmap = hp.synfast(TTs, nside=512)
# hp.mollview(TTmap/1.e3, unit='mK')
# hp.mollview(ilcmap, unit='mK')

# Tmap, Qmap, Umap = hp.synfast([TTs, EEs, BBs, TEs], nside=512, new=True, pol=True)
Tmap, Qmap, Umap = lib.cmb.generate_lcdm_maps(n2r=True)


Sigma is 0.000000 arcmin (0.000000 rad) 
-> fwhm is 0.000000 arcmin

In [182]:
hp.mollview(Tmap, unit='uK', title='T')
hp.mollview(Qmap, unit='uK', title='Q', norm='hist', nest=False)
hp.mollview(Umap, unit='uK', title='U')



In [12]:
#sns.set_style('whitegrid')
sns.set_context("notebook", font_scale=1.5, rc={"lines.linewidth": 2.5})

In [140]:
print ells[:10]
print TTs[:10]


[ 0.  1.  2.  3.  4.  5.  6.  7.  8.  9.]
[   0.            0.          980.78143996  927.74297424  875.7440873
  835.63245842  807.05526219  787.82929472  776.49767959  770.58550408]

In [138]:
#cldict = lib.cmb.get_lcdm_lensed_cls()
cldict = lib.cmb.get_lcdm_cls(True)
ells = cldict['ell']
TTs = cldict['TT']
EEs = cldict['EE']
TEs = cldict['TE']
BBs = cldict['BB']
plot(wmapells, wmapTTs, 'o', label='WMAP TT')
plot(ells, TTs, label='TT')
#plot(ells, EEs, label='EE')
#plot(ells, np.abs(TEs), label='|TE|')
#plot(ells, BBs, label='BB')

xlim([0, 1000])
#ylim(ymin=1.e-3)
yscale('log')
xlabel('$\ell$')
ylabel('$D_\ell^{XX} (\mu K^2)$')
legend(loc='lower right')#bbox_to_anchor=(1.5, 0.7))


Out[138]:
<matplotlib.legend.Legend at 0x116e19750>

In [51]:
import lib.ilc

ilcmap = lib.wmap.load_wmap_maps(['data/wmap_ilc_9yr_v5.fits'])[0]

In [87]:
ilcTTs = hp.anafast(ilcmap/1.e3, lmax=3500, pol=False)

In [82]:
wmapells, wmapTTs, _, _, _ = np.loadtxt('data/wmap_tt_spectrum_9yr_v5.txt').T

In [107]:
wmapells2 = np.hstack([[0, 1], wmapells])
wmapTTs2 = np.hstack([[0., 0.], wmapTTs])

In [115]:
# synfast expects a set of bare Cl's ranging from ell = 0 to ell = lmax. 
# The Cl's are related to the "standard" D_\ell by D_\ell = \ell*(\ell+1)/(2\pi)*C_\ell
# synfast will preserve units, e.g. if the Cl's are in uK^2, then the resulting map will be in uK
# and if the Cl's are in normalized (unitless) form, then the resulting map will be fractional
# anisotropies (unitless). 
wmapTTmap = hp.synfast(wmapTTs*2*np.pi/(wmapells*(wmapells+1)), nside=512)  # uK^2 -> uK
#hp.mollview(wmapTTmap/1.e3, unit='mK', title='Synfast(WMAP $C_\ell^{TT}$)')
#hp.mollview(ilcmap, unit='mK', title='WMAP ILC Temperature')

wmapTTspost = hp.anafast(wmapTTmap)  # uK -> uK^2
wmapellspost = np.arange(len(wmapTTspost))


Sigma is 0.000000 arcmin (0.000000 rad) 
-> fwhm is 0.000000 arcmin
The order of the input cl's will change in a future release.
Use new=True keyword to start using the new order.
See documentation of healpy.synalm.

In [129]:
print wmapTTs[:5] # l(l+1)/(2pi)Cl in uK^2
print (wmapTTspost)[2:5+2]  # uK^2 (???)
plot(wmapells, wmapTTs)
plot(wmapellspost, wmapTTspost*wmapellspost*(wmapellspost+1)/(2*np.pi))
xlim(xmax=10)


[  150.6398   902.1805   730.3626  1467.8517   688.5324]
[ 255.73787426  395.72130696   82.59767457  184.56615428   38.65654762]
Out[129]:
(0.0, 10)

In [133]:
wmapTTspostmap = hp.synfast(wmapTTspost, nside=512)
hp.mollview(wmapTTspostmap)


Sigma is 0.000000 arcmin (0.000000 rad) 
-> fwhm is 0.000000 arcmin
The order of the input cl's will change in a future release.
Use new=True keyword to start using the new order.
See documentation of healpy.synalm.

In [92]:
ilcells = np.arange(len(ilcTTs))
print ilcells[:5]
print (ilcells*(ilcells+1)*ilcTTs*(2.735)**2)[:5]

print TTs[:5]


[0 1 2 3 4]
[  0.00000000e+00   2.32285425e-06   2.68767302e-05   1.20872182e-04
   4.85548724e-05]
[ 980.78143996  927.74297424  875.7440873   835.63245842  807.05526219]

In [343]:
lcdm_cl = np.loadtxt('data/lcdm_cl.dat')
lcdm_cl.shape
ells = lcdm_cl.T[0]
TTs = lcdm_cl.T[1]

plot(ells, TTs)


Out[343]:
[<matplotlib.lines.Line2D at 0x11a172050>]

In [339]:
pdict = lib.cmb.get_params_from_file('lcdm.ini')
cosmo = Class()
cosmo.set(pdict)
cosmo.struct_cleanup()
cosmo.compute()


---------------------------------------------------------------------------
CosmoSevereError                          Traceback (most recent call last)
<ipython-input-339-f0b1a0aefb8a> in <module>()
      3 cosmo.set(pdict)
      4 cosmo.struct_cleanup()
----> 5 cosmo.compute()

/usr/local/lib/python2.7/site-packages/classy.so in classy.Class.compute (classy.c:4312)()

CosmoSevereError: 

Error in Class: Class did not read input parameter(s): selection, P_k_max_h/Mpc, w0_fld, selection_width, bessels_verbose, z_pk, l_max_tensors, selection_mean, bessel file, wa_fld, cs2_fld

In [256]:
from classy import Class

In [311]:
# Define your cosmology (what is not specified will be set to CLASS default parameters)
params = {
    'output': 'tCl lCl pCl',
    'l_max_scalars': 2000,
    'lensing': 'yes',
    'A_s': 2.3e-9,
    'n_s': 0.9624, 
    'h': 0.6711,
    'omega_b': 0.022068,
    'omega_cdm': 0.12029,
    'modes': 's,t',
    'r': 0.01}

In [258]:
# Create an instance of the CLASS wrapper
cosmo = Class()

In [312]:
# Set the parameters to the cosmological code
cosmo.set(params)

# Run the whole code. Depending on your output, it will call the
# CLASS modules more or less fast. For instance, without any
# output asked, CLASS will only compute background quantities,
# thus running almost instantaneously.
# This is equivalent to the beginning of the `main` routine of CLASS,
# with all the struct_init() methods called.
cosmo.struct_cleanup()
cosmo.compute()

In [313]:
# Access the lensed cl until l=2000
lcls = cosmo.lensed_cl(2000)
rcls = cosmo.raw_cl(2000)

# Print on screen to see the output
print lcls.keys()
print rcls.keys()


['pp', 'bb', 'ee', 'tt', 'tp', 'te']
['pp', 'bb', 'ee', 'tt', 'tp', 'te']

In [320]:
ells = np.arange(len(cls['tt']))
prefix = ells*(ells+1)/(2*np.pi)*(cosmo.T_cmb()*1.e6)**2
plot(ells, np.sqrt(prefix*lcls['tt']), label='lensed TT')
plot(ells, np.sqrt(prefix*lcls['ee']), label='lensed EE')
plot(ells, np.sqrt(prefix*lcls['bb']), label='lensed BB')
plot(ells, np.sqrt(prefix*rcls['tt']), label='raw TT')
plot(ells, np.sqrt(prefix*rcls['ee']), label='raw EE')
plot(ells, np.sqrt(prefix*rcls['bb']), label='raw BB')

xlabel('$\ell$')
ylabel('$D_\ell^{XX} (\mu K)$')
xscale('log')
yscale('log')

#xlim([0, 50])
legend(bbox_to_anchor=(1.4, 0.7))


Out[320]:
<matplotlib.legend.Legend at 0x1189a3210>

In [290]:
# Clean CLASS (the equivalent of the struct_free() in the `main`
# of CLASS. This step is primordial when running in a loop over different
# cosmologies, as you will saturate your memory very fast if you ommit
# it.
cosmo.struct_cleanup()

ILC Functions


In [253]:
from lib.wmap import load_wmap_masks, load_wmap_maps

load_wmap_masks('/Users/jlazear/github/cmb/data/wmap_ilc_rgn_defn_9yr_v5.fits')


Out[253]:
(array([[False, False, False, ..., False, False, False],
        [False, False, False, ..., False, False, False],
        [False, False, False, ..., False, False, False],
        ..., 
        [False, False, False, ..., False, False, False],
        [False, False, False, ..., False, False, False],
        [False, False, False, ..., False, False, False]], dtype=bool),
 array([0, 0, 0, ..., 0, 0, 0], dtype=uint16))

In [137]:
def load_wmap_maps(fnames, n2r=True):
    fname0 = fnames[0]
    Nbands = len(fnames)
    
    # Initialize array and populate first band
    with fits.open(fname0) as f0:
        Npix = f0[1].header['NAXIS2']
        a = np.empty([Nbands, Npix])
        a[0] = f0[1].data['TEMPERATURE']
    
    for i, fname in enumerate(fnames[1:]):
        with fits.open(fname) as f:
            a[i+1] = f[1].data['TEMPERATURE']
    
    if n2r:
        for i in range(a.shape[0]):
            a[i] = hp.reorder(a[i], n2r=True)
    
    return a

def load_wmap_masks(fname='data/wmap_ilc_rgn_defn_9yr_v5.fits', n2r=True):
    with fits.open(fname) as f:
        fdata = f[1]
        bitmask = fdata.data['TEMPERATURE'].astype('uint16')
        region = fdata.data['N_OBS']
        if n2r:
            bitmask = hp.reorder(bitmask, n2r=True)
            region = hp.reorder(region, n2r=True)

    masks = np.empty([12, len(bitmask)])
    for i in range(12):
        masks[i] = (bitmask >> i) & 1    
    return masks, region

def compute_ilc_weights(maps, mask=None):
    nd = np.ndim(mask)
    if nd == 0:
        return _compute_ilc_weights(maps)
    elif np.ndim(mask) == 1:
        maps = np.ma.masked_array(maps, np.logical_not(np.outer(np.ones(len(maps)), mask)))
        return _compute_ilc_weights(maps)
    else:
        ws = []
        for i in range(len(mask)):
            ms = np.ma.masked_array(maps, np.logical_not(np.outer(np.ones(len(maps)), mask[i])))
            w = _compute_ilc_weights(ms)
            ws.append(w)
        ws = np.array(ws)
        return ws
    
def _compute_ilc_weights(maps):
    if np.ma.isMaskedArray(maps):
        covfunc = lambda x: np.ma.compress_cols(np.ma.cov(x))
    else:
        covfunc = np.cov
    cov = covfunc(maps)
#    return cov #DELME
    icov = np.linalg.pinv(cov)#, rcond=1.e-4)  # Naive inversion, since cov ~ k x k is small
    sumicov = icov.sum(0)
    totalicov = sumicov.sum()
    w = sumicov/totalicov
    return w

def ilc_map_from_weights(maps, weights, regions_map, 
                         sigma=1.5*np.pi/180, 
                         return_weights_map=False):
    """
    Construct an ILC map from a set of raw temperature maps,
    a set of weights per region, and a map of regions.
    
    `maps` must have shape (Nfreq, Npix). `weights` must 
    have shape (Nregions, Nfreq). `regions_map` must have
    shape (Npix,) and each pixel should contain the integer
    identifier to which region the pixel is assigned. 
    
    `sigma` is the smoothing factor to reduce edge effects.
    It is applied to each region's weight map before
    multiplying into the raw maps. See Bennett 2003 or 
    Hinshaw 2007 for details. 
    
    All of the maps must be in RING format!
    """
    Thats = np.dot(weights, maps)
    That = np.zeros(Thats.shape[1])
    weights_map = np.zeros(Thats.shape[1])
    for i in range(len(Thats)):
        m = np.zeros_like(That)
        m[regions_map == i] = 1
        if sigma is not None:
            m = hp.smoothing(m, sigma=sigma, verbose=False) + m.mean()
        That += Thats[i]*m
        weights_map += m
        print "i={0} m.min, m.max = ".format(i), m.min(), m.max() #DELME
    That = That/weights_map
    if return_weights_map:
        return That, weights_map
    else:
        return That

In [254]:
# jlazear ILC map
bands = ['K', 'Ka', 'Q', 'V', 'W']
fnames = ['data/wmap_band_smth_imap_r9_9yr_{0}_v5.fits'.format(band) for band in bands]
maskdata = 'data/wmap_ilc_rgn_defn_9yr_v5.fits'

maps = load_wmap_maps(fnames)
masks, region = load_wmap_masks(maskdata)
ws = compute_ilc_weights(maps, masks)

T_hat_JL, wm_JL = ilc_map_from_weights(maps, ws, region, return_weights_map=True, sigma=1.5*np.pi/180)

# WMAP 
wmap_ws = np.array([[.1555, -.7572, -.2689, 2.2845, -.4138],
           [.0375, -.5137, .0223, 2.0378, -.5839],
           [.0325, -.3585, -.3103, 1.8521, -.2157],
           [-.0910, .1741, -.6267, 1.5870, -.0433],
           [-.0762, .0907, -.4273, .9707, .4421],
           [.1998, -.7758, -.4295, 2.4684, -.4629],
           [-.0880, .1712, -.5306, 1.0097, 0.4378],
           [.1578, -.8074, -.0923, 2.1966, -.4547],
           [.1992, -.1736, -1.8081, 3.7271, -.9446],
           [-.0813, -.1579, -.0551, 1.2108, .0836],
           [.1717, -.8713, -.1700, 2.8314, -.9618],
           [.2353, -.8325, -.6333, 2.8603, -.6298]])

T_hat_WMAP, wm_WMAP = ilc_map_from_weights(maps, wmap_ws, region, return_weights_map=True, sigma=1.5*np.pi/180)

ilcmap = load_wmap_maps(['data/wmap_ilc_9yr_v5.fits'])[0]


Sigma is 0.000000 arcmin (0.000000 rad) 
-> fwhm is 0.000000 arcmin
Sigma is 0.000000 arcmin (0.000000 rad) 
-> fwhm is 0.000000 arcmin
Sigma is 0.000000 arcmin (0.000000 rad) 
-> fwhm is 0.000000 arcmin
Sigma is 0.000000 arcmin (0.000000 rad) 
-> fwhm is 0.000000 arcmin
Sigma is 0.000000 arcmin (0.000000 rad) 
-> fwhm is 0.000000 arcmin
Sigma is 0.000000 arcmin (0.000000 rad) 
-> fwhm is 0.000000 arcmin
Sigma is 0.000000 arcmin (0.000000 rad) 
-> fwhm is 0.000000 arcmin
Sigma is 0.000000 arcmin (0.000000 rad) 
-> fwhm is 0.000000 arcmin
Sigma is 0.000000 arcmin (0.000000 rad) 
-> fwhm is 0.000000 arcmin
Sigma is 0.000000 arcmin (0.000000 rad) 
-> fwhm is 0.000000 arcmin
Sigma is 0.000000 arcmin (0.000000 rad) 
-> fwhm is 0.000000 arcmin
Sigma is 0.000000 arcmin (0.000000 rad) 
-> fwhm is 0.000000 arcmin
Sigma is 0.000000 arcmin (0.000000 rad) 
-> fwhm is 0.000000 arcmin
Sigma is 0.000000 arcmin (0.000000 rad) 
-> fwhm is 0.000000 arcmin
Sigma is 0.000000 arcmin (0.000000 rad) 
-> fwhm is 0.000000 arcmin
Sigma is 0.000000 arcmin (0.000000 rad) 
-> fwhm is 0.000000 arcmin
Sigma is 0.000000 arcmin (0.000000 rad) 
-> fwhm is 0.000000 arcmin
Sigma is 0.000000 arcmin (0.000000 rad) 
-> fwhm is 0.000000 arcmin
Sigma is 0.000000 arcmin (0.000000 rad) 
-> fwhm is 0.000000 arcmin
Sigma is 0.000000 arcmin (0.000000 rad) 
-> fwhm is 0.000000 arcmin
Sigma is 0.000000 arcmin (0.000000 rad) 
-> fwhm is 0.000000 arcmin
Sigma is 0.000000 arcmin (0.000000 rad) 
-> fwhm is 0.000000 arcmin
Sigma is 0.000000 arcmin (0.000000 rad) 
-> fwhm is 0.000000 arcmin
Sigma is 0.000000 arcmin (0.000000 rad) 
-> fwhm is 0.000000 arcmin

In [140]:
print "jlazear ILC weights"
linestr = "{0:^10} {1:^10} {2:^10} {3:^10} {4:^10} {5:^10}"
print linestr.format(*(['GROUP'] + bands))
for i in range(len(ws)):
    linelist = [i] + [str(w).ljust(8, '0')[:8] for w in ws[i]]
    print linestr.format(*linelist)
    
print "\nWMAP9 ILC weights"
print linestr.format(*(['GROUP'] + bands))
for i in range(len(wmap_ws)):
    linelist = [i] + [str(w).ljust(8, '0')[:8] for w in wmap_ws[i]]
    print linestr.format(*linelist)
    
wmin = min(ws.min(), wmap_ws.min())
wmax = max(ws.max(), wmap_ws.max())
dmin = (ws - wmap_ws).min()
dmax = (ws - wmap_ws).max()
dabs = max(abs(dmin), abs(dmax))


jlazear ILC weights
  GROUP        K          Ka         Q          V          W     
    0       0.155986   -0.76841   -0.24397   2.257494   -0.40108 
    1       0.055803   -0.64892   0.208869   1.945169   -0.56091 
    2       0.048720   -0.43987   -0.23540   1.851762   -0.22519 
    3       -0.09964   0.205965   -0.64152   1.560084   -0.02487 
    4       -0.07787   0.094598   -0.42622   0.967506   0.441998 
    5       0.209118   -0.80956   -0.41398   2.482037   -0.46760 
    6       -0.08915   0.169148   -0.51964   1.000231   0.439416 
    7       0.159301   -0.82604   -0.05642   2.157844   -0.43467 
    8       0.201858   -0.23001   -1.69478   3.628966   -0.90602 
    9       -0.09946   -0.10505   -0.05691   1.156837   0.104587 
    10      0.176710   -0.92922   -0.06826   2.746804   -0.92602 
    11      0.237338   -0.85279   -0.59759   2.827484   -0.61443 

WMAP9 ILC weights
  GROUP        K          Ka         Q          V          W     
    0       0.155500   -0.75720   -0.26890   2.284500   -0.41380 
    1       0.037500   -0.51370   0.022300   2.037800   -0.58390 
    2       0.032500   -0.35850   -0.31030   1.852100   -0.21570 
    3       -0.09100   0.174100   -0.62670   1.587000   -0.04330 
    4       -0.07620   0.090700   -0.42730   0.970700   0.442100 
    5       0.199800   -0.77580   -0.42950   2.468400   -0.46290 
    6       -0.08800   0.171200   -0.53060   1.009700   0.437800 
    7       0.157800   -0.80740   -0.09230   2.196600   -0.45470 
    8       0.199200   -0.17360   -1.80810   3.727100   -0.94460 
    9       -0.08130   -0.15790   -0.05510   1.210800   0.083600 
    10      0.171700   -0.87130   -0.17000   2.831400   -0.96180 
    11      0.235300   -0.83250   -0.63330   2.860300   -0.62980 

In [255]:
#hp.mollview(maps[0], nest=True)
mollret = hp.mollview(ilcmap, unit='mK', title='WMAP ILC',  norm='hist')
hp.mollview(T_hat_JL, unit='mK', title='jlazear ILC', min=-.452123, max=.426628)
hp.mollview((T_hat_JL - ilcmap), unit='mK', title='(jlazear - WMAP) ILC')
hp.mollview(T_hat_WMAP, unit='mK', title='WMAP ILC from coefficients')
hp.mollview(T_hat_WMAP - ilcmap, unit='mK', title='WMAP - WMAP')
delta = T_hat_WMAP - T_hat_JL
sigma_delta = delta.std()
mean_delta = delta.mean()
hp.mollview((T_hat_WMAP - T_hat_JL), unit='mK', title='WMAPr - JL', min=mean_delta-2*sigma_delta, max=mean_delta+2*sigma_delta)
print delta.min(), delta.mean(), delta.std(), delta.max()


-0.0530019504225 -7.11846264709e-05 0.000584945909012 0.0427570859055

In [141]:
nrows = 2
ncols = 3
gs = gridspec.GridSpec(nrows, ncols, height_ratios=[10, 1])
#fig, axes = subplots(nrows=2, ncols=3, height_ratios=[10, 1])
fig = figure()
axes = np.array([[subplot(gs[j, i]) for i in range(ncols)] for j in range(nrows)])

im0 = axes[0, 0].imshow(ws, interpolation='nearest', vmin=wmin, vmax=wmax)
im1 = axes[0, 1].imshow(wmap_ws, interpolation='nearest', vmin=wmin, vmax=wmax)
im2 = axes[0, 2].imshow(ws-wmap_ws, interpolation='nearest', vmin=-dabs, vmax=dabs)

cb0 = fig.colorbar(im0, cax=axes[1, 0], orientation='horizontal')
cb1 = fig.colorbar(im1, cax=axes[1, 1], orientation='horizontal')
cb2 = fig.colorbar(im2, cax=axes[1, 2], orientation='horizontal')

for cb in [cb0, cb1, cb2]:
    tick_locator = ticker.MaxNLocator(nbins=4)
    cb.locator = tick_locator
    cb.update_ticks()
    
for i, title in enumerate(['JL', 'WMAP', 'DIFF']):
    ax = axes[0, i]
    ax.set_title(title)



In [142]:
print wm_JL.max(), wm_JL.min()
hp.mollview(wm_JL, nest=False)
total_wmap = np.zeros_like(wm_JL)
for i in range(len(masks)):
    hp.mollview(masks[i], title='Region {0}'.format(i))
#    smmask0 = hp.smoothing(masks[i], sigma=1.5*np.pi/180., verbose=False)
#    hp.mollview(smmask0, title='smoothed {0}'.format(i))
#    total_wmap += smmask0
#    hp.mollview(total_wmap, title='Total {0}'.format(i))


1.0 1.0

Regions


In [3]:
fnameregn = 'data/wmap_ilc_rgn_defn_9yr_v5.fits'
fregnfits = fits.open(fnameregn)
fregndata = fregnfits[1]
bitmask = fregnfits[1].data['TEMPERATURE'].astype('uint16')
region = fregnfits[1].data['N_OBS']

regions = [(bitmask >> n) & 1 for n in range(12)]

In [127]:
import scipy.constants
h = scipy.constants.h
kB = scipy.constants.k
c = scipy.constants.c

def B(T, nu):
    return (2*h*nu**3/(c**2)) / (np.exp(h*nu/(kB*T)) - 1)

def Binv(I, nu):
    return (h*nu/kB) / np.log(1 + 2*h*nu**3/(I*c*c))

def Binvapprox(I, nu):
    return c*c*I/(2*kB*nu*nu)

bands = ['K', 'Ka', 'Q', 'V', 'W']
fnames = ['data/wmap_band_smth_imap_r9_9yr_{0}_v5.fits'.format(band) for band in bands]
maskdata = 'data/wmap_ilc_rgn_defn_9yr_v5.fits'

maps = load_wmap_maps(fnames, n2r=False)

Dust


In [178]:
reload(lib.foregrounds)
import lib.foregrounds

In [179]:
dust200 = lib.foregrounds.generate_simple_dust_map(200.e9, Nside=512)
hp.mollview(dust200, nest=True, norm='hist')
hp.npix2nside(len(dust200))


Out[179]:
512

In [139]:
nu0 = 353.e9
#print 2*h*nu0**3/(c*c)
#2*h*nu0**3/(dust353map*1.e6*c*c)
hp.mollview(Binv(dust353map/1.e20, nu0), nest=True, norm='hist', unit='K CMB')
hp.mollview(Binv(planck353map/1.e20, 353.e9), nest=True, norm='hist', unit='K CMB', remove_mono=True)
hp.mollview(maps[0], nest=True, norm='hist', unit='uK CMB')


monopole: 1.59035

In [143]:
planck353fname = 'data/HFI_SkyMap_545_2048_R1.10_nominal.fits'
planck353 = fits.open(planck353fname)
planck353map = planck353[1].data['I_STOKES']  # MJy/sr
hp.mollview(Binv(planck353map/1.e20, 353.e9), nest=True, norm='hist', unit='K CMB', remove_mono=True)
planck353map512 = hp.ud_grade(planck353map, nside_out=512, order_in='NESTED', order_out='NESTED')
hp.mollview(Binv(planck353map512/1.e20, 353.e9), nest=True, norm='hist', unit='K CMB', remove_mono=True)


monopole: 1.59035
monopole: 1.58848

In [113]:
dust353[1].header


Out[113]:
XTENSION= 'BINTABLE'           / binary table extension                         
BITPIX  =                    8 / 8-bit bytes                                    
NAXIS   =                    2 / 2-dimensional binary table                     
NAXIS1  =                   24 / width of table in bytes                        
NAXIS2  =               786432 / number of rows in table                        
PCOUNT  =                    0 / size of special data area                      
GCOUNT  =                    1 / one data group (required keyword)              
TFIELDS =                    6 / number of fields in each row                   
TTYPE1  = 'I       '           / label for field   1                            
TFORM1  = 'E       '           / data format of field: 4-byte REAL              
TUNIT1  = 'MJy/sr  '           / physical unit of field                         
TTYPE2  = 'I_stdev '           / label for field   2                            
TFORM2  = 'E       '           / data format of field: 4-byte REAL              
TUNIT2  = 'MJy/sr  '           / physical unit of field                         
TTYPE3  = 'Em      '           / label for field   3                            
TFORM3  = 'E       '           / data format of field: 4-byte REAL              
TUNIT3  = 'uK_CMB  '           / physical unit of field                         
TTYPE4  = 'Em_stdev'           / label for field   4                            
TFORM4  = 'E       '           / data format of field: 4-byte REAL              
TUNIT4  = 'none    '           / physical unit of field                         
TTYPE5  = 'T       '           / label for field   5                            
TFORM5  = 'E       '           / data format of field: 4-byte REAL              
TUNIT5  = 'uK_CMB  '           / physical unit of field                         
TTYPE6  = 'T_stdev '           / label for field   6                            
TFORM6  = 'E       '           / data format of field: 4-byte REAL              
TUNIT6  = 'uK_CMB  '           / physical unit of field                         
EXTNAME = 'COMP-MAP'                                                            
DATE    = '2013-02-13T13:31:15' / file creation date (YYYY-MM-DDThh:mm:ss UT)   
CHECKSUM= '9OAOJO7N9OANGO7N'   / HDU checksum updated 2013-02-13T13:31:15       
DATASUM = '4139938263'         / data unit checksum updated 2013-02-13T13:31:15 
COMMENT                                                                         
COMMENT *** Planck params ***                                                   
COMMENT                                                                         
PIXTYPE = 'HEALPIX '           / HEALPIX pixelisation                           
ORDERING= 'NESTED  '           / Pixel ordering scheme, either RING or NESTED   
NSIDE   =                  256 / Resolution parameter for HEALPIX               
FIRSTPIX=                    0 / First pixel # (0 based)                        
LASTPIX =               786431 / Last pixel # (0 based)                         
INDXSCHM= 'IMPLICIT'           / Indexing: IMPLICIT or EXPLICIT                 
COORDSYS= 'GALACTIC'                                                            
FILENAME= 'COM_CompMap_dust-commrul_0256_R1.00.fits'                            
COMMENT ---------------------------------------------------------------------   
COMMENT The intensity is normalized at 353 GHz                                  
COMMENT Object:                                                                 
COMMENT /sci_planck/lfi_dpc_test/ashdown/repository/compsep_outputs/commander-ru
COMMENT ler/delta_dx9/planck/dx9delta_v1_dust_amp_mean.fits                     
COMMENT /sci_planck/lfi_dpc_test/ashdown/repository/compsep_outputs/commander-ru
COMMENT ler/delta_dx9/planck/dx9delta_v1_dust_amp_stddev.fits                   
COMMENT /sci_planck/lfi_dpc_test/ashdown/repository/compsep_outputs/commander-ru
COMMENT ler/delta_dx9/planck/dx9delta_v1_dust_em_mean.fits                      
COMMENT /sci_planck/lfi_dpc_test/ashdown/repository/compsep_outputs/commander-ru
COMMENT ler/delta_dx9/planck/dx9delta_v1_dust_em_stddev.fits                    
COMMENT /sci_planck/lfi_dpc_test/ashdown/repository/compsep_outputs/commander-ru
COMMENT ler/delta_dx9/planck/dx9delta_v1_dust_T_mean.fits                       
COMMENT /sci_planck/lfi_dpc_test/ashdown/repository/compsep_outputs/commander-ru
COMMENT ler/delta_dx9/planck/dx9delta_v1_dust_T_stddev.fits                     
COMMENT ---------------------------------------------------------------------   

In [222]:
reload(lib.wmap)
import lib.wmap

mcmcw = fits.open('data/wmap_mcmc_w_dust_stk_q_7yr_v4p1.fits')

hp.mollview(mcmcw[1].data['BESTFIT'], nest=True)



In [251]:
fnames = ['docs/img/Q.png', 'docs/img/U.png']
all(map(os.path.isfile, fnames))


Out[251]:
True

In [246]:
qmap = lib.wmap.load_wmap_maps_QU(['data/wmap_mcmc_k_synch_stk_q_7yr_v4p1.fits',
                                   'data/wmap_mcmc_k_synch_stk_u_7yr_v4p1.fits'], 
                                  Nside=256, n2r=True)

qmap[0] = hp.smoothing(qmap[0], sigma=3*np.pi/180, verbose=False)
qmap[1] = hp.smoothing(qmap[1], sigma=3*np.pi/180, verbose=False)


hp.mollview(qmap[0], nest=False, title='WMAP 23 GHz sync Q map', norm='hist')
hp.mollview(qmap[1], nest=False, title='WMAP 23 GHz sync U map', norm='hist')
hp.mollview(0.5*np.arctan2(-qmap[1], qmap[0]), nest=False)


Sigma is 0.000000 arcmin (0.000000 rad) 
-> fwhm is 0.000000 arcmin
Sigma is 0.000000 arcmin (0.000000 rad) 
-> fwhm is 0.000000 arcmin

In [196]:
ls data


COM_CompMap_dust-commrul_0256_R1.00.fits
NOTES.md
get_data.py
wmap_band_smth_imap_r9_9yr_K_v5.fits
wmap_band_smth_imap_r9_9yr_Ka_v5.fits
wmap_band_smth_imap_r9_9yr_Q_v5.fits
wmap_band_smth_imap_r9_9yr_V_v5.fits
wmap_band_smth_imap_r9_9yr_W_v5.fits
wmap_band_smth_imap_r9_nineyear_v5.tar.gz
wmap_ilc_9yr_v5.fits
wmap_ilc_rgn_defn_9yr_v5.fits
wmap_mcmc_w_dust_stk_q_7yr_v4p1.fits
wmap_mcmc_w_dust_stk_u_7yr_v4p1.fits

In [191]:
reload(lib.foregrounds)
import lib.foregrounds

QUmaps = lib.foregrounds.generate_simple_dust_QU_map()

hp.mollview(QUmaps[0], nest=True, norm='hist', title='353 GHz Q map realization, p = 0.2', unit='MJy/sr')
hp.mollview(QUmaps[1], nest=True, norm='hist', title='353 GHz U map realization, p = 0.2', unit='MJy/sr')
hp.mollview(QUmaps[1] - QUmaps[0], nest=True, norm='hist', 
            title='353 GHz Q-U map realization, p = 0.2', unit='MJy/sr')
hp.mollview(np.sqrt(QUmaps[0]**2 + QUmaps[1]**2), nest=True, norm='hist', unit='MJy/sr',
            title='353 GHz I_pol map, p=0.2')
hp.mollview(np.sqrt(QUmaps[0]**2 + QUmaps[1]**2) - 0.2*lib.foregrounds.generate_simple_dust_map(),
            nest=True, norm='hist', title='353 GHz I_pol-I_pol map, p=0.2', unit='MJy/sr')



In [147]:
def clean_dust_map(map):
    """Replaces non-positive pixels with smallest positive pixel value."""
    pos = map[map > 0]
    posmin = pos.min()
    map[map <= 0] = 0.# posmin
    return map

dust353fname = 'data/COM_CompMap_dust-commrul_0256_R1.00.fits'
dust353 = fits.open(dust353fname)
dust353map = dust353[1].data['I']
dust353map = clean_dust_map(dust353map)
dust353em = dust353[1].data['Em']
dust353T = dust353[1].data['T']

hp.mollview(dust353map, nest=True, norm='hist', unit='MJy/sr', title='')
#savefig('planck_dust_353GHz.png')
dust353mapT = Binv(dust353map*dust353em/1.e20, 353.e9)
hp.mollview(dust353mapT, nest=True, norm='hist', unit='K CMB', title='')
hp.mollview(dust353em, nest=True, norm='hist', unit='unitless', title='emissivity')
hp.mollview(dust353em*dust353map, nest=True, norm='hist', unit='MJy/sr', title='emissivity*I')
hp.mollview(dust353T, nest=True, norm='hist', unit='uK CMB', title='T')

#dust353[1].header



In [149]:
dustmodelfname = 'data/HFI_CompMap_ThermalDustModel_2048_R1.20.fits'
dustmodel = fits.open(dustmodelfname)
dustmodeltau = dustmodel[1].data['TAU353']
dustmodeltemp = dustmodel[1].data['TEMP']

dustmodelI = B(dustmodeltemp, 353.e9)*dustmodeltau*1.e20  # MJy/sr
hp.mollview(dustmodeltau, nest=True, unit='unitless', norm='hist')
hp.mollview(dustmodeltemp, nest=True, unit='K CMB', norm='hist')
hp.mollview(dustmodelI, unit='MJy/sr', nest=True, norm='hist')
hp.mollview(Binv(dustmodelI/1.e20, 353.e9), unit='K CMB', norm='hist', nest=True)



In [120]:
hp.mollview(region, nest=True)



In [115]:
summap = sum(regions[i]*(i+1) for i in range(len(regions)))
hp.mollview(summap, nest=True, title='Groups')