Merge pull request #2 from skumblex/dev

Dev

Author: hep-mh

.gitignore

@@ -3,6 +3,8 @@
 /plots/*.pdf
 /plots/*
 
+tmp/*
+
 # Byte-compiled / optimized / DLL files
 __pycache__
 *.py[cod]

CHANGELOG

@@ -3,7 +3,7 @@
 **A generiC fRamework fOr Photodisintegration Of LIght elementS**
 
 ![Language: Python3](https://img.shields.io/badge/language-Python3-blue.svg?style=flat-square)
-![Version: 1.1](https://img.shields.io/badge/current_version-1.1-blue.svg?style=flat-square)
+![Version: 1.2](https://img.shields.io/badge/current_version-1.2-blue.svg?style=flat-square)
 
 When using this code, please cite the following papers
 
@@ -19,6 +19,14 @@ The remarkable agreement between observations of the primordial light element ab
 
 # Changelog
 
+v1.2
+ - Speed improvements when running non-thermal nucleosynthesis (by a factor 7)
+ - Modified the directory structure by moving ./data to ./acropolis/data to transform ACROPOLIS into a proper package, which can be installed via ``python3 setup.py install --user`` (also putting the executables ``decay`` and ``annihilation`` into your ``PATH``)
+ - Added the decay of neutrons and tritium to the calculation
+ - Included a new script 'tools/create_sm_cosmo_file.py' which allows to generate the file cosmo_file.dat for sm.tar.gz and can easily be modified by the user
+ - For AnnihilationModel, it is now possible to freely choose the dark-matter density parameter (default is 0.12)
+
+
 v1.1
  - For the source terms it is now possible to specify arbitrary monochromatic and continuous contributions, meaning that the latter one is no longer limited to only final-state radiation of photons
  - By including additional JIT compilation steps, the runtime without database files was drastically increased (by approximately a factor 15)
@@ -33,7 +41,7 @@ v1.0
 
 # Install the dependencies
 
-ACROPOLIS is written in Python3.7 (remember that Python2 is dead) and depends on the following packages (older versions might work, but have not been thoroughly testes)
+ACROPOLIS is written in Python3.7 (remember that Python2 is dead) and depends on the following packages (older versions might work, but have not been thoroughly tested)
 
  - NumPy (> 1.19.1)
  - SciPy (>1.5.2)
@@ -47,9 +55,19 @@ python3 -m pip install numpy, scipy, numba --user
 
 If these dependencies conflict with those for other programs in your work environment, it is strongly advised to utilise the capabilities of Python's virtual environments.
 
+# Install ACROPOLIS using pip
+
+ACROPOLIS also comes with a ``setup.py`` file, which allows to simply install it via pip. To do this, simply run the following command in the root of the cloned repository
+
+```
+python3 -m pip install .
+```
+
+Afterwards, the different packages of ACROPOLIS can be imported into our own code, just like any other python package. The above command, also copies the executable ``decay`` and ``annihilation`` into your ``PATH`` and makes sure that all dependencies are fulfilled.
+
 # Use the example models
 
-Within the ACROPOLIS main directory there are two executables, ``decay`` and ``annihilation``, which wrap the scenarios discussed in section 4.1 and section 4.2 from the manual, respectively. Both of these files need to be called with six command-line arguments each, a list of which can be obtained by running the command of choice without any arguments at all. On that note, the following command runs the process of photodisintegration for an unstable mediator with a mass of 10MeV and a lifetime of 1e5s that decays exclusively into photons and has an abundance of 1e-10 relative to photons at a reference temperature of 10MeV
+Within the ACROPOLIS main directory there are two executables, ``decay`` and ``annihilation``, which wrap the scenarios discussed in section 4.1 and section 4.2 from the manual, respectively. Both of these files need to be called with six command-line arguments each, a list of which can be obtained by running the command of choice without any arguments at all. On that note, the following command runs the process of photodisintegration for an unstable mediator with a mass of 10MeV and a lifetime of 1e5s that decays exclusively into photons and has an abundance of 1e-10 relative to photons at a reference temperature of 10MeV (**if you installed ACROPOLIS using pip, you can drop the ``./`` at the beginning, since the executables are in your ``PATH``.**)
 
 ```
 ./decay 10 1e5 10 1e-10 0 1

README.md

@@ -0,0 +1,29 @@
+# functools
+from functools import wraps
+
+
+def cached_member(f_uncached):
+    # Define the cache as a dictionary
+    cache = {}
+    cT = {"_": -1.}
+
+    # Define the wrapper function
+    @wraps(f_uncached)
+    def f_cached(*args):
+        T     = args[-1]
+        # Drop the first argument 'self'
+        # ! only for member functions !
+        pargs = args[1:]
+
+        # For each new temperature,
+        # clear the cache and start over
+        if T != cT["_"]:
+            cT["_"] = T
+            cache.clear()
+
+        if pargs not in cache:
+            cache[pargs] = f_uncached(*args)
+
+        return cache[pargs]
+
+    return f_cached

acropolis/cache.py

@@ -1,7 +1,5 @@
 # math
 from math import pi, log, log10, exp, sqrt
-# functools
-from functools import wraps
 # numpy
 import numpy as np
 # scipy
@@ -13,41 +11,18 @@
 import warnings
 
 # db
-from .db import import_data_from_db
-from .db import in_rate_db, interp_rate_db
+from acropolis.db import import_data_from_db
+from acropolis.db import in_rate_db, interp_rate_db
+# cache
+from acropolis.cache import cached_member
 # pprint
-from .pprint import print_warning, print_error
+from acropolis.pprint import print_warning, print_error
 # params
-from .params import me, me2, alpha, re
-from .params import zeta3, pi2
-from .params import Emin
-from .params import approx_zero, eps, Ephb_T_max, E_EC_cut
-from .params import NE_pd, NE_min
-
-
-def _cached(f_uncached):
-    # Define the cache as a dictionary
-    cache = {}
-    cT = {"_": -1.}
-
-    # Define the wrapper function
-    @wraps(f_uncached)
-    def f_cached(*args):
-        T     = args[-1]
-        pargs = args[1:]
-
-        # For each new temperature,
-        # clear the cache and start over
-        if T != cT["_"]:
-            cT["_"] = T
-            cache.clear()
-
-        if pargs not in cache:
-            cache[pargs] = f_uncached(*args)
-
-        return cache[pargs]
-
-    return f_cached
+from acropolis.params import me, me2, alpha, re
+from acropolis.params import zeta3, pi2
+from acropolis.params import Emin
+from acropolis.params import approx_zero, eps, Ephb_T_max, E_EC_cut
+from acropolis.params import NE_pd, NE_min
 
 
 # _ReactionWrapperScaffold ####################################################
@@ -79,8 +54,9 @@ def _JIT_G(Ee, Eph, Ephb):
     # ATTENTION: This range is absent in 'astro-ph/9412055'
     # Here we adopt the original result from
     # 'link.springer.com/content/pdf/10.1007/BF01005624.pdf'
-    Ee_lim_m = ( Eph + Ephb - (Eph - Ephb)*sqrt( 1. - me2/( Eph*Ephb ) ) )/2.
-    Ee_lim_p = ( Eph + Ephb + (Eph - Ephb)*sqrt( 1. - me2/( Eph*Ephb ) ) )/2.
+    dE_sqrt  = (Eph - Ephb)*sqrt( 1. - me2/( Eph*Ephb ) )
+    Ee_lim_m = ( Eph + Ephb - dE_sqrt )/2.
+    Ee_lim_p = ( Eph + Ephb + dE_sqrt )/2.
 
     if not ( me < Ee_lim_m <= Ee <= Ee_lim_p ):
         # CHECKED to never happen, since the intergration
@@ -104,23 +80,23 @@ def _JIT_G(Ee, Eph, Ephb):
 # _PhotonReactionWrapper ######################################################
 
 @nb.jit(cache=True)
-def _JIT_PH_rate_pair_creation(y, x, T):
+def _JIT_ph_rate_pair_creation(y, x, T):
     # Define beta as a function of y
     b = sqrt(1. - 4.*me2/y)
 
     # Define the kernel for the 2d-integral; y = s, x = epsilon_bar
     #                     f/E^2                              s   \sigma_DP
-    return ( 1./(pi**2) )/( exp(x/T) - 1. ) * y * .5*pi*(re**2.)*(1.-b**2.)*( (3.-b**4.)*log( (1.+b)/(1.-b) ) - 2.*b*(2.-b**2.) )
     # ATTENTION: There is an error in 'astro-ph/9412055.pdf'
     # In the integration for \bar{\epsilon}_\gamma the lower
     # limit of integration should be me^2/\epsilon_\gamma
     # (the written limit is unitless, which must be wrong)
     # This limit is a consequence of the constraint on
     # the center-of-mass energy
+    return ( 1./(pi**2) )/( exp(x/T) - 1. ) * y * .5*pi*(re**2.)*(1.-b**2.)*( (3.-b**4.)*log( (1.+b)/(1.-b) ) - 2.*b*(2.-b**2.) )
 
 
 @nb.jit(cache=True)
-def _JIT_PH_kernel_inverse_compton(y, E, Ep, T):
+def _JIT_ph_kernel_inverse_compton(y, E, Ep, T):
     # Return the integrand for the 1d-integral in log-space; x = Ephb
     x = exp(y)
 
@@ -130,21 +106,21 @@ def _JIT_PH_kernel_inverse_compton(y, E, Ep, T):
 # _ElectronReactionWrapper ####################################################
 
 @nb.jit(cache=True)
-def _JIT_EL_rate_inverse_compton(y, x, E, T):
+def _JIT_el_rate_inverse_compton(y, x, E, T):
     # Return the integrand for the 2d-integral; y = Eph, x = Ephb
     return _JIT_F(y, E, x)*x/( (pi**2.)*(exp(x/T) - 1.) )
 
 
 @nb.jit(cache=True)
-def _JIT_EL_kernel_inverse_compton(y, E, Ep, T):
+def _JIT_el_kernel_inverse_compton(y, E, Ep, T):
     # Define the integrand for the 1d-integral in log-space; x = Ephb
     x = exp(y)
 
     return _JIT_F(Ep+x-E, Ep, x)*( x/(pi**2) )/( exp(x/T) - 1. ) * x
 
 
 @nb.jit(cache=True)
-def _JIT_EL_kernel_pair_creation(y, E, Ep, T):
+def _JIT_el_kernel_pair_creation(y, E, Ep, T):
     # Define the integrand for the 1d-integral in log-space; x = Ephb
     x = exp(y)
 
@@ -184,7 +160,7 @@ def _JIT_dsdE_Z2(Ee, Eph):
 
 @nb.jit(cache=True)
 def _JIT_set_spectra(F, i, Fi, cond):
-    F[:,i] = Fi
+    F[:, i] = Fi
     # In the strongly compressed regime, manually
     # set the photon spectrum to zero in order to
     # avoid floating-point errors
@@ -329,7 +305,7 @@ def _rate_pair_creation(self, E, T):
         # In general, the threshold is E ~ me^2/(22*T)
         # However, here we use a slighlty smaller threshold
         # in order to guarantee a smooth transition
-        if E < me2/(30.*T):
+        if E < me2/(50.*T):
             return 0.
 
         # Define the integration limits from the
@@ -340,13 +316,13 @@ def _rate_pair_creation(self, E, T):
         # CHECKED!
 
         # Perform the integration in lin space
-        I_fso_E2 = dblquad(_JIT_PH_rate_pair_creation, llim, ulim, lambda x: 4.*me2, lambda x: 4.*E*x, epsrel=eps, epsabs=0, args=(T,))
+        I_fso_E2 = dblquad(_JIT_ph_rate_pair_creation, llim, ulim, lambda x: 4.*me2, lambda x: 4.*E*x, epsrel=eps, epsabs=0, args=(T,))
 
         return I_fso_E2[0]/( 8.*E**2. )
 
 
     def _rate_pair_creation_db(self, E, T):
-        if E < me2/(30.*T):
+        if E < me2/(50.*T):
             return 0.
 
         E_log, T_log = log10(E), log10(T)
@@ -385,7 +361,7 @@ def _kernel_compton(self, E, Ep, T):
 
 
     # INVERSE COMPTON SCATTERING ##############################################
-    @_cached
+    @cached_member
     def _kernel_inverse_compton(self, E, Ep, T):
         # Incorporate the non-generic integration limit as
         # the algorithm requires Ep > E and not Ep > E + me
@@ -410,7 +386,7 @@ def _kernel_inverse_compton(self, E, Ep, T):
             return 0.
 
         # Perform the integration in log space
-        I_fF_E = quad(_JIT_PH_kernel_inverse_compton, log(llim), log(ulim), epsrel=eps, epsabs=0, args=(E, Ep, T))
+        I_fF_E = quad(_JIT_ph_kernel_inverse_compton, log(llim), log(ulim), epsrel=eps, epsabs=0, args=(E, Ep, T))
 
         # ATTENTION: Kawasaki considers a combined e^+/e^- spectrum
         # Therefore the factor 2 should not be there in our case
@@ -445,7 +421,7 @@ def __init__(self, Y0, eta, db):
     # T is the temperature of the background photons
 
     # INVERSE COMPTON SCATTERING ##############################################
-    @_cached
+    @cached_member
     def _rate_inverse_compton(self, E, T):
         # Define the upper limit for the integration over x
         ulim = min( E - me2/(4.*E), Ephb_T_max*T )
@@ -459,7 +435,7 @@ def _rate_inverse_compton(self, E, T):
         # ATTENTION:
         # The integral over \epsilon_\gamma should start at 0.
         # In fact, for \epsilon_\gamma > \epsilon_e, we have q < 0.
-        I_fF_E = dblquad(_JIT_EL_rate_inverse_compton, 0., ulim, lambda x: x, lambda x: 4.*x*E*E/( me2 + 4.*x*E ), epsrel=eps, epsabs=0, args=(E, T))
+        I_fF_E = dblquad(_JIT_el_rate_inverse_compton, 0., ulim, lambda x: x, lambda x: 4.*x*E*E/( me2 + 4.*x*E ), epsrel=eps, epsabs=0, args=(E, T))
 
         return 2.*pi*(alpha**2.)*I_fF_E[0]/(E**2.)
 
@@ -483,7 +459,7 @@ def total_rate(self, E, T):
     # T  is the temperature of the background photons
 
     # INVERSE COMPTON SCATTERING ##############################################
-    @_cached
+    @cached_member
     def _kernel_inverse_compton(self, E, Ep, T):
         # E == Ep leads to a divergence in
         # the Bose-Einstein distribution
@@ -492,11 +468,12 @@ def _kernel_inverse_compton(self, E, Ep, T):
             return 0.
 
         # Calculate appropriate integration limits
-        pf = .25*me2/Ep - E                     # <= 0.
-        qf = .25*me2*(Ep-E)/Ep                  # >= 0.
+        pf = .25*me2/Ep - E               # <= 0.
+        qf = .25*me2*(Ep-E)/Ep            # >= 0.
 
-        z1 = -pf/2. - sqrt( (pf/2.)**2 - qf )  # smaller
-        z2 = -pf/2. + sqrt( (pf/2.)**2 - qf )  # larger
+        sqrt_d = sqrt( (pf/2.)**2. - qf )
+        z1 = -pf/2. - sqrt_d              # smaller
+        z2 = -pf/2. + sqrt_d              # larger
 
         # Define the integration limits from
         # the range that is specified in '_JIT_F'
@@ -513,7 +490,7 @@ def _kernel_inverse_compton(self, E, Ep, T):
             return 0.
 
         # Perform the integration in log space
-        I_fF_E = quad(_JIT_EL_kernel_inverse_compton, log(llim), log(ulim), epsrel=eps, epsabs=0, args=(E, Ep, T))
+        I_fF_E = quad(_JIT_el_kernel_inverse_compton, log(llim), log(ulim), epsrel=eps, epsabs=0, args=(E, Ep, T))
 
         return 2.*pi*(alpha**2.)*I_fF_E[0]/(Ep**2.)
 
@@ -541,7 +518,7 @@ def _kernel_compton(self, E, Ep, T):
 
 
     # BETHE_HEITLER PAIR CREATION #############################################
-    @_cached
+    @cached_member
     def _kernel_bethe_heitler(self, E, Ep, T):
         # Incorporate the non-generic integration limit as
         # the algorithm requires Ep > E and not Ep > E + me
@@ -553,13 +530,13 @@ def _kernel_bethe_heitler(self, E, Ep, T):
 
 
     # DOUBLE PHOTON PAIR CREATION #############################################
-    @_cached
+    @cached_member
     def _kernel_pair_creation(self, E, Ep, T):
         # In general, the threshold is Ep >~ me^2/(22*T)
         # However, here we use a slighlty smaller threshold
         # in acordance with the implementation we use in
         # '_PhotonReactionWrapper._rate_pair_creation'
-        if Ep < me2/(30.*T):
+        if Ep < me2/(50.*T):
             return 0.
 
         dE, E2 = Ep - E, E**2.
@@ -586,7 +563,7 @@ def _kernel_pair_creation(self, E, Ep, T):
             return 0.
 
         # Perform the integration in log space
-        I_fG_E2 = quad(_JIT_EL_kernel_pair_creation, log(llim), log(ulim), epsrel=eps, epsabs=0, args=(E, Ep, T))
+        I_fG_E2 = quad(_JIT_el_kernel_pair_creation, log(llim), log(ulim), epsrel=eps, epsabs=0, args=(E, Ep, T))
 
         return 0.25*pi*(alpha**2.)*me2*I_fG_E2[0]/(Ep**3.)