Shift-Powermap/Shift-Cls : differences in final angular power spectrum

[dlanzieri]

This issue aims to discuss the possible reasons why the theoretical angular power spectrum does not match the one computed from the log-normal simulated maps when we sample from correlated maps.
After conducting several tests, I start to think that this is not related to how we are sampling from correlated maps, but rather to how we are generating the log-normal power spectrum in this new implementation.

To illustrate this, I have computed the power spectrum for a single redshift bin (the 5th in the following plot) using two different procedures.

In the first case, I generated the power map from the theoretical angular power spectrum first and then shifted the map (map 1). In the second case (map 2), I first shifted the angular power spectrum and then generated the power map.

Below are a few lines of code that can help to understand the differences.

def make_power_map(pk_fn, N, map_size, zero_freq_val=0.0, model_type=None ):
    k = 2 * jnp.pi * jnp.fft.fftfreq(N, d=map_size / N)
    kcoords = jnp.meshgrid(k, k)
    k = jnp.sqrt(kcoords[0]**2 + kcoords[1]**2)
    ps_map = pk_fn(k)
    ps_map = ps_map.at[0, 0].set(zero_freq_val)
    if model_type=='no_correlation':
        power_map = ps_map * (N / map_size)**2
    else:
        power_map = ps_map * (N / map_size)
    return power_map

def make_lognormal_power_map(power_map, shift, zero_freq_val=0.0):
    power_spectrum_for_lognorm = jnp.fft.ifft2(power_map).real
    power_spectrum_for_lognorm = jnp.log(1 +
                                         power_spectrum_for_lognorm / shift**2)
    power_spectrum_for_lognorm = jnp.abs(
        jnp.fft.fft2(power_spectrum_for_lognorm))
    power_spectrum_for_lognorm = power_spectrum_for_lognorm.at[0, 0].set(0.)
    return power_spectrum_for_lognorm


def make_lognormal_cl(cl,shift):
    xsi = jnp.fft.ifft(cl)
    xsi_shift = jnp.log(1 + xsi/(shift)**2)
    cl_shifted = jnp.fft.fft(xsi_shift).real
    return  cl_shifted

N=128            
map_size=5   
pix_area = (map_size * 60 / N)**2 # arcmin2 
map_size = map_size / 180 * jnp.pi    # radians
omega_c = 0.3
sigma_8 =0.8
cosmo = jc.Planck15(Omega_c=omega_c, sigma8=sigma_8)
ell_tab = 2 * jnp.pi * abs(jnp.fft.fftfreq(N, d=map_size / N))
tracer = jc.probes.WeakLensing([nz_shear[-1]])
shift = shift_fn(cosmo.Omega_m, sigma_8)
cell_tab = jc.angular_cl.angular_cl(cosmo, ell_tab, [tracer])[0]

P_1 = lambda k: jc.scipy.interpolate.interp(k.flatten(), ell_tab, cell_tab).reshape(k.shape)
power_map_1 = make_power_map(P_1, N, map_size, model_type='no_correlation') 
power_map_1 = make_lognormal_power_map(power_map_1, shift)


P_2 = lambda k: jc.scipy.interpolate.interp(k.flatten(), ell_tab, cell_tab_shifted).reshape(k.shape)
power_map_2 = make_power_map(P_2, N, map_size, model_type='no_correlation') 


field_1 = jnp.fft.ifft2(jnp.fft.fft2(z) * jnp.sqrt(power_map_1)).real
field_1 = shift * (jnp.exp(field_1 - jnp.var(field_1) / 2) - 1)
field_2 = jnp.fft.ifft2(jnp.fft.fft2(z) * jnp.sqrt(power_map_2)).real
field_2 = shift * (jnp.exp(field_2 - jnp.var(field_2) / 2) - 1)

Below are the results computed from the two maps, averaged over 20 realizations:

If you remember @EiffL, when we compared the two power maps from the two methods, we noted some differences. I am starting to think that these differences are the ones causing the bias in the results.

[EiffL]

Thanks Denise, but I'm not sure I agree. What I remember is that we obtained identical power maps, as it should be, in the case of bin 1, no?

[EiffL]

this is solved now