Skip to content

Commit

Permalink
added docs
Browse files Browse the repository at this point in the history
  • Loading branch information
rsexton2 committed Nov 22, 2024
1 parent 2f11479 commit ee28090
Show file tree
Hide file tree
Showing 2 changed files with 75 additions and 15 deletions.
51 changes: 36 additions & 15 deletions basicrta/gibbs.py
Original file line number Diff line number Diff line change
@@ -1,6 +1,3 @@
"""Analysis functions
"""

import os
import gc
import pickle
Expand All @@ -22,7 +19,16 @@

class ParallelGibbs(object):
"""
A module to take a contact map and run Gibbs samplers for each residue
A module to take a contact map and run Gibbs samplers for each residue.
:param contacts: Contact pickle file (`contacts-{cutoff}.pkl`).
:type contacts: str
:param nproc: Number of processes to use in running Gibbs samplers.
:type nproc: int
:param ncomp: Number of mixture components to use in the Gibbs sampler.
:type ncomp: int
:param niter: Number of iterations of the Gibbs sampler to perform.
:type niter: int
"""

def __init__(self, contacts, nproc=1, ncomp=15, niter=110000):
Expand All @@ -34,6 +40,13 @@ def __init__(self, contacts, nproc=1, ncomp=15, niter=110000):
self.contacts = contacts

def run(self, run_resids=None):
"""
The :meth:`run` method executes the Gibbs samplers for all residues of
`sel1` present in the contact map, or a list of resids can be provided.
:param run_resids: Resid(s) for which to run a Gibbs sampler.
:type run_resids: int or list, optional
"""
from basicrta.util import run_residue

with open(self.contacts, 'r+b') as f:
Expand Down Expand Up @@ -79,7 +92,7 @@ class Gibbs(object):
r"""Gibbs sampler to estimate parameters of an exponential mixture for a set
of data. Results are stored in :class:`gibbs.results`, which uses
:class:`MDAnalysis.analysis.base.Results()`. If 'results=None' the gibbs
sampler has not been executed, which requires calling :meth:`run`
sampler has not been executed, which requires calling :meth:`run`.
:param times: Set of residence times to analyze
:type times: array, optional
Expand All @@ -91,7 +104,7 @@ class Gibbs(object):
:type ncomp: int
:param niter: Number of iterations to run the Gibbs sampler
:type niter: int
:param cutoff: Cutoff calue used in contact analysis, used to determine
:param cutoff: Cutoff value used in contact analysis, used to determine
directory to load/save results. Allows for multiple cutoffs
to be tested in directory containing contacts.
:type cutoff: float
Expand Down Expand Up @@ -162,7 +175,8 @@ def _prepare(self):

def run(self):
r"""
Execute the Gibbs sampler and save the results to :attr:`Gibbs.results`
Execute the Gibbs sampler and save the raw data to the instance of
:class:`Gibbs`.
"""
# initialize weights and rates
self._prepare()
Expand Down Expand Up @@ -260,7 +274,9 @@ def cluster(self, method="GaussianMixture", **kwargs):

def process_gibbs(self):
r"""
Process the data collected from the Gibbs sampler
Process the samples collected from the Gibbs sampler.
:meth:`process_gibbs` can be called multiple times to check the
robustness of the results.
"""
from basicrta.util import mixture_and_plot
from scipy import stats
Expand Down Expand Up @@ -296,7 +312,7 @@ def result_plot(self, remove_noise=False, **kwargs):
Generate the combined result plot with option to change kwargs without
re-clustering.
:param remove_noise: Whether to remove noise
:param remove_noise: Option to remove noise clusters
:type remove_noise: bool
"""
from basicrta.util import mixture_and_plot
Expand All @@ -319,7 +335,7 @@ def _sample_indicator(self):

def save(self):
"""
Save current state of the Gibbs
Save current state of the :class:`Gibbs` instance.
"""
savedir = f'basicrta-{self.cutoff}/{self.residue}/'
filename = f'gibbs_{self.niter}.pkl'
Expand All @@ -335,9 +351,9 @@ def save(self):
@staticmethod
def load(file):
"""
Load Gibbs
Load an instance of :class:`Gibbs`.
:param file: Filename of Gibbs to load
:param file: Path to instance of :class:`Gibbs`
:type file: str
"""
from basicrta.util import get_s
Expand Down Expand Up @@ -366,7 +382,8 @@ def load(file):

def plot_tau_hist(self, scale=1, save=False):
r"""
Plot histogram of tau values
Plot histogram of tau values. The figure aspect ratio is 4:3, and can be
made larger/smaller using the `scale` argument.
:param scale: Increase plot size by this factor
:type scale: float
Expand Down Expand Up @@ -672,11 +689,15 @@ def _estimate_params(self):
setattr(rp, 'intervals', np.array([wbounds, rbounds]))

def estimate_tau(self):
"""
Estimate the posterior maximum for the tau distribution
r"""
Estimate the posterior maximum and confidence interval (CI) for the
:math:`tau` distribution of the slowest process. NOTE: In the future
this will return an array containing :math:`tau` and CI for all
clusters.
:return: An array containing the posterior maximum and bounds of the
95% confidence interval in the format [LB, max, UB].
:rtype: list
"""
rp = self.processed_results

Expand Down
39 changes: 39 additions & 0 deletions basicrta/kinetics.py
Original file line number Diff line number Diff line change
Expand Up @@ -7,6 +7,18 @@


class MapKinetics(object):
r"""
The MapKinetics class takes processed Gibbs sampler data (an instance of the
Gibbs class) and the contact file to create costomized trajectories,
containing all of `sel1` from the initial contact analysis and a single
`sel2` residue.
:param gibbs: Processed instance of :class:`Gibbs` class
:type gibbs: str
:param contacts: Contact pickle file (`contacts-{cutoff}.pkl`)
:type contacts: str
"""
def __init__(self, gibbs, contacts):
self.gibbs = gibbs
self.cutoff = float(contacts.split('/')[-1].strip('.pkl').
Expand Down Expand Up @@ -65,6 +77,20 @@ def _create_data(self):
j += len(tmp)

def create_traj(self, top_n=None):
r"""
Create the customized trajectories for the individual mixture components
of the model. If `top_n` is None, a single trajectory is created
with all of `sel1` and a single `sel2` residue, with every contact
accounted for (ie. a single frame in the original trajectory may be
used multiple times due to multiple contacts formed with the `sel1`
residue of interest at that frame).
:param top_n: Number of frames desired for the individual trajectories
(sorted in order of decreasing classification
probability).
:type top_n: int
"""

if os.path.exists(self.fulltraj) and top_n is None:
raise FileExistsError(f'{self.fulltraj} exists, remove then rerun')

Expand Down Expand Up @@ -99,6 +125,19 @@ def create_traj(self, top_n=None):
W.write(ag1 + ag2.select_atoms(f'resid {int(tmp[i, 1])}'))

def weighted_densities(self, step=1, top_n=None, filterP=0):
"""
Create weighted densities based on the marginal posterior probabilities
of the model component classification.
:param step: Use every `Nth` frame
:type step: int
:param top_n: Use the `N` most likely frames for each component to
compute the weighted densities from.
:type top_n: int
:param filterP: Probabilities less than `filterP` are not used in the
weighted density calculation
:type filterP: float
"""
if not os.path.exists(self.fulltraj):
self.create_traj()

Expand Down

0 comments on commit ee28090

Please sign in to comment.