Merge branch 'main' into nmrpca-fix-egligibility

Author: comcon1

src/fairmd/lipids/analib/maicos.py

@@ -25,7 +25,7 @@
 from maicos.lib.weights import density_weights
 from tqdm import tqdm
 
-from fairmd.lipids.auxiliary.jsonEncoders import CompactJSONEncoder
+from fairmd.lipids.auxiliary import CompactJSONEncoder
 from fairmd.lipids.core import System
 from fairmd.lipids.molecules import lipids_set
 

src/fairmd/lipids/analyze.py

@@ -43,8 +43,7 @@
     traj_centering_for_maicos_mda_parallel,
 )
 from fairmd.lipids.api import UniverseConstructor, mda_gen_selection_mols
-from fairmd.lipids.auxiliary import mollib
-from fairmd.lipids.auxiliary.jsonEncoders import CompactJSONEncoder
+from fairmd.lipids.auxiliary import CompactJSONEncoder, mollib
 from fairmd.lipids.core import System
 from fairmd.lipids.molecules import lipids_set
 

src/fairmd/lipids/api.py

@@ -33,6 +33,7 @@
 import numpy as np
 
 from fairmd.lipids import FMDL_SIMU_PATH
+from fairmd.lipids.auxiliary import block_average_time_series
 from fairmd.lipids.core import System
 from fairmd.lipids.databankio import download_resource_from_uri, resolve_file_url
 from fairmd.lipids.molecules import Molecule, lipids_set
@@ -218,13 +219,14 @@ def get_quality(
     return q
 
 
-def get_mean_ApL(system: System) -> float:  # noqa: N802 (API name)
+def get_ApL_data(system: System, blocksize: float | None = None) -> np.ndarray:  # noqa: N802 (API name)
     """
-    Calculate average area per lipid for a system.
+    Return Area-per-lipid data as a numpy array (block-averaging possible).
 
-    :param system: Simulation object.
+    :param system: Simulation object
+    :param blocksize: Averaged t-series by <blocksize> ps
 
-    :return: area per lipid (Å^2)
+    :return: Array (t, value) with blocksize step.
     """
     path = os.path.join(FMDL_SIMU_PATH, system["path"], "apl.json")
     try:
@@ -236,8 +238,27 @@ def get_mean_ApL(system: System) -> float:  # noqa: N802 (API name)
     except json.JSONDecodeError as e:
         msg = "Area per lipid data for system #{} in {} is invalid.".format(system["ID"], path)
         raise ValueError(msg) from e
-    vals = np.array(list(data.values()))
-    return vals.mean()
+    df = np.vstack(
+        [
+            np.array(list(data.keys()), dtype=float),
+            np.array(list(data.values()), dtype=float),
+        ]
+    ).T
+    if blocksize is not None:
+        df = block_average_time_series(df, blocksize)
+    return df
+
+
+def get_mean_ApL(system: System) -> float:  # noqa: N802 (API name)
+    """
+    Calculate average area per lipid for a system.
+
+    :param system: Simulation object.
+
+    :return: area per lipid (Å^2)
+    """
+    df = get_ApL_data(system)
+    return df[:, 1].mean()
 
 
 def get_total_area(system: System) -> float:

src/fairmd/lipids/auxiliary/__init__.py

@@ -1 +1,9 @@
 # file is required for fairmd.lipids.settings to be accessible as a subpackage
+
+from .jsonEncoders import CompactJSONEncoder
+from .routines import block_average_time_series
+
+__all__ = [
+    "CompactJSONEncoder",
+    "block_average_time_series",
+]

src/fairmd/lipids/auxiliary/jsonEncoders.py

@@ -1,9 +1,8 @@
-#!/usr/bin/env python3
+"""
+Compact JSON encoder class.
 
-#
-# Compact JSON from @jannismain is used here.
-# https://gist.github.com/jannismain/e96666ca4f059c3e5bc28abb711b5c92
-#
+:meta: private
+"""
 
 import json
 import math
@@ -12,7 +11,11 @@
 
 
 class CompactJSONEncoder(json.JSONEncoder):
-    """A JSON Encoder that puts small containers on single lines."""
+    """A JSON Encoder that puts small containers on single lines.
+
+    Originally taken from `the GiST of @jannismain
+    <https://gist.github.com/jannismain/e96666ca4f059c3e5bc28abb711b5c92>`_.
+    """
 
     CONTAINER_TYPES = (list, tuple, dict)
     """Container datatypes include primitives or other containers."""
@@ -100,30 +103,5 @@ def indent_str(self) -> str:
             return " " * (self.indentation_level * self.indent)
         if isinstance(self.indent, str):
             return self.indentation_level * self.indent
-        raise ValueError(
-            f"indent must either be of type int or str (is: {type(self.indent)})",
-        )
-
-
-if __name__ == "__main__":
-    import sys
-
-    if "--example" in sys.argv:
-        data = {
-            "compact_object": {"first": "element", "second": 2},
-            "compact_list": ["first", "second"],
-            "long_list": [
-                "this",
-                "is",
-                "a",
-                "rather",
-                "long\nlist",
-                "and should be broken up because of its width",
-            ],
-            "non_ascii": "汉语",
-            1: 2,
-        }
-        json.dump(data, sys.stdout, cls=CompactJSONEncoder, ensure_ascii=False)
-        exit()
-
-    json.dump(json.load(sys.stdin), sys.stdout, cls=CompactJSONEncoder)
+        msg = f"Indent must either be of type int or str (is: {type(self.indent)})"
+        raise ValueError(msg)

src/fairmd/lipids/auxiliary/routines.py

@@ -0,0 +1,52 @@
+"""
+Algorithmic routines. Pure numpy vector stuff.
+
+:meta: private
+"""
+
+import numpy as np
+
+
+def block_average_time_series(arr: np.ndarray, blocksize: float) -> np.ndarray:
+    """
+    Return block-average of (N,2) array.
+
+    :param arr: numpy array (N,2)
+    :param blocksize: timeframe for slice blocks
+
+    :return: block-averaged array
+    """
+    # arr: (N, 2) -> [time, value]
+    t = arr[:, 0]
+    if t[1] - t[0] > blocksize:
+        msg = f"Blocksize ({blocksize}) must be greater than write timestep ({t[1] - t[0]})"
+        raise ValueError(msg)
+    x = arr[:, 1]
+
+    t0 = t.min()
+    t1 = t.max()
+    # Bin edges: [t0, t0+blocksize, t0+2*blocksize, ...]
+    edges = np.arange(t0, t1 + blocksize + 1e-3, blocksize)  # +1e-3 required to make t1 // blocksize working!
+    nbins = len(edges) - 1
+
+    # Bin index for each sample: 0..len(edges)
+    idx = np.digitize(t, edges) - 1
+
+    # Prepare output arrays
+    sums = np.zeros(nbins)
+    counts = np.zeros(nbins, dtype=int)
+
+    # Accumulate sums and counts per bin
+    np.add.at(sums, idx, x)
+    np.add.at(counts, idx, 1)
+
+    # Avoid division by zero
+    mask = counts > 0
+    avg_vals = np.empty(nbins)
+    avg_vals[:] = np.nan
+    avg_vals[mask] = sums[mask] / counts[mask]
+
+    # Define a representative time for each bin (e.g. center)
+    bin_times = (edges[:-1] + edges[1:]) / 2.0
+
+    return np.column_stack((bin_times, avg_vals))

src/fairmd/lipids/bin/evaluate_quality.py

@@ -22,7 +22,7 @@
 
 import fairmd.lipids.quality as qq
 from fairmd.lipids import FMDL_EXP_PATH, FMDL_SIMU_PATH
-from fairmd.lipids.auxiliary.jsonEncoders import CompactJSONEncoder
+from fairmd.lipids.auxiliary import CompactJSONEncoder
 
 
 def _round_quality_values(obj: dict | list, ndigits: int = 4) -> dict | list:

tests/test_api.py

@@ -119,6 +119,29 @@ def test_get_mean_apl(systems, systemid, result):
     assert apm == pytest.approx(result, abs=6e-4)
 
 
+@pytest.mark.parametrize(
+    "systemid, nlines",
+    [
+        (281, 1001),
+        (566, 401),  # 1
+    ],
+)
+def test_get_apl_data(systems, systemid, nlines):
+    from fairmd.lipids.api import get_ApL_data
+
+    s = systems.loc(systemid)
+    df = get_ApL_data(s)
+    check.is_true(isinstance(df, np.ndarray))
+    check.equal(df.shape[1], 2)
+    check.equal(df.shape[0], nlines)
+    # block-average behavior
+    df1k = get_ApL_data(s, blocksize=1000)
+    df2k = get_ApL_data(s, blocksize=2000)
+    df3k = get_ApL_data(s, blocksize=3000)
+    check.almost_equal(df1k[0:2, 1].mean(), df2k[0, 1], abs=1e-7)
+    check.almost_equal(df1k[0:3, 1].mean(), df3k[0, 1], abs=1e-7)
+
+
 @pytest.mark.parametrize(
     "systemid, result",
     [(281, 4142.234), (566, 3923.568), (787, 4694.191), (243, 2241.920), (86, 3869.417)],

tests/test_misc.py

@@ -11,7 +11,10 @@
 from contextlib import contextmanager
 import os
 import shutil
+import json
 
+import numpy as np
+import numpy.testing as npt
 import pytest
 import pytest_check as check
 
@@ -215,3 +218,56 @@ def test_get_tails():
     check.equal(len(tails), 4)
     check.is_in("sn-1 1", tails)
     check.is_in("sn-2 1", tails)
+
+
+def test_json_encoder(tmpdir):
+    """Fixture for a mock experiment path with no data files."""
+    from fairmd.lipids.auxiliary import CompactJSONEncoder
+
+    exp_dir = tmpdir.mkdir("jsonenc")
+
+    data = {
+        "compact_object": {"first": "element", "second": 2},
+        "compact_list": ["first", "second"],
+        "long_list": [
+            "this",
+            "is",
+            "a",
+            "rather",
+            "long\nlist",
+            "and should be broken up because of its width",
+        ],
+        "non_ascii": "汉语",
+        1: 2,
+    }
+    with open(exp_dir.join("test.yaml"), "w", encoding="utf-8") as fd:
+        json.dump(data, fd, cls=CompactJSONEncoder, ensure_ascii=False)
+
+
+def test_average_block():
+    from fairmd.lipids.auxiliary import block_average_time_series
+
+    # times 0..9, values = times
+    t = np.arange(10, dtype=float)
+    x = t.copy()
+    arr = np.column_stack((t, x))
+
+    out = block_average_time_series(arr, blocksize=2.0)
+    # expected bins: [0,2), [2,4), [4,6), [6,8), [8,10]
+    expected_times = np.array([1, 3, 5, 7, 9], dtype=float)
+    expected_vals = np.array(
+        [
+            (0 + 1) / 2,
+            (2 + 3) / 2,
+            (4 + 5) / 2,
+            (6 + 7) / 2,
+            (8 + 9) / 2,
+        ],
+        dtype=float,
+    )
+
+    npt.assert_allclose(out[:, 0], expected_times)
+    npt.assert_allclose(out[:, 1], expected_vals)
+
+    arr = arr[:-1, :]
+    out = block_average_time_series(arr, blocksize=2.0)