.. _demo_bouts-label:

===============
 Bout analysis
===============

Here is a brief demo on bout analysis with the :py:mod:`bouts` module for
data generated by mixtures of random Poisson processes.

Set up the environment.  Consider loading the :py:mod:`logging` module and
setting up a logger to monitor progress to this section.

.. jupyter-execute::

   # Set up
   import importlib.resources as rsrc
   import pandas as pd
   import matplotlib.pyplot as plt
   from skdiveMove import calibrate
   import skdiveMove.bouts as skbouts

   # Declare figure sizes
   _FIG1X1 = (7, 6)
   _FIG1X2 = (12, 5)
   _FIG3X1 = (11, 11)

.. jupyter-execute::
   :hide-code:
   :hide-output:

   import numpy as np

   pd.set_option("display.precision", 3)
   np.set_printoptions(precision=3, sign="+")
   %matplotlib inline


Calculate postdive duration
===========================

Create a :class:`~skdiveMove.TDR` object to easily calculate the necessary
statistics:

.. jupyter-execute::
   :linenos:

   config_file = (rsrc.files("skdiveMove") / "config_examples" /
                  "ag_mk7_2002_022_config.json")
   tdr_file = (rsrc.files("skdiveMove") / "tests" /
               "data" / "ag_mk7_2002_022.nc")
   tdrX = calibrate(tdr_file, config_file)
   stats = tdrX.dive_stats()
   stamps = tdrX.stamp_dives(ignore_z=True)
   stats_tab = pd.concat((stamps, stats), axis=1)
   stats_tab.info()

Extract postdive duration for further analysis.

.. jupyter-execute::
   :linenos:

   postdives = stats_tab["postdive_dur"][stats_tab["phase_id"] == 4]
   postdives_diff = (postdives.dt.total_seconds()
                     .diff().iloc[1:].abs())
   # Remove isolated dives
   postdives_diff = postdives_diff[postdives_diff < 2000]


Non-linear least squares via "broken-stick" model
=================================================

`skdiveMove` provides the :class:`bouts.BoutsNLS` class for fitting
non-linear least squares (NLS) models to a modified histogram of a given
variable.

The first step is to generate a modified histogram of postdive duration,
and this requires choosing the bin width for the histogram.

.. jupyter-execute::
   :linenos:

   postdives_nlsbouts = skbouts.BoutsNLS(postdives_diff, 0.1)
   print(postdives_nlsbouts)


Two-process model
~~~~~~~~~~~~~~~~~

Assuming a 2-process model, calculate starting values, providing a guess at
50 s interdive interval.

.. jupyter-execute::
   :linenos:

   fig, ax = plt.subplots(figsize=_FIG1X1)
   init_pars2 = postdives_nlsbouts.init_pars([50], plot=True, ax=ax)

Fit the two-process model.

.. jupyter-execute::
   :linenos:

   coefs2, pcov2 = postdives_nlsbouts.fit(init_pars2)
   # Coefficients
   print(coefs2)

.. jupyter-execute::
   :linenos:

   # Covariance between parameters
   print(pcov2)

Calculate bout-ending criterion.

.. jupyter-execute::
   :linenos:

   # `bec` returns ndarray, and we have only one here
   print("bec = {[0]:.2f}".format(postdives_nlsbouts.bec(coefs2)))

Plot the fit.

.. jupyter-execute::
   :linenos:

   fig, ax = plt.subplots(figsize=_FIG1X1)
   postdives_nlsbouts.plot_fit(coefs2, ax=ax);


Three-process model
~~~~~~~~~~~~~~~~~~~

Attempt to discern three processes in the data.

.. jupyter-execute::
   :linenos:

   fig, ax = plt.subplots(figsize=_FIG1X1)
   init_pars3 = postdives_nlsbouts.init_pars([50, 550], plot=True, ax=ax)

Fit three-process model.

.. jupyter-execute::
   :linenos:

   coefs3, pcov3 = postdives_nlsbouts.fit(init_pars3)
   # Coefficients
   print(coefs3)

.. jupyter-execute::
   :linenos:

   # Covariance between parameters
   print(pcov3)

Plot the fit.

.. jupyter-execute::
   :linenos:

   fig, ax = plt.subplots(figsize=_FIG1X1)
   postdives_nlsbouts.plot_fit(coefs3, ax=ax);

Compare the cumulative frequency distributions of two- vs three-process
models.

.. jupyter-execute::
   :linenos:

   fig, axs = plt.subplots(1, 2, figsize=_FIG1X2)
   postdives_nlsbouts.plot_ecdf(coefs2, ax=axs[0])
   postdives_nlsbouts.plot_ecdf(coefs3, ax=axs[1]);

The three-process model does not follow the observed data as well as the
two-process model.


Maximum likelihood estimation
=============================

Another way to model Poisson mixtures that does not rely on the
subjectively created histogram, and involves fewer parameters, requires
fitting via maximum likelihood estimation (MLE). This approach is available
in :class:`bouts.BoutsMLE`.

Set up an instance.

.. jupyter-execute::
   :linenos:

   postdives_mlebouts = skbouts.BoutsMLE(postdives_diff, 0.1)
   print(postdives_mlebouts)

Again, assuming a 2-process model, calculate starting values.

.. jupyter-execute::
   :linenos:

   fig, ax = plt.subplots(figsize=_FIG1X1)
   init_pars = postdives_mlebouts.init_pars([50], plot=True, ax=ax)

Fit the two-process model.  It is important, but optional, to supply
reasonable bounds to help the optimization algorithm.  Otherwise, the
algorithm may fail to converge.  The fitting procedure is done in two
steps: with and without a reparameterized log-likelihood function.
Therefore, there are two sets of bounds required.

.. jupyter-execute::
   :linenos:

   p_bnd = (-2, None)                 # bounds for `p`
   lda1_bnd = (-5, None)              # bounds for `lambda1`
   lda2_bnd = (-10, None)             # bounds for `lambda2`
   bnd1 = (p_bnd, lda1_bnd, lda2_bnd)
   p_bnd = (1e-2, None)
   lda1_bnd = (1e-4, None)
   lda2_bnd = (1e-8, None)
   bnd2 = (p_bnd, lda1_bnd, lda2_bnd)
   fit1, fit2 = postdives_mlebouts.fit(init_pars,
                                       fit1_opts=dict(method="L-BFGS-B",
			                              bounds=bnd1),
			               fit2_opts=dict(method="L-BFGS-B",
			  	                      bounds=bnd2))

.. jupyter-execute::
   :linenos:

   # First fit
   print(fit1)

.. jupyter-execute::
   :linenos:

   # Second fit
   print(fit2)

Calculate bout-ending criterion (BEC).

.. jupyter-execute::
   :linenos:

   print("bec = {:.2f}".format(postdives_mlebouts.bec(fit2)))

Plot the fit.

.. jupyter-execute::
   :linenos:

   fig, ax = plt.subplots(figsize=_FIG1X1)
   postdives_mlebouts.plot_fit(fit2, ax=ax);

Compare the cumulative frequency distribution between NLS and MLE model
estimates.

.. jupyter-execute::
   :linenos:

   fig, axs = plt.subplots(1, 2, figsize=_FIG1X2)
   postdives_nlsbouts.plot_ecdf(coefs2, ax=axs[0])
   axs[0].set_title("NLS")
   postdives_mlebouts.plot_ecdf(fit2, ax=axs[1])
   axs[1].set_title("MLE");

Label bouts based on BEC from the last MLE model.  Note that `Timedelta`
type needs to be converted to total seconds to allow comparison with BEC.

.. jupyter-execute::
   :linenos:

   bec = postdives_mlebouts.bec(fit2)
   skbouts.label_bouts(postdives.dt.total_seconds(), bec, as_diff=True)

Feel free to download a copy of this demo
(:jupyter-download-script:`demo_bouts`).