gistfile1.py

from __future__ import division
import numpy as np
import os
import commander as cm

#
# Parameters (just variables that are reused further down)\
#

output_filename = 'newmonodip-em7-%d.h5' % os.getpid()
lmax = 500
lprecond = 50
eps = 1e-7
cache_path = 'cache'
seed = None # integer, or None for getting random seed from OS


#
# Data
#

observations = []
for name, da_list in [('K', [('K1', '1.437 mK')]),
                      ('Ka', [('Ka1', '1.470 mK')]),
                      ('Q', [('Q1', '2.254 mK'), ('Q2', '2.140 mK')]),
                      ('V', [('V1', '3.319 mK'), ('V2', '2.955 mK')]),
                      ('W', [('W1', '5.906 mK'), ('W2', '6.572 mK'),
                             ('W3', '6.941 mK'), ('W4', '6.778 mK')])]:
    obs = cm.AverageSkyObservations(name, [
        cm.SkyObservation(
            name=da_name,
            description='Raw WMAP r9 data, from http://lambda.gsfc.nasa.gov',
            T_map=('$DATA/wmap/n512/wmap_da_imap_r9_7yr_%s_v4.fits' % name, 1, 'TEMPERATURE'),
            TT_nobs_map=('$DATA/wmap/n512/wmap_da_imap_r9_7yr_%s_v4.fits' % name, 1, 'N_OBS'),
            TT_sigma0=da_sigma0,
            beam_transfer='$DATA/wmap/n512/wmap_%s_ampl_bl_7yr_v4.txt' % name,
            mask_map=('$DATA/wmap/n512/wmap_processing_mask_r9_7yr_v4.fits', 1, 0),
            lmax_beam=lmax)
        for da_name, da_sigma0 in da_list])

K, Ka, Q, V, W = observations


#
# Model
#

def get_mixing_maps(comp_num):
    # Create a mapping from observation to Commander 1 output file
    d = {}
    for idx, obs in enumerate(observations):
        fitsfile = '$DATA/wmap/mixing/mixmat_comp%02d_band%02d_k03999.fits' % (comp_num, idx + 1)
        d[obs] = (fitsfile, 1, 0)
    return d

# Priors for synch/dust
ls = np.arange(lmax + 1)
ls[0] = ls[1]
prior_synch = 3e5 * ls**-2.1
prior_dust = 1e3 * ls**-1.8
prior_synch[:2] = 0
prior_dust[:2] = 0

# Construct model
cmb = cm.IsotropicGaussianCmbSignal(
    name='cmb',
    power_spectrum='$DATA/wmap/wmap_lcdm_sz_lens_wmap5_cl_v3.fits',
    lmax=lmax
    )

dust = cm.MixingMatrixSignal(
    name='dust',
    lmax=lmax,
    power_spectrum=prior_dust,
    mixing_maps=get_mixing_maps(1)
    )

synchrotron = cm.MixingMatrixSignal(
    name='synch',
    lmax=lmax,
    power_spectrum=prior_synch,
    mixing_maps=get_mixing_maps(2)
    )

monodipole = cm.MonoAndDipoleSignal('monodipole', fix_at=[])

signal_components = [cmb, dust, synchrotron, monodipole]

ctx = cm.CommanderContext(cache_path=cache_path, seed=seed)

realizations = cm.app.constrained_realization(
    ctx,
    output_filename,
    observations,
    signal_components,
    wiener_only=False, # want samples
    lprecond=lprecond,
    eps=eps,
    filemode='w')