Add decoders algorithms for contrast and deconvolution based encoders

docs/api/python/decoders/index.rst

@@ -0,0 +1,17 @@
+.. _decoding:
+
+:octicon:`file-directory;0.9em;sd-mr-1 fill-primary` decoders
+=============================================================
+
+The ``decoders`` module within the spikify library contains the essential components for converting spike trains back into continuous signals or meaningful representations. This section is organized to reflect the primary structure of the decoding algorithms included in the library:
+
+- **Temporal Coding**: Encloses algorithms that decode data based on the precise timing of spikes, reconstructing signals from spike timing information.
+
+Each submodule is dedicated to a specific family of decoding techniques, making it easy to navigate and understand the purpose of each algorithm within the library structure.
+
+Below, you will find links to the specific modules library for each decoding method:
+
+.. toctree::
+   :maxdepth: 1
+
+   temporal_decoding/index

docs/api/python/decoders/temporal_decoding/contrast/decoder.rst

@@ -0,0 +1,8 @@
+.. _contrast_decoder_function:
+
+.. title:: Contrast Decoder
+
+.. automodule:: spikify.decoders.temporal.contrast.decoder_algorithm
+   :members: contrast_decoder
+   :undoc-members:
+   :show-inheritance:

docs/api/python/decoders/temporal_decoding/contrast/index.rst

@@ -0,0 +1,17 @@
+.. _contrast_decoder:
+
+:octicon:`file-directory;0.9em;sd-mr-1 fill-primary` contrast
+================================================================
+
+The ``contrast`` folder within the temporal decoding module includes the algorithm responsible for reconstructing the original continuous signal from spike trains generated by any of the Contrast family encoders. Since the encoding methods in this family share a common spike representation based on threshold increments, a single unified decoder is able to handle all of them.
+
+Contents of the ``contrast`` folder:
+
+- **Contrast Decoder**: A unified decoding algorithm compatible with spike trains produced by the Moving Window (MW), Step-Forward (SF), Threshold-Based (TB) encoders ans Zero-Crossing Step Forward (ZCSF).
+
+Below, you will find links to the specific modules for the contrast decoding algorithm:
+
+.. toctree::
+   :maxdepth: 1
+
+   decoder

docs/api/python/decoders/temporal_decoding/deconvolution/decoder.rst

@@ -0,0 +1,8 @@
+.. _deconvolution_decoder_function:
+
+.. title:: Deconvolution Decoder
+
+.. automodule:: spikify.decoders.temporal.deconvolution.decoder_algorithm
+   :members: deconvolution_decoder
+   :undoc-members:
+   :show-inheritance:

docs/api/python/decoders/temporal_decoding/deconvolution/index.rst

@@ -0,0 +1,22 @@
+.. _deconvolution_decoder:
+
+:octicon:`file-directory;0.9em;sd-mr-1 fill-primary` deconvolution
+====================================================================
+
+The ``deconvolution`` module within the temporal decoding module includes the algorithm responsible for
+reconstructing the original continuous signal from spike trains generated by any of the Deconvolution
+family encoders. Since the encoding methods in this family share a common reconstruction approach based
+on causal linear filtering with the same FIR filter used during encoding, a single unified decoder is
+able to handle all of them.
+
+Contents of the ``deconvolution`` folder:
+
+- **Deconvolution Decoder**: A unified decoding algorithm compatible with spike trains produced by the
+  Hough Spiker (HSA), Modified Hough Spiker (MHSA), and Ben's Spiker (BSA) encoders.
+
+Below, you will find links to the specific modules for the deconvolution decoding algorithm:
+
+.. toctree::
+   :maxdepth: 1
+
+   decoder

docs/api/python/decoders/temporal_decoding/index.rst

@@ -0,0 +1,19 @@
+.. _temporal_decoding:
+
+:octicon:`file-directory;0.9em;sd-mr-1 fill-primary` temporal
+================================================================
+
+The ``temporal`` module within the decoding module contains algorithms that reconstruct signals based on the timing of individual spikes. This method is especially effective for scenarios where the timing of events conveys more information than their frequency.
+
+Contents of the ``temporal`` module:
+
+- **Contrast**: Algorithms that decode based on variations in the signal over time.
+- **Deconvolution**: Approaches that reconstruct the original signal from a spike train using deconvolution techniques.
+
+Below, you will find links to the specific modules for each temporal decoding algorithm:
+
+.. toctree::
+   :maxdepth: 1
+
+   contrast/index
+   deconvolution/index

docs/api/python/index.rst

@@ -21,3 +21,4 @@ Below, you will find links to the specific modules within the Python API:
 
    filters/index
    encoders/index
+   decoders/index

spikify/decoders/__init__.py

@@ -0,0 +1 @@
+"""Decoders package."""

spikify/decoders/temporal/__init__.py

@@ -0,0 +1 @@
+"""Temporal package."""

spikify/decoders/temporal/contrast/__init__.py

@@ -0,0 +1,5 @@
+"""Temporal Contrast package."""
+
+from .decoder_algorithm import contrast_decoder
+
+__all__ = ["contrast_decoder"]

spikify/decoders/temporal/contrast/decoder_algorithm.py

@@ -0,0 +1,107 @@
+"""
+.. raw:: html
+
+    <h2>Contrast Decoder</h2>
+"""
+
+import numpy as np
+
+
+def contrast_decoder(
+    spikes: np.ndarray,
+    thresholds: np.ndarray,
+    start_point: float | int | list[float | int] | np.ndarray,
+) -> np.ndarray:
+    """
+    Perform Contrast family decoding on the input spike train.
+
+    This function takes a spike train encoded with any of the Contrast family algorithms,
+    Threshold-Based Representation (TBR), Step-Forward (SF), Moving Window (MW) or Zero-Crossing Step-Forward (ZCSF),
+    and reconstructs the original continuous signal. The reconstruction is performed by starting from the initial
+    signal value and applying the threshold incrementally based on the spike values.
+
+    .. note::
+        - The ``thresholds`` parameter should be the ``thresholds`` array returned directly by the encoder
+          (TBR, SF, or MW). All three encoders return this array as part of their output to ensure the step
+          size used for reconstruction exactly matches the one used during encoding.
+        - The ``start_point`` parameter should be set to the first sample of the original signal prior to
+          encoding (e.g. ``signal[0]`` for a 1D signal, or ``signal[0, :]`` for a multi-feature/channel signal).
+          Since the Contrast family encodes only differences between consecutive samples, the absolute signal
+          level is lost during encoding and must be restored by anchoring reconstruction to this initial value.
+
+
+    **Code Example:**
+
+    .. code-block:: python
+
+        import numpy as np
+        from spikify.encoders.temporal.contrast import moving_window
+        from spikify.decoders.temporal.contrast import contrast_decoder
+
+        signal = np.array([0.1, 0.3, 0.2, 0.5, 0.8, 1.0])
+        window_length = 3
+        spikes, thresholds = moving_window(signal, window_length, threshold=0.2)
+        reconstructed_signal = contrast_decoder(spikes, thresholds, start_point=signal[0])
+
+    .. doctest::
+        :hide:
+
+        >>> import numpy as np
+        >>> from spikify.encoders.temporal.contrast import moving_window
+        >>> from spikify.decoders.temporal.contrast import contrast_decoder
+        >>> signal = np.array([0.1, 0.3, 0.2, 0.5, 0.8, 1.0])
+        >>> spikes, thresholds = moving_window(signal, window_length=3, threshold=0.2)
+        >>> reconstructed_signal = contrast_decoder(spikes, thresholds, start_point=signal[0])
+        >>> reconstructed_signal.flatten()
+        array([0.1, 0.1, 0.1, 0.3, 0.5, 0.7])
+
+    :param spikes: The input spike train to be decoded. This should be a numpy ndarray with values in {-1, 0, +1},
+        as produced by any of the Contrast family encoders (TBR, SF, MW, ZCSF).
+    :type spikes: numpy.ndarray
+    :param threshold: Per-feature or channels threshold values used during encoding, as returned directly by the TBR,
+        SF, MW or ZCSF encoder. Passing the encoder's returned ``thresholds`` array
+        ensures the reconstruction step size exactly matches the one used during encoding.
+    :type threshold: numpy.ndarray
+    :param start_point: Initial signal value(s) for reconstruction (scalar or 1D sequence matching features).
+        Should be set to the first sample of the original signal before encoding (e.g. ``signal[0]``), as the
+        Contrast family encodes only signal differences and the absolute offset must be restored manually.
+    :type start_point: float | int | list[float | int] | numpy.ndarray
+    :return: A numpy array representing the reconstructed continuous signal approximation.
+    :rtype: numpy.ndarray
+    :raises ValueError: If the input spike train is empty or if the start_point dimensions do not match
+        the spike train feature dimensions.
+
+    """
+    # Check for empty spike train
+    if len(spikes) == 0:
+        raise ValueError("Spike train cannot be empty.")
+
+    # Ensure 2D processing (T, F)
+    if spikes.ndim == 1:
+        spikes = spikes.reshape(-1, 1)
+
+    T, F = spikes.shape
+
+    # Handle startpoint
+    if np.isscalar(start_point):
+        start_points = np.full(F, float(start_point))
+    else:
+        start_points = np.asarray(start_point, dtype=float)
+        if start_points.ndim != 1:
+            raise ValueError("Startpoint must be a scalar or a 1D sequence of numbers.")
+        if start_points.size != F:
+            raise ValueError("Startpoint must match the number of features in the spike train.")
+
+    signal = np.zeros((T, F), dtype=float)
+    signal[0] = start_points
+
+    for t in range(1, T):
+        for f in range(F):
+            if spikes[t, f] == 1:
+                signal[t, f] = signal[t - 1, f] + thresholds[f]
+            elif spikes[t, f] == -1:
+                signal[t, f] = signal[t - 1, f] - thresholds[f]
+            else:
+                signal[t, f] = signal[t - 1, f]
+
+    return signal

spikify/decoders/temporal/deconvolution/__init__.py

@@ -0,0 +1,5 @@
+"""Temporal Deconvolution package."""
+
+from .decoder_algorithm import deconvolution_decoder
+
+__all__ = ["deconvolution_decoder"]

spikify/decoders/temporal/deconvolution/decoder_algorithm.py

@@ -0,0 +1,102 @@
+"""
+.. raw:: html
+
+    <h2>Deconvolution Decoder</h2>
+"""
+
+import numpy as np
+from scipy.signal import lfilter
+
+
+def deconvolution_decoder(
+    spikes: np.ndarray,
+    fir_bank: np.ndarray,
+    shift: np.ndarray,
+    norm: np.ndarray | None = None,
+) -> np.ndarray:
+    """
+    Perform signal reconstruction via deconvolution decoding of a spike train.
+
+    The Deconvolution Decoder reverses the encoding performed by spike-based deconvolution algorithms such as
+    the Hough Spiker (HSA), its Modified variant (MHSA) and Bens Spiker (BSA). Given a binary spike train and the FIR
+    filter bank produced during encoding, it reconstructs an approximation of the original continuous signal
+    by applying each feature's filter to its corresponding spike train using causal linear filtering.
+
+    .. note::
+        - The ``fir_bank``, ``shift``, and ``norm`` parameters should be the values returned directly by the
+          encoder algorithm to ensure accurate inversion of all pre-processing steps.
+        - If ``norm`` is ``None``, no amplitude rescaling is applied; only the per-feature shift is restored.
+          Pass ``norm`` explicitly whenever the encoder performed normalization. Only MHSA and HSA have normalization
+          steps, while BSA does not.
+
+    **Code Example:**
+
+    .. code-block:: python
+
+        import numpy as np
+        from spikify.encoders.temporal.deconvolution import modified_hough_spiker
+        from spikify.decoders.temporal.deconvolution import deconvolution_decoder
+
+        signal = np.array([0.1, 0.2, 0.3, 1.0, 0.5, 0.3, 0.1])
+        window_length = 3
+        threshold = 0.5
+        cutoff = 0.1
+        spikes, fir_bank, shift, norm = modified_hough_spiker(signal, window_length, cutoff, threshold)
+        reconstructed = deconvolution_decoder(spikes, fir_bank, shift, norm)
+
+    .. doctest::
+        :hide:
+
+        >>> import numpy as np
+        >>> from spikify.encoders.temporal.deconvolution import modified_hough_spiker
+        >>> from spikify.decoders.temporal.deconvolution import deconvolution_decoder
+        >>> signal = np.array([0.0, 0.0, 0.5, 0.9, 0.7, 0.4, 0.1, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0])
+        >>> window_length = 5
+        >>> threshold = 0.3
+        >>> cutoff = 0.1
+        >>> spikes, fir_bank, shift, norm = modified_hough_spiker(signal, window_length, cutoff, threshold)
+        >>> reconstructed = deconvolution_decoder(spikes, fir_bank, shift, norm)
+        >>> reconstructed.flatten()
+        array([0.        , 0.24793707, 0.50412586, 0.49587414, 0.50412586,
+               0.24793707, 0.        , 0.        , 0.        , 0.        ,
+               0.        , 0.        , 0.        , 0.        , 0.24793707])
+
+    :param spikes: Binary spike train to decode (values in {0, 1}), as produced by
+        any of the deconvolution family encoders (HSA, MHSA, BSA).
+    :type spikes: numpy.ndarray
+    :param fir_bank: FIR filter coefficients used during encoding.
+        Each column corresponds to the filter applied to one feature or channel. This is the ``fir_bank`` value
+        returned directly by the encoder.
+    :type fir_bank: numpy.ndarray
+    :param shift: Per-feature shift values subtracted during encoding to ensure signal non-negativity.
+                  This is the ``shift`` value returned directly by the encoder.
+    :type shift: numpy.ndarray
+    :param norm: Per-feature normalization values used during encoding to scale the signal to [0, 1].
+                 If ``None``, no amplitude rescaling is applied and only the shift is restored.
+                 This is the ``norm`` value returned directly by the encoder.
+    :type norm: numpy.ndarray | None
+    :return: A numpy array representing the reconstructed continuous signal approximation.
+    :rtype: numpy.ndarray
+    :raises ValueError: If the input spike train is empty.
+
+    """
+
+    # Check for empty spike train
+    if len(spikes) == 0:
+        raise ValueError("Spike train cannot be empty.")
+
+    # Ensure 2D processing (T, F)
+    if spikes.ndim == 1:
+        spikes = spikes.reshape(-1, 1)
+
+    T, F = spikes.shape
+    signal = np.zeros((T, F))
+
+    for f in range(F):
+        x = lfilter(fir_bank[:, f], 1, spikes[:, f], axis=0)
+        if norm is not None:
+            signal[:, f] = x * norm[f] + shift[f]
+        else:
+            signal[:, f] = x + shift[f]
+
+    return signal

spikify/encoders/temporal/contrast/moving_window_algorithm.py

@@ -9,7 +9,7 @@
 
 def moving_window(
     signal: np.ndarray, window_length: int, threshold: float | int | list[float | int] | np.ndarray
-) -> np.ndarray:
+) -> tuple[np.ndarray, np.ndarray]:
     """
     Perform Moving Window (MW) encoding on the input signal.
 
@@ -29,8 +29,7 @@ def moving_window(
         signal = np.array([0.1, 0.3, 0.2, 0.5, 0.8, 1.0])
         window_length = 3
         threshold = 0.2
-        encoded_signal = moving_window(signal, window_length, threshold)
-        encoded_signal
+        encoded_signal, thresholds = moving_window(signal, window_length, threshold)
 
     .. doctest::
         :hide:
@@ -40,20 +39,23 @@ def moving_window(
         >>> signal = np.array([0.1, 0.3, 0.2, 0.5, 0.8, 1.0])
         >>> window_length = 3
         >>> threshold = 0.2
-        >>> encoded_signal = moving_window(signal, window_length, threshold)
-        >>> encoded_signal
+        >>> encoded_signal, _ = moving_window(signal, window_length, threshold)
+        >>> encoded_signal.flatten()
         array([0, 0, 0, 1, 1, 1], dtype=int8)
 
-    :param signal: The input signal to be encoded. This should be a numpy ndarray.
+    :param signal: Input signal to encode (1D or 2D: time × features or channels).
     :type signal: numpy.ndarray
     :param window_length: The size of the sliding window for calculating the signal base mean.
     :type window_length: int
     :param threshold: Threshold(s) for spike generation; scalar or 1D sequence matching features.
     :type threshold: float | int | list[float | int] | numpy.ndarray
-    :return: A numpy array representing the encoded spike train.
-    :rtype: numpy.ndarray
-    :raises ValueError: If the input signal is empty or if the threshold dimensions do not match the signal features.
-    :raises TypeError: If the threshold parameter is of invalid dimension.
+    :return:
+        - spikes: A numpy array representing the encoded spike train. (values in {-1, 0, +1})
+        - thresholds: Per-feature or channel thresholds used for encoding, returned for use in decoding,
+          shape (features or channels,).
+    :rtype: tuple[numpy.ndarray, numpy.ndarray]
+    :raises ValueError: If the input signal is empty or if the threshold dimensions do not match the signal f
+                        eature dimensions.
 
     """
 
@@ -73,7 +75,7 @@ def moving_window(
     else:
         thresholds = np.asarray(threshold, dtype=float)
         if thresholds.ndim != 1:
-            raise TypeError("Threshold must be a scalar or a 1D sequence of numbers.")
+            raise ValueError("Threshold must be a scalar or a 1D sequence of numbers.")
         if thresholds.size != F:
             raise ValueError("Threshold must match the number of features in the signal.")
 
@@ -100,8 +102,4 @@ def moving_window(
             elif signal[t, f] < base - thresholds[f]:
                 spikes[t, f] = -1
 
-    # Flatten if input was 1D
-    if F == 1:
-        spikes = spikes.flatten()
-
-    return spikes
+    return spikes, thresholds