Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
8 changes: 6 additions & 2 deletions CHANGELOG.md
Original file line number Diff line number Diff line change
Expand Up @@ -15,7 +15,7 @@ The rules for this file:
* YYYY-MM-DD date format (following ISO 8601)
* accompany each entry with github issue/PR number (Issue #xyz)
-->
## [1.1.3] - UNRELEASES
## [1.1.3] - 2025-09-11

### Authors
* @orbeckst
Expand All @@ -28,11 +28,15 @@ The rules for this file:
(Issue #48)

### Changed
The following changes do not functionally change the code but users relying
on default values need to be aware of the changes to gibbs.Gibbs.

* Default kwargs for the skipping in the Gibbs sampler are now
gibbs.Gibbs(g=100, gskip=1) (used to be g=50, gskip=2) but for most users
gskip for processing data is not important and it makes more sense to focus
on g as the stride at which we sample AND process data (#48, PR #49)

* internal code refactor and clean-up: moved util.run_residue() to
gibbs.run_residue() and use direct assignment instead of setattr() (PR #50)

## [1.1.2] - 2025-07-22

Expand Down
16 changes: 8 additions & 8 deletions basicrta/cluster.py
Original file line number Diff line number Diff line change
Expand Up @@ -119,10 +119,10 @@ def reprocess(self, nproc=1):

taus = np.array(taus)
bars = get_bars(taus)
setattr(self, 'taus', taus[:, 1])
setattr(self, 'bars', bars)
setattr(self, 'residues', np.array(residues))
setattr(self, 'files', np.array(results))
self.taus = taus[:, 1]
self.bars = bars
self.residues = np.array(residues)
self.files = np.array(results)

def get_taus(self, nproc=1):
r"""Get :math:`\tau` and 95\% confidence interval bounds for the slowest
Expand Down Expand Up @@ -158,10 +158,10 @@ def get_taus(self, nproc=1):

taus = np.array(taus)
bars = get_bars(taus)
setattr(self, 'taus', taus[:, 1])
setattr(self, 'bars', bars)
setattr(self, 'residues', np.array(residues))
setattr(self, 'files', np.array(results))
self.taus = taus[:, 1]
self.bars = bars
self.residues = np.array(residues)
self.files = np.array(results)
return taus[:, 1], bars

def write_data(self, fname='tausout'):
Expand Down
50 changes: 40 additions & 10 deletions basicrta/gibbs.py
Original file line number Diff line number Diff line change
Expand Up @@ -8,6 +8,7 @@
from tqdm import tqdm
from MDAnalysis.analysis.base import Results
from basicrta.util import confidence_interval
import multiprocessing
from multiprocessing import Pool, Lock
import MDAnalysis as mda
from basicrta import istarmap
Expand Down Expand Up @@ -47,8 +48,6 @@ def run(self, run_resids=None, g=100):
: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:
contacts = pickle.load(f)

Expand Down Expand Up @@ -96,6 +95,38 @@ def run(self, run_resids=None, g=100):
pass


def run_residue(residue, time, proc, ncomp, niter, cutoff, g, from_combined=False):
"""Run Gibbs sampler for a single residue.

:param residue: Residue name
:type residue: str
:param time: Residence times data
:type time: array-like
:param proc: Process number for progress bar positioning
:type proc: int
:param ncomp: Number of mixture components
:type ncomp: int
:param niter: Number of iterations
:type niter: int
:param cutoff: Cutoff value used in contact analysis
:type cutoff: float
:param g: Gibbs skip parameter
:type g: int
:param from_combined: Whether data comes from combined contacts
:type from_combined: bool
"""
x = np.array(time)
if len(x) != 0:
try:
proc = int(multiprocessing.current_process().name.split('-')[-1])
except ValueError:
proc = 1

gib = Gibbs(times=x, residue=residue, loc=proc, ncomp=ncomp, niter=niter, cutoff=cutoff, g=g)
gib._from_combined_contacts = from_combined
gib.run()


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
Expand Down Expand Up @@ -301,8 +332,8 @@ def cluster(self, method="GaussianMixture", **kwargs):
pindicator[tmpind, mapinds[i]] += 1

pindicator = (pindicator.T / pindicator.sum(axis=1)).T
setattr(self.processed_results, 'indicator', pindicator)
setattr(self.processed_results, 'labels', all_labels)
self.processed_results.indicator = pindicator
self.processed_results.labels = all_labels

def process_gibbs(self, show=True):
r"""
Expand All @@ -329,9 +360,8 @@ def process_gibbs(self, show=True):

self.cluster(n_init=117, n_components=lmode)
labels, presorts = mixture_and_plot(self, show=show)
setattr(self.processed_results, 'labels', labels)
setattr(self.processed_results, 'indicator',
self.processed_results.indicator[:, presorts])
self.processed_results.labels = labels
self.processed_results.indicator = self.processed_results.indicator[:, presorts]

attrs = ["weights", "rates", "ncomp", "residue", "iteration", "niter"]
values = [fweights, frates, lmode, self.residue, indices, self.niter]
Expand Down Expand Up @@ -364,7 +394,7 @@ def _sample_indicator(self):
# sample indicator
s = np.argmax(rng.multinomial(1, z), axis=1)
indicator[i] = s
setattr(self, 'indicator', indicator)
self.indicator = indicator
return indicator[burnin_ind::self.gskip]

def save(self):
Expand Down Expand Up @@ -719,8 +749,8 @@ def _estimate_params(self):
params = np.array([[wh[1][np.argmax(wh[0])], rh[1][np.argmax(rh[0])]]
for wh, rh in zip(whists, rhists)])

setattr(rp, 'parameters', params)
setattr(rp, 'intervals', np.array([wbounds, rbounds]))
rp.parameters = params
rp.intervals = np.array([wbounds, rbounds])

def estimate_tau(self):
r"""
Expand Down
13 changes: 0 additions & 13 deletions basicrta/util.py
Original file line number Diff line number Diff line change
Expand Up @@ -472,19 +472,6 @@ def plot_protein(residues, t_slow, bars, prot=None, label_cutoff=3, ylim=None,
# gib.run()


def run_residue(residue, time, proc, ncomp, niter, cutoff, g, from_combined=False):
from basicrta.gibbs import Gibbs
x = np.array(time)
if len(x) != 0:
try:
proc = int(multiprocessing.current_process().name.split('-')[-1])
except ValueError:
proc = 1

gib = Gibbs(x, residue, proc, ncomp=ncomp, niter=niter, cutoff=cutoff)
gib._from_combined_contacts = from_combined
gib.run()


def check_results(residues, times, ts):
if not os.path.exists('result_check'):
Expand Down
Loading