MADCAP 3.0

MADCAP 3.0 has been completely re-written to enable the analysis of

(i) polarization data, in the form of i- and/or q- and u-maps
(ii) data from multiple uncorrelated detectors, including different beam sizes

In addition, it now exploits an additional level of parallelism to decrease wallclock runtime, and includes simple user interface codes to launch successful massively parallel jobs on NERSC's IBM SP seaborg .

MADCAP 3.0 implements maximum likelihood algorithms for both map making & noise correlating and power spectrum estimation.

Map Making & Noise Correlating

Since maps are an important diagnostic for systematic effects in a dataset, and often need to be recalculated several times before a power spectrum analysis is undertaken, we recommend that they are initially made using MADmap - a fast preconditioned conjugate gradient map-maker. In general, only once such tests have been completed should the much more computationally expensive step of calculating the full pixel-pixel noise correlation matrix be performed. MADCAP also offers its own conjugate-gradient map-making routing, but uses the explicitly-constructed inverse pixel-pixel noise correlation matrix rather than the faster Fourier domain methods.

The map-making algorithms used by MADCAP and MADmap use the same data inputs. In both cases we assume that the time-ordered data comprises piecewise stationary instrument noise and one or more signal components each of which is discretized in a unique, non-temporal, basis. For example

d_t = n_t + A_tp s_p + B_ti x_i

where	s_p	is the sky-signal (intensity and/or polarization) discretized into pixels.
	A_tp	is the associated pointing matrix, indicating which pixels were observed, and with what weights, at each time.
	x_i	is some other signal in the timestream (neither time- nor sky-synchronous) discretized in its i-basis.
	B_ti	is the associated pointing matrix, indicating which i-elements were observed, and with what weights, at each time.

Combining the pixel- and other element-numbering into a single metapixelization and composite metapointing matrix, the map-making inputs are:

(i) one or more unformatted binary time-ordered data files containing

initial observation number (8 byte integer)
final observation number (8 byte integer)
for each observation
       data (8 byte float)

(ii) one or more unformatted binary pointing files containing
initial observation number (8 byte integer)
final observation number (8 byte integer)
number of pixels in the pointing class (8 byte integer)
number of signal components (of all types) in any observation (8 byte integer)
for each observation
       for each signal component in the observation
                     pixel weight (4 byte float)
                     pixel index (4 byte integer)
(iii) one or more unformatted binary piecewise stationary inverse time-time noise correlation function files each containing

              initial applicable observation number (8 byte integer)
              final applicable observation number (8 byte integer)
            maximum time-time correlation length (white noise = 0) (8 byte integer)
              for each separation
                       correlation (8 byte float)

Once these data files are in place, run the code

noise.x $XTT_MAX $BLOCKSIZE $PCG $IMAX $EPSILON

where	XTT_MAX	is the maximum time-time correlation length to allow (overriding the file value if necessary).
	QBLOCKSIZE	is the optimal ScaLAPACK blocksize for the particular computer.
	PCG	is a flag for using the PCG solver for the map (and not inverting the inverse pixel-pixel noise matrix).
	IMAX	is the maximum number of iterations to use in if the PCG option is set.
	EPSILON	is the required PCG precision.

On seaborg this can most easily be achieved by running ~borrill/MADCAP/sea_launch_nc.x which will check the input data file consistency and resource availability, determine the command-line parameters, estimate the wallclock runtime, and generate and submit an appropriate parallel job batch script.

This will generate firstly unformatted binary files (i) of 4-byte floats inv_qq_noise, inv_qq_noise.diag & inv_qq_noise.q_data and of 4-byte ints qhits, and (ii) either of 4-byte floats q_data.pcg if the PCG option is chosen or of 4-byte floats qq_noise & q_data otherwise.

Note that all the pixel-domain output files except qhits only include data for pixels that are included in the observation, although the pixel numbering can (and should) be global.

Power Spectrum Estimation

Once the quality of the map(s) is satisfactory, and approriate cuts have been taken (for example to excise bright point sources, or galactic foregrounds, or very low signal-to-noise pixels), MADCAP can estimate up to 6 power spectra - TT, EE, BB, TE, TB & EB, which are labelled T, E, B, X, Y & Z respectively - using Newton-Raphson iteration to locate the peak of the Gaussian likelihood function.

For each dataset the power spectrum estimation inputs are (assuming N_p total pixels, with N_i i-pixels, N_q q-pixels & N_u u-pixels)

(i) an unformatted binary file of N_p 4-byte floats p_data.{n} containing the concatenated i- , q- and u-maps in dataset n.
(ii) an unformatted binary file of (N_p x N_p) 4-byte floats pp_noise.{n} containing the {i,q,u} pixel-pixel noise correlations in dataset n.

A schematic representation of the p_data and pp_noise files for a dataset with i, q & u pixels. Note that the number of pixels of each type need not be the same - in this example the higher signal-to-noise intensity data is pixelized twice as finely as the polarization.
This is the data format generated by MADmap and MADCAP if the i, q & u pixels are numbered consecutively in the t_data metapixelization.

(iii) unformatted binary files of 2 x N_i , N_q & N_u 4-byte floats pixels_i.{n} , pixels_q.{n} & pixels_u.{n} containing the (RA, Dec) coordinates in degrees of the i-, q- & u-pixels respectively in dataset n.

(iv) optional: an unformatted binary file of (N_t x N_p) 4-byte floats p_templates.{n} containing N_t pixel-templates to be marginalized over.

(v) ascii files of N_l floats mwindow_i.{n} , mwindow_q.{n} & mwindow_u.{n} containing the convolved beam+pixel multipole window functions for the i-, q- & u-maps respectively in dataset n.

Given the impact of incomplete sky coverage and finite pixel resolution, the spectra are estimated in bins of multipoles, where for each spectral type

C_l = C_b C_shape(l)

The remaining input files, common to all the datasets, are
(vi) ascii files of N_l floats C{T,E,B,X,Y,Z}_shape containing the T, E, B, X, Y & Z spectral multipole shape functions respectively.

(vii) ascii files of 2 x N_b integers bin{T,E,B,X,Y,Z} containing the minimum and maximum multipoles in each bin of the T,E,B,X,Y & Z spectra respectively.

(viii) optional: an ascii file of N_b floats C_0 giving the initial values of the bin powers for all of the spectra concatenated in T,E,B,X,Y,Z order. If this file is omitted then the initial bin powers are all assumed to be zero, enabling the calculation of the offset lognormal approximation to the bin error bars.

Many of the above files may not be needed for a particular analysis, in which case they are simply omitted. For example, a dataset with CMB temperature only would not use any of the polarization files. If a class of spectral bin file is omitted, that spectrum is taken to be fixed at zero.

Once these data files are in place, run the code

spectrum.x NO_GANG BLOCKSIZE CONVERGENCE STEP

twice for each Newton-Raphson iteration.

where	NO_GANG	is the number of independent work gangs to divide the processors into.
	BLOCKSIZE	is the optimal ScaLAPACK blocksize for the particular computer.
	CONVERGENCE	is the convergence threshold, where convergence is deemed to have occured when, for every bin, the change in the bin power is less than CONVERGENCE times the standard deviation of the bin power estimated from the diagonal elements of the inverse fisher matrix.
	STEP	is either 0 or 1, denoting which phase of the code is to be executed; step 0 sould be completed (marked by the generation of the file done_step0) before step 1 is initiated.

On seaborg this can most easily be achieved by running ~borrill/MADCAP/sea_launch_ps.x which will check the input data file consistency and resource availability, determine the command-line parameters, estimate the wallclock runtime, and generate and submit an appropriate parallel job batch script.