Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
Observable vector
Returns
-------
val: float
Equilibrium expectation value of the given observable
Notes
-----
The equilibrium expectation value of an observable a is defined as follows
.. math::
\mathbb{E}_{\mu}[a] = \sum_i \mu_i a_i
:math:`\mu=(\mu_i)` is the stationary vector of the transition matrix :math:`T`.
"""
# check input and go
a = _types.ensure_ndarray(a, ndim=1, size=self.nstates, kind='numeric')
return _np.dot(a, self.stationary_distribution)
def fingerprint_relaxation(self, p0, a, k=None, ncv=None):
# basic checks for a and b
p0 = _types.ensure_ndarray(p0, ndim=1, kind='numeric')
a = _types.ensure_ndarray(a, ndim=1, kind='numeric', size=len(p0))
# are we on microstates space?
if len(a) == self.nstates_obs:
p0 = _np.dot(self.observation_probabilities, p0)
a = _np.dot(self.observation_probabilities, a)
# now we are on macrostate space, or something is wrong
if len(a) == self.nstates:
return _MSM.fingerprint_relaxation(self, p0, a)
else:
raise ValueError('observable vectors have size ' + len(a) + ' which is incompatible with both hidden (' +
self.nstates + ') and observed states (' + self.nstates_obs + ')')
J. Stat. Phys. 123: 503-523 (2006)
.. [2] P. Metzner, C. Schuette and E. Vanden-Eijnden.
Transition Path Theory for Markov Jump Processes.
Multiscale Model Simul 7: 1192-1219 (2009)
.. [3] F. Noe, Ch. Schuette, E. Vanden-Eijnden, L. Reich and T. Weikl:
Constructing the Full Ensemble of Folding Pathways from Short Off-Equilibrium Simulations.
Proc. Natl. Acad. Sci. USA, 106, 19011-19016 (2009)
"""
from msmtools.flux import flux_matrix, to_netflux
import msmtools.analysis as msmana
T = msmobj.transition_matrix
mu = msmobj.stationary_distribution
A = _types.ensure_ndarray(A, kind='i')
B = _types.ensure_ndarray(B, kind='i')
if len(A) == 0 or len(B) == 0:
raise ValueError('set A or B is empty')
n = T.shape[0]
if len(A) > n or len(B) > n or max(A) > n or max(B) > n:
raise ValueError('set A or B defines more states, than given transition matrix.')
# forward committor
qplus = msmana.committor(T, A, B, forward=True)
# backward committor
if msmana.is_reversible(T, mu=mu):
qminus = 1.0 - qplus
else:
qminus = msmana.committor(T, A, B, forward=False, mu=mu)
# gross flux
grossflux = flux_matrix(T, mu, qminus, qplus, netflux=False)
val: float
Equilibrium expectation value fo the given observable
Notes
-----
The equilibrium expectation value of an observable :math:`a` is defined as follows
.. math::
\mathbb{E}_{\mu}[a] = \sum_i \pi_i a_i
:math:`\pi=(\pi_i)` is the stationary vector of the transition matrix :math:`P`.
"""
# check input and go
a = _types.ensure_ndarray(a, ndim=1, size=self.nstates, kind='numeric')
return _np.dot(a, self.stationary_distribution)
>>> wham.meval('stationary_distribution') # doctest: +ELLIPSIS +REPORT_NDIFF
[array([ 0.5..., 0.4...]), array([ 0.6..., 0.3...])]
References
----------
.. [1] Ferrenberg, A.M. and Swensen, R.H. 1988.
New Monte Carlo Technique for Studying Phase Transitions.
Phys. Rev. Lett. 23, 2635--2638
.. [2] Kumar, S. et al 1992.
The Weighted Histogram Analysis Method for Free-Energy Calculations on Biomolecules. I. The Method.
J. Comp. Chem. 13, 1011--1021
"""
self.bias_energies_full = _types.ensure_ndarray(bias_energies_full, ndim=2, kind='numeric')
self.stride = stride
self.dt_traj = dt_traj
self.maxiter = maxiter
self.maxerr = maxerr
self.save_convergence_info = save_convergence_info
# set derived quantities
self.nthermo, self.nstates_full = bias_energies_full.shape
# set iteration variables
self.therm_energies = None
self.conf_energies = None
def relaxation(self, p0, a, maxtime=None, k=None, ncv=None):
# basic checks for a and b
p0 = _types.ensure_ndarray(p0, ndim=1, kind='numeric')
a = _types.ensure_ndarray(a, ndim=1, kind='numeric', size=len(p0))
# are we on microstates space?
if len(a) == self.nstates_obs:
p0 = _np.dot(self.observation_probabilities, p0)
a = _np.dot(self.observation_probabilities, a)
# now we are on macrostate space, or something is wrong
if len(a) == self.nstates:
return _MSM.relaxation(self, p0, a, maxtime=maxtime)
else:
raise ValueError('observable vectors have size ' + len(a) + ' which is incompatible with both hidden (' +
self.nstates + ') and observed states (' + self.nstates_obs + ')')
def test_estimator(self, test_estimator):
self._test_estimator = test_estimator
self.active_set = types.ensure_ndarray(np.array(test_estimator.active_set), kind='i') # create a copy
# map from the full set (here defined by the largest state index in active set) to active
self._full2active = np.zeros(np.max(self.active_set)+1, dtype=int)
self._full2active[self.active_set] = np.arange(self.nstates)
def fingerprint_relaxation(self, p0, a, k=None, ncv=None):
# basic checks for a and b
p0 = _types.ensure_ndarray(p0, ndim=1, kind='numeric')
a = _types.ensure_ndarray(a, ndim=1, kind='numeric', size=len(p0))
# are we on microstates space?
if len(a) == self.nstates_obs:
p0 = _np.dot(self.observation_probabilities, p0)
a = _np.dot(self.observation_probabilities, a)
# now we are on macrostate space, or something is wrong
if len(a) == self.nstates:
return _MSM.fingerprint_relaxation(self, p0, a)
else:
raise ValueError('observable vectors have size %s which is incompatible with both hidden (%s)'
' and observed states (%s)' % (len(a), self.nstates, self.nstates_obs))
conf : float, optional, default = 0.95
confidence interval
Return
------
lower : ndarray(shape)
element-wise lower bounds
upper : ndarray(shape)
element-wise upper bounds
"""
if conf < 0 or conf > 1:
raise ValueError('Not a meaningful confidence level: '+str(conf))
try:
data = types.ensure_ndarray(data, kind='numeric')
except:
# if 1D array of arrays try to fuse it
if isinstance(data, np.ndarray) and np.ndim(data) == 1:
newshape = tuple([len(data)] + list(data[0].shape))
newdata = np.zeros(newshape)
for i in range(len(data)):
newdata[i, :] = data[i]
data = newdata
types.assert_array(data, kind='numeric')
if np.ndim(data) == 1:
m, lower, upper = _confidence_interval_1d(data, conf)
return lower, upper
else:
I = _indexes(data[0])