Computing the power spectrum¶

In this notebook, I will show how to use the Mexican-hat filtering from Arévalo & al. to compute the power spectra of images with arbitrary masks such as excluded point sources. The cell below is the mandatory data loading.

In [1]:

import astropy.units as u

from xsb_fluc.data.cluster import Cluster

cluster = Cluster(
    imglink='data/A2142/mosaic_a2142.fits.gz',
    explink='data/A2142/mosaic_a2142_expo.fits.gz',
    bkglink='data/A2142/mosaic_a2142_bkg.fits.gz',
    reglink='data/A2142/src_ps.reg',
    nhlink='data/A2142/A2142_nh.fits',
    ra=239.58615,
    dec=27.229433,
    r_500=1.403*u.Mpc, 
    redshift=0.09,
)

cluster = cluster.reduce_to_r500(1).rebin(5)
cluster.plot_cluster()

import astropy.units as u from xsb_fluc.data.cluster import Cluster cluster = Cluster( imglink='data/A2142/mosaic_a2142.fits.gz', explink='data/A2142/mosaic_a2142_expo.fits.gz', bkglink='data/A2142/mosaic_a2142_bkg.fits.gz', reglink='data/A2142/src_ps.reg', nhlink='data/A2142/A2142_nh.fits', ra=239.58615, dec=27.229433, r_500=1.403*u.Mpc, redshift=0.09, ) cluster = cluster.reduce_to_r500(1).rebin(5) cluster.plot_cluster()

WARNING: FITSFixedWarning: RADECSYS= 'FK5 ' / Stellar reference frame 
the RADECSYS keyword is deprecated, use RADESYSa. [astropy.wcs.wcs]
WARNING: FITSFixedWarning: EQUINOX = '2000.0 ' / Coordinate system equinox 
a floating-point value was expected. [astropy.wcs.wcs]

In this cell, I load the previously determined best-fit parameters from the MCMC notebook. There is also a bit of formatting to extract the Maximum A Posteriori and build an appropriate dictionary that can be fed to haiku functions.

In [2]:

import haiku as hk
import json
import jax.numpy as jnp
from jax.tree_util import tree_map
from xsb_fluc.simulation.mock_image import MockXrayCountsBetaModel

# This is the function we used to fit the X-ray surface brightness model
images_simulator = hk.without_apply_rng(hk.transform(lambda : MockXrayCountsBetaModel(cluster)()))

# We reload the parameters and parse them to the appropriate format
with open('data/A2142/posterior_parameters.json', 'r') as file:
    posterior_parameters = {key:jnp.asarray(value) for key, value in json.load(file).items()}
    # I have to manually set the background here because json encoding is buggy 
    posterior_parameters['log_bkg'] = jnp.asarray(-50.)

# We compute the median of the posterior samples, which is coincident with the 
# MAP estimate when the posterior distributions are constrained
posterior_median = tree_map(lambda x: jnp.median(x), posterior_parameters)

# Here we change the structure of the dictionary so it is compliant with haiku
hk_posterior_median = hk.data_structures.to_haiku_dict(images_simulator.init(None))

for module, parameter, prior in hk.data_structures.traverse(hk_posterior_median):
    
    if parameter in list(posterior_parameters.keys()):
        hk_posterior_median[module][parameter] = posterior_median[parameter]
        
posterior_image = images_simulator.apply(hk_posterior_median)

import haiku as hk import json import jax.numpy as jnp from jax.tree_util import tree_map from xsb_fluc.simulation.mock_image import MockXrayCountsBetaModel # This is the function we used to fit the X-ray surface brightness model images_simulator = hk.without_apply_rng(hk.transform(lambda : MockXrayCountsBetaModel(cluster)())) # We reload the parameters and parse them to the appropriate format with open('data/A2142/posterior_parameters.json', 'r') as file: posterior_parameters = {key:jnp.asarray(value) for key, value in json.load(file).items()} # I have to manually set the background here because json encoding is buggy posterior_parameters['log_bkg'] = jnp.asarray(-50.) # We compute the median of the posterior samples, which is coincident with the # MAP estimate when the posterior distributions are constrained posterior_median = tree_map(lambda x: jnp.median(x), posterior_parameters) # Here we change the structure of the dictionary so it is compliant with haiku hk_posterior_median = hk.data_structures.to_haiku_dict(images_simulator.init(None)) for module, parameter, prior in hk.data_structures.traverse(hk_posterior_median): if parameter in list(posterior_parameters.keys()): hk_posterior_median[module][parameter] = posterior_median[parameter] posterior_image = images_simulator.apply(hk_posterior_median)

Using these best-fit parameters, we can compute a single fluctuation map as done before, and the associated mask which is a combination between the excluded regions, the CCD gaps and the $R_{500}$ region.

In [3]:

import matplotlib.pyplot as plt 
import cmasher as cmr 
from matplotlib.colors import SymLogNorm

mask = (cluster.exp>0) & (cluster.center.separation(cluster.coords) < cluster.t_500)
fluc = jnp.nan_to_num((cluster.img - posterior_image)/(2*cluster.exp))

fig, axs = plt.subplots(nrows=1, ncols=2, subplot_kw={'projection':cluster.wcs})
map_fluc = axs[0].imshow(jnp.where(mask, fluc, jnp.nan), norm=SymLogNorm(vmin=-5e-6, vmax=5e-6, linthresh=1e-7), cmap=cmr.guppy)
map_mask = axs[1].imshow(mask)

plt.colorbar(map_fluc, ax=axs[0], location='bottom', label='Fluctuations')
plt.colorbar(map_mask, ax=axs[1], location='bottom', label='Mask');

import matplotlib.pyplot as plt import cmasher as cmr from matplotlib.colors import SymLogNorm mask = (cluster.exp>0) & (cluster.center.separation(cluster.coords) < cluster.t_500) fluc = jnp.nan_to_num((cluster.img - posterior_image)/(2*cluster.exp)) fig, axs = plt.subplots(nrows=1, ncols=2, subplot_kw={'projection':cluster.wcs}) map_fluc = axs[0].imshow(jnp.where(mask, fluc, jnp.nan), norm=SymLogNorm(vmin=-5e-6, vmax=5e-6, linthresh=1e-7), cmap=cmr.guppy) map_mask = axs[1].imshow(mask) plt.colorbar(map_fluc, ax=axs[0], location='bottom', label='Fluctuations') plt.colorbar(map_mask, ax=axs[1], location='bottom', label='Mask');

With this mask and image, we can define some scales of interest (in units of $R_{500}$) and compute the associated power spectrum using the Mexican-Hat filtering implemented in this package. This method is bound to the PowerSpectrum operator. Below I show the power spectrum associated to the absolute fluctuation map of A2142.

In [4]:

import jax.numpy as jnp
import haiku as hk
from xsb_fluc.simulation.power_spectrum import PowerSpectrum

scales = jnp.geomspace(0.05, 1, 20)

# Note that you can also fix the scales in the lambda function
power_spectrum = hk.without_apply_rng(
    hk.transform(
        lambda img, s: PowerSpectrum(cluster, mask=mask)(img, s)
    )
)

ps_fluc = power_spectrum.apply(None, fluc, scales)

plt.plot(scales, ps_fluc)
plt.xlabel('Scales [$R_{500}$]')
plt.ylabel('Power spectrum')
plt.loglog();

import jax.numpy as jnp import haiku as hk from xsb_fluc.simulation.power_spectrum import PowerSpectrum scales = jnp.geomspace(0.05, 1, 20) # Note that you can also fix the scales in the lambda function power_spectrum = hk.without_apply_rng( hk.transform( lambda img, s: PowerSpectrum(cluster, mask=mask)(img, s) ) ) ps_fluc = power_spectrum.apply(None, fluc, scales) plt.plot(scales, ps_fluc) plt.xlabel('Scales [$R_{500}$]') plt.ylabel('Power spectrum') plt.loglog();

Note that this mask can be used to compute power spectra in various regions of the fluctuation map. In the following cell, I compute the power spectrum in 4 concentric annuli.

In [5]:

import jax.numpy as jnp
import haiku as hk
import cmasher as cmr 
from xsb_fluc.simulation.power_spectrum import PowerSpectrum

scales = jnp.geomspace(0.05, 1, 20)
# Separation from the Planck center in units of R500
distance = jnp.asarray((cluster.center.separation(cluster.coords)/cluster.t_500).to(u.dimensionless_unscaled))

fig, axs = plt.subplots(nrows=1, ncols=2, figsize=(12, 5))

mask_img = jnp.ones_like(distance)*jnp.nan

for low_rad, high_rad in [(0, 0.25), (0.25, 0.5), (0.5, 0.75), (0.75, 1.)]:

    # Note that you can also fix the scales in the lambda function
    local_mask = (low_rad <= distance) & (distance< high_rad) & mask
    mask_img = mask_img.at[local_mask].set(high_rad) 
    
    power_spectrum = hk.without_apply_rng(
        hk.transform(
            lambda img, s: PowerSpectrum(cluster, mask=local_mask)(img, s)
        )
    )

    ps_fluc = power_spectrum.apply(None, fluc, scales)
    
    axs[0].plot(scales, ps_fluc, color=cmr.amber_r(high_rad))

axs[0].loglog()
axs[0].set_xlabel('Scales [$R_{500}$]')
axs[0].set_ylabel('Power spectrum')
axs[1].imshow(mask_img, cmap=cmr.ember_r)
axs[1].set_axis_off()
plt.show();

import jax.numpy as jnp import haiku as hk import cmasher as cmr from xsb_fluc.simulation.power_spectrum import PowerSpectrum scales = jnp.geomspace(0.05, 1, 20) # Separation from the Planck center in units of R500 distance = jnp.asarray((cluster.center.separation(cluster.coords)/cluster.t_500).to(u.dimensionless_unscaled)) fig, axs = plt.subplots(nrows=1, ncols=2, figsize=(12, 5)) mask_img = jnp.ones_like(distance)*jnp.nan for low_rad, high_rad in [(0, 0.25), (0.25, 0.5), (0.5, 0.75), (0.75, 1.)]: # Note that you can also fix the scales in the lambda function local_mask = (low_rad <= distance) & (distance< high_rad) & mask mask_img = mask_img.at[local_mask].set(high_rad) power_spectrum = hk.without_apply_rng( hk.transform( lambda img, s: PowerSpectrum(cluster, mask=local_mask)(img, s) ) ) ps_fluc = power_spectrum.apply(None, fluc, scales) axs[0].plot(scales, ps_fluc, color=cmr.amber_r(high_rad)) axs[0].loglog() axs[0].set_xlabel('Scales [$R_{500}$]') axs[0].set_ylabel('Power spectrum') axs[1].imshow(mask_img, cmap=cmr.ember_r) axs[1].set_axis_off() plt.show();