Complex Stellar Populations (#64)

This submodule introduces the possibility to create Complex Stellar Populations with BPASS data.

The two main modules are `csp.eventrate` and `csp.spectra`
`csp.eventrate` contains the `CSPEventRate` class, which can be used to calculate event rates from a given Stellar Formation History and Metallicity evolution.

The following methods are available in this class:

    `at_time`: calculates the event rates at a given lookback time.
    `over_time`: calculates the event rates over a lookback time (binned).
    `grid_at_time`: alternative input. Takes a grid of points instead of a function. Does the same as at_time.
    `grid_over_time`: alternative input. Takes a grid of points instead of a function. Does the same as over_time.

`csp.spectra` contains the `CSPSpectra` class, which can be used to calculate synthetic spectra for the given population.
The same methods are available in this class as `CSPEventRate`, but provide spectra rather than even rates.

`csp.sfh` contains a class with which a Stellar Formation History (SFH) can be made. A collection of common parametric SFH are available (see documentation), but a custom one can also be given.

As support for these high level functions, a collection of utility and helper functions has been added. These can be found in `csp.utils`. Some these are optimised for performance (`numba` JIT compiled) and are not user friendly.

* Created CSP archistecture and basic sfh class
* Added variables to constants
* implements utility functions & new test files
* Implements event rate calculation over time
* BPASS data_files for testing file_loading
* Adds `at_time` function to the event rate class
* Adds docstrings and extra tests
* Implements calculate spectrum at_time and over_time functions
* Add normalisation function for the BPASS spectrum
* fixed hrd plot
* adding csp folder
* minor changes
* Moved tests to CSP folder + fix in function name _type_check
* added a hokitype error
* Added some code review for Max and minor stylistic changes
* Added some code review for Max and minor stylistic changes
* Updated the tests
* bla
* Adjust tests and docstrings
* Adjust eventrate functions to take lists of functions
* Disabled spectra test
* Added parametric star formation histories to SFH and refactored some variable names
* fixing CI
* Auto stash before merge of "csp" and "official_hoki/csp"
* upload files for remote work
* Improve testing and functions
* Add SpectraCompiler test
* Adds extensive tests for the CSP package
* Update SpectraCompiler
* Documentation update + SFH input
* Removed CSP from utils due to cyclical import
* Updated setup.cfg to include csp_test_data subfolder
* Removes leftover print statement + prints loading precompiled failed
* Correct metallicity mid-points
* Added the parametric star formation histories as public functions
* at_time new default sample rate + load functions moved to hoki.load
* Fixed typo
* Removes cyclical import
* Adjusted tests that use SpectraCompiler
* Getting tox to work
* private optimised mass_per_time_bin calculation
* Refactored Event Rate calculation from array data
* Removed dependency and added more doc information
* BUG: pass on kwarg to spectra compiler.
* Fix SpectraCompiler
* BUG: use np instead of pd in SpectraCompiler.
* BUG: fix index and make test pass.
* BUG: fix another regression bug.
* Fix further regression bug.
* Clean up loading module
* Updates mass per bin calculation for vectorized funciton input
* Renames CSPEventRate functions and adds grid functions
* np.zeros to np.empty where possible
* Fixed incorrect test of at_time and over_time
* Set default sample rate for at_time
* Adds `cache` parameter to over_time_spectrum
* Adds grid functions for the spectra calculation
* Adds return_time_edges to grid_over_time & bug fix sfh vectorisation
* Bug fix in nr_time_bins (grid_over_time functions)
* Docs string cleanup
* Solve cyclical import error
* Fix tox and remove corner plot
* File check, combined data_compiler tests, renames func+variables

Co-authored-by: Max Briel <max.briel@auckland.ac.nz>
Co-authored-by: heloises <heloisefw@gmail.com>
Co-authored-by: Heloise <HeloiseS@users.noreply.github.com>
Co-authored-by: findesgh <findessp@yandex.ru>
Co-authored-by: Martin Glatzle <mglatzle@mpa-garching.mpg.de>

Author: 4 people

.github/workflows/testing.yml

@@ -9,7 +9,6 @@ on:
 
   pull_request:
 
-
 # A workflow run is made up of one or more jobs that can run sequentially or in parallel
 jobs:
   # This workflow contains a single job called "build"

README.rst

@@ -58,8 +58,7 @@ This way it'll show you each test as they pass or FAIL. In the pip and github ve
 
 .. code-block:: none
 
-
-   `astropy`, `numpy`, `pandas`, `matplotlib`, `pyyaml`, `wheel`, `emcee`, `corner`
+   `astropy`, `numpy`, `pandas`, `matplotlib`, `pyyaml`, `wheel`, `emcee`, `corner`, `numba`
 
 **Note:** Python 2 is not supported
 

hoki/__init__.py

@@ -33,4 +33,3 @@ class UnsupportedPythonError(Exception):
     #__all__ += ['module_input', 'module_output']
     # or you can keep everything from the subpackage with the following instead
     # __all__ += example_mod.__all__
-

hoki/age_utils.py

@@ -480,6 +480,7 @@ def calculate_distributions_dt_divided(obs_df, model):
     # We add the time bins in there because it can make plotting extra convenient.
     # distributions_df['time_bins'] = hoki.constants.BPASS_TIME_BINS
 
+
     return likelihoods_df
 
 
@@ -490,7 +491,7 @@ def calculate_distributions_dt_divided(obs_df, model):
 
 def calculate_sample_pdf(distributions_df, not_you=None):
     """
-    Multiplies together all the columns in given in DataFrame apart from the "time_bins" column
+    Multiplies together all the columns in given in DataFrame apart from the "time_axis" column
 
     Parameters
     ----------

hoki/constants.py

@@ -23,7 +23,23 @@
 OUTPUTS_PATH = settings['outputs_path'] # This constant for dev purposes.
 BPASS_TIME_BINS = np.arange(6.0, 11.1, 0.1)
 BPASS_TIME_INTERVALS = np.array([10**(t+0.05) - 10**(t-0.05) for t in BPASS_TIME_BINS])
-BPASS_TIME_WEIGHT_GRID = np.array([np.zeros((100,100)) + dt for dt in BPASS_TIME_INTERVALS])
+BPASS_TIME_WEIGHT_GRID = np.array([np.zeros((100, 100)) + dt for dt in BPASS_TIME_INTERVALS])
+
+BPASS_LINEAR_TIME_EDGES = np.append([0.0], 10**np.arange(6.05, 11.15, 0.1))
+BPASS_LINEAR_TIME_INTERVALS = np.diff(BPASS_LINEAR_TIME_EDGES)
+
+BPASS_METALLICITIES = ["zem5", "zem4", "z001","z002", "z003", "z004", "z006", "z008", "z010", "z014", "z020", "z030", "z040"]
+BPASS_NUM_METALLICITIES = np.array([0.00001, 0.0001, 0.001, 0.002, 0.003, 0.004, 0.006, 0.008, 0.010,
+                                    0.014, 0.020, 0.030, 0.040])
+BPASS_METALLICITY_MID_POINTS = (BPASS_NUM_METALLICITIES[1:] + BPASS_NUM_METALLICITIES[:-1])/2
+
+BPASS_EVENT_TYPES = ["Ia", "IIP", "II", "Ib", "Ic", "LGRB", "PISNe", "low_mass"]
+HOKI_NOW = 13.799e9
+
+BPASS_IMFS = ["imf_chab100", "imf_chab300", "imf100_100", "imf100_300",
+              "imf135_100", "imf135_300", "imf135all_100", "imf170_100",
+              "imf170_300"]
+
 
 # Create a deprecation warning when using dummy_dict
 dummy_dict = {'timestep': 0, 'age': 1, 'log(R1)': 2, 'log(T1)': 3, 'log(L1)': 4, 'M1': 5, 'He_core1': 6, 'CO_core1': 7,
@@ -202,4 +218,3 @@ def set_default_bpass_version(version):
 
     print('Looks like everything went well! You can check the path was correctly updated by looking at this file:'
           '\n'+path_to_settings)
-

hoki/csp/__init__.py

@@ -0,0 +1,110 @@
+"""
+Author: Max Briel
+
+Parent class of a complex stellar population
+"""
+
+from hoki.csp.sfh import SFH
+from hoki.utils.exceptions import HokiTypeError, HokiFormatError
+from hoki.constants import HOKI_NOW
+
+
+############################################
+# Complex Stellar Populations Parent Class #
+############################################
+
+
+class CSP(object):
+    """
+    Complex Stellar Population class
+
+    Notes
+    -----
+    Parent class for `CSPEventRate` and `CSPSpectra`
+
+    Attributes
+    ----------
+    now : `float`
+        The age of the universe.
+    """
+    now = HOKI_NOW
+
+    def __init__(self):
+        pass
+
+    def _type_check_histories(self, sfh, zeh):
+        """
+        Function to check sfh and zeh are the correct type and transform them
+        into a consistent format.
+
+        Notes
+        -----
+        sfh and zeh can either be a callable python function or a list of
+        callables.
+
+        Input
+        -----
+        sfh
+            Stellar formation history
+        zeh
+            Z evolution history
+
+        Returns
+        -------
+        `tuple` ([sfh callables,], [zeh callables,])
+            A tuple containing the sfh callables and zeh callables as arrays.
+        """
+        # Check sfh list
+        if isinstance(sfh, list):
+
+            # check if zeh is also a list
+            if isinstance(zeh, list):
+                # have to be equal lengths
+                if len(sfh) == len(zeh):
+                    # have to be all callables
+                    if (all(callable(val) for val in sfh) and all(callable(val) for val in zeh)):
+                        return (sfh, zeh)
+                    # A non-callable is present
+                    else:
+                        raise HokiTypeError(
+                            "SFH or ZEH contains an object that's not a SFH object or function.")
+                # sfh and zeh are not equal length
+                else:
+                    raise HokiFormatError(
+                        "sfh_functions and Z_functions must have the same length.")
+            # zeh is not a list
+            else:
+                # sfh has to be 1 length and a callable, zeh has to be a callable
+                if len(sfh) == 1:
+                    if (callable(sfh[0]) and callable(zeh)):
+
+                        return (sfh, [zeh])
+                    else:
+                        raise HokiTypeError(
+                            "SFH or ZEH contains an object that's not a SFH object or python callable.")
+                else:
+                    raise HokiFormatError(
+                        "SFH must have length 1, be a python callable, or a SFH object.")
+        # sfh is a callable
+        elif callable(sfh):
+            # zeh is a list
+            if isinstance(zeh, list):
+                # list has to be 1 long, because sfh is callable
+                if len(zeh) == 1:
+                    return ([sfh], zeh)
+                # list it too long
+                else:
+                    raise HokiFormatError(
+                        "ZEH must be either length 1 or a python callable")
+
+            # zeh is also a callable return
+            elif callable(zeh):
+                return ([sfh], [zeh])
+            else:
+                raise HokiTypeError(
+                    "ZEH is not a python callable or a list of callables.")
+        # sfh cannot be identified
+        else:
+            raise HokiTypeError(
+                "SFH type is not a python callable or a SFH object."
+            )