Python MSMLib.get_counts_from_traj方法代码示例

# 需要导入模块: from msmbuilder import MSMLib [as 别名]
# 或者: from msmbuilder.MSMLib import get_counts_from_traj [as 别名]
def ndgrid_msm_likelihood_score(estimator, sequences):
    """Log-likelihood score function for an (NDGrid, MarkovStateModel) pipeline

    Parameters
    ----------
    estimator : sklearn.pipeline.Pipeline
        A pipeline estimator containing an NDGrid followed by a MarkovStateModel
    sequences: list of array-like, each of shape (n_samples_i, n_features)
        Data sequences, where n_samples_i in the number of samples
        in sequence i and n_features is the number of features.

    Returns
    -------
    log_likelihood : float
        Mean log-likelihood per data point.

    Examples
    --------
    >>> pipeline = Pipeline([
    >>>    ('grid', NDGrid()),
    >>>    ('msm', MarkovStateModel())
    >>> ])
    >>> grid = GridSearchCV(pipeline, param_grid={
    >>>    'grid__n_bins_per_feature': [10, 20, 30, 40]
    >>> }, scoring=ndgrid_msm_likelihood_score)
    >>> grid.fit(dataset)
    >>> print grid.grid_scores_

    References
    ----------
    .. [1] McGibbon, R. T., C. R. Schwantes, and V. S. Pande. "Statistical
       Model Selection for Markov Models of Biomolecular Dynamics." J. Phys.
       Chem B. (2014)
    """
    import msmbuilder.MSMLib as msmlib
    from mixtape import cluster

    grid = [model for (name, model) in estimator.steps if isinstance(model, cluster.NDGrid)][0]
    msm = [model for (name, model) in estimator.steps if isinstance(model, MarkovStateModel)][0]

    # NDGrid supports min/max being different along different directions, which
    # means that the bin widths are coordinate dependent. But I haven't
    # implemented that because I've only been using this for 1D data
    if grid.n_features != 1:
        raise NotImplementedError("file an issue on github :)")

    transition_log_likelihood = 0
    emission_log_likelihood = 0
    logtransmat = np.nan_to_num(np.log(np.asarray(msm.transmat_.todense())))
    width = grid.grid[0, 1] - grid.grid[0, 0]

    for X in grid.transform(sequences):
        counts = np.asarray(
            _apply_mapping_to_matrix(msmlib.get_counts_from_traj(X, n_states=grid.n_bins), msm.mapping_).todense()
        )
        transition_log_likelihood += np.multiply(counts, logtransmat).sum()
        emission_log_likelihood += -1 * np.log(width) * len(X)

    return (transition_log_likelihood + emission_log_likelihood) / sum(len(x) for x in sequences)

# 需要导入模块: from msmbuilder import MSMLib [as 别名]
# 或者: from msmbuilder.MSMLib import get_counts_from_traj [as 别名]
def build_classic_from_memberships(memberships, lag_time=1):
    """Build a classic msm by turning a membership array into a state list.

    This function uses msmbuilder code to calculate the count matrix. Use this
    for compairing quantized versions of the fuzzy count matrix building
    for consistency.
    """
    states = np.zeros(memberships.shape[0], dtype='int')
    n_states = memberships.shape[1]

    for i in xrange(memberships.shape[0]):
        memb = memberships[i]
        state = np.argmax(memb)
        states[i] = state

    counts = msm.get_counts_from_traj(states, n_states, lag_time)
    rev_counts, t_matrix, populations, mapping = msm.build_msm(counts)
    return rev_counts, t_matrix, populations, mapping

# 需要导入模块: from msmbuilder import MSMLib [as 别名]
# 或者: from msmbuilder.MSMLib import get_counts_from_traj [as 别名]
import os
import sys
import numpy as np
import scipy.sparse

from msmbuilder import io, msm_analysis, MSMLib
from bayesmutant import SimpleMutantSampler

P = np.loadtxt('base_transition_matrix.dat')

mutant_transition_matrix = P + 0.2*scipy.sparse.rand(P.shape[0], P.shape[1], density=0.1).todense()
mutant_transition_matrix /= np.sum(mutant_transition_matrix, axis=1)

trajectory =  np.array(msm_analysis.sample(P, 0, 5000))
base_counts = MSMLib.get_counts_from_traj(trajectory).todense()


print 'base counts'
print base_counts

ms = SimpleMutantSampler(base_counts, mutant_transition_matrix)
ms.step(5000)

print 'observed counts'
print ms.counts

io.saveh('sampling.h5', base_counts=base_counts, samples=ms.samples,
                        observed_counts=ms.counts, scores=ms.scores,
                        transition_matrix=mutant_transition_matrix)