gistfile1.txt

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

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

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


#
# Data
#

observations_by_name = {}
observations = []
for name, sigma0 in [('K1', '1.437 mK'),
                     ('Ka1', '1.470 mK'),
                     ('Q1', '2.254 mK'),
                     ('Q2', '2.140 mK'),
                     ('V1', '3.319 mK'),
                     ('V2', '2.955 mK'),
                     ('W1', '5.906 mK'),
                     ('W2', '6.572 mK'),
                     ('W3', '6.941 mK'),
                     ('W4', '6.778 mK')]:
    obs = cm.SkyObservation(
        name=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=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)
    observations_by_name[name] = obs

    if name.endswith('1'):
        # TODO: Add averaging over DAs; for now just take the first one
        observations.append(obs)

K, Ka, Q, V, W = observations


#
# Model
#

# map map to index of mixing matrix map to use
def get_mixing_maps(comp_num):
    obs_to_idx = {'K1' : 1,
                  'Ka1' : 2,
                  'Q1' : 3, 'Q2' : 3,
                  'V1' : 4, 'V2' : 4,
                  'W1' : 5, 'W2' : 5, 'W4' : 5, 'W4' : 5}
    return dict((observations_by_name[name],
                 ('$DATA/wmap/mixing/mixmat_comp%02d_band%02d_k03999.fits' % (comp_num, obsidx),
                  1, 0))
                for name, obsidx in obs_to_idx.iteritems())

# 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=[K, Q])

signal_components = [cmb, dust, synchrotron]

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')