Skip to content

Commit

Permalink
fix style
Browse files Browse the repository at this point in the history
  • Loading branch information
santiagocasas committed Dec 1, 2024
1 parent c43e2c9 commit 4b73bd2
Show file tree
Hide file tree
Showing 14 changed files with 153 additions and 138 deletions.
2 changes: 1 addition & 1 deletion cosmicfishpie/CMBsurvey/CMB_cov.py
Original file line number Diff line number Diff line change
Expand Up @@ -97,7 +97,7 @@ def getclsnoise(self, cls):
]

Bell = [
np.exp(ang**2.0 * noisy_cls["ells"] * (noisy_cls["ells"] + 1) / 2.0) for ang in thetab
np.exp(ang ** 2.0 * noisy_cls["ells"] * (noisy_cls["ells"] + 1) / 2.0) for ang in thetab
]

wtemp = [
Expand Down
2 changes: 1 addition & 1 deletion cosmicfishpie/LSSsurvey/photo_cov.py
Original file line number Diff line number Diff line change
Expand Up @@ -142,7 +142,7 @@ def getclsnoise(self, cls):
)
elif obs == "WL":
noisy_cls[obs + " " + str(ind) + "x" + obs + " " + str(ind)] += (
self.ellipt_error**2.0
self.ellipt_error ** 2.0
) / self.ngalbin[ind - 1]

return noisy_cls
Expand Down
2 changes: 1 addition & 1 deletion cosmicfishpie/LSSsurvey/photo_window.py
Original file line number Diff line number Diff line change
Expand Up @@ -78,7 +78,7 @@ def dNdz(self, z):
pref = z / self.z0_p
expo = z / self.z0

return pref**2 * np.exp(-(expo**self.ngamma))
return pref ** 2 * np.exp(-(expo ** self.ngamma))

def n_i(self, z, i):
"""Function to compute the unnormalized dN/dz(z) with a window picking function applied to it
Expand Down
60 changes: 33 additions & 27 deletions cosmicfishpie/LSSsurvey/spectro_cov.py
Original file line number Diff line number Diff line change
Expand Up @@ -18,8 +18,7 @@

class SpectroCov:
def __init__(
self, fiducialpars, fiducial_specobs=None, bias_samples=["g", "g"],
configuration=None
self, fiducialpars, fiducial_specobs=None, bias_samples=["g", "g"], configuration=None
):
"""
Initializes an object with specified fiducial parameters and computes
Expand Down Expand Up @@ -75,11 +74,11 @@ def __init__(
self.fsky_spectro = self.area_survey / upm.areasky()
if fiducial_specobs is None:
self.pk_obs = spec_obs.ComputeGalSpectro(
fiducialpars,
fiducial_cosmopars=fiducialpars,
bias_samples=bias_samples,
configuration=self.config
)
fiducialpars,
fiducial_cosmopars=fiducialpars,
bias_samples=bias_samples,
configuration=self.config,
)
# if no other parameters are provided, the method will use the fiducials from config
else:
self.pk_obs = fiducial_specobs
Expand Down Expand Up @@ -216,7 +215,7 @@ def veff(self, zi, k, mu):
The effective volume for a given wavenumber, angle and redshift
"""
npobs = self.n_density(zi) * self.pk_obs.observed_Pgg(zi, k, mu)
prefactor = 1 / (8 * (np.pi**2))
prefactor = 1 / (8 * (np.pi ** 2))
covterm = prefactor * (npobs / (1 + npobs)) ** 2
if zi < self.inter_z_bin_mids[0] or zi > self.inter_z_bin_mids[-1]:
covterm = np.zeros_like(covterm)
Expand Down Expand Up @@ -272,14 +271,16 @@ def P_noise_21(self, z, k, mu, temp_dim=True, beam_term=False):
elif temp_dim:
temp = 1
pref = (2 * np.pi * self.pk_obs.fsky_IM) / (self.pk_obs.f_21 * self.pk_obs.t_tot)
cosmo = ((1 + z) ** 2 * self.pk_obs.fiducialcosmo.comoving(z) ** 2) / self.pk_obs.fiducialcosmo.Hubble(z)
cosmo = (
(1 + z) ** 2 * self.pk_obs.fiducialcosmo.comoving(z) ** 2
) / self.pk_obs.fiducialcosmo.Hubble(z)
T_term = (self.Tsys_func(z) / temp) ** 2 # in K
alpha = self.pk_obs.alpha_SD()
if beam_term:
beta = self.pk_obs.beta_SD(z, k, mu)
else:
beta = np.ones_like(k)
noise = pref * cosmo * T_term * (alpha / beta**2)
noise = pref * cosmo * T_term * (alpha / beta ** 2)
return noise

def veff_II(self, zi, k, mu):
Expand All @@ -300,7 +301,7 @@ def veff_II(self, zi, k, mu):
"""
pobs = self.pk_obs.observed_P_ij(zi, k, mu, si="I", sj="I")
pnoisy = self.noisy_P_ij(zi, k, mu, si="I", sj="I")
prefactor = 1 / (8 * (np.pi**2))
prefactor = 1 / (8 * (np.pi ** 2))
covterm = prefactor * (pobs / pnoisy) ** 2
if zi < self.inter_z_bin_mids[0] or zi > self.inter_z_bin_mids[-1]:
covterm = np.zeros_like(covterm)
Expand Down Expand Up @@ -329,24 +330,25 @@ def veff_Ig(self, zi, k, mu):
pnoisy_Ig = self.noisy_P_ij(zi, k, mu, si="I", sj="g")
pnoisy_II = self.noisy_P_ij(zi, k, mu, si="I", sj="I")
pnoisy_gg = self.noisy_P_ij(zi, k, mu, si="g", sj="g")
covterm = pobs_Ig**2 / (pnoisy_gg * pnoisy_II + pnoisy_Ig * pnoisy_Ig)
prefactor = 1 / (4 * (np.pi**2))
covterm = pobs_Ig ** 2 / (pnoisy_gg * pnoisy_II + pnoisy_Ig * pnoisy_Ig)
prefactor = 1 / (4 * (np.pi ** 2))
covterm = prefactor * covterm
if zi < self.inter_z_bin_mids[0] or zi > self.inter_z_bin_mids[-1]:
covterm = np.zeros_like(covterm)
return covterm

def noisy_P_ij(self, z, k, mu, si="I", sj="g"):
if si == "I" and sj == "I":
noiseterm = self.P_noise_21(z, k, mu, temp_dim=True)
elif si == "g" and sj == "g":
noiseterm = 1/self.n_density(z)
noiseterm = 1 / self.n_density(z)
else:
noiseterm = 0
pobs_ij = self.pk_obs.observed_P_ij(z, k, mu, si=si, sj=sj)
pnoisy_ij = pobs_ij + noiseterm
return pnoisy_ij


class SpectroDerivs:
def __init__(
self,
Expand Down Expand Up @@ -440,16 +442,16 @@ def initialize_obs(self, allpars):
else:
IMbiaspars = None
self.pobs = spec_obs.ComputeGalSpectro(
cosmopars=cosmopars,
fiducial_cosmopars=self.fiducial_cosmopars,
spectrobiaspars=spectrobiaspars,
spectrononlinearpars=spectrononlinearpars,
PShotpars=PShotpars,
IMbiaspars=IMbiaspars,
fiducial_cosmo=self.fiducial_cosmo,
bias_samples=self.bias_samples,
configuration=self.config,
)
cosmopars=cosmopars,
fiducial_cosmopars=self.fiducial_cosmopars,
spectrobiaspars=spectrobiaspars,
spectrononlinearpars=spectrononlinearpars,
PShotpars=PShotpars,
IMbiaspars=IMbiaspars,
fiducial_cosmo=self.fiducial_cosmo,
bias_samples=self.bias_samples,
configuration=self.config,
)
strdic = str(sorted(cosmopars.items()))
hh = hash(strdic)
self.cosmology_variations_dict[hh] = self.pobs.cosmo
Expand All @@ -473,11 +475,15 @@ def get_obs(self, allpars):
result_array["z_bins"] = self.z_array
for ii, zzi in enumerate(self.z_array):
if self.bias_samples == ["I", "I"]:
result_array[ii] = self.pobs.lnpobs_ij(zzi, self.pk_kmesh, self.pk_mumesh, si="I", sj="I")
result_array[ii] = self.pobs.lnpobs_ij(
zzi, self.pk_kmesh, self.pk_mumesh, si="I", sj="I"
)
elif self.bias_samples == ["g", "g"]:
result_array[ii] = self.pobs.lnpobs_gg(zzi, self.pk_kmesh, self.pk_mumesh)
elif self.bias_samples == ["I", "g"] or self.bias_samples == ["g", "I"]:
result_array[ii] = self.pobs.lnpobs_ij(zzi, self.pk_kmesh, self.pk_mumesh, si="I", sj="g")
result_array[ii] = self.pobs.lnpobs_ij(
zzi, self.pk_kmesh, self.pk_mumesh, si="I", sj="g"
)
return result_array

def exact_derivs(self, par):
Expand Down
32 changes: 17 additions & 15 deletions cosmicfishpie/LSSsurvey/spectro_obs.py
Original file line number Diff line number Diff line change
Expand Up @@ -166,7 +166,7 @@ def __init__(
self.IMbiaspars = self.fiducial_IMbiaspars
else:
self.IMbiaspars = IMbiaspars

self.nuisance = nuisance.Nuisance(
configuration=self.config,
spectrobiasparams=self.spectrobiaspars,
Expand Down Expand Up @@ -211,7 +211,7 @@ def __init__(
if "IM" in self.observables:
self.set_IM_specs()
self.set_IM_bias_specs()

tend = time()
upt.time_print(
feedback_level=self.feed_lvl,
Expand Down Expand Up @@ -376,7 +376,7 @@ def kper(self, z, k, mu):
Observed perpendicular projection of wavevector onto the line of sight with AP-effect corrected for
"""

k_per = k * np.sqrt(1 - mu**2) * (1 / self.qperpendicular(z))
k_per = k * np.sqrt(1 - mu ** 2) * (1 / self.qperpendicular(z))
return k_per

def k_units_change(self, k):
Expand Down Expand Up @@ -471,7 +471,7 @@ def spec_err_z(self, z, k, mu):
elif self.dz_type == "z-dependent":
spec_dz_err = self.dz_err * (1 + z)
err = spec_dz_err * (1 / self.cosmo.Hubble(z)) * self.kpar(z, k, mu)
return np.exp(-(1 / 2) * err**2)
return np.exp(-(1 / 2) * err ** 2)

def BAO_term(self, z):
"""Calculates the BAO term. This is the rescaling of the Fourier volume by the AP-effect
Expand Down Expand Up @@ -538,7 +538,7 @@ def bterm_fid(self, z, k=None, bias_sample="g"):
raise ValueError(
f"Bias sample {bias_sample} not found. "
f"Please use {self.IM_bias_sample} bias sample."
)
)
bterm_z = self.nuisance.IM_bias_at_z(z)
bterm_k = self.nuisance.gcsp_bias_kscale(k)
bterm_zk = bterm_z * bterm_k
Expand Down Expand Up @@ -593,9 +593,9 @@ def kaiserTerm(self, z, k, mu, b_i=None, just_rsd=False, bias_sample="g"):
fterm = self.cosmo.f_growthrate(z, k, tracer=self.tracer)

if not just_rsd:
kaiser = bterm + fterm * mu**2
kaiser = bterm + fterm * mu ** 2
elif just_rsd:
kaiser = 1 + (fterm / bterm) * mu**2
kaiser = 1 + (fterm / bterm) * mu ** 2

return kaiser

Expand Down Expand Up @@ -679,7 +679,7 @@ def sigmavNL(self, zz, mu):
f0 = self.P_ThetaTheta_Moments(zz, 0)
f1 = self.P_ThetaTheta_Moments(zz, 1)
f2 = self.P_ThetaTheta_Moments(zz, 2)
sv = np.sqrt(f0 + 2 * mu**2 * f1 + mu**2 * f2)
sv = np.sqrt(f0 + 2 * mu ** 2 * f1 + mu ** 2 * f2)
if self.vary_sigmav:
sv *= self.nuisance.vectorized_gcsp_rescale_sigmapv_at_z(zz, sigma_key="sigmav")
return sv
Expand Down Expand Up @@ -714,7 +714,7 @@ def f_mom(k):
pp = cosmoF.matpow(zz, self.k_grid).flatten()
integrand = pp * ff
Int = integrate.trapezoid(integrand, x=self.k_grid)
ptt = (1 / (6 * np.pi**2)) * Int
ptt = (1 / (6 * np.pi ** 2)) * Int
return ptt

def normalized_pdd(self, z, k):
Expand Down Expand Up @@ -800,8 +800,8 @@ def dewiggled_pdd(self, z, k, mu):

self.p_dd = self.normalized_pdd(z, k)
self.p_dd_NW = self.normalized_pnw(z, k)
self.p_dd_DW = self.p_dd * np.exp(-gmudamping * k**2) + self.p_dd_NW * (
1 - np.exp(-gmudamping * k**2)
self.p_dd_DW = self.p_dd * np.exp(-gmudamping * k ** 2) + self.p_dd_NW * (
1 - np.exp(-gmudamping * k ** 2)
)
return self.p_dd_DW

Expand Down Expand Up @@ -860,7 +860,7 @@ def observed_Pgg(self, z, k, mu, b_i=None):
lorentzFoG = self.FingersOfGod(z, k, mu, mode="Lorentz")
p_dd_DW = self.dewiggled_pdd(z, k, mu)

pgg_obs = baoterm * (kaiser**2) * p_dd_DW * lorentzFoG * (error_z**2) + extra_shotnoise
pgg_obs = baoterm * (kaiser ** 2) * p_dd_DW * lorentzFoG * (error_z ** 2) + extra_shotnoise

tend = time()
upt.time_print(
Expand Down Expand Up @@ -940,7 +940,7 @@ def beta_SD(self, z, k, mu):
tol = 1.0e-12
k = np.atleast_1d(k)
mu = np.atleast_1d(mu)
expo = k**2 * (1 - mu**2) * self.fiducialcosmo.comoving(z) ** 2 * self.theta_b(z) ** 2
expo = k ** 2 * (1 - mu ** 2) * self.fiducialcosmo.comoving(z) ** 2 * self.theta_b(z) ** 2
bet = np.exp(-expo / (16.0 * np.log(2.0)))
bet[np.abs(bet) < tol] = tol
return bet
Expand All @@ -949,9 +949,11 @@ def observed_P_ij(self, z, k, mu, bsi_z=None, bsj_z=None, si="I", sj="g"):
error_z = self.spec_err_z(z, k, mu)
beam_damping_term_si = self.beta_SD(z, k, mu) if si == "I" else 1
beam_damping_term_sj = self.beta_SD(z, k, mu) if sj == "I" else 1
k = self.k_units_change(k) # h-bug set to False by default, leaving here for cross-check of old cases
k = self.k_units_change(
k
) # h-bug set to False by default, leaving here for cross-check of old cases
k, mu = self.kmu_alc_pac(z, k, mu)
#if self.bias_samples is not None:
# if self.bias_samples is not None:
# si = self.bias_samples[0]
# sj = self.bias_samples[1]
baoterm = self.BAO_term(z)
Expand Down
43 changes: 18 additions & 25 deletions cosmicfishpie/analysis/fishconsumer.py
Original file line number Diff line number Diff line change
Expand Up @@ -101,7 +101,7 @@ def display_colors(colors, figsize=(6, 6)):
# rotation_angle = angle
x1, y1 = wedges[i].center
x2, y2 = np.cos(np.deg2rad(angle)), np.sin(np.deg2rad(angle))
dx, dy = 1.2 * np.array([x2, y2]) / np.sqrt(x2**2 + y2**2)
dx, dy = 1.2 * np.array([x2, y2]) / np.sqrt(x2 ** 2 + y2 ** 2)
ax.annotate(
color,
xy=(x1, y1),
Expand Down Expand Up @@ -627,8 +627,7 @@ def prepare_settings_plot(


def chainfishplot(
return_dictionary,
**cckwargs,
return_dictionary, **cckwargs,
):
"""
Chain fish plot function
Expand Down Expand Up @@ -776,6 +775,7 @@ def chainfishplot(
print("Plot saved to: ", plotfilename)
return fig


def simple_fisher_plot(
fisher_list,
params_to_plot,
Expand All @@ -784,7 +784,7 @@ def simple_fisher_plot(
save_plot=False,
legend=True,
n_samples=10000,
output_file="fisher_plot.pdf"
output_file="fisher_plot.pdf",
):
"""Create a triangle plot from Fisher matrices using ChainConsumer.
Expand All @@ -810,34 +810,27 @@ def simple_fisher_plot(
"""
# Initialize ChainConsumer
c = ChainConsumer()

# Default colors if none provided
if colors is None:
colors = ['#3a86ff', '#fb5607', '#8338ec', '#ffbe0b', '#d11149']
colors = colors[:len(fisher_list)] # Truncate to needed length
colors = ["#3a86ff", "#fb5607", "#8338ec", "#ffbe0b", "#d11149"]
colors = colors[: len(fisher_list)] # Truncate to needed length

# Default labels if none provided
if labels is None:
labels = [f"Fisher {i+1}" for i in range(len(fisher_list))]

# Generate samples for each Fisher matrix
n_samples = 100000
for i, fisher in enumerate(fisher_list):
# Get samples from multivariate normal using Fisher matrix
samples = multivariate_normal(
fisher.param_fiducial,
fisher.inverse_fisher_matrix(),
size=n_samples
fisher.param_fiducial, fisher.inverse_fisher_matrix(), size=n_samples
)

# Add chain to plot
c.add_chain(
samples,
parameters=fisher.get_param_names(),
name=labels[i],
color=colors[i]
)

c.add_chain(samples, parameters=fisher.get_param_names(), name=labels[i], color=colors[i])

# Configure plot settings
c.configure(
plot_hists=True,
Expand All @@ -848,15 +841,15 @@ def simple_fisher_plot(
shade_alpha=0.3,
bar_shade=True,
linewidths=2,
legend_kwargs={"fontsize": 12}
legend_kwargs={"fontsize": 12},
)

# Create the plot
fig = c.plotter.plot(parameters=params_to_plot, legend=legend)

# Save if requested
if save_plot:
fig.savefig(output_file, bbox_inches='tight', dpi=200)
fig.savefig(output_file, bbox_inches="tight", dpi=200)
print(f"Plot saved to: {output_file}")

return fig
Loading

0 comments on commit 4b73bd2

Please sign in to comment.