Skip to content

Commit

Permalink
Merge branch 'master' into xcor_cmp
Browse files Browse the repository at this point in the history
  • Loading branch information
@vitenti
vitenti committed Apr 13, 2024
2 parents 7a8eb33 + 1efaffe commit efc29fa
Show file tree
Hide file tree
Showing 43 changed files with 4,691 additions and 369 deletions.
2 changes: 0 additions & 2 deletions .flake8
View file
Original file line number Diff line number Diff line change
@@ -1,5 +1,3 @@
[flake8]
max-line-length = 88
select = C,E,F,W,B,B950
extend-ignore = E203, E501

1 change: 1 addition & 0 deletions devel_environment.yml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -17,6 +17,7 @@ dependencies:
- gobject-introspection <=1.78
- gsl
- hdf5
- healpy
- lcov
- libflint
- libfyaml
Expand Down
1 change: 1 addition & 0 deletions docs/numcosmo-docs.sgml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -434,6 +434,7 @@
<xi:include href="xml/nc_data_dist_mu.xml"/>
<xi:include href="xml/nc_data_snia_cov.xml"/>
<xi:include href="xml/nc_data_cluster_ncount.xml"/>
<xi:include href="xml/nc_data_cluster_ncounts_gauss.xml"/>
<xi:include href="xml/nc_data_cluster_pseudo_counts.xml"/>
<xi:include href="xml/nc_data_cluster_wl.xml"/>
<xi:include href="xml/nc_data_reduced_shear_cluster_mass.xml"/>
Expand Down
106 changes: 47 additions & 59 deletions examples/example_sphere_proj.py
100755 → 100644
View file
Original file line number Diff line number Diff line change
Expand Up @@ -24,26 +24,25 @@

"""Simple example computing spherical projections."""

import math
import numpy as np
import matplotlib
import matplotlib.pyplot as plt

from matplotlib.colors import SymLogNorm

from numcosmo_py import Nc, Ncm


matplotlib.rcParams.update({"font.size": 11})

#
# Initializing the library objects, this must be called before
# any other library function.
#
Ncm.cfg_init()
def compute_spherical_projection():
"""Compute spherical projections of the power spectrum."""

matplotlib.rcParams.update({"font.size": 11})

def test_spherical_projection() -> None:
"""Example computing spherical projection."""
#
# Initializing the library objects, this must be called before
# any other library function.
#
Ncm.cfg_init()

#
# New homogeneous and isotropic cosmological model NcHICosmoDEXcdm
Expand All @@ -52,17 +51,17 @@ def test_spherical_projection() -> None:
cosmo.omega_x2omega_k()
cosmo.param_set_by_name("Omegak", 0.0)
cosmo.param_set_by_name("w", -1.0)
cosmo.param_set_by_name("Omegab", 0.04909244421)
cosmo.param_set_by_name("Omegac", 0.26580755578)
cosmo.param_set_by_name("Omegab", 0.045)
cosmo.param_set_by_name("Omegac", 0.255)
cosmo.param_set_by_name("massnu_0", 0.06)
cosmo.param_set_by_name("ENnu", 2.0328)

reion = Nc.HIReionCamb.new()
prim = Nc.HIPrimPowerLaw.new()

cosmo.param_set_by_name("H0", 67.31)
cosmo.param_set_by_name("H0", 67.0)

prim.param_set_by_name("n_SA", 0.9658)
prim.param_set_by_name("n_SA", 0.96)
prim.param_set_by_name("ln10e10ASA", 3.0904)

reion.param_set_by_name("z_re", 9.9999)
Expand All @@ -73,7 +72,7 @@ def test_spherical_projection() -> None:
#
# Printing the parameters used.
#
print("# Model parameters: ", end="")
print("# Model parameters: ", end=" ")
cosmo.params_log_all()
print(f"# Omega_X0: {cosmo.E2Omega_de(0.0): 22.15g}")

Expand All @@ -82,11 +81,10 @@ def test_spherical_projection() -> None:
# Redshift bounds
z_min = 0.0
z_max = 2.0
# zdiv = 0.49999999999

# Mode bounds
k_min = 1.0e-6
k_max = 1.0e5
k_min = 1.0e-8
k_max = 1.0e6

ps_lin.set_kmin(k_min)
ps_lin.set_kmax(k_max)
Expand All @@ -95,22 +93,28 @@ def test_spherical_projection() -> None:

ps_lin.set_intern_k_min(k_min)
ps_lin.set_intern_k_max(50.0)
psf = Ncm.PowspecFilter.new(ps_lin, Ncm.PowspecFilterType.TOPHAT)
psf.set_best_lnr0()

ps = Nc.PowspecMNLHaloFit.new(ps_lin, 2.0, 1.0e-5)
ps = Nc.PowspecMNLHaloFit.new(ps_lin, 2.0, 1.0e-8)

ps.set_kmin(k_min)
ps.set_kmax(k_max)
ps.require_zi(z_min)
ps.require_zf(z_max)

ell_min = 250
ell_max = 250
old_amplitude = math.exp(prim.props.ln10e10ASA)
prim.props.ln10e10ASA = math.log((0.81 / cosmo.sigma8(psf)) ** 2 * old_amplitude)

ell_min = 0
ell_max = 1

sproj = Ncm.PowspecSphereProj.new(ps, ell_min, ell_max)
sproj.props.reltol = 1.0e-4

xi_i = 1.0e1
xi_f = 1.0e4
sproj.props.reltol = 1.0e-8
sproj.props.reltol_z = 1.0e-8
xi_i = 1.0e-1
xi_f = 1.0e7

sproj.set_xi_i(xi_i)
sproj.set_xi_f(xi_f)
Expand All @@ -135,36 +139,21 @@ def test_spherical_projection() -> None:

xi_a = np.geomspace(xi_i, xi_f, num=5000)
z_a = np.array([dist.inv_comoving(cosmo, xi / RH_Mpc) for xi in xi_a])
# index = np.logical_and((z_a > 0.01), (z_a < 10.0))
index = np.logical_and((z_a > 0.0), (z_a < 10.0))

z_a = np.linspace(0.05, 0.5, num=2000)
z_a = np.linspace(0.2, 2, num=2001)[1:]

# print (xi_a)
# print (z_a)

for ell in range(ell_min, ell_min + 1):
# limber = np.array(
# [
# ps.eval(cosmo, z, (0.5 + ell) / xi) / (xi * xi * RH_Mpc / cosmo.E(z))
# for (xi, z) in zip(xi_a, z_a)
# ]
# )

# plt.plot (xi_a[index], limber[index], label = r'$\ell$ = %d, limber' % ell)

# for w_i in range (20):
# w = sproj.get_w (w_i)
# Cell_s = sproj.peek_ell_spline (w_i, ell)
# Cell = [Cell_s.eval (math.log (xi)) for xi in xi_a]

# for xi, Cell_i in zip (xi_a, Cell):
# sCell_i = sproj.eval_Cell_xi1_xi2 (ell, xi / w, xi * w)
# Cell_int_i = ps.sproj (cosmo, 1.0e-5, ell, 0.0, 0.0, xi / w, xi * w)

# z_a = np.array (z_a)
# Cell = np.array ([(gf.eval (cosmo, z)**2 * Cell_i) for (Cell_i, z) in zip (Cell, z_a)])
limber = np.array(
[
ps.eval(cosmo, z, (0.5 + ell) / xi) / (xi * xi * RH_Mpc / cosmo.E(z))
for (xi, z) in zip(xi_a, z_a)
]
)

# plt.plot (xi_a, Cell, label = r'$\ell$ = %d, $w$ = %f, full' % (ell, w))
m = [
[
sproj.eval_Cell_xi1_xi2(
Expand All @@ -180,25 +169,23 @@ def test_spherical_projection() -> None:
for z2 in z_a
]
cov = np.asanyarray(m)
# print (np.diag (cov))
std_ = np.sqrt(np.diag(cov))
corr = cov / np.outer(std_, std_)
# print (m)
# print (corr)
# plt.matshow (np.abs (corr), norm = LogNorm (vmin=1.0e-8, vmax = 1.0))
# plt.matshow (corr, extent = np.array ([z_a[0], z_a[-1], z_a[0], z_a[-1]]),
# origin = "lower")
print(len(cov))
print(len(cov[0]))
plt.matshow(
corr,
cov,
extent=np.array([z_a[0], z_a[-1], z_a[0], z_a[-1]]),
origin="lower",
norm=SymLogNorm(1.0e-8),
norm=SymLogNorm(1.0e-6),
)
plt.xlabel(r"$z_2$")
plt.ylabel(r"$z_1$")
# plt.matshow (corr)
# plt.matshow (corr, norm = SymLogNorm (1.0e-4))
plt.title(rf"$C_{{{ell}}}(z_1, z_2)$")
plt.title(r"$C_{%d}(z_1, z_2)$" % (ell))
plt.colorbar()

plt.show()
# plt.xticks(xi_a)

# plt.xlabel (r'$\xi \; [\mathrm{Mpc}]$')
Expand All @@ -208,8 +195,9 @@ def test_spherical_projection() -> None:
# plt.yscale ('symlog', linthreshy=1.0e-8)
# plt.xlim ([200.0, 4000.0])

plt.show()
# plt.show ()



if __name__ == "__main__":
test_spherical_projection()
compute_spherical_projection()
46 changes: 34 additions & 12 deletions examples/twofluids_wkb_spec.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -44,7 +44,7 @@

def test_two_fluids_wkb_spec() -> None:
"""Compute WKB approximation for the two-fluids model spectrum."""

#
# New homogeneous and isotropic cosmological model NcHICosmoQGRW
#
Expand Down Expand Up @@ -80,8 +80,8 @@ def test_two_fluids_wkb_spec() -> None:
Ps_S2 = []
Ps_Pzeta1 = []
Ps_PS1 = []
# Ps_Pzeta2 = []
# Ps_PS2 = []
Ps_Pzeta2 = []
Ps_PS2 = []

out_file = open("twofluids_spectrum_{w}.dat", "w", encoding="utf-8")

Expand All @@ -93,18 +93,16 @@ def test_two_fluids_wkb_spec() -> None:
pert.set_mode_k(k)
k_a.append(k)

alphaf = cosmo.abs_alpha(1.0e20)
alphaf = cosmo.abs_alpha(1.0e10)

# print ("# Evolving mode %e from %f to %f" % (k, alphai, alphaf))

alphai = -cosmo.abs_alpha(start_alpha1 * k**2)
pert.get_init_cond_zetaS(cosmo, alphai, 1, 0.25 * math.pi, ci)
pert.set_init_cond(cosmo, alphai, 1, False, ci)

print(f"# Mode 1 k {k: 21.15e}, state module {pert.get_state_mod():f}")

pert.evolve(cosmo, alphaf)
pert.evolve(cosmo, alphaf)
v, _alphac = pert.peek_state(cosmo)
print ("# Evolving mode %e from %f to %f" % (k, alphai, alphaf))

Delta_zeta1 = (
k**3
Expand Down Expand Up @@ -148,9 +146,9 @@ def test_two_fluids_wkb_spec() -> None:

alphai = -cosmo.abs_alpha(start_alpha2 * k**2)
pert.get_init_cond_zetaS(cosmo, alphai, 2, 0.25 * math.pi, ci)
pert.set_init_cond(cosmo, alphai, 0, False, ci)
pert.set_init_cond(cosmo, alphai, 2, False, ci)

print("# Mode 2 k {k: 21.15e}, state module {pert.get_state_mod():f}")
print(f"# Mode 2 k {k: 21.15e}, state module {pert.get_state_mod():f}")

pert.evolve(cosmo, alphaf)
v, _alphac = pert.peek_state(cosmo)
Expand Down Expand Up @@ -192,8 +190,8 @@ def test_two_fluids_wkb_spec() -> None:

Ps_zeta2.append(Delta_zeta2)
Ps_S2.append(Delta_S2)
Ps_zeta2.append(Delta_Pzeta2)
Ps_S2.append(Delta_PS2)
Ps_Pzeta2.append(Delta_Pzeta2)
Ps_PS2.append(Delta_PS2)

out_file.write(
f"{k: 20.15e} {Delta_zeta1: 20.15e} {Delta_zeta2: 20.15e} {Delta_S1: 20.15e} "
Expand All @@ -216,3 +214,27 @@ def test_two_fluids_wkb_spec() -> None:

plt.show()
plt.clf()

test_two_fluids_wkb_spec()


'''
import pandas as pd
data = pd.read_csv('twofluids_spectrum_{w}.dat', header=None, delim_whitespace= True)
lnk = data[0]
Ps1zeta = data[1]
Ps1S = data[2]
Ps2zeta = data[3]
Ps2S = data[4]
ns1zeta = np.polyfit(lnk,Ps1zeta,1)
ns1S = np.polyfit(lnk,Ps1S,1)
ns2zeta = np.polyfit(lnk,Ps2zeta,1)
ns2S = np.polyfit(lnk,Ps2S,1)
print(ns1zeta)
print(ns1S)
print(ns2zeta)
print(ns2S)
'''

0 comments on commit efc29fa

Please sign in to comment.