Simple component separation

%load_ext autoreload
%autoreload 2
%matplotlib widget

import pandas

import healpy as hp
import pysm

from fgbuster import (CMB, Dust, Synchrotron,
                      basic_comp_sep,
                      get_observation, get_noise_realization, get_instrument, get_sky)
from fgbuster.segmentation import define_ids

import matplotlib.pyplot as plt

Simulate your frequency maps

You have some frequecy maps to clean, they can be either data or simulations. For example, you can use pysm to simulate microwave skys with different complexities. ForeGroundBuster adds a couple of functions that make the process even easier.

instrument = get_instrument('LiteBIRD')
instrument.depth_p

0     37.5
1     24.0
2     19.9
3     16.2
4     13.5
5     11.7
6      9.2
7      7.6
8      5.9
9      6.5
10     5.8
11     7.7
12    13.2
13    19.5
14    37.5
Name: depth_p, dtype: float64

nside = 16
sky = get_sky(nside, 'd1')

sky.components[0].mbb_temperature.value[:] = 20

noiseless_maps = get_observation(instrument, sky)[:, 1:]
noisy_maps = get_observation(instrument, sky, noise=True)[:, 1:]

Define what you fit for

Create your instrument model as a list of components

components = [CMB(), Synchrotron(40., running=None, nu_pivot=70.)]

components = [CMB(), Dust(150.)]

components = [CMB(), Dust(150.), Synchrotron(40., running=None, nu_pivot=70.)]

ids, nsides, df = define_ids(components, instrument, noisy_maps, 25,
                             #bounds = [(None, None)]*2,
                             method='TNC',
                             bounds=[(1.0, 1.9), (13, 30)],
                             options={ 'disp':False, 'gtol': 1e-12, 'eps': 1e-12,
                                        'maxiter': 10000, 'ftol': 1e-12 },
                             tol=1e-18
                             #, (-3.5, -2.5), (-0.01, 0.01)])
)

1
New round
New round
New round
2
New round
New round
New round
4
New round
New round
New round
8
New round
New round
New round
16
New round
New round
New round

idf = df.set_index(['nside', 'round'])

fig, axs = plt.subplots(4, 1, sharex=True)
nside = 32
i = 0
par = 2
while True:
    print(i)
    try:
        res = idf.loc[(nside, i)]
        axs[0].plot(res.ids[par])
        axs[1].plot(res.nsides[par])
        axs[2].plot(res.significance[par])
        axs[3].plot(res.significance[par] > 1)
    except TypeError:
        break
    i += 1
axs[2].set_ylim(0.1, None)
axs[2].set_yscale('log')
axs[2].axhline(1, color='k')

0
1
2
3
4
5

<matplotlib.lines.Line2D at 0x11dab2d50>

from fgbuster import MixingMatrix

MixingMatrix(*components).params

['Dust.beta_d', 'Dust.temp']

res = idf.iloc[-1]

import numpy as np

fig = plt.figure()
n_row = (len(nsides)  - 1) / 2 + 1
n_col = 2 if (len(nsides) > 1) else 1
for i, (nsides_i, par) in enumerate(zip(nsides, MixingMatrix(*components).params), 1):
    hp.mollview(np.log10(res.chi2), sub=(n_row, n_col, i), title=par)

idf.to_pickle('nside16_d1n_constTd_bd_td_tr25')

ids_nl, nsides_nl, df_nl = define_ids(components, instrument, noiseless_maps, 5,
                             #bounds = [(None, None)]*2,
                             method='pso',
                             bounds=[(1.0, 1.9), (13, 30)]#, (-3.5, -2.5), (-0.01, 0.01)])
)

1
New round
New round
New round
2
New round
New round
New round
New round
New round
New round
New round
New round
4
New round
New round
New round
New round
New round
New round
New round
New round
New round
New round
New round
8
New round
New round
New round
New round
New round
New round
New round
New round
New round
New round
New round
New round
New round
16
New round
New round
New round
New round
New round
New round
New round
New round
New round
New round
New round
New round
New round
New round
New round
New round
New round
New round
New round

ids_ny, nsides_ny, df_ny = define_ids(components, instrument, noisy_maps, 5,
                             #bounds = [(None, None)]*2,
                             method='pso',
                             bounds=[(1.0, 1.9), (13, 30)]#, (-3.5, -2.5), (-0.01, 0.01)])
)

res.significance.dtype

dtype('int64')

Component separation

The tools inside ForeGroundBuster allow for very flexible and diverse component separation approaches. However, we also provide a set of predefined function that perform component separation out of the box. They suit most common use cases.

result = basic_comp_sep(components, instrument, freq_maps)

Explore the results

You have just solved for both the spectral parameters of your components and their amplitudes.

Get the spectral parameters name and values with

print(result.params)
print(result.x)

You can also have a look at their semi-analytic covariance, but remember that it is accurate only the high signal-to-noise regime

corner_norm(result.x, result.Sigma, labels=result.params)

The amplitudes of the components are stacked in the s attribute and they are in the same format of the input frequency maps: Q and U healpix maps, in this case.

hp.mollview(result.s[0,1], title='CMB')
hp.mollview(result.s[1,1], title='Dust', norm='hist')
hp.mollview(result.s[2,1], title='Synchrotron', norm='hist')