first iteration of documentation

.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

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"
 

docs/index.rst

@@ -29,12 +29,12 @@ Contents
 .. toctree::
    :maxdepth: 2
 
+   Module Reference <modules>
    Overview <readme>
    Contributions & Help <contributing>
    License <license>
    Authors <authors>
    Changelog <changelog>
-   Module Reference <api/modules>
 
 
 Indices and tables

docs/modules.rst

@@ -0,0 +1,7 @@
+stratagemc
+==========
+
+.. toctree::
+   :maxdepth: 4
+
+   stratagemc

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

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

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

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

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]

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,29 +144,40 @@ 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)
 
     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