From 0ffce2f17a9405d069789f23b5176039d31d0b74 Mon Sep 17 00:00:00 2001 From: robrobo Date: Thu, 17 Jul 2025 18:45:11 +0200 Subject: [PATCH] added full_output option to non gaussian parameter so that correct statistics can be done afterwards --- src/mdevaluate/correlation.py | 36 +++++++++++++++++++++++------------ 1 file changed, 24 insertions(+), 12 deletions(-) diff --git a/src/mdevaluate/correlation.py b/src/mdevaluate/correlation.py index 11be7da..07d667c 100644 --- a/src/mdevaluate/correlation.py +++ b/src/mdevaluate/correlation.py @@ -18,7 +18,7 @@ def log_indices(first: int, last: int, num: int = 100) -> np.ndarray: return np.unique(np.int_(ls) - 1 + first) -@autosave_data(2) +@autosave_data(nargs=2, kwargs_keys=('selector', 'segments', 'skip', 'window', 'average', 'points',), version=1.0) def shifted_correlation( function: Callable, frames: Coordinates, @@ -460,6 +460,7 @@ def non_gaussian_parameter( end_frame: CoordinateFrame, trajectory: Coordinates = None, axis: str = "all", + full_output = False, ) -> float: r""" Calculate the non-Gaussian parameter. @@ -472,30 +473,41 @@ def non_gaussian_parameter( else: vectors = displacements_without_drift(start_frame, end_frame, trajectory) if axis == "all": - r = (vectors**2).sum(axis=1) + r2 = (vectors**2).sum(axis=1) dimensions = 3 elif axis == "xy" or axis == "yx": - r = (vectors[:, [0, 1]]**2).sum(axis=1) + r2 = (vectors[:, [0, 1]]**2).sum(axis=1) dimensions = 2 elif axis == "xz" or axis == "zx": - r = (vectors[:, [0, 2]]**2).sum(axis=1) + r2 = (vectors[:, [0, 2]]**2).sum(axis=1) dimensions = 2 elif axis == "yz" or axis == "zy": - r = (vectors[:, [1, 2]]**2).sum(axis=1) + r2 = (vectors[:, [1, 2]]**2).sum(axis=1) dimensions = 2 elif axis == "x": - r = vectors[:, 0] ** 2 + r2 = vectors[:, 0] ** 2 dimensions = 1 elif axis == "y": - r = vectors[:, 1] ** 2 + r2 = vectors[:, 1] ** 2 dimensions = 1 elif axis == "z": - r = vectors[:, 2] ** 2 + r2 = vectors[:, 2] ** 2 dimensions = 1 else: raise ValueError('Parameter axis has to be ether "all", "x", "y", or "z"!') - mean_r = np.mean(r) - if mean_r == 0.0: - return 0.0 - return (np.mean(r**2) / ((1 + 2 / dimensions) * (mean_r ** 2))) - 1 + m2 = np.mean(r2) + m4 = np.mean(r2**2) + if m2 == 0.0: + if full_output: + return 0.0, 0.0, 0.0 + else: + return 0.0 + + alpha_2 = (m4 / ((1 + 2 / dimensions) * m2**2)) - 1 + if full_output: + return alpha_2, m2, m4 + else: + return alpha_2 + +