diff --git a/cmeutils/sampling.py b/cmeutils/sampling.py index 4f1162b..c43455b 100644 --- a/cmeutils/sampling.py +++ b/cmeutils/sampling.py @@ -103,3 +103,23 @@ def is_equilibrated(data, threshold_fraction=0.50, threshold_neff=50, nskip=1): return [True, t0, g, Neff] else: return [False, None, None, None] + + +def autocorr1D(array): + """Takes in a linear np array, performs autocorrelation + function and returns normalized array with half the length + of the input. + + Parameters + ---------- + data : numpy.typing.Arraylike, required + 1-D series of data to perform autocorrelation on. + + Returns + ------- + 1D np.array + + """ + ft = np.fft.rfft(array - np.average(array)) + acorr = np.fft.irfft(ft * np.conjugate(ft)) / (len(array) * np.var(array)) + return acorr[0 : len(acorr) // 2] # noqa: E203 diff --git a/cmeutils/tests/test_gsd.py b/cmeutils/tests/test_gsd.py index 46cf35a..3be87d9 100644 --- a/cmeutils/tests/test_gsd.py +++ b/cmeutils/tests/test_gsd.py @@ -50,7 +50,7 @@ def test_ellipsoid_gsd(self, butane_gsd): def test_create_rigid_snapshot(self): benzene = mb.load("c1ccccc1", smiles=True) benzene.name = "Benzene" - box = mb.fill_box(benzene, 5, box=[1, 1, 1]) + box = mb.fill_box(benzene, 5, box=[3, 3, 3], seed=42) box.label_rigid_bodies(discrete_bodies="Benzene") rigid_init_snap = create_rigid_snapshot(box) @@ -60,7 +60,7 @@ def test_create_rigid_snapshot(self): def test_update_rigid_snapshot(self): benzene = mb.load("c1ccccc1", smiles=True) benzene.name = "Benzene" - box = mb.fill_box(benzene, 5, box=[1, 1, 1]) + box = mb.fill_box(benzene, 5, box=[3, 3, 3], seed=42) box.label_rigid_bodies(discrete_bodies="Benzene") rigid_init_snap = create_rigid_snapshot(box) diff --git a/cmeutils/tests/test_sampling.py b/cmeutils/tests/test_sampling.py index 8fbe6bc..f6a004b 100644 --- a/cmeutils/tests/test_sampling.py +++ b/cmeutils/tests/test_sampling.py @@ -1,7 +1,8 @@ import numpy as np import pytest +import scipy -from cmeutils.sampling import equil_sample, is_equilibrated +from cmeutils.sampling import autocorr1D, equil_sample, is_equilibrated from cmeutils.tests.base_test import BaseTest @@ -82,3 +83,28 @@ def test_trim_high_threshold(self, correlated_data_tau100_n10000): [equil_data, uncorr_indices, prod_start, Neff] = equil_sample( data, threshold_fraction=0.75, threshold_neff=10000 ) + + def test_autocorr1D(self): + randoms = np.random.randint(100, size=(100)) + integers = np.array(np.arange(100)) + + autocorrelation_randoms = autocorr1D(randoms) + autocorrelation_integers = autocorr1D(integers) + + def expfunc(x, a): + return np.exp(-x / a) + + x_values = np.array(range(len(autocorrelation_integers))) + + exp_coeff_randoms = scipy.optimize.curve_fit( + expfunc, x_values, autocorrelation_randoms + )[0][0] + exp_coeff_int = scipy.optimize.curve_fit( + expfunc, x_values, autocorrelation_integers + )[0][0] + + assert ( + round(autocorrelation_integers[0], 10) + and round(autocorrelation_randoms[0], 10) == 1 + ) + assert exp_coeff_int > 5 and exp_coeff_randoms < 1 diff --git a/environment-dev.yml b/environment-dev.yml index 650908f..65c1156 100644 --- a/environment-dev.yml +++ b/environment-dev.yml @@ -10,7 +10,7 @@ dependencies: - mbuild >=0.16.4 - numpy - pip - - python >= 3.10 + - python >= 3.10,<3.12 - matplotlib - pre-commit - pymbar >= 4.0 diff --git a/environment.yml b/environment.yml index 1be941e..eafb255 100644 --- a/environment.yml +++ b/environment.yml @@ -10,7 +10,7 @@ dependencies: - mbuild >=0.16.4 - numpy - pip - - python >= 3.10 + - python >= 3.10,<3.12 - matplotlib - pymbar >= 4.0 - rowan