Notebook

Large Scale Structure Statistics¶

Material used to prepare my 2 Feb 2015 LSST DESC Dark Energy School lecture. You are welcome to modify and improve this material. Contact me at dkirkby@uci.edu with questions or corrections. This material is released under the MIT License.

Some of the data files used for the plots are created with external programs, using command lines provided below. There is also an accompanying Mathematica notebook used to generate of the files we read below.

Initialization¶

In [1]:

%pylab inline

Populating the interactive namespace from numpy and matplotlib

In [2]:

import math

In [3]:

import astropy.units as u
import astropy.constants as const

In [4]:

import scipy.interpolate
import scipy.fftpack

In [5]:

# yt is not included in anaconda but installs easily with "pip install yt".
# It is only used to make one plot below, so you skip this if you don't need that plot.
import yt

Fiducial Cosmology¶

In [5]:

from astropy.cosmology import Planck13 as fiducial

In [6]:

h0 = fiducial.H0/(100*u.km/u.s/u.Mpc)
print 'Fiducial h0 =',h0

Fiducial h0 = 0.6777

The median redshift of the LSST "gold" weak-lensing sample (i < 25.3, SNR > 20 for point sources). Redshift errors in this sample should be $\sigma(z)/(1+z)$ = (1-2)% with 90% of galaxies better than 4% (from Section 4.1 of http://arxiv.org/abs/0805.2366).

Load the Planck13 power spectrum at z=1.2 which has sigma8 = 0.465 at z=1.2 and 0.835 at z=0. We are using the accompanying CAMB input file camb/Planck13.camb and running the July 2013 version of CAMB since that was the version used for the Planck13 likelihood calculations. The CAMB input file was generated using the CosmoTools mathematica package (using the accompanying camb/Planck13.nb Mathematica notebook) and would need some minor modifications to run with more recent versions of CAMB due to changes in how the neutrino sector is parameterized.

In [7]:

kvec,planck13_Pk = np.loadtxt('camb/Planck13_matterpower_1.dat').transpose()

Create an interpolating function that is linear in log(k),log(P(k)) and returns zero for k=0, where units are (h/Mpc):

In [8]:

pk_interpolator = scipy.interpolate.InterpolatedUnivariateSpline(np.log(kvec),np.log(planck13_Pk),k=1)

In [9]:

def pk_func(kvec_hMpc):
    pk = np.zeros_like(kvec_hMpc)
    nonzero = (kvec_hMpc>0)
    pk[nonzero] = np.exp(pk_interpolator(np.log(kvec_hMpc[nonzero])))
    # Units of P(k) are (h/Mpc)**2
    return pk

Tabulate the k-space window function for evaluating sigma8:

In [10]:

kR = 8.*kvec
sigma8_window = (3/kR**3*(np.sin(kR)-kR*np.cos(kR)))**2*kvec**3/(2*np.pi**2)

Plot the sigma8 window function superimposed on $P(k)$ (not $k^3 P(k)/(2\pi^2)$):

In [11]:

fig = plt.figure(figsize=(6,4))
plt.fill_between(kvec,planck13_Pk,edgecolor='blue',lw=1,facecolor=(0.,0.,1.,0.25))
plt.xlim(kvec[0],kvec[-1])
plt.xlabel('Wavenumber k (h/Mpc)')
plt.ylabel('Power P(k) (h/Mpc)^3',color='blue')
plt.xscale('log')
ax2 = plt.gca().twinx()
ax2.fill_between(kvec,sigma8_window,edgecolor='red',lw=0.5,facecolor=(1.,0.,0.,0.25))
ax2.set_xlim(kvec[0],kvec[-1])
ax2.set_yticklabels([])
ax2.set_yticks([]);
ax2.set_ylabel('Sigma8 Window Function',color='red');
plt.savefig('plots/sigma8window.pdf')

Alternative Cosmologies¶

Define some alternative power spectra tabulated on the same k-vector grid:

In [13]:

k0 = 0.05 # h/Mpc
r0 = 100 # Mpc/h
powerlaw_Pk = 1e6*kvec/(1+(kvec/k0)**3.5)
wiggles_Pk = powerlaw_Pk*(1 + np.sin(kvec*r0)*np.exp(-(kvec*r0)**2/1e4))
bump1_Pk = 4e4*np.exp(-np.log(kvec/k0)**2)
bump2_Pk = 4e4*np.exp(-np.log(kvec/(2*k0))**2)
noscale_Pk = kvec**-2

In [14]:

plt.plot(kvec/k0,planck13_Pk,'k-')
plt.plot(kvec/k0,powerlaw_Pk,'r-')
plt.plot(kvec/k0,wiggles_Pk,'b-')
plt.plot(kvec/k0,bump1_Pk,'g-')
plt.plot(kvec/k0,bump2_Pk,'g--')
plt.plot(kvec/k0,noscale_Pk,'magenta')
plt.xlim(kvec[0]/k0,kvec[-1]/k0)
plt.ylim(1e-3,1e5)
plt.xscale('log')
plt.yscale('log')
plt.xlabel('Wavenumber k/k0')
plt.ylabel('Power P(k) (h/Mpc)^3')
plt.grid();

Save these tabulated power spectra:

In [15]:

np.savetxt('alt/powerlaw_Pk.dat',np.vstack((kvec,powerlaw_Pk)).T)
np.savetxt('alt/wiggles_Pk.dat',np.vstack((kvec,wiggles_Pk)).T)
np.savetxt('alt/bump1_Pk.dat',np.vstack((kvec,bump1_Pk)).T)
np.savetxt('alt/bump2_Pk.dat',np.vstack((kvec,bump2_Pk)).T)
np.savetxt('alt/noscale_Pk.dat',np.vstack((kvec,noscale_Pk)).T)

Calculate the corresponding r-space correlation functions (times 2*pi^2) using the following command lines:

cosmotrans -i camb/Planck13_matterpower_1.dat -o camb/Planck13_xi.dat --min 1 --max 200
cosmotrans -i alt/powerlaw_Pk.dat -o alt/powerlaw_xi.dat --min 1 --max 200
cosmotrans -i alt/wiggles_Pk.dat -o alt/wiggles_xi.dat --min 1 --max 200
cosmotrans -i alt/bump1_Pk.dat -o alt/bump1_xi.dat --min 1 --max 200
cosmotrans -i alt/bump2_Pk.dat -o alt/bump2_xi.dat --min 1 --max 200
cosmotrans -i alt/noscale_Pk.dat -o alt/noscale_xi.dat --min 1 --max 200

The cosmotrans program is part of the cosmo package used by baofit.

Load and plot these correlation functions:

In [16]:

rvec,planck13_xi = np.loadtxt('camb/Planck13_xi.dat').transpose()
rvec,powerlaw_xi = np.loadtxt('alt/powerlaw_xi.dat').transpose()
rvec,wiggles_xi = np.loadtxt('alt/wiggles_xi.dat').transpose()
rvec,bump1_xi = np.loadtxt('alt/bump1_xi.dat').transpose()
rvec,bump2_xi = np.loadtxt('alt/bump2_xi.dat').transpose()
rvec,noscale_xi = np.loadtxt('alt/noscale_xi.dat').transpose()
planck13_xi /= 2*np.pi**2
powerlaw_xi /= 2*np.pi**2
wiggles_xi /= 2*np.pi**2
bump1_xi /= 2*np.pi**2
bump2_xi /= 2*np.pi**2
noscale_xi /= 2*np.pi**2

In [17]:

plt.plot(rvec*k0,rvec**2*planck13_xi,'k-')
plt.plot(rvec*k0,rvec**2*powerlaw_xi,'r-')
plt.plot(rvec*k0,rvec**2*wiggles_xi,'b-')
plt.plot(rvec*k0,rvec**2*bump1_xi,'g-')
plt.plot(rvec*k0,rvec**2*bump2_xi,'g--')
plt.plot(rvec*k0,rvec**2*noscale_xi,'magenta')
plt.xlim(0.,rvec[-1]*k0)
plt.ylim(-150.,500.)
plt.xlabel('Separation r*k0')
plt.ylabel('Correlation function r^2*xi(r) (Mpc/h)^2');

Run the following commands to generate (100 Mpc/h)^3 Gaussian random fields from these power spectra and save 2D slices:

cosmogrf --verbose --spacing 0.5 --nx 200 --load-power camb/Planck13_matterpower_1.dat --save-delta-slice camb/planck13_slice.dat
cosmogrf --verbose --spacing 0.5 --nx 200 --load-power alt/powerlaw_Pk.dat --save-delta-slice alt/powerlaw_slice.dat
cosmogrf --verbose --spacing 0.5 --nx 200 --load-power alt/wiggles_Pk.dat --save-delta-slice alt/wiggles_slice.dat
cosmogrf --verbose --spacing 0.5 --nx 200 --load-power alt/bump1_Pk.dat --save-delta-slice alt/bump1_slice.dat
cosmogrf --verbose --spacing 0.5 --nx 200 --load-power alt/bump2_Pk.dat --save-delta-slice alt/bump2_slice.dat
cosmogrf --verbose --spacing 0.5 --nx 200 --load-power alt/noscale_Pk.dat --save-delta-slice alt/noscale_slice.dat

The cosmogrf program is also in the cosmo package.

Load these realization slices:

In [18]:

planck13_slice = np.loadtxt('camb/planck13_slice.dat')
powerlaw_slice = np.loadtxt('alt/powerlaw_slice.dat')
wiggles_slice = np.loadtxt('alt/wiggles_slice.dat')
bump1_slice = np.loadtxt('alt/bump1_slice.dat')
bump2_slice = np.loadtxt('alt/bump2_slice.dat')
noscale_slice = np.loadtxt('alt/noscale_slice.dat')

Create images where students have to match P(k), xi(r) and a realization.

In [19]:

def plot_Pk(data,label):
    plt.plot(kvec/k0,data,'k-')
    visible = (1e-2 < kvec/k0) & (kvec/k0 < 1e2)
    pmax = np.max(data[visible])
    plt.xlim(1e-2,1e2)
    plt.ylim(2e-4*pmax,2*pmax)
    plt.xscale('log')
    plt.yscale('log')
    plt.xlabel('Wavenumber $k/k_{0}$')
    plt.ylabel('Power $P(k)$')
    plt.gca().get_yaxis().set_ticklabels([])
    plt.grid()
    plt.annotate(label,xy=(0.85,0.85),xycoords='axes fraction',fontsize='xx-large')
    
def plot_xi(data,label):
    plt.plot(k0*rvec,rvec**2*data,'k-')
    plt.xlim(0.,rvec[-1]*k0)
    if np.min(data) >= 0:
        plt.ylim(-0.05*np.max(rvec**2*data),None)
    plt.xlabel('Separation $k_{0} r$')
    plt.ylabel('Correlation $r^{2} \\xi(r)$')
    plt.gca().get_yaxis().set_ticklabels([])
    plt.gca().axhline(0., linestyle='--', color='k')
    plt.annotate(label,xy=(0.85,0.85),xycoords='axes fraction',fontsize='xx-large')
    
# Display a 2D density field using a color map where 0 = white, >0 is red, <0 is blue.
# The input data is assumed to have a mean near zero.
def slice_plot(data,clip_percent=1.,**kwargs):
    lo,hi = np.percentile(data,(clip_percent,100.-clip_percent))
    assert lo < 0. and hi > 0.
    cut = 0.5*(hi-lo)
    plt.imshow(data,cmap='RdBu_r',vmin=lo,vmax=hi,**kwargs)
    
def plot_slice(data,label):
    slice_plot(data,extent=(0.,5.,0.,5.))
    plt.xlabel('$k_{0}x$')
    plt.ylabel('$k_{0}y$')
    plt.annotate(label,xy=(0.85,0.85),xycoords='axes fraction',fontsize='xx-large')

In [20]:

def match1(save_name=None):
    fig = plt.figure(figsize=(7.5,7.5))
    # planck13 P(k)
    plt.subplot(3,3,1)
    plot_Pk(planck13_Pk,'A')
    # bump2 P(k)
    plt.subplot(3,3,4)
    plot_Pk(bump2_Pk,'B')
    # bump1 P(k)
    plt.subplot(3,3,7)
    plot_Pk(bump1_Pk,'C')
    # bump2 xi(r)
    plt.subplot(3,3,2)
    plot_xi(bump2_xi,'A')
    # bump1 xi(r)
    plt.subplot(3,3,5)
    plot_xi(bump1_xi,'B')
    # planck13 xi(r)
    plt.subplot(3,3,8)
    plot_xi(planck13_xi,'C')
    # bump1 realization
    plt.subplot(3,3,3)
    plot_slice(bump1_slice,'A')
    # planck13 realization
    plt.subplot(3,3,6)
    plot_slice(planck13_slice,'B')
    # bump2 realization
    plt.subplot(3,3,9)
    plot_slice(bump2_slice,'C')
    #
    plt.tight_layout()
    if save_name:
        plt.savefig(save_name)
match1('plots/match1.pdf')

In [21]:

def match2(save_name=None):
    fig = plt.figure(figsize=(7.5,7.5))
    # noscale P(k)
    plt.subplot(3,3,1)
    plot_Pk(noscale_Pk,'A')
    # powerlaw P(k)
    plt.subplot(3,3,4)
    plot_Pk(powerlaw_Pk,'B')
    # wiggles P(k)
    plt.subplot(3,3,7)
    plot_Pk(wiggles_Pk,'C')
    # powerlaw xi(r)
    plt.subplot(3,3,2)
    plot_xi(powerlaw_xi,'A')
    # noscale xi(r)
    plt.subplot(3,3,5)
    plot_xi(noscale_xi,'B')
    # wiggles xi(r)
    plt.subplot(3,3,8)
    plot_xi(wiggles_xi,'C')
    # noscale realization
    plt.subplot(3,3,3)
    plot_slice(noscale_slice,'A')
    # powerlaw realization
    plt.subplot(3,3,6)
    plot_slice(powerlaw_slice,'B')
    # wiggles realization
    plt.subplot(3,3,9)
    plot_slice(wiggles_slice,'C')
    #
    plt.tight_layout()
    if save_name:
        plt.savefig(save_name)
match2('plots/match2.pdf')