Input/Output Power spectra not matching

[florpi]

Thanks for the repo :) It's very clear and easy to use.

I am having issues generating a density field with nbodykit's power spectrum and powerbox, maybe you can help me. This is what I'm doing,

import matplotlib.pyplot as plt
import numpy as np
from astropy.cosmology import Planck18
from nbodykit.cosmology.power.halofit import HalofitPower
import powerbox as pbox
from powerbox import get_power

redshift = 0.
pk = HalofitPower(Planck18,redshift=redshift)
boxsize = 512.

pb = pbox.PowerBox(
    N=256,                     
    dim=2,                     
    pk = pk, 
    boxlength = boxsize,           
    ensure_physical=True       
)

delta_x = pb.delta_x()

samples = pb.create_discrete_sample(
    delta_x=delta_x,
    nbar=0.1,    
    min_at_zero=True  
)
p_k_field, bins_field = get_power(delta_x, pb.boxlength)
p_k_samples, bins_samples = get_power(samples, pb.boxlength,N=pb.N)
plt.plot(bins_field,pk(bins_field), label='Input')
plt.plot(bins_field, p_k_field,label="Output Field Power")
plt.plot(bins_samples, p_k_samples,label="Output Sample Power")
plt.legend()
plt.xscale('log')
plt.yscale('log')

And this what I find for the output power spectra:

As you can see, the amplitude of each of one them is different. Any idea of why they don't match up?

Also, I have noticed that when running the code in 3D it becomes quite slow. Do you have any suggestions to speed it up?

Thank you in advance!

[steven-murray]

Thanks @florpi for the detailed issue. This is confusing indeed. Can you make sure that using a power spectrum of say lambda k: k**-2 gives the correct normalization for you?

[florpi]

That one is closer

[steven-murray]

I wonder if this has something to do with the slope of the power spectrum? I don't recall this being an issue before though.

In terms of speed: you should get some speed increase by using pyfftw if you have it. Other than that, I don't think you can really get any faster. 3D cubes are just large :-)

[florpi]

It seems to be coming from the ensure_physical flag

But then it doesn't work if you sample from the non physical field, I get this error message,

ValueError: lam < 0 or lam contains NaNs

[steven-murray]

Ahhhh I see. You are definitely messing up the power spectrum by clipping the density values to be "physical". This will work "ok" if your density field only just dips below -1, but will get worse if it goes down to -50 or more.

You might want to try the LogNormal density field?

[florpi]

Yes, I tried that too but it is still problematic

This is the code I ran,

import matplotlib.pyplot as plt
import numpy as np
from astropy.cosmology import Planck18
from nbodykit.cosmology.power.halofit import HalofitPower
import powerbox as pbox
from powerbox import get_power

redshift = 0.
pk = HalofitPower(Planck18,redshift=redshift)
boxsize = 512.
pb = pbox.LogNormalPowerBox(
    N=256,                     
    dim=2,                     
    pk = pk, 
    boxlength = boxsize,           
)

delta_x = pb.delta_x()
p_k_field, bins_field = get_power(delta_x, pb.boxlength)

samples = pb.create_discrete_sample(
    delta_x=delta_x,
    nbar=1., 
    min_at_zero=True,
)
p_k_samples, bins_samples = get_power(samples, pb.boxlength,N=pb.N)
plt.plot(bins_field,pk(bins_field), label='Input')
plt.plot(bins_field, p_k_field,label="Log-normal field")
plt.plot(bins_samples, p_k_samples,label="Log-normal Samples")
plt.legend()
plt.xscale('log')
plt.yscale('log')
plt.show()

The measured power spectrum seems to vary a lot with each random sample though (orders of magnitude?) This is an example of the same code but different sample,

[steven-murray]

Interesting. If it's fast enough, can you run ~100 samples and plot the mean and std?

[florpi]

This is for 500 samples

Std is large

[steven-murray]

Indeed it is. I'm not sure if I expect this or not. It seems like if you treated the bins as independent, the overall significance of the deviation would be more than 1sigma.

I can't think of what else to consider right now... maybe making a similar plot for k^-2 (for lognormal)

[florpi]

It seems to be a matter of amplitude. When P(k) = 1/k, the variance is low,

But when P(k) = 5.e2/k (to match the amplitude of halofit's power spectrum),

Conversely, this is what happens when I use P(k) = P_halofit(k)/5.e2

[steven-murray]

Yes, this makes sense -- the variance should be proportional to the power spectrum, but that should look roughly constant on a log plot, I'd have thought. I'm wondering if there is an issue with the units of the box length or something -- physically the cosmological power spectrum shouldn't give crazy fields.

Maybe check the little h's?