diff --git a/examples/example_usage.py b/examples/example_usage.py index be7d947..2c3c1ef 100644 --- a/examples/example_usage.py +++ b/examples/example_usage.py @@ -92,6 +92,36 @@ def run_examples(): ) print("Example 4 completed.\n") + # Example 5: E1 topology with modified initial power spectrum (PS_mod=True, E18_mod=True) + print("Running Example 5: E1 Topology with modified power spectrum and E18_mod=True") + + run_topology( + topology='E1', + l_max=5, + Lx=1.0, + Ly=1.0, + Lz=1.0, + beta=90, + alpha=90, + gamma=0, + do_polarization=False, + normalize=True, + l_range=np.array([[2, 5]]), + lp_range=np.array([[2, 5]]), + + # --- Power spectrum modification parameters --- + PS_mod=True, # Enable modification of the initial power spectrum + E18_mod=True, # Apply the same modification to E18 in KL calculations + powerspec='wavepacket', # Type of modification + amp=1.0, # Overall amplitude + width=1.5, # Width parameter (used by selected models) + freq=10, # Frequency parameter (used by oscillatory models) + x_cutoff=1.0 # Cutoff scale (within recommended range) + ) + + print("Example 5 completed.\n") + + if __name__ == '__main__': run_examples() print("The results for these examples have been saved in ./CMBtopology/runs/. " \ diff --git a/tests/test_run_topology.py b/tests/test_run_topology.py index 74a1318..440e915 100644 --- a/tests/test_run_topology.py +++ b/tests/test_run_topology.py @@ -75,5 +75,39 @@ def test_missing_required_parameter(self): with self.assertRaises(TypeError): run_topology(topology='E1') # Missing l_max + def test_ps_mod_allowed_topology(self): + """PS_mod=True with an allowed topology should not raise errors.""" + try: + run_topology( + topology='E1', + l_max=15, + Lx=1.0, + Ly=1.0, + Lz=1.0, + do_polarization=False, + l_range=np.array([[2, 15]]), + lp_range=np.array([[2, 15]]), + PS_mod=True, + powerspec='cutoff', + amp=1.0, + width=1.5, + freq=10, + x_cutoff=1.0 + ) + except Exception as e: + self.fail(f"Unexpected error with PS_mod on allowed topology: {e}") + + def test_ps_mod_disallowed_topology_raises(self): + """PS_mod=True with a disallowed topology should raise ValueError.""" + with self.assertRaises(ValueError): + run_topology( + topology='E8', #not yet implemented + l_max=20, + PS_mod=True, + powerspec='cutoff' + ) + + + if __name__ == '__main__': unittest.main() \ No newline at end of file diff --git a/topology/parameter_files/default_PS.py b/topology/parameter_files/default_PS.py new file mode 100644 index 0000000..16a120e --- /dev/null +++ b/topology/parameter_files/default_PS.py @@ -0,0 +1,52 @@ +# Parameter file +# Specify the initial power spectrum of the non-trivial topology (E18 uses the standard power law) +import numpy as np + +power_parameter = { + # POWERSPECTRUM PARAMETERS + 'PS_mod': False, # Default: no modification to the initial power spectrum + 'E18_mod': False, # False/True means a given modification is applied to the non-trivial topology only/also E18 in the KL calculation + 'powerspec': 'powlaw', + 'amp': 1.0, + 'width': 1.5, + 'freq': 10, + 'x_cutoff': 1, +} + +''' +Available power spectra, their name and free parameters: + +Standard power law (default) + 'powlaw' + As*(k/0.05)**(ns-1) + +Wavepacket (made to fit with the shape of the Planck posterior) + 'wavepacket' + As*(k/0.05)**(ns-1)*(1+ np.sin(x*freq)*amp*np.exp(-(x/xc))) + parameters: + 'amp' Amplitude + 'freq' Frecuency of the oscillation + 'width' The xc, controls the exponential surpression of the oscillation + +Exponential cutoff + 'cutoff' + As*(k/0.05)**(ns-1)*(1-amp*np.exp(-(x/xc))) + parameters: + 'x_cutoff' Position of the cutoff + 'amp' Larger value means a steeper cut-off + +Enhancement + 'enhance' + As*(k/0.05)**(ns-1)*(1+amp*np.exp(-(x/xc))) + parameters: + 'x_cutoff' Position of the rise + 'amp' Larger value means a steeper increase + +Logarithmic oscillation model + 'logosci' + As*(k/0.05)**(ns-1)*(1+amp*np.cos(np.log(k/0.05)*freq)) + parameters: + 'amp' Amplitude + 'freq' Frequency + +''' \ No newline at end of file diff --git a/topology/run_topology.py b/topology/run_topology.py index 8dfd5b6..e6f70e7 100644 --- a/topology/run_topology.py +++ b/topology/run_topology.py @@ -8,6 +8,7 @@ import argparse import importlib import numpy as np +import warnings from .src.E1 import E1 from .src.E2 import E2 from .src.E3 import E3 @@ -27,6 +28,8 @@ def run_topology( normalize=False, l_range=None, lp_range=None, + PS_mod=False, + x_cutoff=None, **topology_params ): """Run covariance matrix computation for a specified topology. @@ -64,9 +67,23 @@ def run_topology( lp_range = np.array([[l_min, l_max]]) # Dynamically import the parameter file + # --- Load global PS defaults --- try: - parameter_module = importlib.import_module(f".parameter_files.default_{topology}", package="topology") - param = parameter_module.parameter + ps_module = importlib.import_module( + ".parameter_files.default_PS", + package="topology", + ) + param = ps_module.power_parameter.copy() + except ImportError: + raise ValueError( + "Power spectrum parameter file not found. " + "Expected: topology/parameter_files/default_PS.py" + ) + + # --- Load topology-specific defaults --- + try: + topology_module = importlib.import_module(f".parameter_files.default_{topology}", package="topology") + param.update(topology_module.parameter) except ImportError: raise ValueError(f"Parameter file for topology {topology} not found. Expected: topology/parameter_files/default_{topology}.py") @@ -74,6 +91,26 @@ def run_topology( for key in topology_params: if key not in param: raise ValueError(f"Invalid parameter '{key}' for topology '{topology}'. Valid parameters are: {list(param.keys())}") + + PS_ALLOWED_TOPOLOGIES = {"E1", "E2", "E3", "E4", "E5", "E6"} + # --- Validate PS_mod vs topology --- + if PS_mod and topology not in PS_ALLOWED_TOPOLOGIES: + raise ValueError( + f"PS_mod=True is only supported for topologies {sorted(PS_ALLOWED_TOPOLOGIES)}. " + f"Received topology='{topology}'." + ) + + X_CUTOFF_MIN = 0.0 + X_CUTOFF_MAX = 4.0 + # --- Validate x_cutoff range --- + if x_cutoff is not None: + if not (X_CUTOFF_MIN <= x_cutoff <= X_CUTOFF_MAX): + warnings.warn( + f"x_cutoff={x_cutoff} is outside the recommended range " + f"[{X_CUTOFF_MIN}, {X_CUTOFF_MAX}]. Normalization and k_max will probably not be accurate." + "Proceeding anyway.", + UserWarning + ) # Update parameters param.update(topology_params) @@ -169,14 +206,59 @@ def main(): help="Observer position as three floats (default: [0.0, 0.0, 0.0])" ) + # Allow primordial power spectrum specific parameters as optional arguments + parser.add_argument("--PS_mod", + action="store_true", + help="Include a non-standard primordial power spectrum (implemented for E1-E6 only)" + ) + + parser.add_argument("--E18_mod", + action="store_true", + help="Use the non-standard primordial power spectrum also for E18 in the KL divergence computation" + ) + + parser.add_argument( + "--powerspec", + type=str, + choices=["powlaw", "wavepacket", "cutoff", "enhance"], + help="Type of power spectrum modification (powlaw, wavepacket, cutoff, enhance)" + ) + + parser.add_argument("--amp", type=float, help="Amplitude of the modification") + parser.add_argument("--width", type=float, help="Width parameter for modified power spectrum") + parser.add_argument("--freq",type=float,help="Frequency parameter for modified power spectrum") + parser.add_argument("--x_cutoff",type=float,help="Modification cutoff scale in terms of x=kL/2pi") + args = parser.parse_args() - # Collect topology-specific parameters + PS_ALLOWED_TOPOLOGIES = {"E1", "E2", "E3", "E4", "E5", "E6"} + # --- Validate PS_mod vs topology --- + if args.PS_mod and args.topology not in PS_ALLOWED_TOPOLOGIES: + raise ValueError( + f"--PS_mod is only supported for topologies {sorted(PS_ALLOWED_TOPOLOGIES)}. " + f"Received topology='{args.topology}'." + ) + + X_CUTOFF_MIN = 0.0 + X_CUTOFF_MAX = 4.0 + # --- Validate x_cutoff range --- + if args.x_cutoff is not None: + if not (X_CUTOFF_MIN <= args.x_cutoff <= X_CUTOFF_MAX): + warnings.warn( + f"x_cutoff={args.x_cutoff} is outside the recommended range " + f"[{X_CUTOFF_MIN}, {X_CUTOFF_MAX}]. Normalization and k_max will probably not be accurate." + "Proceeding anyway.", + UserWarning + ) + + + # Collect topology-specific parameters and primordial power spectrum parameters topology_params = {} for param in [ 'Lx', 'Ly', 'Lz', 'beta', 'alpha', 'gamma', 'r_x', 'r_y', 'r_z', - 'LAx', 'LAy', 'L1y', 'L2x', 'L2z', 'LBx', 'LBz', 'LCy', 'x0' + 'LAx', 'LAy', 'L1y', 'L2x', 'L2z', 'LBx', 'LBz', 'LCy', 'x0', + 'powerspec', 'amp', 'width', 'freq', 'x_cutoff' ]: if hasattr(args, param) and getattr(args, param) is not None: topology_params[param] = getattr(args, param) @@ -192,6 +274,7 @@ def main(): l_min=args.l_min, do_polarization=args.do_polarization, normalize=args.normalize, + PS_mod=args.PS_mod, **topology_params ) diff --git a/topology/src/config.py b/topology/src/config.py index 816fe46..01c56e2 100644 --- a/topology/src/config.py +++ b/topology/src/config.py @@ -3,4 +3,8 @@ """ # Physical constants -L_LSS = 13824.9 * 2 # Last scattering surface diameter (Mpc) \ No newline at end of file +L_LSS = 13824.9 * 2 # Last scattering surface diameter (Mpc) + +# Primordial power spectrum +A_s = 2e-9 # Amplitude +n_s = 0.965 # Spectral index \ No newline at end of file diff --git a/topology/src/topology.py b/topology/src/topology.py index e90ac19..82bc41c 100755 --- a/topology/src/topology.py +++ b/topology/src/topology.py @@ -7,6 +7,7 @@ """ import camb +import copy import os import numpy as np from numpy import pi, sqrt @@ -19,18 +20,22 @@ import time from .tools import * from sys import getsizeof +from .config import L_LSS, A_s, n_s class Topology: """Base class for CMB covariance matrix computation in non-trivial topologies. Attributes: - param (dict): Configuration parameters (e.g., topology, l_max, c_l_accuracy). + param (dict): Configuration parameters (e.g., topology, l_max, c_l_accuracy, PS_mod, powerspec). topology (str): Topology type (e.g., 'E1'--'E10'). l_max (int): Maximum multipole. l_min (int): Minimum multipole. c_l_accuracy (float): Accuracy for wavevector cutoff (e.g., 0.99). C_l_type_array (np.ndarray): Array of correlation types ['TT', 'EE', 'TE']. do_polarization (bool): Whether to compute polarization (EE, TE). + PS_mod (bool): Whether to use a non-standard primordial power spectrum + E18_mod (bool): Whether the non-standard prim. PS is also used for E18 in the KL div. computation + powerspec (str): The type of prim. PS modification fig_name (str): Base name for output figures. debug (bool): Enable debug output. make_run_folder (bool): Create output directories. @@ -71,6 +76,13 @@ def __init__(self, param, debug=True, make_run_folder = True): self.l_max = param['l_max'] self.l_min = param['l_min'] self.c_l_accuracy = param['c_l_accuracy'] + self.PS_mod = param['PS_mod'] + self.E18_mod = param['E18_mod'] + self.powerspec = param['powerspec'] + self.amp = param['amp'] + self.width = param['width'] + self.freq = param['freq'] + self.x_cutoff = param['x_cutoff'] C_l_type_array= np.array(['TT','EE','TE']) self.C_l_type_array = C_l_type_array self.do_polarization = param['do_polarization'] @@ -107,7 +119,7 @@ def do_pre_processing(self): # Configure CAMB parameters pars = camb.CAMBparams() pars.set_cosmology(H0=67.5, ombh2=0.022, omch2=0.122, mnu=0.06, omk=0, tau=0.06) - pars.InitPower.set_params(As=2e-9, ns=0.965, r=0) + pars.InitPower.set_params(As=A_s, ns=n_s, r=0) pars.set_for_lmax(self.l_max) # More accurate transfer functions. Takes longer time to run @@ -115,14 +127,57 @@ def do_pre_processing(self): pars.Accuracy.IntkAccuracyBoost = 2 pars.Accuracy.SourcekAccuracyBoost = 2 - # Compute transfer functions and power spectra + # Compute transfer functions and power spectra for the standard power law data = camb.get_transfer_functions(pars) results = camb.get_results(pars) self.powers = results.get_cmb_power_spectra(pars, raw_cl=True, CMB_unit='muK')['unlensed_scalar'] transfer_function = data.get_cmb_transfer_data(tp='scalar') + # Compute the non-standard primordial power spectrum + # reference scale is the parameter L_mod, the minimum distance between any point in the manifold and its nearest clone + l_min = min(self.param['Lx'], self.param['Ly'], self.param['Lz']) + if self.topology in ['E1', 'E2', 'E3', 'E4', 'E5']: + L_PS = l_min + elif self.topology=='E6': + L_PS = np.sqrt(2)*l_min + else: + raise ValueError(f"Topology '{self.topology}' is not implemented. Supported types are E1 through E6.") + + #custom power spectrum function ( power law with one wavepacket) + def PK(k, L, amp, freq, xc): + x=k*L*L_LSS/(2*np.pi) + return A_s*(k/0.05)**(n_s-1)*(1+ np.sin(x*freq)*amp*np.exp(-(x/xc))) + + #exponential cutoff + def PK_exp(k, L, xc, amp): + x=k*L*L_LSS/(2*np.pi) + return A_s*(k/0.05)**(n_s-1)*(1-amp*np.exp(-(x/xc))) + + def PK_rising(k, L, xc, amp): + x=k*L*L_LSS/(2*np.pi) + return A_s*(k/0.05)**(n_s-1)*(1+amp*np.exp(-(x/xc))) + + #standard power law + def PK_pow(k): + return A_s*(k/0.05)**(n_s-1) + + # Obtain the C_ls (power spectrum) again but for a modified initial power spectrum + pars_mod = copy.deepcopy(pars) + if self.powerspec == 'powlaw': + pars_mod.set_initial_power_function(PK_pow) + elif self.powerspec == 'wavepacket': + pars_mod.set_initial_power_function(PK, args=(L_PS, self.amp, self.freq, self.width)) + elif self.powerspec == 'cutoff': + pars_mod.set_initial_power_function(PK_exp, args=(L_PS, self.x_cutoff, self.amp)) + elif self.powerspec == 'enhance': + pars_mod.set_initial_power_function(PK_rising, args=(L_PS, self.x_cutoff, self.amp)) + + results_mod = camb.get_results(pars_mod) + self.powers_mod = results_mod.get_cmb_power_spectra(pars_mod, raw_cl=True, CMB_unit='muK')['unlensed_scalar'] + # To get C_\ell in units of umK, we multiply by 1e6 (K to micro K) and the temperature of the CMB in K + # The transfer functions are not affected by a modified primordial power spectrum transfer_data = np.array(transfer_function.delta_p_l_k) * 1e6 * 2.7255 # Convert to muK print(f'Shape of transfer function from CAMB: {transfer_data.shape}') @@ -150,6 +205,7 @@ def do_pre_processing(self): self.transfer_E_interpolate_k_l_list[l] = scipy.interpolate.interp1d(self.k_list, transfer_data[1, l-2, :], kind='cubic') # Compute wavevector cutoffs + # keep computing them with the standard PS to enable easier comparison, to do: check differences self.get_kmax_as_function_of_ell(pars.scalar_power) # We find all allowed |k|, phi, theta and put them in big lists @@ -167,7 +223,11 @@ def do_pre_processing(self): # Get P(k) / k^3 for all unique |k| values - self.scalar_pk_k3 = pars.scalar_power(self.k_amp_unique) / self.k_amp_unique**3 + # Using the modified PS for the non-trivial topology if specified + if self.PS_mod: + self.scalar_pk_k3 = pars_mod.scalar_power(self.k_amp_unique) / self.k_amp_unique**3 + else: + self.scalar_pk_k3 = pars.scalar_power(self.k_amp_unique) / self.k_amp_unique**3 # Get the transfer function for all unique |k| values self.transfer_T_delta_kl = self.get_transfer_functions_multi(self.transfer_T_interpolate_k_l_list) @@ -552,4 +612,27 @@ def get_c_lmlpmp_multiprocessing(self, ell_range, ell_p_range): lm_index = ell * (ell+1) + m - l_min * l_min c_lmlpmp[lm_p_index, lm_index] = np.conjugate(c_lmlpmp[lm_index, lm_p_index]) - return c_lmlpmp \ No newline at end of file + return c_lmlpmp + + def calculate_exact_kl_divergence(self): #normalize cov matrix with the one from E18 to get matrix for KL div calculation + print('Calculating KL divergence') + + c_lmlpmp_ordered = self.calculate_c_lmlpmp(only_diag=False, normalize=False, plotting = False) + #normalize with self.powers_mod to have E18 also with a mod. initial power spec# + if self.E18_mod: + A_ssp = normalize_c_lmlpmp(c_lmlpmp_ordered, self.powers_mod[:, 0], cl_accuracy = self.c_l_accuracy, l_min=self.l_min, lp_min=self.l_min, l_max=self.l_max, lp_max=self.l_max) + else: + A_ssp = normalize_c_lmlpmp(c_lmlpmp_ordered, self.powers[:, 0], cl_accuracy = self.c_l_accuracy, l_min=self.l_min, lp_min=self.l_min, l_max=self.l_max, lp_max=self.l_max) + + w, _ = np.linalg.eig(A_ssp) + kl_P_assuming_Q = 0 + kl_Q_assuming_P = 0 + for eig in w: + kl_P_assuming_Q += (np.log(np.abs(eig)) + 1/eig - 1)/2 + kl_Q_assuming_P += (-np.log(np.abs(eig)) + eig - 1)/2 + + np.fill_diagonal(A_ssp, 0) + + a_t = np.sqrt(np.sum(np.abs(A_ssp)**2)) + + return np.real(kl_P_assuming_Q), np.real(kl_Q_assuming_P), np.real(a_t) \ No newline at end of file