[NSETM-2281] Fix calculation of dynamic offset (#19)

.github/workflows/publish-sdist.yml

@@ -12,7 +12,7 @@ jobs:
     steps:
       - uses: actions/checkout@v4
       - name: Set up Python 3.11
-        uses: actions/setup-python@v4
+        uses: actions/setup-python@v5
         with:
           python-version: 3.11
 

.github/workflows/run-tox.yml

@@ -16,7 +16,7 @@ jobs:
     steps:
     - uses: actions/checkout@v4
     - name: Set up Python ${{ matrix.python-version }}
-      uses: actions/setup-python@v4
+      uses: actions/setup-python@v5
       with:
         python-version: ${{ matrix.python-version }}
 

CHANGELOG.rst

@@ -1,6 +1,15 @@
 Changelog
 =========
 
+Version 0.7.0
+-------------
+
+New Features
+~~~~~~~~~~~~
+
+- Allow to specify ``trial_steps_label`` to calculate the dynamic offset of trial steps [NSETM-2281]
+
+
 Version 0.6.0
 -------------
 

src/blueetl/config/analysis.py

@@ -5,7 +5,7 @@
 from copy import deepcopy
 from itertools import chain
 from pathlib import Path
-from typing import NamedTuple, Optional, Union
+from typing import NamedTuple, Union
 
 from blueetl.config.analysis_model import (
     FeaturesConfig,
@@ -20,13 +20,23 @@
 L = logging.getLogger(__name__)
 
 
-def _resolve_paths(global_config: MultiAnalysisConfig, base_path: Optional[Path] = None) -> None:
+def _resolve_paths(global_config: MultiAnalysisConfig, base_path: Path) -> None:
     """Resolve any relative path."""
     base_path = base_path or Path()
     global_config.output = base_path / global_config.output
     global_config.simulation_campaign = base_path / global_config.simulation_campaign
 
 
+def _resolve_trial_steps(global_config: MultiAnalysisConfig):
+    """Set trial_steps_config.base_path to the same value as global_config.output.
+
+    In this way, the custom function can use it as the base path to save any figure.
+    """
+    for config in global_config.analysis.values():
+        for trial_steps_config in config.extraction.trial_steps.values():
+            trial_steps_config.base_path = str(global_config.output)
+
+
 def _resolve_windows(global_config: MultiAnalysisConfig) -> None:
     """Calculate the hash of any referenced windows in each single analysis configuration.
 
@@ -140,13 +150,12 @@ def _resolve_analysis_configs(global_config: MultiAnalysisConfig) -> None:
         config.features = _resolve_features(config.features)
 
 
-def init_multi_analysis_configuration(
-    global_config: dict, base_path: Optional[Path] = None
-) -> MultiAnalysisConfig:
+def init_multi_analysis_configuration(global_config: dict, base_path: Path) -> MultiAnalysisConfig:
     """Return a config object from a config dict."""
     validate_config(global_config, schema=read_schema("analysis_config"))
     config = MultiAnalysisConfig(**global_config)
     _resolve_paths(config, base_path=base_path)
+    _resolve_trial_steps(config)
     _resolve_windows(config)
     _resolve_analysis_configs(config)
     return config

src/blueetl/config/analysis_model.py

@@ -68,8 +68,6 @@ class WindowConfig(BaseModel):
     @model_validator(mode="after")
     def validate_values(self):
         """Validate the values after loading them."""
-        if self.trial_steps_label:
-            raise ValueError("trial_steps_label cannot be used yet, see NSETM-2281")
         if self.trial_steps_list and (self.n_trials or self.trial_steps_value):
             raise ValueError("trial_steps_list cannot be set with n_trials or trial_steps_value")
         if self.n_trials > 1 and not self.trial_steps_value:
@@ -84,13 +82,22 @@ class TrialStepsConfig(BaseModel):
         **BaseModel.model_config,
         "extra": "allow",
     }
+    _forbidden_extra_fields: set[str] = {
+        "initial_offset",
+    }
     function: str
-    initial_offset: float = 0.0
     bounds: tuple[float, float]
     population: Optional[str] = None
     node_set: Optional[str] = None
     limit: Optional[int] = None
 
+    @model_validator(mode="after")
+    def forbid_fields(self):
+        """Verify that the forbidden extra fields have not been specified."""
+        if found := self._forbidden_extra_fields.intersection(self.model_extra):
+            raise ValueError(f"Forbidden extra fields: {found}")
+        return self
+
 
 class NeuronClassConfig(BaseModel):
     """NeuronClassConfig Model."""

src/blueetl/constants.py

@@ -21,8 +21,6 @@
 DURATION = "duration"
 WINDOW_TYPE = "window_type"
 COUNT = "count"
-TRIAL_STEPS_LABEL = "trial_steps_label"
-TRIAL_STEPS_VALUE = "trial_steps_value"
 TIMES = "times"
 BIN = "bin"
 VALUE = "value"

src/blueetl/external/bnac/calculate_trial_step.py

@@ -1,41 +1,50 @@
 """Trial steps functions adapted from BlueNetworkActivityComparison/bnac/onset.py."""
 
+import logging
+from pathlib import Path
+
 import numpy as np
 from scipy.ndimage import gaussian_filter
 
+L = logging.getLogger(__name__)
+
+
+def _get_bounds(params):
+    lower_bound, upper_bound = params["bounds"]
+    assert lower_bound <= 0
+    assert upper_bound >= 0
+    return lower_bound, upper_bound
+
 
 def _histogram_from_spikes(spikes, params):
-    onset_test_window_length = params["post_window"][1] - params["pre_window"][0]
+    lower_bound, upper_bound = _get_bounds(params)
+    onset_test_window_length = upper_bound - lower_bound
     histogram, _ = np.histogram(
         spikes,
-        range=[params["pre_window"][0], params["post_window"][1]],
+        range=(lower_bound, upper_bound),
         bins=int(onset_test_window_length * params["histo_bins_per_ms"]),
     )
     return histogram
 
 
 def _onset_from_histogram(histogram, params):
-    onset_dict = {}
+    lower_bound, _ = _get_bounds(params)
     smoothed_histogram = gaussian_filter(histogram, sigma=params["smoothing_width"])
 
-    onset_pre_window_length = params["pre_window"][1] - params["pre_window"][0]
-    onset_zeroed_post_start = params["post_window"][0] - params["pre_window"][0]
-    pre_smoothed_histogram = smoothed_histogram[
-        : onset_pre_window_length * params["histo_bins_per_ms"]
-    ]
-    post_smoothed_histogram = smoothed_histogram[
-        onset_zeroed_post_start * params["histo_bins_per_ms"] :
-    ]
-
-    onset_dict["pre_mean"] = np.mean(pre_smoothed_histogram)
-    onset_dict["pre_std"] = np.std(pre_smoothed_histogram)
-    onset_dict["post_max"] = np.max(post_smoothed_histogram)
+    pre_window_length = -lower_bound
+    pre_window_bins = int(pre_window_length * params["histo_bins_per_ms"])
+    pre_smoothed_histogram = smoothed_histogram[:pre_window_bins]
+    post_smoothed_histogram = smoothed_histogram[pre_window_bins:]
+
+    onset_dict = {
+        "pre_mean": np.mean(pre_smoothed_histogram),
+        "pre_std": np.std(pre_smoothed_histogram),
+        "post_max": np.max(post_smoothed_histogram),
+    }
     onset_dict["pre_mean_post_max_ratio"] = onset_dict["pre_mean"] / onset_dict["post_max"]
 
-    where_above_thresh = np.where(
-        post_smoothed_histogram
-        > (onset_dict["pre_mean"] + params["threshold_std_multiple"] * onset_dict["pre_std"])
-    )[0]
+    threshold = onset_dict["pre_mean"] + params["threshold_std_multiple"] * onset_dict["pre_std"]
+    where_above_thresh = np.where(post_smoothed_histogram > threshold)[0]
     index_above_std = 0
     if len(where_above_thresh) > 0:
         index_above_std = where_above_thresh[0]
@@ -45,13 +54,10 @@ def _onset_from_histogram(histogram, params):
     cortical_onset_index = index_above_std
 
     onset_dict["cortical_onset"] = (
-        float(cortical_onset_index) / float(params["histo_bins_per_ms"])
-        + float(params["post_window"][0])
-        + params["ms_post_offset"]
+        float(cortical_onset_index) / float(params["histo_bins_per_ms"]) + params["ms_post_offset"]
     )
-    if params["fig_paths"]:
+    if params.get("figures_path"):
         _plot(smoothed_histogram, params, onset_dict)
-    onset_dict["trial_steps_value"] = onset_dict["cortical_onset"]
     return onset_dict
 
 
@@ -63,27 +69,44 @@ def _plot(smoothed_histogram, params, onset_dict):
     import matplotlib.pyplot as plt
     import seaborn as sns
 
+    lower_bound, upper_bound = _get_bounds(params)
     plt.figure()
     sns.set()
     sns.set_style("ticks")
-    x_vals = list(np.arange(params["pre_window"][0], params["post_window"][1], 0.2))
+    x_vals = list(np.arange(lower_bound, upper_bound, 0.2))
     plt.plot(x_vals, smoothed_histogram)
     plt.scatter(
         onset_dict["cortical_onset"],
         [onset_dict["pre_mean"] + params["threshold_std_multiple"] * onset_dict["pre_std"]],
     )
-    plt.gca().set_xlim([params["pre_window"][0], params["post_window"][1]])
+    plt.gca().set_xlim([lower_bound, upper_bound])
     plt.gca().set_xlabel("Time (ms)")
     plt.gca().set_ylabel("Number of spikes")
     plt.gca().spines["right"].set_visible(False)
     plt.gca().spines["top"].set_visible(False)
-    for fig_path in params["fig_paths"]:
-        plt.savefig(fig_path)
+    filepath = Path(params["base_path"], params["figures_path"], "plot.pdf")
+    L.info("Figures path: %s", filepath)
+    filepath.parent.mkdir(exist_ok=True)
+    plt.savefig(filepath)
     plt.close()
 
 
-def onset_from_spikes(spikes, params):
-    """Calculate trial steps from spikes."""
+def onset_from_spikes(spikes_list, params):
+    """Calculate the cortical onset from a list of spikes, one for each trial.
+
+    Args:
+        spikes_list: list of spikes as numpy arrays.
+        params: dictionary of parameters from the trial steps configuration.
+
+    Returns:
+        float representing the dynamic offset to be added to the initial offset of each trial step.
+    """
+    L.info(
+        "onset_from_spikes: processing %s arrays of spikes using params=%r",
+        len(spikes_list),
+        params,
+    )
+    spikes = np.concatenate(spikes_list)
     histogram = _histogram_from_spikes(spikes, params)
     onset_dict = _onset_from_histogram(histogram, params)
-    return onset_dict
+    return onset_dict["cortical_onset"]