The testsystems Module: `pymbar.testsystems`¶

The pymbar.testsystems module contains a number of test systems with analytically or numerically computable expectations or free energies we use to validate its implementation. These test systems are also convenient to use if you want to easily generate synthetic data to experiment with the capabilities of ``pymbar`.

class pymbar.testsystems.harmonic_oscillators.HarmonicOscillatorsTestCase(O_k=(0, 1, 2, 3, 4), K_k=(1, 2, 4, 8, 16), beta=1.0)¶

Test cases using harmonic oscillators.

Examples

Generate energy samples with default parameters.

>>> testcase = HarmonicOscillatorsTestCase()
>>> [x_kn, u_kln, N_k, s_n] = testcase.sample()

Retrieve analytical properties.

>>> analytical_means = testcase.analytical_means()
>>> analytical_variances = testcase.analytical_variances()
>>> analytical_standard_deviations = testcase.analytical_standard_deviations()
>>> analytical_free_energies = testcase.analytical_free_energies()
>>> analytical_x_squared = testcase.analytical_observable('position^2')

Generate energy samples with default parameters in one line.

>>> (x_kn, u_kln, N_k, s_n) = HarmonicOscillatorsTestCase().sample()

Generate energy samples with specified parameters.

>>> testcase = HarmonicOscillatorsTestCase(O_k=[0, 1, 2, 3, 4], K_k=[1, 2, 4, 8, 16])
>>> (x_kn, u_kln, N_k, s_n) = testcase.sample(N_k=[10, 20, 30, 40, 50])

Test sampling in different output modes.

>>> (x_kn, u_kln, N_k) = testcase.sample(N_k=[10, 20, 30, 40, 50], mode='u_kln')
>>> (x_n, u_kn, N_k, s_n) = testcase.sample(N_k=[10, 20, 30, 40, 50], mode='u_kn')

Generate test case with exponential distributions.

Parameters:

O_k : np.ndarray, float, shape=(n_states)

Offset parameters for each state.

K_k : np.ndarray, float, shape=(n_states)

Force constants for each state.

beta : float, optional, default=1.0

Inverse temperature.

Notes

We assume potentials of the form U(x) = (k / 2) * (x - o)^2 Here, k and o are the corresponding entries of O_k and K_k. The equilibrium distribution is given analytically by p(x;beta,K) = sqrt[(beta K) / (2 pi)] exp[-beta K (x-x_0)**2 / 2] The dimensionless free energy is therefore f(beta,K) = - (1/2) * ln[ (2 pi) / (beta K) ]

classmethod evenly_spaced_oscillators(n_states, n_samples_per_state, lower_O_k=1.0, upper_O_k=5.0, lower_k_k=1.0, upper_k_k=3.0)¶

Generate samples from evenly spaced harmonic oscillators.

Parameters:

n_states : np.ndarray, int

number of states

n_samples_per_state : np.ndarray, int

number of samples per state. The total number of samples n_samples will be equal to n_states * n_samples_per_state

lower_O_k : float, optional, default=1.0

Lower bound of O_k values

upper_O_k : float, optional, default=5.0

Upper bound of O_k values

lower_k_k : float, optional, default=1.0

Lower bound of O_k values

upper_k_k : float, optional, default=3.0

Upper bound of k_k values

Returns:

name: str

Name of testsystem

testsystem : TestSystem

The testsystem object

x_n : np.ndarray, shape=(n_samples)

Coordinates of the samples

u_kn : np.ndarray, shape=(n_states, n_samples)

Reduced potential energies

N_k : np.ndarray, shape=(n_states)

Number of samples drawn from each state

s_n : np.ndarray, shape=(n_samples)

State of origin of each sample

sample(N_k=[10, 20, 30, 40, 50], mode='u_kn', seed=None)¶

Draw samples from the distribution.

Parameters:

N_k : np.ndarray, int

number of samples per state

mode : str, optional, default=’u_kn’

If ‘u_kln’, return K x K x N_max matrix where u_kln[k,l,n] is reduced potential of sample n from state k evaluated at state l. If ‘u_kn’, return K x N_tot matrix where u_kn[k,n] is reduced potential of sample n (in concatenated indexing) evaluated at state k.

seed: int, optional, default=None. Provides control over the random seed for replicability.

Returns:

if mode == ‘u_kn’:

x_n : np.ndarray, shape=(n_states*n_samples), dtype=float

x_n[n] is sample n (in concatenated indexing)

u_kn : np.ndarray, shape=(n_states, n_states*n_samples), dtype=float

u_kn[k,n] is reduced potential of sample n (in concatenated indexing) evaluated at state k.

N_k : np.ndarray, shape=(n_states), dtype=float

N_k[k] is the number of samples generated from state k

s_n : np.ndarray, shape=(n_samples), dtype=’int’

s_n is the state of origin of x_n

x_kn : np.ndarray, shape=(n_states, n_samples), dtype=float

1D harmonic oscillator positions

u_kln : np.ndarray, shape=(n_states, n_states, n_samples), dytpe=float, only if mode=’u_kln’

u_kln[k,l,n] is reduced potential of sample n from state k evaluated at state l.

N_k : np.ndarray, shape=(n_states), dtype=int32

N_k[k] is the number of samples generated from state k

pymbar.testsystems.timeseries.correlated_timeseries_example(N=10000, tau=5.0, seed=None)¶

Generate synthetic timeseries data with known correlation time.

Parameters:

N : int, optional

length (in number of samples) of timeseries to generate

tau : float, optional

correlation time (in number of samples) for timeseries

seed : int, optional

If not None, specify the numpy random number seed.

Returns:

dih : np.ndarray, shape=(num_dihedrals), dtype=float

dih[i,j] gives the dihedral angle at traj[i] correponding to indices[j].

Notes

Synthetic timeseries generated using bivariate Gaussian process described by Janke (Eq. 41 of Ref. [1]).

As noted in Eq. 45-46 of Ref. [1], the true integrated autocorrelation time will be given by tau_int = (1/2) coth(1 / 2 tau) = (1/2) (1+rho)/(1-rho) which, for tau >> 1, is approximated by tau_int = tau + 1/(12 tau) + O(1/tau^3) So for tau >> 1, tau_int is approximately the given exponential tau.

References

[R11]

Janke W. Statistical analysis of simulations: Data correlations and error estimation. In ‘Quantum Simulations of Complex Many-Body Systems: From Theory to Algorithms’. NIC Series, VOl. 10, pages 423-445, 2002.

Examples

Generate a timeseries of length 10000 with correlation time of 10.

>>> A_t = correlated_timeseries_example(N=10000, tau=10.0)

Generate an uncorrelated timeseries of length 1000.

>>> A_t = correlated_timeseries_example(N=1000, tau=1.0)

Generate a correlated timeseries with correlation time longer than the length.

>>> A_t = correlated_timeseries_example(N=1000, tau=2000.0)

class pymbar.testsystems.exponential_distributions.ExponentialTestCase(rates=[1, 2, 3, 4, 5], beta=1.0)¶

Test cases using exponential distributions.

Examples

Generate energy samples with default parameters.

>>> testcase = ExponentialTestCase()
>>> [x_kn, u_kln, N_k] = testcase.sample()

Retrieve analytical properties.

>>> analytical_means = testcase.analytical_means()
>>> analytical_variances = testcase.analytical_variances()
>>> analytical_standard_deviations = testcase.analytical_standard_deviations()
>>> analytical_free_energies = testcase.analytical_free_energies()
>>> analytical_x_squared = testcase.analytical_x_squared()

Generate energy samples with default parameters in one line.

>>> [x_kn, u_kln, N_k] = ExponentialTestCase().sample()

Generate energy samples with specified parameters.

>>> testcase = ExponentialTestCase(rates=[1., 2., 3., 4., 5.])
>>> [x_kn, u_kln, N_k] = testcase.sample(N_k=[10, 20, 30, 40, 50])

Test sampling in different output modes.

>>> [x_kn, u_kln, N_k] = testcase.sample(N_k=[10, 20, 30, 40, 50], mode='u_kln')
>>> [x_n, u_kn, N_k, s_n] = testcase.sample(N_k=[10, 20, 30, 40, 50], mode='u_kn')

Generate test case with exponential distributions.

Parameters:

rates : np.ndarray, float, shape=(n_states)

Rate parameters (e.g. lambda) for each state.

beta : float, optional, default=1.0

Inverse temperature.

Notes

We assume potentials of the form U(x) = lambda x.

analytical_free_energies()¶: Return the FE: -log(Z)

classmethod evenly_spaced_exponentials(n_states, n_samples_per_state, lower_rate=1.0, upper_rate=3.0)¶

Generate samples from evenly spaced exponential distributions.

Parameters:

n_states : np.ndarray, int

number of states

n_samples_per_state : np.ndarray, int

number of samples per state. The total number of samples n_samples will be equal to n_states * n_samples_per_state

lower_O_k : float, optional, default=1.0

Lower bound of O_k values

upper_O_k : float, optional, default=5.0

Upper bound of O_k values

lower_k_k : float, optional, default=1.0

Lower bound of O_k values

upper_k_k : float, optional, default=3.0

Upper bound of k_k values

Returns:

name: str

Name of testsystem

testsystem : TestSystem

The testsystem object

x_n : np.ndarray, shape=(n_samples)

Coordinates of the samples

u_kn : np.ndarray, shape=(n_states, n_samples)

Reduced potential energies

N_k : np.ndarray, shape=(n_states)

Number of samples drawn from each state

s_n : np.ndarray, shape=(n_samples)

State of origin of each sample

sample(N_k=(10, 20, 30, 40, 50), mode='u_kln', seed=None)¶

Draw samples from the distribution.

Parameters:

N_k : np.ndarray, int

number of samples per state

mode : str, optional, default=’u_kln’

If ‘u_kln’, return K x K x N_max matrix where u_kln[k,l,n] is reduced potential of sample n from state k evaluated at state l. If ‘u_kn’, return K x N_tot matrix where u_kn[k,n] is reduced potential of sample n (in concatenated indexing) evaluated at state k.

seed: int, optional, default=None. Provides control over the random seed for replicability.

Returns:

if mode == ‘u_kn’:

x_n : np.ndarray, shape=(n_states*n_samples), dtype=float

x_n[n] is sample n (in concatenated indexing)

u_kn : np.ndarray, shape=(n_states, n_states*n_samples), dtype=float

u_kn[k,n] is reduced potential of sample n (in concatenated indexing) evaluated at state k.

N_k : np.ndarray, shape=(n_states), dtype=float

N_k[k] is the number of samples generated from state k

s_n : np.ndarray, shape=(n_samples), dtype=’int’

s_n is the state of origin of x_n

x_kn : np.ndarray, shape=(n_states, n_samples), dtype=float

1D harmonic oscillator positions

u_kln : np.ndarray, shape=(n_states, n_states, n_samples), dytpe=float, only if mode=’u_kln’

u_kln[k,l,n] is reduced potential of sample n from state k evaluated at state l.

N_k : np.ndarray, shape=(n_states), dtype=float

N_k[k] is the number of samples generated from state k

pymbar.testsystems.gaussian_work.gaussian_work_example(N_F=200, N_R=200, mu_F=2.0, DeltaF=None, sigma_F=1.0, seed=None)¶

Generate samples from forward and reverse Gaussian work distributions.

Parameters:

N_F : int, optional

number of forward measurements (default: 200)

N_R : float, optional

number of reverse measurements (default: 200)

mu_F : float, optional

mean of forward work distribution (default: 2.0)

DeltaF : float, optional

the free energy difference, which can be specified instead of mu_F (default: None)

sigma_F : float, optional

variance of the forward work distribution (default: 1.0)

seed : int, optional

If not None, specify the numpy random number seed. Old state is restored after completion.

Returns:

w_F : np.ndarray, dtype=float

forward work values

w_R : np.ndarray, dtype=float

reverse work values

Notes

By the Crooks fluctuation theorem (CFT), the forward and backward work distributions are related by

P_R(-w) = P_F(w) exp[DeltaF - w]

If the forward distribution is Gaussian with mean mu_F and std dev sigma_F, then

P_F(w) = (2 pi)^{-1/2} sigma_F^{-1} exp[-(w - mu_F)^2 / (2 sigma_F^2)]

With some algebra, we then find the corresponding mean and std dev of the reverse distribution are

mu_R = - mu_F + sigma_F^2 sigma_R = sigma_F exp[mu_F - sigma_F^2 / 2 + Delta F]

where all quantities are in reduced units (e.g. divided by kT).

Note that mu_F and Delta F are not independent! By the Zwanzig relation,

E_F[exp(-w)] = int dw exp(-w) P_F(w) = exp[-Delta F]

which, with some integration, gives

Delta F = mu_F + sigma_F^2/2

which can be used to determine either mu_F or DeltaF.

Examples

Generate work values with default parameters.

>>> [w_F, w_R] = gaussian_work_example()

Generate 50 forward work values and 70 reverse work values.

>>> [w_F, w_R] = gaussian_work_example(N_F=50, N_R=70)

Generate work values specifying the work distribution parameters.

>>> [w_F, w_R] = gaussian_work_example(mu_F=3.0, sigma_F=2.0)

Generate work values specifying the work distribution parameters, specifying free energy difference instead of mu_F.

>>> [w_F, w_R] = gaussian_work_example(mu_F=None, DeltaF=3.0, sigma_F=2.0)

Generate work values with known seed to ensure reproducibility for testing.

>>> [w_F, w_R] = gaussian_work_example(seed=0)

The testsystems Module: pymbar.testsystems¶

The testsystems Module: `pymbar.testsystems`¶