Merge branch 'dev'

docs/api-invisible.rst

@@ -14,6 +14,6 @@
    ProxyRecord
 
    read_14c
-   read_dates
+   read_chron
    read_proxy
 

docs/api.rst

@@ -62,7 +62,7 @@ Creating a ChronRecord
    :toctree: generated/
 
    ChronRecord
-   read_dates
+   read_chron
 
 
 Attributes

docs/whats_new.rst

@@ -10,7 +10,9 @@ v0.0.2
 
 Enhancements
 ~~~~~~~~~~~~
-- DateRecord class has been renamed ChronRecord to avoid confusion with DatedProxyRecord (`Issue #3 <https://github.com/brews/snakebacon/issues/3>`_).
+- DateRecord class has been renamed ChronRecord to avoid confusion with DatedProxyRecord (`Issue #3 <https://github.com/brews/snakebacon/issues/3>`_). The read_chron() is also now read_chron().
+
+- AgeDepthModels now can plot prior distributions.
 
 - Can now query AgeDepthModels to get prior distribution for chronology dates, sediment rates, and sediment rate memory (`Issue #5 <https://github.com/brews/snakebacon/issues/5>`_).
 

snakebacon/__init__.py

@@ -1,6 +1,6 @@
 from .agedepth import AgeDepthModel
 from .records import CalibCurve, ChronRecord, ProxyRecord, DatedProxyRecord
-from .records import read_14c, read_dates, read_proxy
+from .records import read_14c, read_chron, read_proxy
 from .utils import suggest_accumulation_rate
 
 try:

snakebacon/agedepth.py

@@ -1,6 +1,9 @@
 import logging as logging
 
+import scipy
 import matplotlib.pylab as plt
+from matplotlib.collections import PatchCollection
+from matplotlib.patches import Polygon
 import numpy as np
 
 from .mcmc import McmcSetup
@@ -109,7 +112,7 @@ def plot(self, agebins=50, p=(2.5, 97.5), ax=None):
         ax.step(self.depth, self.age_percentile(p[1]), where='mid', color='red', linestyle=':')
         ax.set_ylabel('Age (cal yr BP)')
         ax.set_xlabel('Depth (cm)')
-        ax.grid()
+        ax.grid(True)
         return ax
 
     def agedepth(self, d):
@@ -146,13 +149,89 @@ def agedepth(self, d):
         return out
 
     def prior_dates(self):
-        self.mcmcsetup.prior_dates()
+        return self.mcmcsetup.prior_dates()
+
+    def plot_prior_dates(self, dwidth=30, ax=None):
+        """Plot prior chronology dates in age-depth plot"""
+        if ax is None:
+            ax = plt.gca()
+        depth, probs = self.prior_dates()
+        pat = []
+        for i, d in enumerate(depth):
+            p = probs[i]
+            z = np.array([p[:, 0], dwidth * p[:, 1] / np.sum(p[:, 1])])  # Normalize
+            z = z[:, z[0].argsort(kind='mergesort')]  # np.interp requires `xp` arg to be sorted
+            zy = np.linspace(np.min(z[0]), np.max(z[0]), num=200)
+            zp = np.interp(x=zy, xp=z[0], fp=z[1])
+            pol = np.vstack([np.concatenate([d + zp, d - zp[::-1]]),
+                             np.concatenate([zy, zy[::-1]])])
+            pat.append(Polygon(pol.T))
+        p = PatchCollection(pat)
+        p.set_labels('Prior dates')
+        ax.add_collection(p)
+        ax.autoscale_view()
+        ax.set_ylabel('Age (cal yr BP)')
+        ax.set_xlabel('Depth (cm)')
+        ax.grid(True)
+        return ax
 
     def prior_sediment_rate(self):
-        self.mcmcsetup.prior_sediment_rate()
+        return self.mcmcsetup.prior_sediment_rate()
+
+    def plot_sediment_rate(self, ax=None):
+        """Plot sediment accumulation rate prior and posterior distributions"""
+        if ax is None:
+            ax = plt.gca()
+
+        y_prior, x_prior = self.prior_sediment_rate()
+        ax.plot(x_prior, y_prior, label='Prior')
+
+        y_posterior = self.mcmcfit.sediment_rate
+        density = scipy.stats.gaussian_kde(y_posterior.flat)
+        density.covariance_factor = lambda: 0.25
+        density._compute_covariance()
+        ax.plot(x_prior, density(x_prior), label='Posterior')
+
+        acc_shape = self.mcmcsetup.mcmc_kwargs['acc_shape']
+        acc_mean = self.mcmcsetup.mcmc_kwargs['acc_mean']
+        annotstr_template = 'acc_shape: {0}\nacc_mean: {1}'
+        annotstr = annotstr_template.format(acc_shape, acc_mean)
+        ax.annotate(annotstr, xy=(0.9, 0.9), xycoords='axes fraction',
+                    horizontalalignment='right', verticalalignment='top')
+
+        ax.set_ylabel('Density')
+        ax.set_xlabel('Acc. rate (yr/cm)')
+        ax.grid(True)
+        return ax
 
     def prior_sediment_memory(self):
-        self.mcmcsetup.prior_sediment_memory()
+        return self.mcmcsetup.prior_sediment_memory()
+
+    def plot_sediment_memory(self, ax=None):
+        """Plot sediment memory prior and posterior distributions"""
+        if ax is None:
+            ax = plt.gca()
+
+        y_prior, x_prior = self.prior_sediment_memory()
+        ax.plot(x_prior, y_prior, label='Prior')
+
+        y_posterior = self.mcmcfit.sediment_memory
+        density = scipy.stats.gaussian_kde(y_posterior ** (1/self.thick))
+        density.covariance_factor = lambda: 0.25
+        density._compute_covariance()
+        ax.plot(x_prior, density(x_prior), label='Posterior')
+
+        mem_mean = self.mcmcsetup.mcmc_kwargs['mem_mean']
+        mem_strength = self.mcmcsetup.mcmc_kwargs['mem_strength']
+        annotstr_template = 'mem_strength: {0}\nmem_mean: {1}\nK: {2}'
+        annotstr = annotstr_template.format(mem_strength, mem_mean, self.thick)
+        ax.annotate(annotstr, xy=(0.9, 0.9), xycoords='axes fraction',
+                    horizontalalignment='right', verticalalignment='top')
+
+        ax.set_ylabel('Density')
+        ax.set_xlabel('Memory (ratio)')
+        ax.grid(True)
+        return ax
 
 
 class NeedFitError(Exception):

snakebacon/mcmcbackends/__init__.py

@@ -64,7 +64,7 @@ def prior_sediment_rate(*args, **kwargs):
         # TODO(brews): Check that these stats are correctly translated to scipy.stats distribs.
         acc_mean = kwargs['acc_mean']
         acc_shape = kwargs['acc_shape']
-        x = np.linspace(0, 3 * np.max(acc_mean), 100)
+        x = np.linspace(0, 6 * np.max(acc_mean), 100)
         y = stats.gamma.pdf(x, a=acc_shape,
                             scale=1 / (acc_shape/acc_mean))
         return y, x

snakebacon/records.py

@@ -24,8 +24,8 @@ def read_14c(fl):
     return outcurve
 
 
-def read_dates(fl):
-    """Create proxy instance from Bacon proxy file
+def read_chron(fl):
+    """Create ChronRecord instance from Bacon file
     """
     indata = pd.read_table(fl, sep=r'\s*\,\s*', index_col=None, engine='python')
     outcore = ChronRecord(age=indata['age'],

snakebacon/tests/test_agedepth.py

@@ -5,7 +5,7 @@
 import numpy as np
 import pandas as pd
 
-from snakebacon import read_dates, ProxyRecord
+from snakebacon import read_chron, ProxyRecord
 from snakebacon.agedepth import AgeDepthModel
 
 
@@ -16,12 +16,13 @@
                    d_r=[0], d_std=[0], t_a=[3], t_b=[4], k=20,
                    minyr=-1000, maxyr=1e6, th01=4147, th02=4145,
                    acc_mean=20, acc_shape=1.5, mem_strength=4, mem_mean=0.7)
-fullrun_agemodel = AgeDepthModel(read_dates(path.join(here, 'MSB2K.csv')),
+fullrun_agemodel = AgeDepthModel(read_chron(path.join(here, 'MSB2K.csv')),
                                  mcmc_kwargs=mcmc_kwargs)
 
 
 class TestAgeDepth(unittest.TestCase):
     def setUp(self):
+        np.random.seed(123)
         self.testdummy = deepcopy(fullrun_agemodel)
 
     def test_init(self):
@@ -78,6 +79,29 @@ def test_agedepth(self):
         np.testing.assert_allclose(victim.mean(), goal_mean, atol=10)
         np.testing.assert_allclose(victim.var(), goal_var, atol=550)
 
+    def test_prior_sediment_memory(self):
+        goal_mean = 0.98457848264590286
+        goal_std = 0.71613816177236256
+        goal_n = 100
+
+        victim, x = self.testdummy.prior_sediment_memory()
+
+        np.testing.assert_equal(len(victim), goal_n)
+        np.testing.assert_allclose(victim.mean(), goal_mean, atol=1e-3)
+        np.testing.assert_allclose(victim.std(), goal_std, atol=1e-3)
+
+
+    def test_prior_sediment_rate(self):
+        goal_mean = 0.008194080517375957
+        goal_std = 0.01172658825323754
+        goal_n = 100
+
+        victim, x = self.testdummy.prior_sediment_rate()
+
+        np.testing.assert_equal(len(victim), goal_n)
+        np.testing.assert_allclose(victim.mean(), goal_mean, atol=1e-3)
+        np.testing.assert_allclose(victim.std(), goal_std, atol=1e-3)
+
 
 if __name__ == '__main__':
     unittest.main()

snakebacon/tests/test_bacon_utils.py

@@ -5,8 +5,7 @@
 
 from snakebacon.mcmcbackends.bacon import fetch_calibcurve
 from snakebacon.mcmcbackends.bacon.utils import d_cal, calibrate_dates
-# from snakebacon.mcmcbackends import Bacon
-from snakebacon import read_dates
+from snakebacon import read_chron
 import snakebacon.records as curve
 
 
@@ -60,7 +59,7 @@ def test_calibrate_dates(self):
         dgoal_n = 40
         pgoal_n = dgoal_n
 
-        chron = read_dates(os.path.join(here, 'MSB2K.csv'))
+        chron = read_chron(os.path.join(here, 'MSB2K.csv'))
 
         d_r = [0]
         d_std = [0]

snakebacon/tests/test_baconwrap.py

@@ -35,7 +35,7 @@ def test__baconin_str(self):
 
     def test_run_baconmcmc(self):
         testcore_path = path.join(here, 'MSB2K.csv')
-        c = snek.read_dates(testcore_path)
+        c = snek.read_chron(testcore_path)
         fullrun_victim = baconwrap.run_baconmcmc(core_labid=c.labid, core_age=c.age, core_error=c.error,
                                                  core_depth=c.depth, depth_min=1.5, depth_max=99.5, cc=[1],
                                                  cc1='IntCal13', cc2='Marine13', cc3='SHCal13', cc4='ConstCal',

snakebacon/tests/test_mcmc.py

@@ -4,7 +4,7 @@
 
 import numpy as np
 
-from snakebacon import read_dates
+from snakebacon import read_chron
 from snakebacon.mcmc import McmcResults, McmcSetup
 
 
@@ -16,7 +16,7 @@
                    d_r=[0], d_std=[0], t_a=[3], t_b=[4], k=20,
                    minyr=-1000, maxyr=1e6, th01=4147, th02=4145,
                    acc_mean=20, acc_shape=1.5, mem_strength=4, mem_mean=0.7)
-fullrun_setup = McmcSetup(read_dates(path.join(here, 'MSB2K.csv')), **mcmc_kwargs)
+fullrun_setup = McmcSetup(read_chron(path.join(here, 'MSB2K.csv')), **mcmc_kwargs)
 fullrun_victim = McmcResults(fullrun_setup)
 
 

snakebacon/tests/test_mcmcbackend_bacon.py

@@ -4,7 +4,7 @@
 import numpy as np
 
 from snakebacon.mcmcbackends import Bacon
-from snakebacon import read_dates
+from snakebacon import read_chron
 
 
 here = os.path.abspath(os.path.dirname(__file__))
@@ -30,7 +30,7 @@ def test_prior_dates(self):
         dgoal_n = 40
         pgoal_n = dgoal_n
 
-        chron = read_dates(os.path.join(here, 'MSB2K.csv'))
+        chron = read_chron(os.path.join(here, 'MSB2K.csv'))
         d_target, p_target = Bacon.prior_dates(chron, **self.mcmc_kwargs)
         np.testing.assert_allclose(d_target[-3:], dgoal)
         np.testing.assert_equal(len(d_target), dgoal_n)
@@ -45,7 +45,7 @@ def test_prior_dates(self):
     def test_prior_sediment_rate(self):
         np.random.seed(123)
 
-        goal_mean = 0.015989306883717701
+        goal_mean = 0.008194080517375957
         goal_std = 0.01172658825323754
         goal_n = 100
 
@@ -68,31 +68,6 @@ def test_prior_sediment_memory(self):
         np.testing.assert_allclose(victim.mean(), goal_mean, atol=1e-3)
         np.testing.assert_allclose(victim.std(), goal_std, atol=1e-3)
 
-    def test_calibrate_dates(self):
-        pgoal_age_mean = 6702.5
-        pgoal_age_std = 121.23496470353207
-
-        pgoal_density_mean = 0.011642671907680679
-        pgoal_density_std = 0.010893804880746285
-
-        dgoal = np.array([77.5, 79.5, 99.5])
-        dgoal_n = 40
-        pgoal_n = dgoal_n
-
-        chron = read_dates(os.path.join(here, 'MSB2K.csv'))
-
-        d_target, p_target = Bacon.prior_dates(chron, **self.mcmc_kwargs)
-        np.testing.assert_allclose(d_target[-3:], dgoal)
-        np.testing.assert_equal(len(d_target), dgoal_n)
-
-        np.testing.assert_allclose(np.mean(p_target[-1][:, 0]), pgoal_age_mean, atol=7)
-        np.testing.assert_allclose(np.std(p_target[-1][:, 0]), pgoal_age_std, atol=2)
-        np.testing.assert_equal(len(p_target), pgoal_n)
-
-
-        np.testing.assert_allclose(np.mean(p_target[-1][:, 1]), pgoal_density_mean, atol=1e-2)
-        np.testing.assert_allclose(np.std(p_target[-1][:, 1]), pgoal_density_std, atol=1e-2)
-
 
 if __name__ == '__main__':
     unittest.main()