Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
"""
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])
lower = np.zeros(data[0].shape)
upper = np.zeros(data[0].shape)
for i in I:
col = _column(data, i)
m, lower[i], upper[i] = _confidence_interval_1d(col, conf)
# return
return lower, upper
def _estimate(self, X):
ttrajs, dtrajs_full, btrajs = X
# shape and type checks
assert len(ttrajs) == len(dtrajs_full) == len(btrajs)
for t in ttrajs:
_types.assert_array(t, ndim=1, kind='i')
for d in dtrajs_full:
_types.assert_array(d, ndim=1, kind='i')
for b in btrajs:
_types.assert_array(b, ndim=2, kind='f')
# find dimensions
self.nstates_full = max(_np.max(d) for d in dtrajs_full) + 1
self.nthermo = max(_np.max(t) for t in ttrajs) + 1
# dimensionality checks
for t, d, b, in zip(ttrajs, dtrajs_full, btrajs):
assert t.shape[0] == d.shape[0] == b.shape[0]
assert b.shape[1] == self.nthermo
# cast types and change axis order if needed
ttrajs = [_np.require(t, dtype=_np.intc, requirements='C') for t in ttrajs]
dtrajs_full = [_np.require(d, dtype=_np.intc, requirements='C') for d in dtrajs_full]
btrajs = [_np.require(b, dtype=_np.float64, requirements='C') for b in btrajs]
# find state visits
self.state_counts_full = _util.state_counts(ttrajs, dtrajs_full)
self.therm_state_counts_full = self.state_counts_full.sum(axis=1)
def set_model_params(self, models=None, f_therm=None, pi=None, f=None, label='ground state'):
# don't normalize f, because in a multiensemble the relative energy levels matter
_StationaryModel.set_model_params(self, pi=pi, f=f, normalize_f=False)
# check and set other parameters
_types.assert_array(f_therm, ndim=1, kind='numeric')
f_therm = _np.array(f_therm, dtype=float)
for m in models:
assert issubclass(m.__class__, _Model)
self.update_model_params(models=models, f_therm=f_therm)
def _estimate(self, trajs):
# check input
assert isinstance(trajs, (tuple, list))
assert len(trajs) == 2
ttrajs = trajs[0]
dtrajs = trajs[1]
# validate input
for ttraj, dtraj in zip(ttrajs, dtrajs):
_types.assert_array(ttraj, ndim=1, kind='numeric')
_types.assert_array(dtraj, ndim=1, kind='numeric')
assert _np.shape(ttraj)[0] == _np.shape(dtraj)[0]
# harvest state counts
self.state_counts_full = _util.state_counts(
ttrajs, dtrajs, nthermo=self.nthermo, nstates=self.nstates_full)
# active set
self.active_set = _np.where(self.state_counts_full.sum(axis=0) > 0)[0]
self.state_counts = _np.ascontiguousarray(
self.state_counts_full[:, self.active_set].astype(_np.intc))
self.bias_energies = _np.ascontiguousarray(
self.bias_energies_full[:, self.active_set], dtype=_np.float64)
# run estimator
pg = _ProgressReporter()
def _estimate(self, X):
ttrajs, dtrajs_full, btrajs = X
# shape and type checks
assert len(ttrajs) == len(dtrajs_full) == len(btrajs)
for t in ttrajs:
_types.assert_array(t, ndim=1, kind='i')
for d in dtrajs_full:
_types.assert_array(d, ndim=1, kind='i')
for b in btrajs:
_types.assert_array(b, ndim=2, kind='f')
# find dimensions
nstates_full = max(_np.max(d) for d in dtrajs_full) + 1
if self.nstates_full is None:
self.nstates_full = nstates_full
elif self.nstates_full < nstates_full:
raise RuntimeError("Found more states (%d) than specified by nstates_full (%d)" % (
nstates_full, self.nstates_full))
self.nthermo = max(_np.max(t) for t in ttrajs) + 1
# dimensionality checks
for t, d, b, in zip(ttrajs, dtrajs_full, btrajs):
assert t.shape[0] == d.shape[0] == b.shape[0]
assert b.shape[1] == self.nthermo
# cast types and change axis order if needed
ttrajs = [_np.require(t, dtype=_np.intc, requirements='C') for t in ttrajs]
dtrajs_full = [_np.require(d, dtype=_np.intc, requirements='C') for d in dtrajs_full]
def _estimate(self, trajs):
# check input
assert isinstance(trajs, (tuple, list))
assert len(trajs) == 2
ttrajs = trajs[0]
dtrajs = trajs[1]
# validate input
for ttraj, dtraj in zip(ttrajs, dtrajs):
_types.assert_array(ttraj, ndim=1, kind='numeric')
_types.assert_array(dtraj, ndim=1, kind='numeric')
assert _np.shape(ttraj)[0] == _np.shape(dtraj)[0]
# harvest state counts
self.state_counts_full = _util.state_counts(
ttrajs, dtrajs, nthermo=self.nthermo, nstates=self.nstates_full)
# active set
self.active_set = _np.where(self.state_counts_full.sum(axis=0) > 0)[0]
self.state_counts = _np.ascontiguousarray(
self.state_counts_full[:, self.active_set].astype(_np.intc))
self.bias_energies = _np.ascontiguousarray(
self.bias_energies_full[:, self.active_set], dtype=_np.float64)
# run estimator
pg = _ProgressReporter()
stage = 'WHAM'
normalize_energy : bool, default=True
If parametrized by free energy f, normalize them such that
:math:`\sum_i \pi_i = 1`, which is achieved by :math:`\log \sum_i \exp(-f_i) = 0`.
label : str, default=None
Human-readable description for the thermodynamic state of this
model. May contain a temperature description, such as '300 K' or
a description of bias energy such as 'unbiased' or 'Umbrella 1'.
"""
if f is not None:
_types.assert_array(f, ndim=1, kind='numeric')
f = _np.array(f, dtype=float)
if normalize_f:
f += _logsumexp(-f) # normalize on the level on energies to achieve sum_i pi_i = 1
pi = _np.exp(-f)
elif pi is not None: # if f is not given, use pi. pi can't be None at this point
_types.assert_array(pi, ndim=1, kind='numeric')
pi = _np.array(pi, dtype=float)
f = -_np.log(pi)
f += _logsumexp(-f) # always shift f when set by pi
else:
raise ValueError(
"Trying to initialize model without parameters: both pi (stationary distribution)" \
" and f (free energy) are None. At least one of them needs to be set.")
# set parameters (None does not overwrite)
self.update_model_params(pi=pi, f=f, normalize_energy=normalize_f, label=label)