Source code for wale.RateFunction

import numpy as np
from scipy.optimize import root
from numpy import newaxis


[docs] def get_tau(rho): """ Compute the large-deviation theory mapping τ(ρ) for a given density contrast ρ. Parameters ---------- rho : float or array_like Local density ρ = 1 + δ. Returns ------- tau : float or ndarray The corresponding τ(ρ) value based on the ν parameter. """ nu = 1.4 return nu * (1.0 - rho ** (-1.0 / nu))
[docs] def get_psi_2cell(variance, chi, recal, z, delta1, delta2, theta1, theta2): """ Compute the 2-cell action ψ(δ₁, δ₂) using large-deviation theory. Parameters ---------- variance : object Object providing nonlinear_sigma2 method for computing variances. chi : float or ndarray Comoving distance(s). recal : float Empirical recalibration factor. z : float Redshift. delta1, delta2 : float Density contrasts in each cell. theta1, theta2 : float Angular scales of the two cells. Returns ------- psi : float Value of the large-deviation action ψ(δ₁, δ₂). """ delta1 = np.array([delta1]) delta2 = np.array([delta2]) tau1sq = get_tau(1.0 + delta1) * get_tau(1.0 + delta1) tau2sq = get_tau(1.0 + delta2) * get_tau(1.0 + delta2) tau12sq = get_tau(1.0 + delta1) * get_tau(1.0 + delta2) # sig2lR11 = 1.0 sig2lr12 = variance.nonlinear_sigma2( redshift=z, R1=chi * ((1 + delta1[:, newaxis]) ** 0.5) * theta1, R2=chi * ((1 + delta2[:, newaxis]) ** 0.5) * theta2, ) sig2lr11 = variance.nonlinear_sigma2( redshift=z, R1=chi * ((1 + delta1[:, newaxis]) ** 0.5) * theta1, R2=chi * ((1 + delta1[:, newaxis]) ** 0.5) * theta1, ) sig2lr22 = variance.nonlinear_sigma2( redshift=z, R1=chi * ((1 + delta2[:, newaxis]) ** 0.5) * theta2, R2=chi * ((1 + delta2[:, newaxis]) ** 0.5) * theta2, ) det = (sig2lr11 * sig2lr22) - (sig2lr12 * sig2lr12) psi = ( (sig2lr11 * tau2sq - 2.0 * sig2lr12 * tau12sq + sig2lr22 * tau1sq) * recal / (det * 2.0) ) return psi
[docs] def get_phi_projec_2cell( theta1, theta2, zarr, chis, dchis, w, y, recal, variance, **kwargs ): """ Compute projected 2-cell φ(y) by solving saddle-point equations. Parameters ---------- theta1, theta2 : float Angular scales of the two cells. zarr : ndarray Redshifts at which to evaluate. chis : ndarray Comoving distances corresponding to redshifts. dchis : float Integration step size in χ. w : ndarray Lensing weight function W(χ). y : ndarray Slope values for SCGF (y-grid). recal : float Empirical recalibration factor. variance : object Provides nonlinear_sigma2(χ) interface. deld : float, optional Step size for finite differences (default: 1e-8). Returns ------- phi_proj : ndarray Projected φ(y) values computed via saddle-point approximation. """ deld = kwargs.get("deld", 1e-8) nchi = len(chis) ny = len(y) def to_solve2(delta, A, recal, chi, z): delta_1, delta_2 = delta psi_at_delta = get_psi_2cell( variance, chi, recal=recal, z=z, delta1=delta_1, delta2=delta_2, theta1=theta1, theta2=theta2, ) psi_at_delta1_plus_epsilon = get_psi_2cell( variance, chi=chi, recal=recal, z=z, delta1=delta_1 + deld, delta2=delta_2, theta1=theta1, theta2=theta2, ) psi_at_delta2_plus_epsilon = get_psi_2cell( variance, chi=chi, recal=recal, z=z, delta1=delta_1, delta2=delta_2 + deld, theta1=theta1, theta2=theta2, ) psi_derivative_delta1 = psi_at_delta1_plus_epsilon - psi_at_delta psi_derivative_delta2 = psi_at_delta2_plus_epsilon - psi_at_delta equation1 = A + psi_derivative_delta1 equation2 = -A + psi_derivative_delta2 return [equation1[0], equation2[0]] delta_zeros = np.zeros((nchi, ny, 2)) for i in range(nchi): print(f"Iteration {i*100/nchi} %", end="\r") for j in range(ny): A = y[j] * w[i] * deld if j == 0: x0 = 0.0 y0 = 0.0 else: x0 = delta_zeros[i, j - 1][0] y0 = delta_zeros[i, j - 1][1] delta_zeros[i, j] = root( to_solve2, [x0, y0], args=(A, recal, np.array([chis[i]]), zarr[i]), method="hybr", ).x phi_proj = [] for i in range(ny): phi_ = 0.0 for j in range(nchi): phi_ += ( y[i] * w[j] * (-delta_zeros[j, i, 0] + delta_zeros[j, i, 1]) - get_psi_2cell( variance, chi=np.array([chis[j]]), recal=recal, z=zarr[j], delta1=delta_zeros[j, i, 0], delta2=delta_zeros[j, i, 1], theta1=theta1, theta2=theta2, ) ) * dchis # [j] phi_proj.append(phi_) phi_proj = np.array(phi_proj) return phi_proj
[docs] def get_scaled_cgf( theta1, theta2, zarr, chis, dchis, lensing_weight, y, recal, variance ): """ Wrapper to compute the scaled cumulant generating function (SCGF). Parameters ---------- theta1, theta2 : float Angular scales. zarr : ndarray Redshift values. chis : ndarray Comoving distances. dchis : float Integration step size. lensing_weight : ndarray Lensing kernel W(χ). y : ndarray SCGF slope values. recal : float Recalibration factor. variance : object Provides nonlinear sigma²(R₁, R₂, z). Returns ------- scgf : ndarray The scaled cumulant generating function φ(y). """ scgf = get_phi_projec_2cell( theta1, theta2, zarr, chis, dchis, lensing_weight, y, recal, variance ) return scgf
[docs] def get_psi_derivative_delta1( deld, variance, chi, recal, z, delta1, delta2, theta1, theta2 ): """ Compute ∂ψ/∂δ₁ using central finite differences. Returns ------- derivative : float Numerical partial derivative of ψ with respect to δ₁. """ delh = deld # Small step size for numerical differentiation delta1_plus_h = delta1 + delh delta1_minus_h = delta1 - delh psi_plus_h = get_psi_2cell( variance, chi, recal, z, delta1_plus_h, delta2, theta1, theta2 ) psi_minus_h = get_psi_2cell( variance, chi, recal, z, delta1_minus_h, delta2, theta1, theta2 ) derivative = (psi_plus_h - psi_minus_h) / (2 * delh) return derivative
[docs] def get_psi_derivative_delta2( deld, variance, chi, recal, z, delta1, delta2, theta1, theta2 ): """ Compute ∂ψ/∂δ₂ using central finite differences. Returns ------- derivative : float Numerical partial derivative of ψ with respect to δ₂. """ delh = deld # Small step size for numerical differentiation delta2_plus_h = delta2 + delh delta2_minus_h = delta2 - delh psi_plus_h = get_psi_2cell( variance, chi, recal, z, delta1, delta2_plus_h, theta1, theta2 ) psi_minus_h = get_psi_2cell( variance, chi, recal, z, delta1, delta2_minus_h, theta1, theta2 ) derivative = (psi_plus_h - psi_minus_h) / (2 * delh) return derivative
[docs] def get_psi_2nd_derivative_delta1( deld, variance, chi, recal, z, delta1, delta2, theta1, theta2 ): """ Compute ∂²ψ/∂δ₁² using second-order central finite differences. Returns ------- second_derivative : float Second partial derivative of ψ with respect to δ₁. """ delh = deld delta1_plus_h = delta1 + delh delta1_minus_h = delta1 - delh psi_plus_h = get_psi_2cell( variance, chi, recal, z, delta1_plus_h, delta2, theta1, theta2 ) psi_minus_h = get_psi_2cell( variance, chi, recal, z, delta1_minus_h, delta2, theta1, theta2 ) psi_at_delta1 = get_psi_2cell( variance, chi, recal, z, delta1, delta2, theta1, theta2 ) second_derivative = (psi_plus_h - (2.0 * psi_at_delta1) + psi_minus_h) / (delh**2.0) return second_derivative
[docs] def get_psi_2nd_derivative_delta2( deld, variance, chi, recal, z, delta1, delta2, theta1, theta2 ): """ Compute ∂²ψ/∂δ₂² using second-order central finite differences. Returns ------- second_derivative : float Second partial derivative of ψ with respect to δ₂. """ delh = deld delta2_plus_h = delta2 + delh delta2_minus_h = delta2 - delh psi_plus_h = get_psi_2cell( variance, chi, recal, z, delta1, delta2_plus_h, theta1, theta2 ) psi_minus_h = get_psi_2cell( variance, chi, recal, z, delta1, delta2_minus_h, theta1, theta2 ) psi_at_delta2 = get_psi_2cell( variance, chi, recal, z, delta1, delta2, theta1, theta2 ) second_derivative = (psi_plus_h - (2.0 * psi_at_delta2) + psi_minus_h) / (delh**2.0) return second_derivative
[docs] def get_psi_mixed_derivative_delta1_delta2( deld, variance, chi, recal, z, delta1, delta2, theta1, theta2 ): """ Compute the mixed partial derivative ∂²ψ/∂δ₁∂δ₂ using central differences. Returns ------- mixed_derivative : float Mixed second-order partial derivative of ψ. """ delh = deld delta1_plus_h = delta1 + delh delta1_minus_h = delta1 - delh delta2_plus_h = delta2 + delh delta2_minus_h = delta2 - delh psi_delta1_plus_h = get_psi_2cell( variance, chi, recal, z, delta1_plus_h, delta2_plus_h, theta1, theta2 ) psi_delta1_minus_h = get_psi_2cell( variance, chi, recal, z, delta1_minus_h, delta2_plus_h, theta1, theta2 ) psi_delta2_plus_h = get_psi_2cell( variance, chi, recal, z, delta1_plus_h, delta2_minus_h, theta1, theta2 ) psi_delta2_minus_h = get_psi_2cell( variance, chi, recal, z, delta1_minus_h, delta2_minus_h, theta1, theta2 ) mixed_derivative = ( psi_delta1_plus_h - psi_delta1_minus_h - psi_delta2_plus_h + psi_delta2_minus_h ) / (4.0 * delh * delh) return mixed_derivative
[docs] def psi_derivative_determinant( deld, delta1, delta2, z, variance, chi, recal, theta1, theta2 ): """ Compute the determinant of the Hessian matrix of ψ(δ₁, δ₂): det(H) = ∂²ψ/∂δ₁² * ∂²ψ/∂δ₂² − (∂²ψ/∂δ₁∂δ₂)² This determinant is used for computing the normalization in saddle-point approximations. Returns ------- result : float Determinant of the ψ Hessian matrix. """ psi_11 = get_psi_2nd_derivative_delta1( deld, variance, chi, recal, z, delta1, delta2, theta1, theta2 ) psi_22 = get_psi_2nd_derivative_delta2( deld, variance, chi, recal, z, delta1, delta2, theta1, theta2 ) psi_12 = get_psi_mixed_derivative_delta1_delta2( deld, variance, chi, recal, z, delta1, delta2, theta1, theta2 ) result = (psi_11 * psi_22) - (psi_12**2.0) return result