Mod prim ps

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/. " \

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()

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
+
+'''

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,16 +67,50 @@ 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")
 
     # Validate topology-specific parameters
     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
     )
 

topology/src/config.py

@@ -3,4 +3,8 @@
 """
 
 # Physical constants
-L_LSS = 13824.9 * 2  # Last scattering surface diameter (Mpc)
+L_LSS = 13824.9 * 2  # Last scattering surface diameter (Mpc)
+
+# Primordial power spectrum
+A_s = 2e-9      # Amplitude
+n_s = 0.965     # Spectral index