geff.gef
This module defines the BaseGEF class and its class factory compile_model.
1import numpy as np 2import inspect 3from ._docs import generate_docs, docs_gef 4from .bgtypes import BGSystem, Val, Func 5from .models import load_model, pai 6 7from numbers import Number 8from typing import Callable 9 10class BaseGEF: 11 12 GEFSolver = pai.solver 13 """The solver used to solve the GEF equations in `run`.""" 14 ModeSolver = pai.MbM 15 """The mode solver used for mode-by-mode cross checks.""" 16 _input_signature = pai.model_input 17 define_units = staticmethod(pai.define_units) 18 19 def __init__(self, **kwargs): 20 """ 21 Initialize the GEF from user-specified initial data, constants and functions. 22 23 Parameters 24 ---------- 25 kwargs 26 To access a list of kwargs for this model use `print_input`. 27 28 Raises 29 ------ 30 KeyError 31 if a necessary input is missing. 32 TypeError 33 if the input is of the wrong type. 34 """ 35 36 for val in self._input_signature: 37 if val.name not in kwargs.keys(): 38 raise KeyError(f"Missing input '{val.name}'") 39 else: 40 key = val.name 41 item = kwargs[key] 42 if issubclass(val, Val) and not(isinstance(item, Number)): 43 raise TypeError(f"Input for '{key}' should be 'Number' type.") 44 elif issubclass(val, Func) and not(callable(item)): 45 raise TypeError(f"Input for '{key}' must be callable.") 46 47 #check which keys are needed for define_units 48 expected_keys = inspect.getfullargspec(self.define_units).args 49 omega, mu = self.define_units(**{key:kwargs[key] for key in expected_keys}) 50 51 known_objects = set().union(*[item for key, item in self.GEFSolver.known_variables.items() if key!="gauge"] ) 52 53 initial_data = BGSystem(known_objects, omega, mu) 54 55 for key, item in kwargs.items(): 56 initial_data.initialise(key)(item) 57 58 #initialise the other values with dummy variables. 59 for obj in initial_data.quantity_set(): 60 if not(hasattr(initial_data, obj.name)): 61 if issubclass(obj, Val): 62 initial_data.initialise(obj.name)(0) 63 elif issubclass(obj, Func): 64 initial_data.initialise(obj.name)(lambda *x: 0) 65 66 self.initial_data = initial_data 67 68 # test if ODE solver works 69 try: 70 self.GEFSolver(initial_data) 71 except Exception: 72 raise ModelError("Initializing the GEFSolver returned an error.") 73 74 75 @classmethod 76 def print_input(cls): 77 """ 78 Print the input required to initialize the class. 79 """ 80 print("This GEF model expects the following input:\n") 81 for i in cls._input_signature: 82 print(f" * {i.get_description()}") 83 print() 84 return 85 86 @classmethod 87 def pring_ingredients(cls): 88 """ 89 Print a list of known variables, functions and constants for this model 90 """ 91 print("This GEF model knows the following variables:\n") 92 for key, item in cls.GEFSolver.known_variables.items(): 93 if len(item) > 0: 94 print(f"{key}:") 95 for i in item: 96 print(f" * {i.get_description()}") 97 print() 98 return 99 100 def run(self, ntr:int=150, tend:float=120, nmodes:int=500, mbm_attempts:int=5, resume_mbm:bool=True, 101 err_tol:float=0.1, err_thr:float=0.025, binning:int=5, integrator:str="simpson", print_stats:bool=True, **solver_kwargs): 102 """ 103 Solve the ODE's of the GEF using `GEFSolver`. Cross check the solution using `ModeSolver`. 104 105 The `GEFSolver` is initialized using the initial conditions defined by the class. 106 If the solver returns a successful solution, `ModeSolver.compute_spectrum` computes a gauge field spectrum 107 to perform a mode-by-mode cross check (unless `nmodes=None`). 108 If the mode-by-mode cross check is a success, the solution is stored in the underlying `BGSystem` of the class. 109 Otherwise, the `GEFSolver` tries to self correct using the gauge field spectrum. This is attempted for `mbm_attempts` or until successful. 110 111 Parameters 112 ---------- 113 ntr : int 114 initial truncation number `GEFSolver.ntr` 115 tend : float 116 initial target time for `GEFSolver.tend` 117 nmodes : int or None 118 the number of modes computed by `ModeSolver`. If None, no cross-check is performed. 119 mbm_attempts : int 120 number of mode-by-mode self correction attempts 121 resume_mbm : bool 122 if `True` use `ModeSolver.update_spectrum` in case successive mode-by-mode comparisons are needed. 123 err_tol : float 124 passed to `mbm_crosscheck`. 125 err_thr : float 126 passed to `mbm_crosscheck`. 127 binning : int 128 passed to `mbm_crosscheck`. 129 integrator : str 130 integrator for `mbm_crosscheck` ('simpson' is advised) 131 print_stats : bool 132 if `True`, a summary report is printed for the returned solution. 133 solver_kwargs 134 the `settings` of `GEFSolver` (see `GEFSolver.settings`) 135 136 Returns 137 ------- 138 vals : BGSystem 139 contains the background evolution obtained from `GEFSolver` 140 spec : GaugeSpec or None 141 the result of `ModeSolver.compute_spectrum` 142 sol 143 the full result of `GEFSolver.compute_GEF_solution` 144 145 146 Raises 147 ------ 148 RuntimeError 149 if no successful solution was obtained. 150 """ 151 solver = self.GEFSolver(self.initial_data) 152 153 #Configuring GEFSolver 154 solver.ntr=ntr 155 solver.tend=tend 156 solver.update_settings(**solver_kwargs) 157 158 integrator_kwargs = {"integrator":integrator, "epsabs":solver.settings["atol"], "epsrel":solver.settings["rtol"]} 159 160 done=False 161 vals = BGSystem.from_system(self.initial_data, copy=True) 162 vals.units = False 163 sol = None 164 spec = None 165 t_reinit = 0. 166 attempt=0 167 168 print() 169 170 while not(done) and attempt<mbm_attempts: 171 attempt +=1 172 #This can be taken care of internally. The GEF should not need to get sol objects... 173 sol_new = solver.compute_GEF_solution() 174 sol = self._update_sol(sol, sol_new) 175 solver.parse_arr_to_sys(sol.t, sol.y, vals) 176 177 if nmodes is not None: 178 print("Using last successful GEF solution to compute gauge-field mode functions.") 179 MbM = self.ModeSolver(vals) 180 181 rtol = solver.settings["rtol"] 182 183 if resume_mbm and attempt > 1: 184 spec = MbM.update_spectrum(spec, t_reinit, rtol=rtol) 185 else: 186 spec = MbM.compute_spectrum(nmodes, rtol=rtol) 187 print("Performing mode-by-mode comparison with GEF results.") 188 189 agreement, reinit_spec = self.mbm_crosscheck(spec, vals, err_tol= err_tol, err_thr=err_thr, binning=binning, 190 **integrator_kwargs) 191 192 if agreement: 193 print("The mode-by-mode comparison indicates a convergent GEF run.\n") 194 done=True 195 196 else: 197 t_reinit = reinit_spec["t"] 198 print(f"Attempting to solve GEF using self-correction starting from t={reinit_spec['t']:.1f}, N={reinit_spec['N']:.1f}.\n") 199 200 solver.set_initial_conditions_to_MbM(sol, reinit_spec) 201 202 else: 203 done=True 204 205 if done: 206 if print_stats: 207 self._print_summary(sol) 208 if sol.success: 209 vals.units = True 210 return vals, spec, sol 211 else: 212 print("The run terminated on with an error, check output for details.") 213 else: 214 raise RuntimeError(f"GEF did not complete after {attempt} attempts.") 215 216 @staticmethod 217 def _print_summary(sol): 218 print("GEF run completed with the following statistics") 219 for attr in sol.keys(): 220 if attr not in ["y", "t", "y_events", "t_events", "sol", "events"]: 221 print(rf" - {attr} : {getattr(sol, attr)}") 222 events = sol.events 223 if np.array([(len(event)==0) for event in events.values()]).all(): 224 print("No events occurred during the final run.") 225 else: 226 print("The following events occurred during the final run:") 227 for event, time in events.items(): 228 if len(time) > 0: 229 if len(time) > 1: 230 tstr = [f"{t:.1f}" for t in time] 231 else: 232 tstr = f"{time[0]:.1f}" 233 print(f" - {event} at t={tstr}") 234 return 235 236 @staticmethod 237 def _update_sol(sol_old, sol_new): 238 """ 239 Update an old GEF solution with a new one, overwriting the overlap. 240 """ 241 if sol_old is None: 242 return sol_new 243 else: 244 sol = sol_old 245 ind_overlap = np.searchsorted(sol_old.t, sol_new.t[0], "left") 246 sol.t = np.concatenate([sol_old.t[:ind_overlap], sol_new.t]) 247 248 if sol_old.y.shape[0] < sol_new.y.shape[0]: 249 #if ntr increased from one solution to the next, fill up sol_old with zeros to match sol_new 250 fillshape = (sol_new.y.shape[0] - sol_old.y.shape[0], sol_old.y.shape[1]) 251 yfill = np.zeros( fillshape ) 252 sol_old.y = np.concatenate([sol_old.y, yfill], axis=0) 253 254 sol.y = np.concatenate([sol_old.y[:,:ind_overlap], sol_new.y], axis=1) 255 sol.events.update(sol_new.events) 256 for attr in ["nfev", "njev", "nlu"]: 257 setattr(sol, attr, getattr(sol_old, attr) + getattr(sol_new, attr)) 258 for attr in ["message", "success", "status"]: 259 setattr(sol, attr, getattr(sol_new, attr)) 260 return sol 261 262 def mbm_crosscheck(self, spec, vals:BGSystem, err_tol:float, err_thr:float, binning:int, **integrator_kwargs): 263 """ 264 Estimate the error of a GEF solution using `.mbm.GaugeSpec.estimate_GEF_error`. 265 266 If either the RMS error or the final error exceeds `err_tol`, the solution is rejected. 267 268 Parameters 269 ---------- 270 vals : BGSystem 271 contains the GEF solution. 272 err_tol : float 273 the tolerance on the RMS and final error. 274 err_thr : float 275 passed to `estimate_GEF_error`. 276 binning : int 277 passed to `estimate_GEF_error`. 278 integratorkwargs: 279 passed to kwargs of `estimate_GEF_error`. 280 281 Returns 282 ------- 283 agreement : bool 284 indicates if the solution is accepted or rejected. 285 reinit_slice : SpecSlice 286 the spectrum with which the GEF solver is re-initialized. 287 """ 288 GFs = self.GEFSolver.known_variables["gauge"] 289 agreement=True 290 291 for GF in GFs: 292 ref = [a.name for a in GF.zeros] 293 cut = GF.cutoff.name 294 295 errs, terr, _ = spec.estimate_GEF_error(vals, references=ref, cutoff=cut, 296 err_thr=err_thr, binning=binning, **integrator_kwargs) 297 298 reinit_inds = [] 299 300 for err in errs: 301 rmserr = np.sqrt(np.sum(err**2)/len(err)) 302 if max(err[-1], rmserr) > err_tol: 303 agreement=False 304 #find where the error is above 5%, take the earliest occurrence, reduce by 1 305 err_ind = np.where(err > err_thr)[0][0]-1 306 else: 307 err_ind = len(terr)-1 308 reinit_inds.append( err_ind ) 309 310 t0 = terr[min(reinit_inds)] 311 312 ind = np.searchsorted(spec["t"],t0, "left") 313 314 reinit_slice = spec.tslice(ind) 315 316 return agreement, reinit_slice 317 318 def load_GEFdata(self, path : str): 319 """ 320 Load data and return a BGSystem with the results. 321 322 Note, data is always loaded assuming numerical units. 323 324 Parameters 325 ---------- 326 path : None or str 327 Path to load data from 328 329 Raises 330 ------ 331 FileNotFoundError 332 if no file is found at `path`. 333 AttributeError 334 if the file contains a column labeled by a key which does not match any GEF-value name. 335 """ 336 newsys = BGSystem.from_system(self.initial_data, copy=True) 337 338 newsys.load_variables(path) 339 340 return newsys 341 342def compile_model(modelname:str, settings:dict={}): 343 """ 344 Define a custom subclass of BaseGEF adapted to a new GEF model. 345 346 Parameters 347 ---------- 348 modelname : str 349 The name of the GEF model or a full dotted import path (e.g., "path.to.module"). 350 settings : dict 351 a dictionary of settings used by the model 352 353 Returns 354 ------- 355 CustomGEF 356 a custom subclass of BaseGEF 357 358 """ 359 model = load_model(modelname, settings) 360 signature = [inspect.Parameter("self", inspect.Parameter.POSITIONAL_ONLY)] 361 for val in model.model_input: 362 if issubclass(val, Func): 363 annotation = Callable 364 else: 365 annotation = float 366 signature.append(inspect.Parameter(val.name, inspect.Parameter.KEYWORD_ONLY, annotation=annotation)) 367 368 class CustomGEF(BaseGEF): 369 GEFSolver = model.solver 370 """The solver used to solve the GEF equations in `run`.""" 371 ModeSolver = model.MbM 372 """The mode solver used for mode-by-mode cross checks.""" 373 _input_signature = model.model_input 374 define_units = staticmethod(model.define_units) 375 376 CustomGEF.__init__.__signature__ = inspect.Signature(signature) 377 378 CustomGEF.__qualname__ = model.name 379 CustomGEF.__module__ = __name__ 380 381 return CustomGEF 382 383class ModelError(Exception): 384 """ 385 Exception for errors when compiling a GEF model. 386 """ 387 pass 388 389generate_docs(docs_gef.DOCS)
class
BaseGEF:
11class BaseGEF: 12 13 GEFSolver = pai.solver 14 """The solver used to solve the GEF equations in `run`.""" 15 ModeSolver = pai.MbM 16 """The mode solver used for mode-by-mode cross checks.""" 17 _input_signature = pai.model_input 18 define_units = staticmethod(pai.define_units) 19 20 def __init__(self, **kwargs): 21 """ 22 Initialize the GEF from user-specified initial data, constants and functions. 23 24 Parameters 25 ---------- 26 kwargs 27 To access a list of kwargs for this model use `print_input`. 28 29 Raises 30 ------ 31 KeyError 32 if a necessary input is missing. 33 TypeError 34 if the input is of the wrong type. 35 """ 36 37 for val in self._input_signature: 38 if val.name not in kwargs.keys(): 39 raise KeyError(f"Missing input '{val.name}'") 40 else: 41 key = val.name 42 item = kwargs[key] 43 if issubclass(val, Val) and not(isinstance(item, Number)): 44 raise TypeError(f"Input for '{key}' should be 'Number' type.") 45 elif issubclass(val, Func) and not(callable(item)): 46 raise TypeError(f"Input for '{key}' must be callable.") 47 48 #check which keys are needed for define_units 49 expected_keys = inspect.getfullargspec(self.define_units).args 50 omega, mu = self.define_units(**{key:kwargs[key] for key in expected_keys}) 51 52 known_objects = set().union(*[item for key, item in self.GEFSolver.known_variables.items() if key!="gauge"] ) 53 54 initial_data = BGSystem(known_objects, omega, mu) 55 56 for key, item in kwargs.items(): 57 initial_data.initialise(key)(item) 58 59 #initialise the other values with dummy variables. 60 for obj in initial_data.quantity_set(): 61 if not(hasattr(initial_data, obj.name)): 62 if issubclass(obj, Val): 63 initial_data.initialise(obj.name)(0) 64 elif issubclass(obj, Func): 65 initial_data.initialise(obj.name)(lambda *x: 0) 66 67 self.initial_data = initial_data 68 69 # test if ODE solver works 70 try: 71 self.GEFSolver(initial_data) 72 except Exception: 73 raise ModelError("Initializing the GEFSolver returned an error.") 74 75 76 @classmethod 77 def print_input(cls): 78 """ 79 Print the input required to initialize the class. 80 """ 81 print("This GEF model expects the following input:\n") 82 for i in cls._input_signature: 83 print(f" * {i.get_description()}") 84 print() 85 return 86 87 @classmethod 88 def pring_ingredients(cls): 89 """ 90 Print a list of known variables, functions and constants for this model 91 """ 92 print("This GEF model knows the following variables:\n") 93 for key, item in cls.GEFSolver.known_variables.items(): 94 if len(item) > 0: 95 print(f"{key}:") 96 for i in item: 97 print(f" * {i.get_description()}") 98 print() 99 return 100 101 def run(self, ntr:int=150, tend:float=120, nmodes:int=500, mbm_attempts:int=5, resume_mbm:bool=True, 102 err_tol:float=0.1, err_thr:float=0.025, binning:int=5, integrator:str="simpson", print_stats:bool=True, **solver_kwargs): 103 """ 104 Solve the ODE's of the GEF using `GEFSolver`. Cross check the solution using `ModeSolver`. 105 106 The `GEFSolver` is initialized using the initial conditions defined by the class. 107 If the solver returns a successful solution, `ModeSolver.compute_spectrum` computes a gauge field spectrum 108 to perform a mode-by-mode cross check (unless `nmodes=None`). 109 If the mode-by-mode cross check is a success, the solution is stored in the underlying `BGSystem` of the class. 110 Otherwise, the `GEFSolver` tries to self correct using the gauge field spectrum. This is attempted for `mbm_attempts` or until successful. 111 112 Parameters 113 ---------- 114 ntr : int 115 initial truncation number `GEFSolver.ntr` 116 tend : float 117 initial target time for `GEFSolver.tend` 118 nmodes : int or None 119 the number of modes computed by `ModeSolver`. If None, no cross-check is performed. 120 mbm_attempts : int 121 number of mode-by-mode self correction attempts 122 resume_mbm : bool 123 if `True` use `ModeSolver.update_spectrum` in case successive mode-by-mode comparisons are needed. 124 err_tol : float 125 passed to `mbm_crosscheck`. 126 err_thr : float 127 passed to `mbm_crosscheck`. 128 binning : int 129 passed to `mbm_crosscheck`. 130 integrator : str 131 integrator for `mbm_crosscheck` ('simpson' is advised) 132 print_stats : bool 133 if `True`, a summary report is printed for the returned solution. 134 solver_kwargs 135 the `settings` of `GEFSolver` (see `GEFSolver.settings`) 136 137 Returns 138 ------- 139 vals : BGSystem 140 contains the background evolution obtained from `GEFSolver` 141 spec : GaugeSpec or None 142 the result of `ModeSolver.compute_spectrum` 143 sol 144 the full result of `GEFSolver.compute_GEF_solution` 145 146 147 Raises 148 ------ 149 RuntimeError 150 if no successful solution was obtained. 151 """ 152 solver = self.GEFSolver(self.initial_data) 153 154 #Configuring GEFSolver 155 solver.ntr=ntr 156 solver.tend=tend 157 solver.update_settings(**solver_kwargs) 158 159 integrator_kwargs = {"integrator":integrator, "epsabs":solver.settings["atol"], "epsrel":solver.settings["rtol"]} 160 161 done=False 162 vals = BGSystem.from_system(self.initial_data, copy=True) 163 vals.units = False 164 sol = None 165 spec = None 166 t_reinit = 0. 167 attempt=0 168 169 print() 170 171 while not(done) and attempt<mbm_attempts: 172 attempt +=1 173 #This can be taken care of internally. The GEF should not need to get sol objects... 174 sol_new = solver.compute_GEF_solution() 175 sol = self._update_sol(sol, sol_new) 176 solver.parse_arr_to_sys(sol.t, sol.y, vals) 177 178 if nmodes is not None: 179 print("Using last successful GEF solution to compute gauge-field mode functions.") 180 MbM = self.ModeSolver(vals) 181 182 rtol = solver.settings["rtol"] 183 184 if resume_mbm and attempt > 1: 185 spec = MbM.update_spectrum(spec, t_reinit, rtol=rtol) 186 else: 187 spec = MbM.compute_spectrum(nmodes, rtol=rtol) 188 print("Performing mode-by-mode comparison with GEF results.") 189 190 agreement, reinit_spec = self.mbm_crosscheck(spec, vals, err_tol= err_tol, err_thr=err_thr, binning=binning, 191 **integrator_kwargs) 192 193 if agreement: 194 print("The mode-by-mode comparison indicates a convergent GEF run.\n") 195 done=True 196 197 else: 198 t_reinit = reinit_spec["t"] 199 print(f"Attempting to solve GEF using self-correction starting from t={reinit_spec['t']:.1f}, N={reinit_spec['N']:.1f}.\n") 200 201 solver.set_initial_conditions_to_MbM(sol, reinit_spec) 202 203 else: 204 done=True 205 206 if done: 207 if print_stats: 208 self._print_summary(sol) 209 if sol.success: 210 vals.units = True 211 return vals, spec, sol 212 else: 213 print("The run terminated on with an error, check output for details.") 214 else: 215 raise RuntimeError(f"GEF did not complete after {attempt} attempts.") 216 217 @staticmethod 218 def _print_summary(sol): 219 print("GEF run completed with the following statistics") 220 for attr in sol.keys(): 221 if attr not in ["y", "t", "y_events", "t_events", "sol", "events"]: 222 print(rf" - {attr} : {getattr(sol, attr)}") 223 events = sol.events 224 if np.array([(len(event)==0) for event in events.values()]).all(): 225 print("No events occurred during the final run.") 226 else: 227 print("The following events occurred during the final run:") 228 for event, time in events.items(): 229 if len(time) > 0: 230 if len(time) > 1: 231 tstr = [f"{t:.1f}" for t in time] 232 else: 233 tstr = f"{time[0]:.1f}" 234 print(f" - {event} at t={tstr}") 235 return 236 237 @staticmethod 238 def _update_sol(sol_old, sol_new): 239 """ 240 Update an old GEF solution with a new one, overwriting the overlap. 241 """ 242 if sol_old is None: 243 return sol_new 244 else: 245 sol = sol_old 246 ind_overlap = np.searchsorted(sol_old.t, sol_new.t[0], "left") 247 sol.t = np.concatenate([sol_old.t[:ind_overlap], sol_new.t]) 248 249 if sol_old.y.shape[0] < sol_new.y.shape[0]: 250 #if ntr increased from one solution to the next, fill up sol_old with zeros to match sol_new 251 fillshape = (sol_new.y.shape[0] - sol_old.y.shape[0], sol_old.y.shape[1]) 252 yfill = np.zeros( fillshape ) 253 sol_old.y = np.concatenate([sol_old.y, yfill], axis=0) 254 255 sol.y = np.concatenate([sol_old.y[:,:ind_overlap], sol_new.y], axis=1) 256 sol.events.update(sol_new.events) 257 for attr in ["nfev", "njev", "nlu"]: 258 setattr(sol, attr, getattr(sol_old, attr) + getattr(sol_new, attr)) 259 for attr in ["message", "success", "status"]: 260 setattr(sol, attr, getattr(sol_new, attr)) 261 return sol 262 263 def mbm_crosscheck(self, spec, vals:BGSystem, err_tol:float, err_thr:float, binning:int, **integrator_kwargs): 264 """ 265 Estimate the error of a GEF solution using `.mbm.GaugeSpec.estimate_GEF_error`. 266 267 If either the RMS error or the final error exceeds `err_tol`, the solution is rejected. 268 269 Parameters 270 ---------- 271 vals : BGSystem 272 contains the GEF solution. 273 err_tol : float 274 the tolerance on the RMS and final error. 275 err_thr : float 276 passed to `estimate_GEF_error`. 277 binning : int 278 passed to `estimate_GEF_error`. 279 integratorkwargs: 280 passed to kwargs of `estimate_GEF_error`. 281 282 Returns 283 ------- 284 agreement : bool 285 indicates if the solution is accepted or rejected. 286 reinit_slice : SpecSlice 287 the spectrum with which the GEF solver is re-initialized. 288 """ 289 GFs = self.GEFSolver.known_variables["gauge"] 290 agreement=True 291 292 for GF in GFs: 293 ref = [a.name for a in GF.zeros] 294 cut = GF.cutoff.name 295 296 errs, terr, _ = spec.estimate_GEF_error(vals, references=ref, cutoff=cut, 297 err_thr=err_thr, binning=binning, **integrator_kwargs) 298 299 reinit_inds = [] 300 301 for err in errs: 302 rmserr = np.sqrt(np.sum(err**2)/len(err)) 303 if max(err[-1], rmserr) > err_tol: 304 agreement=False 305 #find where the error is above 5%, take the earliest occurrence, reduce by 1 306 err_ind = np.where(err > err_thr)[0][0]-1 307 else: 308 err_ind = len(terr)-1 309 reinit_inds.append( err_ind ) 310 311 t0 = terr[min(reinit_inds)] 312 313 ind = np.searchsorted(spec["t"],t0, "left") 314 315 reinit_slice = spec.tslice(ind) 316 317 return agreement, reinit_slice 318 319 def load_GEFdata(self, path : str): 320 """ 321 Load data and return a BGSystem with the results. 322 323 Note, data is always loaded assuming numerical units. 324 325 Parameters 326 ---------- 327 path : None or str 328 Path to load data from 329 330 Raises 331 ------ 332 FileNotFoundError 333 if no file is found at `path`. 334 AttributeError 335 if the file contains a column labeled by a key which does not match any GEF-value name. 336 """ 337 newsys = BGSystem.from_system(self.initial_data, copy=True) 338 339 newsys.load_variables(path) 340 341 return newsys
This class defines the basic properties of a GEF model.
Each GEF model contains a GEFSolver and a ModeSolver that are used to solve GEF equations with the run method.
Instances of a GEF model can also be used to load data using load_GEFdata.
Loading data this way creates a BGSystem which also knows the appropriate units, constants and functions for this model.
The BaseGEF is a compiled version of geff.models.pai. To compile other models, use the class factory compile_model.
BaseGEF(**kwargs)
20 def __init__(self, **kwargs): 21 """ 22 Initialize the GEF from user-specified initial data, constants and functions. 23 24 Parameters 25 ---------- 26 kwargs 27 To access a list of kwargs for this model use `print_input`. 28 29 Raises 30 ------ 31 KeyError 32 if a necessary input is missing. 33 TypeError 34 if the input is of the wrong type. 35 """ 36 37 for val in self._input_signature: 38 if val.name not in kwargs.keys(): 39 raise KeyError(f"Missing input '{val.name}'") 40 else: 41 key = val.name 42 item = kwargs[key] 43 if issubclass(val, Val) and not(isinstance(item, Number)): 44 raise TypeError(f"Input for '{key}' should be 'Number' type.") 45 elif issubclass(val, Func) and not(callable(item)): 46 raise TypeError(f"Input for '{key}' must be callable.") 47 48 #check which keys are needed for define_units 49 expected_keys = inspect.getfullargspec(self.define_units).args 50 omega, mu = self.define_units(**{key:kwargs[key] for key in expected_keys}) 51 52 known_objects = set().union(*[item for key, item in self.GEFSolver.known_variables.items() if key!="gauge"] ) 53 54 initial_data = BGSystem(known_objects, omega, mu) 55 56 for key, item in kwargs.items(): 57 initial_data.initialise(key)(item) 58 59 #initialise the other values with dummy variables. 60 for obj in initial_data.quantity_set(): 61 if not(hasattr(initial_data, obj.name)): 62 if issubclass(obj, Val): 63 initial_data.initialise(obj.name)(0) 64 elif issubclass(obj, Func): 65 initial_data.initialise(obj.name)(lambda *x: 0) 66 67 self.initial_data = initial_data 68 69 # test if ODE solver works 70 try: 71 self.GEFSolver(initial_data) 72 except Exception: 73 raise ModelError("Initializing the GEFSolver returned an error.")
Initialize the GEF from user-specified initial data, constants and functions.
Parameters
- kwargs: To access a list of kwargs for this model use
print_input.
Raises
- KeyError: if a necessary input is missing.
- TypeError: if the input is of the wrong type.
def
define_units(phi, dphi, V):
68def define_units(phi, dphi, V): 69 #compute Hubble rate at t0 70 rhoK = 0.5*dphi**2 71 rhoV = V(phi) 72 H0 = friedmann( rhoK, rhoV ) 73 74 omega = H0 #Characteristic frequency is the initial Hubble rate in Planck units 75 mu = 1. #Charatcterisic amplitude is the Planck mass (in Planck units) 76 77 return omega, mu
Define how initial data is used to set the reference frequency, $\omega$ and energy scale, $\mu$.
- energy scale: $M_{\rm P}$ in Planck units
- frequency scale: initial Hubble rate, $H_0$ (in Planck units)
@classmethod
def
print_input(cls):
76 @classmethod 77 def print_input(cls): 78 """ 79 Print the input required to initialize the class. 80 """ 81 print("This GEF model expects the following input:\n") 82 for i in cls._input_signature: 83 print(f" * {i.get_description()}") 84 print() 85 return
Print the input required to initialize the class.
@classmethod
def
pring_ingredients(cls):
87 @classmethod 88 def pring_ingredients(cls): 89 """ 90 Print a list of known variables, functions and constants for this model 91 """ 92 print("This GEF model knows the following variables:\n") 93 for key, item in cls.GEFSolver.known_variables.items(): 94 if len(item) > 0: 95 print(f"{key}:") 96 for i in item: 97 print(f" * {i.get_description()}") 98 print() 99 return
Print a list of known variables, functions and constants for this model
def
run( self, ntr: int = 150, tend: float = 120, nmodes: int = 500, mbm_attempts: int = 5, resume_mbm: bool = True, err_tol: float = 0.1, err_thr: float = 0.025, binning: int = 5, integrator: str = 'simpson', print_stats: bool = True, **solver_kwargs):
101 def run(self, ntr:int=150, tend:float=120, nmodes:int=500, mbm_attempts:int=5, resume_mbm:bool=True, 102 err_tol:float=0.1, err_thr:float=0.025, binning:int=5, integrator:str="simpson", print_stats:bool=True, **solver_kwargs): 103 """ 104 Solve the ODE's of the GEF using `GEFSolver`. Cross check the solution using `ModeSolver`. 105 106 The `GEFSolver` is initialized using the initial conditions defined by the class. 107 If the solver returns a successful solution, `ModeSolver.compute_spectrum` computes a gauge field spectrum 108 to perform a mode-by-mode cross check (unless `nmodes=None`). 109 If the mode-by-mode cross check is a success, the solution is stored in the underlying `BGSystem` of the class. 110 Otherwise, the `GEFSolver` tries to self correct using the gauge field spectrum. This is attempted for `mbm_attempts` or until successful. 111 112 Parameters 113 ---------- 114 ntr : int 115 initial truncation number `GEFSolver.ntr` 116 tend : float 117 initial target time for `GEFSolver.tend` 118 nmodes : int or None 119 the number of modes computed by `ModeSolver`. If None, no cross-check is performed. 120 mbm_attempts : int 121 number of mode-by-mode self correction attempts 122 resume_mbm : bool 123 if `True` use `ModeSolver.update_spectrum` in case successive mode-by-mode comparisons are needed. 124 err_tol : float 125 passed to `mbm_crosscheck`. 126 err_thr : float 127 passed to `mbm_crosscheck`. 128 binning : int 129 passed to `mbm_crosscheck`. 130 integrator : str 131 integrator for `mbm_crosscheck` ('simpson' is advised) 132 print_stats : bool 133 if `True`, a summary report is printed for the returned solution. 134 solver_kwargs 135 the `settings` of `GEFSolver` (see `GEFSolver.settings`) 136 137 Returns 138 ------- 139 vals : BGSystem 140 contains the background evolution obtained from `GEFSolver` 141 spec : GaugeSpec or None 142 the result of `ModeSolver.compute_spectrum` 143 sol 144 the full result of `GEFSolver.compute_GEF_solution` 145 146 147 Raises 148 ------ 149 RuntimeError 150 if no successful solution was obtained. 151 """ 152 solver = self.GEFSolver(self.initial_data) 153 154 #Configuring GEFSolver 155 solver.ntr=ntr 156 solver.tend=tend 157 solver.update_settings(**solver_kwargs) 158 159 integrator_kwargs = {"integrator":integrator, "epsabs":solver.settings["atol"], "epsrel":solver.settings["rtol"]} 160 161 done=False 162 vals = BGSystem.from_system(self.initial_data, copy=True) 163 vals.units = False 164 sol = None 165 spec = None 166 t_reinit = 0. 167 attempt=0 168 169 print() 170 171 while not(done) and attempt<mbm_attempts: 172 attempt +=1 173 #This can be taken care of internally. The GEF should not need to get sol objects... 174 sol_new = solver.compute_GEF_solution() 175 sol = self._update_sol(sol, sol_new) 176 solver.parse_arr_to_sys(sol.t, sol.y, vals) 177 178 if nmodes is not None: 179 print("Using last successful GEF solution to compute gauge-field mode functions.") 180 MbM = self.ModeSolver(vals) 181 182 rtol = solver.settings["rtol"] 183 184 if resume_mbm and attempt > 1: 185 spec = MbM.update_spectrum(spec, t_reinit, rtol=rtol) 186 else: 187 spec = MbM.compute_spectrum(nmodes, rtol=rtol) 188 print("Performing mode-by-mode comparison with GEF results.") 189 190 agreement, reinit_spec = self.mbm_crosscheck(spec, vals, err_tol= err_tol, err_thr=err_thr, binning=binning, 191 **integrator_kwargs) 192 193 if agreement: 194 print("The mode-by-mode comparison indicates a convergent GEF run.\n") 195 done=True 196 197 else: 198 t_reinit = reinit_spec["t"] 199 print(f"Attempting to solve GEF using self-correction starting from t={reinit_spec['t']:.1f}, N={reinit_spec['N']:.1f}.\n") 200 201 solver.set_initial_conditions_to_MbM(sol, reinit_spec) 202 203 else: 204 done=True 205 206 if done: 207 if print_stats: 208 self._print_summary(sol) 209 if sol.success: 210 vals.units = True 211 return vals, spec, sol 212 else: 213 print("The run terminated on with an error, check output for details.") 214 else: 215 raise RuntimeError(f"GEF did not complete after {attempt} attempts.")
Solve the ODE's of the GEF using GEFSolver. Cross check the solution using ModeSolver.
The GEFSolver is initialized using the initial conditions defined by the class.
If the solver returns a successful solution, ModeSolver.compute_spectrum computes a gauge field spectrum
to perform a mode-by-mode cross check (unless nmodes=None).
If the mode-by-mode cross check is a success, the solution is stored in the underlying BGSystem of the class.
Otherwise, the GEFSolver tries to self correct using the gauge field spectrum. This is attempted for mbm_attempts or until successful.
Parameters
- ntr (int):
initial truncation number
GEFSolver.ntr - tend (float):
initial target time for
GEFSolver.tend - nmodes (int or None):
the number of modes computed by
ModeSolver. If None, no cross-check is performed. - mbm_attempts (int): number of mode-by-mode self correction attempts
- resume_mbm (bool):
if
TrueuseModeSolver.update_spectrumin case successive mode-by-mode comparisons are needed. - err_tol (float):
passed to
mbm_crosscheck. - err_thr (float):
passed to
mbm_crosscheck. - binning (int):
passed to
mbm_crosscheck. - integrator (str):
integrator for
mbm_crosscheck('simpson' is advised) - print_stats (bool):
if
True, a summary report is printed for the returned solution. - solver_kwargs: the
settingsofGEFSolver(seeGEFSolver.settings)
Returns
- vals (BGSystem):
contains the background evolution obtained from
GEFSolver - spec (GaugeSpec or None):
the result of
ModeSolver.compute_spectrum - sol: the full result of
GEFSolver.compute_GEF_solution
Raises
- RuntimeError: if no successful solution was obtained.
def
mbm_crosscheck( self, spec, vals: geff.bgtypes.BGSystem, err_tol: float, err_thr: float, binning: int, **integrator_kwargs):
263 def mbm_crosscheck(self, spec, vals:BGSystem, err_tol:float, err_thr:float, binning:int, **integrator_kwargs): 264 """ 265 Estimate the error of a GEF solution using `.mbm.GaugeSpec.estimate_GEF_error`. 266 267 If either the RMS error or the final error exceeds `err_tol`, the solution is rejected. 268 269 Parameters 270 ---------- 271 vals : BGSystem 272 contains the GEF solution. 273 err_tol : float 274 the tolerance on the RMS and final error. 275 err_thr : float 276 passed to `estimate_GEF_error`. 277 binning : int 278 passed to `estimate_GEF_error`. 279 integratorkwargs: 280 passed to kwargs of `estimate_GEF_error`. 281 282 Returns 283 ------- 284 agreement : bool 285 indicates if the solution is accepted or rejected. 286 reinit_slice : SpecSlice 287 the spectrum with which the GEF solver is re-initialized. 288 """ 289 GFs = self.GEFSolver.known_variables["gauge"] 290 agreement=True 291 292 for GF in GFs: 293 ref = [a.name for a in GF.zeros] 294 cut = GF.cutoff.name 295 296 errs, terr, _ = spec.estimate_GEF_error(vals, references=ref, cutoff=cut, 297 err_thr=err_thr, binning=binning, **integrator_kwargs) 298 299 reinit_inds = [] 300 301 for err in errs: 302 rmserr = np.sqrt(np.sum(err**2)/len(err)) 303 if max(err[-1], rmserr) > err_tol: 304 agreement=False 305 #find where the error is above 5%, take the earliest occurrence, reduce by 1 306 err_ind = np.where(err > err_thr)[0][0]-1 307 else: 308 err_ind = len(terr)-1 309 reinit_inds.append( err_ind ) 310 311 t0 = terr[min(reinit_inds)] 312 313 ind = np.searchsorted(spec["t"],t0, "left") 314 315 reinit_slice = spec.tslice(ind) 316 317 return agreement, reinit_slice
Estimate the error of a GEF solution using .mbm.GaugeSpec.estimate_GEF_error.
If either the RMS error or the final error exceeds err_tol, the solution is rejected.
Parameters
- vals (BGSystem): contains the GEF solution.
- err_tol (float): the tolerance on the RMS and final error.
- err_thr (float):
passed to
estimate_GEF_error. - binning (int):
passed to
estimate_GEF_error. - integratorkwargs:: passed to kwargs of
estimate_GEF_error.
Returns
- agreement (bool): indicates if the solution is accepted or rejected.
- reinit_slice (SpecSlice): the spectrum with which the GEF solver is re-initialized.
def
load_GEFdata(self, path: str):
319 def load_GEFdata(self, path : str): 320 """ 321 Load data and return a BGSystem with the results. 322 323 Note, data is always loaded assuming numerical units. 324 325 Parameters 326 ---------- 327 path : None or str 328 Path to load data from 329 330 Raises 331 ------ 332 FileNotFoundError 333 if no file is found at `path`. 334 AttributeError 335 if the file contains a column labeled by a key which does not match any GEF-value name. 336 """ 337 newsys = BGSystem.from_system(self.initial_data, copy=True) 338 339 newsys.load_variables(path) 340 341 return newsys
Load data and return a BGSystem with the results.
Note, data is always loaded assuming numerical units.
Parameters
- path (None or str): Path to load data from
Raises
- FileNotFoundError: if no file is found at
path. - AttributeError: if the file contains a column labeled by a key which does not match any GEF-value name.
The solver used to solve the GEF equations in run.
Inherited Members
The mode solver used for mode-by-mode cross checks.
def
compile_model(modelname: str, settings: dict = {}):
343def compile_model(modelname:str, settings:dict={}): 344 """ 345 Define a custom subclass of BaseGEF adapted to a new GEF model. 346 347 Parameters 348 ---------- 349 modelname : str 350 The name of the GEF model or a full dotted import path (e.g., "path.to.module"). 351 settings : dict 352 a dictionary of settings used by the model 353 354 Returns 355 ------- 356 CustomGEF 357 a custom subclass of BaseGEF 358 359 """ 360 model = load_model(modelname, settings) 361 signature = [inspect.Parameter("self", inspect.Parameter.POSITIONAL_ONLY)] 362 for val in model.model_input: 363 if issubclass(val, Func): 364 annotation = Callable 365 else: 366 annotation = float 367 signature.append(inspect.Parameter(val.name, inspect.Parameter.KEYWORD_ONLY, annotation=annotation)) 368 369 class CustomGEF(BaseGEF): 370 GEFSolver = model.solver 371 """The solver used to solve the GEF equations in `run`.""" 372 ModeSolver = model.MbM 373 """The mode solver used for mode-by-mode cross checks.""" 374 _input_signature = model.model_input 375 define_units = staticmethod(model.define_units) 376 377 CustomGEF.__init__.__signature__ = inspect.Signature(signature) 378 379 CustomGEF.__qualname__ = model.name 380 CustomGEF.__module__ = __name__ 381 382 return CustomGEF
Define a custom subclass of BaseGEF adapted to a new GEF model.
Parameters
- modelname (str): The name of the GEF model or a full dotted import path (e.g., "path.to.module").
- settings (dict): a dictionary of settings used by the model
Returns
- CustomGEF: a custom subclass of BaseGEF
class
ModelError(builtins.Exception):
384class ModelError(Exception): 385 """ 386 Exception for errors when compiling a GEF model. 387 """ 388 pass
Exception for errors when compiling a GEF model.