added testing

docs/examples/plot_glm_demo.py

@@ -252,17 +252,20 @@
 excit_b = np.random.uniform(1.1, 5, size=(n_neurons, n_neurons))
 
 # define 2x2 coupling filters of the specific with create_temporal_filter
-coupling_filter_bank = np.zeros((n_neurons, n_neurons, coupling_filter_duration))
+coupling_filter_bank = np.zeros((coupling_filter_duration, n_neurons, n_neurons))
 for unit_i in range(n_neurons):
     for unit_j in range(n_neurons):
-        coupling_filter_bank[unit_i, unit_j, :] = nmo.simulation.difference_of_gammas(
+        coupling_filter_bank[:, unit_i, unit_j] = nmo.simulation.difference_of_gammas(
             coupling_filter_duration,
             inhib_a=inhib_a,
             excit_a=excit_a[unit_i, unit_j],
             inhib_b=inhib_b,
             excit_b=excit_b[unit_i, unit_j],
         )
 
+# shrink the filters for simulation stability
+coupling_filter_bank *= 0.8
+
 # %%
 # If we represent our filters in terms of basis functions, we can simulate our network by
 # directly calling the `simulate` method of the `nmo.glm.GLMRecurrent` class.
@@ -272,7 +275,8 @@
 basis = nmo.basis.RaisedCosineBasisLog(n_basis_funcs)
 
 # approximate the coupling filters in terms of the basis function
-coupling_basis, coupling_coeff = simulation.regress_filter(coupling_filter_bank, basis)
+_, coupling_basis = basis.evaluate_on_grid(coupling_filter_bank.shape[0])
+coupling_coeff = simulation.regress_filter(coupling_filter_bank, coupling_basis)
 intercept = -4 * np.ones(n_neurons)
 
 # %%
@@ -287,7 +291,7 @@
     for unit_j in range(n_neurons):
         axs[unit_i, unit_j].set_title(f"unit {unit_j} -> unit {unit_i}")
         coeff = coupling_coeff[unit_i, unit_j]
-        axs[unit_i, unit_j].plot(coupling_filter_bank[unit_i, unit_j], label="gamma difference")
+        axs[unit_i, unit_j].plot(coupling_filter_bank[:, unit_i, unit_j], label="gamma difference")
         axs[unit_i, unit_j].plot(np.dot(coupling_basis, coeff), ls="--", color="k", label="basis function")
 axs[0, 0].legend()
 plt.tight_layout()

src/nemos/simulation.py

@@ -1,20 +1,17 @@
 """Utility functions for coupling filter definition."""
 
-from typing import Tuple
 
 import numpy as np
 import scipy.stats as sts
 from numpy.typing import NDArray
 
-from .basis import Basis
-
 
 def difference_of_gammas(
     ws: int,
     upper_percentile: float = 0.99,
     inhib_a: float = 1.0,
     excit_a: float = 2.0,
-    inhib_b: float = 2.0,
+    inhib_b: float = 1.0,
     excit_b: float = 2.0,
 ) -> NDArray:
     r"""Generate coupling filter as a Gamma pdf difference.
@@ -49,11 +46,34 @@ def difference_of_gammas(
     filter:
         The coupling filter.
 
+    Raises
+    ------
+    ValueError:
+        - If any of the Gamma parameters is lesser or equal to 0.
+        - If the upper_percentile is not in [0, 1).
+
     References
     ----------
     1. [SciPy Docs - "scipy.stats.gamma"](https://docs.scipy.org/doc/
     scipy/reference/generated/scipy.stats.gamma.html)
     """
+    # check that the gamma parameters are positive (scipy returns
+    # nans otherwise but no exception is raised)
+    variables = {
+        "excit_a": excit_a,
+        "inhib_a": inhib_a,
+        "excit_b": excit_b,
+        "inhib_b": inhib_b,
+    }
+    for name, value in variables.items():
+        if value <= 0:
+            raise ValueError(f"Gamma parameter {name} must be >0.")
+    # check for valid pecentile
+    if upper_percentile < 0 or upper_percentile >= 1:
+        raise ValueError(
+            f"upper_percentile should lie in the [0, 1) interval. {upper_percentile} provided instead!"
+        )
+
     gm_inhibition = sts.gamma(a=inhib_a, scale=1 / inhib_b)
     gm_excitation = sts.gamma(a=excit_a, scale=1 / excit_b)
 
@@ -69,39 +89,58 @@ def difference_of_gammas(
     return gamma_diff
 
 
-def regress_filter(
-    coupling_filter_bank: NDArray, basis: Basis
-) -> Tuple[NDArray, NDArray]:
+def regress_filter(coupling_filters: NDArray, eval_basis: NDArray) -> NDArray:
     """Approximate scipy.stats.gamma based filters with basis function.
 
     Find the ols weights for representing the filters in terms of basis functions.
     This is done to re-use the nsl.glm.simulate method.
 
     Parameters
     ----------
-    coupling_filter_bank:
-        The coupling filters. Shape (n_neurons, n_neurons, window_size)
-    basis:
-        The basis function to instantiate.
+    coupling_filters:
+        The coupling filters. Shape (window_size, n_neurons_receiver, n_neurons_sender)
+    eval_basis:
+        The evaluated basis function, shape (window_size, n_basis_funcs)
 
     Returns
     -------
-    eval_basis:
-        The basis matrix, shape (window_size, n_basis_funcs)
     weights:
-        The weights for each neuron. Shape (n_neurons, n_neurons, n_basis_funcs)
+        The weights for each neuron. Shape (n_neurons_receiver, n_neurons_sender, n_basis_funcs)
     """
-    n_neurons, _, ws = coupling_filter_bank.shape
-    eval_basis = basis.evaluate(np.linspace(0, 1, ws))
-
-    # Reshape the coupling_filter_bank for vectorized least-squares
-    filters_reshaped = coupling_filter_bank.reshape(-1, ws)
+    # check shapes
+    if eval_basis.ndim != 2:
+        raise ValueError(
+            "eval_basis must be a 2 dimensional array, "
+            "shape (window_size, n_basis_funcs). "
+            f"{eval_basis.ndim} dimensional array provided instead!"
+        )
+    if coupling_filters.ndim != 3:
+        raise ValueError(
+            "coupling_filters must be a 3 dimensional array, "
+            "shape (window_size, n_neurons, n_neurons). "
+            f"{coupling_filters.ndim} dimensional array provided instead!"
+        )
+
+    ws, n_neurons_receiver, n_neurons_sender = coupling_filters.shape
+
+    # check that window size matches
+    if eval_basis.shape[0] != ws:
+        raise ValueError(
+            "window_size mismatch. The window size of coupling_filters and eval_basis "
+            f"does not match. coupling_filters has a window size of {ws}; "
+            f"eval_basis has a window size of {eval_basis.shape[0]}."
+        )
+
+    # Reshape the coupling_filters for vectorized least-squares
+    filters_reshaped = coupling_filters.reshape(ws, -1)
 
     # Solve the least squares problem for all filters at once
     # (vecotrizing the features)
-    weights = np.linalg.lstsq(eval_basis, filters_reshaped.T, rcond=None)[0]
+    weights = np.linalg.lstsq(eval_basis, filters_reshaped, rcond=None)[0]
 
     # Reshape back to the original dimensions
-    weights = weights.T.reshape(n_neurons, n_neurons, -1)
+    weights = np.transpose(
+        weights.reshape(-1, n_neurons_receiver, n_neurons_sender), axes=(1, 2, 0)
+    )
 
-    return eval_basis, weights
+    return weights

tests/test_simulation.py

@@ -0,0 +1,147 @@
+import itertools
+from contextlib import nullcontext as does_not_raise
+
+import numpy as np
+import pytest
+
+import nemos.simulation as simulation
+
+
+@pytest.mark.parametrize(
+        "inhib_a, expectation",
+        [
+            (-1, pytest.raises(ValueError, match="Gamma parameter [a-z]+_[a,b] must be >0.")),
+            (0, pytest.raises(ValueError, match="Gamma parameter [a-z]+_[a,b] must be >0.")),
+            (1, does_not_raise()),
+        ],
+    )
+def test_difference_of_gammas_inhib_a(inhib_a, expectation):
+    with expectation:
+        simulation.difference_of_gammas(10, inhib_a=inhib_a)
+
+
+@pytest.mark.parametrize(
+        "excit_a, expectation",
+        [
+            (-1, pytest.raises(ValueError, match="Gamma parameter [a-z]+_[a,b] must be >0.")),
+            (0, pytest.raises(ValueError, match="Gamma parameter [a-z]+_[a,b] must be >0.")),
+            (1, does_not_raise()),
+        ],
+    )
+def test_difference_of_gammas_excit_a(excit_a, expectation):
+    with expectation:
+        simulation.difference_of_gammas(10, excit_a=excit_a)
+
+
+@pytest.mark.parametrize(
+        "inhib_b, expectation",
+        [
+            (-1, pytest.raises(ValueError, match="Gamma parameter [a-z]+_[a,b] must be >0.")),
+            (0, pytest.raises(ValueError, match="Gamma parameter [a-z]+_[a,b] must be >0.")),
+            (1, does_not_raise()),
+        ],
+    )
+def test_difference_of_gammas_excit_a(inhib_b, expectation):
+    with expectation:
+        simulation.difference_of_gammas(10, inhib_b=inhib_b)
+
+
+@pytest.mark.parametrize(
+        "excit_b, expectation",
+        [
+            (-1, pytest.raises(ValueError, match="Gamma parameter [a-z]+_[a,b] must be >0.")),
+            (0, pytest.raises(ValueError, match="Gamma parameter [a-z]+_[a,b] must be >0.")),
+            (1, does_not_raise()),
+        ],
+    )
+def test_difference_of_gammas_excit_a(excit_b, expectation):
+    with expectation:
+        simulation.difference_of_gammas(10, excit_b=excit_b)
+
+
+@pytest.mark.parametrize(
+        "upper_percentile, expectation",
+        [
+            (-0.1, pytest.raises(ValueError, match=r"upper_percentile should lie in the \[0, 1\) interval.")),
+            (0, does_not_raise()),
+            (0.1, does_not_raise()),
+            (1, pytest.raises(ValueError, match=r"upper_percentile should lie in the \[0, 1\) interval.")),
+            (10, pytest.raises(ValueError, match=r"upper_percentile should lie in the \[0, 1\) interval.")),
+        ],
+    )
+def test_difference_of_gammas_percentile_params(upper_percentile, expectation):
+    with expectation:
+        simulation.difference_of_gammas(10, upper_percentile)
+
+
+@pytest.mark.parametrize("window_size", [0, 1, 2])
+def test_difference_of_gammas_output_shape(window_size):
+    result_size = simulation.difference_of_gammas(window_size).size
+    assert result_size == window_size, f"Expected output size {window_size}, but got {result_size}"
+
+
+@pytest.mark.parametrize("window_size", [1, 2, 10])
+def test_difference_of_gammas_output_norm(window_size):
+    result = simulation.difference_of_gammas(window_size)
+    assert np.allclose(np.linalg.norm(result, ord=2),1), "The output of difference_of_gammas is not unit norm."
+
+
+@pytest.mark.parametrize(
+        "coupling_filters, expectation",
+        [
+            (np.zeros((10, )), pytest.raises(ValueError, match=r"coupling_filters must be a 3 dimensional array")),
+            (np.zeros((10, 2)), pytest.raises(ValueError, match=r"coupling_filters must be a 3 dimensional array")),
+            (np.zeros((10, 2, 2)), does_not_raise()),
+            (np.zeros((10, 2, 2, 2)), pytest.raises(ValueError, match=r"coupling_filters must be a 3 dimensional array"))
+        ],
+    )
+def test_regress_filter_coupling_filters_dim(coupling_filters, expectation):
+    ws = coupling_filters.shape[0]
+    with expectation:
+        simulation.regress_filter(coupling_filters, np.zeros((ws, 3)))
+
+
+@pytest.mark.parametrize(
+        "eval_basis, expectation",
+        [
+            (np.zeros((10, )), pytest.raises(ValueError, match=r"eval_basis must be a 2 dimensional array")),
+            (np.zeros((10, 2)), does_not_raise()),
+            (np.zeros((10, 2, 2)), pytest.raises(ValueError, match=r"eval_basis must be a 2 dimensional array")),
+            (np.zeros((10, 2, 2, 2)), pytest.raises(ValueError, match=r"eval_basis must be a 2 dimensional array"))
+        ],
+    )
+def test_regress_filter_eval_basis_dim(eval_basis, expectation):
+    ws = eval_basis.shape[0]
+    with expectation:
+        simulation.regress_filter(np.zeros((ws, 1, 1)), eval_basis)
+
+
+@pytest.mark.parametrize(
+        "delta_ws, expectation",
+        [
+            (-1, pytest.raises(ValueError, match=r"window_size mismatch\. The window size of ")),
+            (0, does_not_raise()),
+            (1, pytest.raises(ValueError, match=r"window_size mismatch\. The window size of ")),
+        ],
+    )
+def test_regress_filter_window_size_matching(delta_ws, expectation):
+    ws = 2
+    with expectation:
+        simulation.regress_filter(np.zeros((ws, 1, 1)), np.zeros((ws + delta_ws, 1)))
+
+
+@pytest.mark.parametrize(
+        "window_size, n_neurons_sender, n_neurons_receiver, n_basis_funcs",
+        [x for x in itertools.product([1, 2], [1, 2], [1, 2], [1, 2])],
+    )
+def test_regress_filter_weights_size(window_size, n_neurons_sender, n_neurons_receiver, n_basis_funcs):
+    weights = simulation.regress_filter(
+        np.zeros((window_size, n_neurons_sender, n_neurons_receiver)),
+        np.zeros((window_size, n_basis_funcs))
+    )
+    assert weights.shape[0] == n_neurons_sender, (f"First dimension of weights (n_neurons_receiver) does not "
+                                                  f"match the second dimension of coupling_filters.")
+    assert weights.shape[1] == n_neurons_receiver, (f"Second dimension of weights (n_neuron_sender) does not "
+                                                    f"match the third dimension of coupling_filters.")
+    assert weights.shape[2] == n_basis_funcs, (f"Third dimension of weights (n_basis_funcs) does not "
+                                               f"match the second dimension of eval_basis.")