CompSep.jl

A julia implementation of parametrized component separation algorithm for CMB experiment, following fgbuster which is a python implementation – sorry I wish it is actually about a super octo rotary phone.

Currently, the code is able to produce the same result as fgbuster, despite being significantly shorter.

Now that I have had enough fun writing obscure codes, it’s time to rewrite it for performance and maintainability, so help is welcomed! For a single pixel search, it is roughly 100 times faster than fgbuster. For a multi-resolution search with nside=8, it is ~300-400 times faster than fgbuster.

@btime fgbuster.basic_comp_sep($components, $instrument, $obs; tol=1);

# result: 5.679 s (699 allocations: 23.70 MiB)

@btime res = compsep($comps, $freqs, $Nmat.^-2, $obs, x₀=[-3,1.5,15.]);

# result: 60.788 ms (420728 allocations: 32.63 MiB)

Examples

Basic

using PyCall

# obtain simulated data using fgbuster
@pyimport fgbuster
sky = fgbuster.observation_helpers.get_sky(nside=32, tag="c1d0s0")
instrument = fgbuster.observation_helpers.get_instrument("LiteBIRD");
obs = fgbuster.observation_helpers.get_observation(instrument, sky, noise=true);
freqs = instrument.frequency.values
Nmat = hcat(instrument.depth_i.values,instrument.depth_p.values,instrument.depth_p.values).^-2

# define components of interests, each of which is a function of frequency and other predefined parameters
comps = [cmb, sync, dust]

# perform simple component separation
res = compsep(comps, freqs, Nmat, obs, x₀=[-3,1.54,20.])

Multi-resolution removal

# same as before
nside = 8
res = compsep(comps, freqs, Nmat, obs, nside; x₀=[-3,1.54,20.])

bandpass integration

# same as before
bands = SimplePassband.(instrument.frequency.values, instrument.bandwidth.values)
res = compsep(comps, bands, Nmat, obs, x₀=[-3,1.54,20.])

Next steps

• [X] allow passing in mask
• [X] allow multi-resolution
• [X] bandpass integration
• [X] add functions to postprocess results
• [X] add some tests
• [ ] more documentation