diff --git a/.github/workflows/documentation.yml b/.github/workflows/documentation.yml new file mode 100644 index 0000000..4e6762e --- /dev/null +++ b/.github/workflows/documentation.yml @@ -0,0 +1,27 @@ +name: documentation + +on: push + +permissions: + contents: write + +jobs: + build: + runs-on: ubuntu-latest + steps: + - uses: actions/checkout@v4 + - uses: actions/setup-python@v5 + - name: Install dependencies + run: | + pip install -r docs/requirements.txt + - name: Sphinx build + run: | + sphinx-build docs _build + - name: Deploy to GitHub Pages + uses: peaceiris/actions-gh-pages@v3 + if: ${{ github.event_name == 'push' && github.ref == 'refs/heads/main' }} + with: + publish_branch: gh-pages + github_token: ${{ secrets.GITHUB_TOKEN }} + publish_dir: _build/ + force_orphan: true \ No newline at end of file diff --git a/docs/conf.py b/docs/conf.py index df9a99b..dd99399 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -28,31 +28,31 @@ # setup.py install" in the RTD Advanced Settings. # Additionally it helps us to avoid running apidoc manually -try: # for Sphinx >= 1.7 - from sphinx.ext import apidoc -except ImportError: - from sphinx import apidoc +# try: # for Sphinx >= 1.7 +# from sphinx.ext import apidoc +# except ImportError: +# from sphinx import apidoc -output_dir = os.path.join(__location__, "api") -module_dir = os.path.join(__location__, "../src/stratagemc") -try: - shutil.rmtree(output_dir) -except FileNotFoundError: - pass +# output_dir = os.path.join(__location__, "api") +# module_dir = os.path.join(__location__, "../src/stratagemc") +# try: +# shutil.rmtree(output_dir) +# except FileNotFoundError: +# pass -try: - import sphinx +# try: +# import sphinx - cmd_line = f"sphinx-apidoc --implicit-namespaces -f -o {output_dir} {module_dir}" +# cmd_line = f"sphinx-apidoc --implicit-namespaces -f -o {output_dir} {module_dir}" - args = cmd_line.split(" ") - if tuple(sphinx.__version__.split(".")) >= ("1", "7"): - # This is a rudimentary parse_version to avoid external dependencies - args = args[1:] +# args = cmd_line.split(" ") +# if tuple(sphinx.__version__.split(".")) >= ("1", "7"): +# # This is a rudimentary parse_version to avoid external dependencies +# args = args[1:] - apidoc.main(args) -except Exception as e: - print("Running `sphinx-apidoc` failed!\n{}".format(e)) +# apidoc.main(args) +# except Exception as e: +# print("Running `sphinx-apidoc` failed!\n{}".format(e)) # -- General configuration --------------------------------------------------- @@ -77,6 +77,8 @@ # Add any paths that contain templates here, relative to this directory. templates_path = ["_templates"] +napoleon_custom_sections = [('Returns', 'params_style')] + # The suffix of source filenames. source_suffix = ".rst" diff --git a/docs/index.rst b/docs/index.rst index b542fd7..62b73ab 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -29,12 +29,12 @@ Contents .. toctree:: :maxdepth: 2 + Module Reference Overview Contributions & Help License Authors Changelog - Module Reference Indices and tables diff --git a/docs/modules.rst b/docs/modules.rst new file mode 100644 index 0000000..34bc1f0 --- /dev/null +++ b/docs/modules.rst @@ -0,0 +1,7 @@ +stratagemc +========== + +.. toctree:: + :maxdepth: 4 + + stratagemc \ No newline at end of file diff --git a/docs/requirements.txt b/docs/requirements.txt index 2ddf98a..b4f9bbc 100644 --- a/docs/requirements.txt +++ b/docs/requirements.txt @@ -3,3 +3,9 @@ # under `install_requires` in `setup.cfg` is also listed here! sphinx>=3.2.1 # sphinx_rtd_theme +numba +numpy +pymc +scipy +tqdm +matplotlib \ No newline at end of file diff --git a/docs/stratagemc.rst b/docs/stratagemc.rst new file mode 100644 index 0000000..4d4b9c9 --- /dev/null +++ b/docs/stratagemc.rst @@ -0,0 +1,28 @@ +stratagemc package +================== + +Module contents +--------------- + +.. automodule:: stratagemc + :imported-members: + :members: + :undoc-members: + :show-inheritance: + :special-members: __init__ + +Submodules +---------- + +stratagemc.geochron module +-------------------------- + +.. automodule:: stratagemc.geochron + :members: + :undoc-members: + :show-inheritance: + :special-members: __init__ + :private-members: + + + diff --git a/setup.cfg b/setup.cfg index 37e6762..7891762 100644 --- a/setup.cfg +++ b/setup.cfg @@ -5,7 +5,7 @@ [metadata] name = stratagemc -description = Add a short description here! +description = Bayesian age modeling of stratigraphy. author = Adrian Tasistro-Hart author_email = adrian@tasistro-hart.com license = GPL-3.0-only @@ -49,7 +49,12 @@ package_dir = # For more information, check out https://semver.org/. install_requires = importlib-metadata; python_version<"3.8" - + numba + numpy + pymc + scipy + tqdm + matplotlib [options.packages.find] where = src diff --git a/src/stratagemc/__init__.py b/src/stratagemc/__init__.py index 2f4d93b..142d7be 100644 --- a/src/stratagemc/__init__.py +++ b/src/stratagemc/__init__.py @@ -16,3 +16,18 @@ del version, PackageNotFoundError from .stratagemc import * + +__all__ = ['geochron_height_check', + 'trim_units', + 'get_times', + 'sigmoid', + 'DT_logp_l_gen', + 'randlike_gen', + 'loglike_gen', + 'floating_age', + 'age', + 'fit_floating_model', + 'age_depth', + 'model_ls', + 'model', + 'model2ages'] diff --git a/src/stratagemc/geochron.py b/src/stratagemc/geochron.py index 6e31a32..39b4268 100644 --- a/src/stratagemc/geochron.py +++ b/src/stratagemc/geochron.py @@ -8,14 +8,17 @@ class Geochron: """Class for handling geochronologic constraints in stratigraphic sections.""" - def __init__(self, h, rv, dt, prob_threshold=1e-16): + def __init__(self, h, rv, dt, prob_threshold=1e-8): """ Initializes the Geochron class. + Args: - h (array-like): List of stratigraphic heights at which temporal constraints are measured. - The coordinate increases up section. + h (array-like): List of stratigraphic heights at which temporal constraints + are measured. The coordinate increases up section. rv (array-like): List of random variables representing the temporal constraints at each height. The members of this list must be objects similar to stats.dist objects with the following methods: pdf, ppf (inverse of cdf). dt (float): Time spacing for computing the time increment pdfs. Units are same as rvs + prob_threshold (float, optional): Probability threshold for determining the temporal range over which to grid. Defaults to 1e-8. + Attributes: h (list): Sorted list of stratigraphic heights. rv (list): Sorted list of random variables representing the temporal constraints. @@ -88,16 +91,21 @@ def model_weights(self, units): """ Calculate the weights for units and contacts based on their bounding by geochron constraints. + Args: - units: 2d array-like (nx2) of bottom and top heights for each of n units in the section + units (numpy.ndarray): Bottom and top heights of units in section. + 2d array (nx2) of bottom and top heights for each of n units + Returns: Tuple[List[List[float]], List[List[float]]]: - alpha: 2D list (NxL) of weights for each of N units and L pairs of geochron constraints. - beta: 2D list (MxL) of weights for each of M contacts and L pairs of geochron constraints. + Raises: AssertionError: If any rows in alpha or beta are entirely zero, indicating unconstrained model parameters. + Notes: - - Run after trim_units. + - Run after :py:func:`stratagemc.trim_units()`. """ unit_thicks = np.diff(units, axis=1).flatten() contacts = units[0:-1, 1] diff --git a/src/stratagemc/stratagemc.py b/src/stratagemc/stratagemc.py index 8209200..9cc7b6b 100644 --- a/src/stratagemc/stratagemc.py +++ b/src/stratagemc/stratagemc.py @@ -25,9 +25,12 @@ def geochron_height_check(unit_heights, geochron_heights): """ Ensure that geochron heights are within units and not at contacts. Also ensure that geochron heights are within the section. Run before trimming section to the top and bottom of the geochron constraints. + Args: - unit_heights: Bottom and top heights of units in section. The heights are given as a 2d array-like (nx2) of bottom and top heights for each of n units, increasing downwards in array - geochron_heights: Heights of geochron constraints in the section, array-like + unit_heights (numpy.ndarray): Bottom and top heights of units in section. + 2d array (nx2) of bottom and top heights for each of n units + geochron_heights (arraylike): Heights of geochron constraints in the section + Raises: AssertionError: If any of the geochron heights are below the section or above the section. AssertionError: If any of the geochron heights are at contacts. @@ -46,12 +49,15 @@ def geochron_height_check(unit_heights, geochron_heights): def trim_units(unit_heights, geochron_heights): """ Trim the top and bottom of the section to the top and bottom of the geochron constraints. Run after geochron_height_check. + Args: - unit_heights: Heights of contacts in section, including top and bottom. + unit_heights (numpy.ndarray): Bottom and top heights of units in section. 2d array-like (nx2) of bottom and top heights for each of n units - geochron_heights: Heights of geochron constraints in section + geochron_heights (arraylike): Heights of geochron constraints in section + Returns: - unit_heights_trim: Trimmed unit heights after adjusting for the top and bottom units + numpy.ndarray: Trimmed unit heights after adjusting for the top and bottom units + Raises: None """ @@ -76,13 +82,15 @@ def trim_units(unit_heights, geochron_heights): def get_times(sed_rates, hiatuses, units): """ Floating times array for units. + Args: - sed_rates: Sedimentation rates for each unit, array-like - hiatuses: Hiatuses between units, array-like - units: Heights of contacts in section, including top and bottom. - 2d array-like (nx2) of bottom and top heights for each of n units + sed_rates (arraylike): Sedimentation rates for each unit + hiatuses (arraylike): Hiatuses between units + units (numpy.ndarray): Bottom and top heights of units in section. + 2d array (nx2) of bottom and top heights for each of n units + Returns: - times: Floating times array for units, array-like + times (numpy.ndarray): Floating times array for units, same shape as units """ n_units = units.shape[0] assert len(sed_rates) == n_units, 'must have sed rate for each unit' @@ -118,14 +126,16 @@ def sigmoid(t, scale=0.001): def DT_logp_l_gen(pdt): """ Generates a log-probability function for a numerical time increment distribution. Returned function interpolates a log probability and takes as inputs time increment(s) at which to evaluate as well as (required) time increment coordinates corresponding to the numerical pdf associated with the function. Sped up with numba. + Args: pdt (array-like): The probability density function values for the time increments. + Returns: - function: A function that computes the log-probability of a given time increment. - The returned function, DT_logp, takes the following parameters: - dt_query (float): The time increment for which the log-probability is to be computed. - dt (array-like): The array of time increments corresponding to the probability density function values. Must be same length as pdt. - The DT_logp function computes the log-probability by interpolating the probability density function values and applying sigmoid functions to ensure the values are within the valid range. + numba.core.registry.CPUDispatcher: A function that computes the log-probability of a given time increment. + The returned function, DT_logp, takes the following parameters: + dt_query (float): The time increment for which the log-probability is to be computed. + dt (array-like): The array of time increments corresponding to the probability density function values. Must be same length as pdt. + The DT_logp function computes the log-probability by interpolating the probability density function values and applying sigmoid functions to ensure the values are within the valid range. """ @njit def DT_logp(dt_query, dt): @@ -134,14 +144,21 @@ def DT_logp(dt_query, dt): def randlike_gen(geochron, units): - """Generate function for generating random draws from likelihood. Permits generating prior_predictive samples from a model with the CustomDist likelihood. + """ + Generate function for generating random draws from likelihood. Permits generating prior_predictive samples from a model with the CustomDist likelihood. Args: geochron (geochron.Geochron): Geochron object containing geochronologic constraints. - units (ndarray): nx2 array of unit bottom and top heights for n units. + units (numpy.ndarray): Bottom and top heights of units in section. + 2d array (nx2) of bottom and top heights for each of n units Returns: function: Function for generating random draws from likelihood. + The returned function, eps_rand, takes the following parameters: + | params (array-like): Array of parameters; concatenation of sedimentation rates and hiatuses. + | rng (numpy.random.Generator): Random number generator. + | size (int or tuple of ints, optional): Output shape. Default is None. Last dimension must be the number of pairs of geochron constraints. + The function returns a numpy.ndarray of random draws. """ n_units = units.shape[0] alpha, beta = geochron.model_weights(units) @@ -149,14 +166,18 @@ def randlike_gen(geochron, units): def eps_rand(params, rng, size): """ Generate random samples from likelihood. Params is the concatenation of sedimentation rates and hiatuses. + Args: params (array-like): Array of parameters; concatenation of sedimentation rates and hiatuses. rng (numpy.random.Generator): Random number generator. size (int or tuple of ints, optional): Output shape. Default is None. Last dimension must be the number of pairs of geochron constraints. + Returns: numpy.ndarray: Random draws. + Raises: AssertionError: If size is incompatible with the number of pairs. + Notes: - The function calculates the random draws based on the given parameters and geochron constraints. - If size is not provided or is None, it will default to the number of pairs. @@ -360,6 +381,7 @@ def time_offset_cost(offset): def age_depth(units, times): """ Convert units, times arrays to age-depth curve + Args: units (numpy.ndarray): An nx2 array of bottom and top elevations for each of n units. times (numpy.ndarray): An nx2 array of bottom and top times for each of n units. @@ -383,9 +405,10 @@ def model_ls(units, geochron, geochron (geochron.Geochron): Geochron object containing geochron constraints. sed_rate_bounds (list, optional): Bounds on sedimentation rates. Defaults to None. hiatus_bounds (list, optional): Bounds on hiatuses. Defaults to None. + Returns: - ndarray: Sedimentation rates for each unit. - ndarray: Hiatuses between units. + sed_rates (numpy.ndarray): Sedimentation rates for each unit. + hiatuses (numpy.ndarray): Hiatuses between units. """ alpha, beta = geochron.model_weights(units) # model matrix @@ -433,6 +456,7 @@ def model(units, geochron, sed_rates_prior, hiatuses_prior, hiatuses_prior (function): Prior distribution for hiatuses. Must be valid as dist argument to pymc.CustomDist(dist=dist). Signature is hiatus_prior(size=size). draws (int, optional): Number of MCMC draws. Defaults to 1000. **kwargs: Additional keyword arguments to pass to the MCMC sampler. + Returns: arviz.InferenceData: ArviZ InferenceData object containing the MCMC trace. """ @@ -483,4 +507,4 @@ def model2ages(trace, n_posterior=None): # if n_posterior is None, use effective sample size to dictate number of posterior samples if n_posterior is None: # n_posterior = trace.n_eff - return + return