From 5e58d747fbb0afa5f945e6438710e95a4b2da7e6 Mon Sep 17 00:00:00 2001
From: BalzaniEdoardo <edoardo.balzani87@gmail.com>
Date: Fri, 13 Oct 2023 11:16:23 -0400
Subject: [PATCH 01/12] added fourier and tests

---
 src/neurostatslib/basis.py | 145 ++++++++++++++++++++++++++++++++-
 tests/test_basis.py        | 163 ++++++++++++++++++++++++++++++++++++-
 2 files changed, 304 insertions(+), 4 deletions(-)

diff --git a/src/neurostatslib/basis.py b/src/neurostatslib/basis.py
index ee7ff73e..b5c2fad4 100644
--- a/src/neurostatslib/basis.py
+++ b/src/neurostatslib/basis.py
@@ -23,6 +23,7 @@
     "OrthExponentialBasis",
     "AdditiveBasis",
     "MultiplicativeBasis",
+    "FourierBasis"
 ]
 
 
@@ -613,7 +614,6 @@ def _evaluate(self, sample_pts: NDArray) -> NDArray:
         The evaluation is performed by looping over each element and using `splev`
         from SciPy to compute the basis values.
         """
-
         # add knots
         knot_locs = self._generate_knots(sample_pts, 0.0, 1.0)
 
@@ -677,7 +677,6 @@ def _evaluate(self, sample_pts: NDArray) -> NDArray:
         The evaluation is performed by looping over each element and using `splev` from
         SciPy to compute the basis values.
         """
-
         knot_locs = self._generate_knots(sample_pts, 0.0, 1.0, is_cyclic=True)
 
         # for cyclic, do not repeat knots
@@ -941,7 +940,8 @@ def _check_rates(self):
                 "linearly dependent set of function for the basis."
             )
 
-    def _check_sample_range(self, sample_pts: NDArray):
+    @staticmethod
+    def _check_sample_range(sample_pts: NDArray):
         """
         Check if the sample points are all positive.
 
@@ -1009,6 +1009,145 @@ def _evaluate(self, sample_pts: NDArray) -> NDArray:
         )
 
 
+class FourierBasis(Basis):
+    """Set of 1D Fourier basis.
+
+    Parameters
+    ----------
+    n_freqs
+            Number of frequencies. The number of basis function will be 2*n_freqs - 1.
+    """
+
+    def __init__(self, n_freqs: int):
+        super().__init__(n_basis_funcs=2 * n_freqs - 1)
+
+        self._frequencies = np.arange(n_freqs, dtype=np.float32)
+        self._n_input_dimensionality = 1
+
+    def _check_n_basis_min(self) -> None:
+        """Check that the user required enough basis elements.
+
+        Checks that the number of basis is at least 1.
+
+        Raises
+        ------
+        ValueError
+            If an insufficient number of basis element is requested for the basis type
+        """
+        if self.n_basis_funcs < 1:
+            raise ValueError(
+                f"Object class {self.__class__.__name__} requires >= 1 basis elements. "
+                f"{self.n_basis_funcs} basis elements specified instead"
+            )
+
+    def _evaluate(self, sample_pts: NDArray) -> NDArray:
+        """Generate basis functions with given spacing.
+
+        Parameters
+        ----------
+        sample_pts
+            Spacing for basis functions.
+
+        Returns
+        -------
+        basis_funcs
+            Evaluated Fourier basis, shape (n_samples, n_basis_funcs).
+
+        Notes
+        -----
+            If the frequencies provided are np.arange(n_freq), convolving a signal
+            of length n_samples with this basis is equivalent, but slower,
+            then computing the FFT truncated to the first n_freq components.
+
+            Therefore, convolving a signal with this basis is equivalent
+            to compute the FFT over sliding window.
+
+            >>> import neurostatslib as nsl
+            >>> import numpy as np
+            >>> n_samples, n_freqs = 1000, 10
+            >>> basis = nsl.basis.FourierBasis(n_freqs*2)
+            >>> eval_basis = basis.evaluate(np.linspace(0, 1, n_samples))
+            >>> sinusoid = np.cos(3 * np.arange(0, 1000) * np.pi * 2 / 1000.)
+            >>> conv = [np.convolve(eval_basis[::-1, k], sinusoid, mode='valid')[0] for k in range(2*n_freqs-1)]
+            >>> fft = np.fft.fft(sinusoid)
+            >>> print('FFT power:   ', np.round(np.real(fft[:10]), 4))
+            >>> print('Convolution: ', np.round(conv[:10], 4))
+        """
+        # assumes equi-spaced samples.
+        if sample_pts.shape[0] / np.max(self._frequencies) < 2:
+            raise ValueError("Not enough samples, aliasing likely to occur!")
+
+        # rescale to [0, 2pi)
+        mn, mx = np.nanmin(sample_pts), np.nanmax(sample_pts)
+        # first sample in 0, last sample in 2 pi - 2 pi / n_samples.
+        sample_pts = (
+            2
+            * np.pi
+            * (sample_pts - mn)
+            / (mx - mn)
+            * (1.0 - 1.0 / sample_pts.shape[0])
+        )
+        # create the basis
+        angles = np.einsum("i,j->ij", sample_pts, self._frequencies)
+        return np.concatenate([np.cos(angles), -np.sin(angles[:, 1:])], axis=1)
+
+
+class BernsteinPolyBasis(Basis):
+    """Set of 1D Bernstein Polynomial basis.
+
+    Parameters
+    ----------
+    n_basis_funcs : int
+        Number of basis functions.
+    """
+
+    def __init__(self, n_basis_funcs: int):
+        super().__init__(n_basis_funcs=n_basis_funcs)
+        self._n_input_dimensionality = 1
+
+    def _check_n_basis_min(self) -> None:
+        """Check that the user required enough basis elements.
+
+        Checks that the number of basis is at least 1.
+
+        Raises
+        ------
+        ValueError
+            If an insufficient number of basis element is requested for the basis type
+        """
+        if self.n_basis_funcs < 1:
+            raise ValueError(
+                f"Object class {self.__class__.__name__} requires >= 1 basis elements. "
+                f"{self.n_basis_funcs} basis elements specified instead"
+            )
+
+    def _evaluate(self, sample_pts: NDArray) -> NDArray:
+        """Generate basis functions with given spacing.
+
+        Parameters
+        ----------
+        sample_pts
+            Spacing for basis functions.
+
+        Returns
+        -------
+        :
+            Evaluated Bernstein polynomial basis, shape (n_samples, n_basis_funcs).
+
+        Examples
+        --------
+        >>> import neurostatslib as nsl
+        >>> import numpy as np
+        >>> bernstein_basis = BernsteinPolyBasis(n_basis_funcs=5)
+        >>> sample_pts = np.linspace(0, 1, 10)
+        >>> basis_funcs = bernstein_basis._evaluate(sample_pts)
+        >>> basis_funcs.shape
+        (10, 5)
+
+        """
+        return bernstein_poly(sample_pts, self.n_basis_funcs, der=0)
+
+
 def mspline(x: NDArray, k: int, i: int, T: NDArray):
     """Compute M-spline basis function.
 
diff --git a/tests/test_basis.py b/tests/test_basis.py
index 7f31fe2e..f43ee173 100644
--- a/tests/test_basis.py
+++ b/tests/test_basis.py
@@ -749,6 +749,165 @@ def test_decay_rate_size_match_n_basis_func(self, decay_rates, n_basis_func):
             self.cls(n_basis_funcs=n_basis_func, decay_rates=decay_rates)
 
 
+class TestFourierBasis(BasisFuncsTesting):
+    cls = basis.FourierBasis
+
+    @pytest.mark.parametrize("n_freqs", [2, 4, 8])
+    @pytest.mark.parametrize("sample_size", [20, 1000])
+    def test_evaluate_returns_expected_number_of_basis(
+            self, n_freqs, sample_size
+    ):
+        """Tests whether the evaluate method returns the expected number of basis functions."""
+        basis_obj = self.cls(n_freqs=n_freqs)
+        eval_basis = basis_obj.evaluate(np.linspace(0, 1, sample_size))
+        assert(eval_basis.shape[1] == 2*n_freqs-1)
+
+    @pytest.mark.parametrize("samples", [[], np.zeros(11)])
+    def test_non_empty_samples(self, samples):
+        if len(samples) == 0:
+            with pytest.raises(
+                ValueError, match="All sample provided must be non empty"
+            ):
+                self.cls(5).evaluate(samples)
+        else:
+            self.cls(5).evaluate(samples)
+
+    @pytest.mark.parametrize(
+        "arraylike", [0, [0]*11, (0,)*11, np.array([0]*11), jax.numpy.array([0]*11)]
+    )
+    def test_input_to_evaluate_is_arraylike(self, arraylike):
+        """
+        Checks that the sample size of the output from the evaluate() method matches the input sample size.
+        """
+        basis_obj = self.cls(n_freqs=5)
+        raise_exception = not isinstance(
+            arraylike, (tuple, list, np.ndarray, jax.numpy.ndarray)
+        )
+        if raise_exception:
+            with pytest.raises(TypeError, match="Input samples must be array-like"):
+                basis_obj.evaluate(arraylike)
+        else:
+            basis_obj.evaluate(arraylike)
+
+    @pytest.mark.parametrize("sample_size", [100, 1000])
+    @pytest.mark.parametrize("n_freqs", [4, 10])
+    def test_sample_size_of_evaluate_matches_that_of_input(
+            self, n_freqs, sample_size
+    ):
+        """
+        Checks that the sample size of the output from the evaluate() method matches the input sample size.
+        """
+        basis_obj = self.cls(n_freqs)
+        eval_basis = basis_obj.evaluate(np.linspace(0, 1, sample_size))
+        if eval_basis.shape[0] != sample_size:
+            raise ValueError(
+                f"Dimensions do not agree: The window size should match the second dimension of the evaluated basis."
+                f"The window size is {sample_size}",
+                f"The second dimension of the evaluated basis is {eval_basis.shape[0]}",
+            )
+
+    @pytest.mark.parametrize("n_freqs", [-1, 0, 1, 3])
+    def test_minimum_number_of_basis_required_is_matched(self, n_freqs):
+        """
+        Verifies that the minimum number of basis functions and order required (i.e., at least 1) and
+        order < #basis are enforced.
+        """
+        n_samples = 10
+        raise_exception = n_freqs < 1
+        if raise_exception:
+            with pytest.raises(
+                    ValueError,
+                    match=rf"Object class {self.cls.__name__} requires >= 1 basis elements",
+            ):
+                basis_obj = self.cls(n_freqs=n_freqs)
+                basis_obj.evaluate(np.linspace(0, 1, n_samples))
+        else:
+            basis_obj = self.cls(n_freqs=n_freqs)
+            basis_obj.evaluate(np.linspace(0, 1, n_samples))
+
+    @pytest.mark.parametrize("n_freqs", [3, 6, 7, 10])
+    def test_minimum_aliasing_detection(self, n_freqs):
+        """
+        Verifies that the minimum number of basis functions and order required (i.e., at least 1) and
+        order < #basis are enforced.
+        """
+        n_samples = 10
+        # highest freq * 2 <= n_samples to avoid aliasing.
+        raise_exception = (n_freqs - 1) * 2 > n_samples
+        if raise_exception:
+            with pytest.raises(
+                    ValueError,
+                    match=rf"Not enough samples, aliasing likely to occur",
+            ):
+                basis_obj = self.cls(n_freqs=n_freqs)
+                basis_obj.evaluate(np.linspace(0, 1, n_samples))
+        else:
+            basis_obj = self.cls(n_freqs=n_freqs)
+            basis_obj.evaluate(np.linspace(0, 1, n_samples))
+
+    @pytest.mark.parametrize(
+        "sample_range", [(0, 1), (0.1, 0.9), (-0.5, 1), (0, 1.5), (-0.5, 1.5)]
+    )
+    def test_samples_range_matches_evaluate_requirements(self, sample_range: tuple):
+        """
+        Verifies that the evaluate() method can handle input range.
+        """
+        basis_obj = self.cls(n_freqs=5)
+        basis_obj.evaluate(np.linspace(*sample_range, 100))
+
+    @pytest.mark.parametrize("n_input", [0, 1, 2, 3])
+    def test_number_of_required_inputs_evaluate(self, n_input):
+        """
+        Confirms that the evaluate() method correctly handles the number of input samples that are provided.
+        """
+        basis_obj = self.cls(n_freqs=5)
+        raise_exception = n_input != basis_obj._n_input_dimensionality
+        inputs = [np.linspace(0, 1, 20)] * n_input
+        if raise_exception:
+            with pytest.raises(
+                    ValueError,
+                    match="Input dimensionality mismatch. This basis evaluation requires [0-9]+ inputs,",
+            ):
+                basis_obj.evaluate(*inputs)
+        else:
+            basis_obj.evaluate(*inputs)
+
+    @pytest.mark.parametrize("sample_size", [-1, 0, 10, 11, 100])
+    def test_evaluate_on_grid_meshgrid_size(self, sample_size):
+        """
+        Checks that the evaluate_on_grid() method returns a grid of the expected size.
+        """
+        basis_obj = self.cls(n_freqs=5)
+        raise_exception = sample_size <= 0
+        if raise_exception:
+            with pytest.raises(
+                    ValueError,
+                    match=r"Invalid input data|"
+                          rf"All sample counts provided must be greater",
+            ):
+                basis_obj.evaluate_on_grid(sample_size)
+        else:
+            grid, _ = basis_obj.evaluate_on_grid(sample_size)
+            assert grid.shape[0] == sample_size
+
+    @pytest.mark.parametrize("n_input", [0, 1, 2])
+    def test_evaluate_on_grid_input_number(self, n_input):
+        """
+        Validates that the evaluate_on_grid() method correctly handles the number of input samples that are provided.
+        """
+        basis_obj = self.cls(n_freqs=5)
+        inputs = [10] * n_input
+        raise_exception = n_input != basis_obj._n_input_dimensionality
+        if raise_exception:
+            with pytest.raises(
+                    ValueError,
+                    match=r"Input dimensionality mismatch\. This basis evaluation requires [0-9]+ inputs",
+            ):
+                basis_obj.evaluate_on_grid(*inputs)
+        else:
+            basis_obj.evaluate_on_grid(*inputs)
+
+
 class TestBSplineBasis(BasisFuncsTesting):
     cls = basis.BSplineBasis
 
@@ -1152,6 +1311,8 @@ def instantiate_basis(n_basis, basis_class):
             basis_obj = basis_class(
                 n_basis_funcs=n_basis, decay_rates=np.arange(1, 1 + n_basis)
             )
+        elif basis_class == basis.FourierBasis:
+            basis_obj = basis_class(n_freqs=n_basis)
         elif basis_class == basis.BSplineBasis:
             basis_obj = basis_class(n_basis_funcs=n_basis, order=3)
         elif basis_class == basis.CyclicBSplineBasis:
@@ -1472,7 +1633,7 @@ def test_evaluate_returns_expected_number_of_basis(
                 f"The first dimension of the evaluated basis is {eval_basis.shape[1]}",
             )
 
-    @pytest.mark.parametrize("sample_size", [6, 30, 35])
+    @pytest.mark.parametrize("sample_size", [12, 30, 35])
     @pytest.mark.parametrize("n_basis_a", [5, 6])
     @pytest.mark.parametrize("n_basis_b", [5, 6])
     @pytest.mark.parametrize(

From 290f223c0a6753c710ef556a0c1deb986172b3b9 Mon Sep 17 00:00:00 2001
From: BalzaniEdoardo <edoardo.balzani87@gmail.com>
Date: Fri, 13 Oct 2023 11:28:06 -0400
Subject: [PATCH 02/12] fixed basis note

---
 docs/developers_notes/01-basis_module.md | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/docs/developers_notes/01-basis_module.md b/docs/developers_notes/01-basis_module.md
index dc3e43ce..49f9cd8f 100644
--- a/docs/developers_notes/01-basis_module.md
+++ b/docs/developers_notes/01-basis_module.md
@@ -25,7 +25,9 @@ Abstract Class Basis
 │   │
 │   └─ Concrete Subclass RaisedCosineBasisLog
 │
-└─ Concrete Subclass OrthExponentialBasis
+├─ Concrete Subclass OrthExponentialBasis
+│
+└─ Concrete Subclass FourierBasis
 ```
 
 The super-class `Basis` provides two public methods, [`evaluate`](#the-public-method-evaluate) and [`evaluate_on_grid`](#the-public-method-evaluate_on_grid). These methods perform checks on both the input provided by the user and the output of the evaluation to ensure correctness, and are thus considered "safe". They both make use of the private abstract method `_evaluate` that is specific for each concrete class. See below for more details.

From 0a576740dfc66a428ee4be83226013c4ffa83786 Mon Sep 17 00:00:00 2001
From: BalzaniEdoardo <edoardo.balzani87@gmail.com>
Date: Fri, 13 Oct 2023 11:29:04 -0400
Subject: [PATCH 03/12] removed bp

---
 src/neurostatslib/basis.py | 58 +-------------------------------------
 1 file changed, 1 insertion(+), 57 deletions(-)

diff --git a/src/neurostatslib/basis.py b/src/neurostatslib/basis.py
index b5c2fad4..9e55b155 100644
--- a/src/neurostatslib/basis.py
+++ b/src/neurostatslib/basis.py
@@ -23,7 +23,7 @@
     "OrthExponentialBasis",
     "AdditiveBasis",
     "MultiplicativeBasis",
-    "FourierBasis"
+    "FourierBasis",
 ]
 
 
@@ -1092,62 +1092,6 @@ def _evaluate(self, sample_pts: NDArray) -> NDArray:
         return np.concatenate([np.cos(angles), -np.sin(angles[:, 1:])], axis=1)
 
 
-class BernsteinPolyBasis(Basis):
-    """Set of 1D Bernstein Polynomial basis.
-
-    Parameters
-    ----------
-    n_basis_funcs : int
-        Number of basis functions.
-    """
-
-    def __init__(self, n_basis_funcs: int):
-        super().__init__(n_basis_funcs=n_basis_funcs)
-        self._n_input_dimensionality = 1
-
-    def _check_n_basis_min(self) -> None:
-        """Check that the user required enough basis elements.
-
-        Checks that the number of basis is at least 1.
-
-        Raises
-        ------
-        ValueError
-            If an insufficient number of basis element is requested for the basis type
-        """
-        if self.n_basis_funcs < 1:
-            raise ValueError(
-                f"Object class {self.__class__.__name__} requires >= 1 basis elements. "
-                f"{self.n_basis_funcs} basis elements specified instead"
-            )
-
-    def _evaluate(self, sample_pts: NDArray) -> NDArray:
-        """Generate basis functions with given spacing.
-
-        Parameters
-        ----------
-        sample_pts
-            Spacing for basis functions.
-
-        Returns
-        -------
-        :
-            Evaluated Bernstein polynomial basis, shape (n_samples, n_basis_funcs).
-
-        Examples
-        --------
-        >>> import neurostatslib as nsl
-        >>> import numpy as np
-        >>> bernstein_basis = BernsteinPolyBasis(n_basis_funcs=5)
-        >>> sample_pts = np.linspace(0, 1, 10)
-        >>> basis_funcs = bernstein_basis._evaluate(sample_pts)
-        >>> basis_funcs.shape
-        (10, 5)
-
-        """
-        return bernstein_poly(sample_pts, self.n_basis_funcs, der=0)
-
-
 def mspline(x: NDArray, k: int, i: int, T: NDArray):
     """Compute M-spline basis function.
 

From 8ebf2dc430c736330f9b3ec380401d41a1978bdf Mon Sep 17 00:00:00 2001
From: BalzaniEdoardo <edoardo.balzani87@gmail.com>
Date: Fri, 13 Oct 2023 11:30:02 -0400
Subject: [PATCH 04/12] linted

---
 docs/examples/plot_1D_basis_function.py | 87 ++++++++++++++++++++++++-
 1 file changed, 85 insertions(+), 2 deletions(-)

diff --git a/docs/examples/plot_1D_basis_function.py b/docs/examples/plot_1D_basis_function.py
index 2344d0ea..4af70b40 100644
--- a/docs/examples/plot_1D_basis_function.py
+++ b/docs/examples/plot_1D_basis_function.py
@@ -65,8 +65,10 @@
 # -----------------
 # Each basis type may necessitate specific hyperparameters for instantiation. For a comprehensive description,
 # please refer to the  [Code References](../../../reference/neurostatslib/basis). After instantiation, all classes
-# share the same syntax for basis evaluation. The following is an example of how to instantiate and
-# evaluate a log-spaced cosine raised function basis.
+# share the same syntax for basis evaluation.
+#
+# ### The Log-Spaced Raised Cosine Basis
+# The following is an example of how to instantiate and  evaluate a log-spaced cosine raised function basis.
 
 # Instantiate the basis noting that the `RaisedCosineBasisLog` does not require an `order` parameter
 raised_cosine_log = nsl.basis.RaisedCosineBasisLog(n_basis_funcs=10)
@@ -81,3 +83,84 @@
 plt.plot(samples, eval_basis)
 plt.show()
 
+# %%
+# ### The Fourier Basis
+# Another type of basis available is the Fourier Basis. Fourier basis are ideal to capture periodic and
+# quasi-periodic patterns. Such oscillatory, rhythmic behavior is a common signature of many neural signals.
+# Additionally, the Fourier basis has the advantage of being orthogonal, which simplifies the estimation and
+# interpretation of the model parameters, each of which will represent the relative contribution of a specific
+# oscillation frequency to the overall signal.
+
+
+# A Fourier basis can be instantiated with the usual syntax.
+# The user can pass the desired frequencies for the basis or
+# the frequencies will be set to np.arange(n_basis_funcs//2).
+# The number of basis function is required to be even.
+fourier_basis = nsl.basis.FourierBasis(n_freqs=4)
+
+# evaluate on equi-spaced samples
+samples, eval_basis = fourier_basis.evaluate_on_grid(1000)
+
+# plot the `sin` and `cos` separately
+plt.figure(figsize=(6, 3))
+plt.subplot(121)
+plt.title("Cos")
+plt.plot(samples, eval_basis[:, :4])
+plt.subplot(122)
+plt.title("Sin")
+plt.plot(samples, eval_basis[:, 4:])
+plt.tight_layout()
+
+# %%
+# !!! note "Fourier basis convolution and Fourier transform"
+#     The Fourier transform of a signal $ s(t) $ restricted to a temporal window $ [t_0,\;t_1] $ is
+#     $$ \\hat{x}(\\omega) = \\int_{t_0}^{t_1} s(\\tau) e^{-j\\omega \\tau} d\\tau. $$
+#     where $ e^{-j\\omega \\tau} = \\cos(\\omega \\tau) - j \\sin (\\omega \\tau) $.
+#
+#     When computing the cross-correlation of a signal with the Fourier basis functions,
+#     we essentially measure how well the signal correlates with sinusoids of different frequencies,
+#     within a specified temporal window. This process mirrors the operation performed by the Fourier transform.
+#     Therefore, it becomes clear that computing the cross-correlation of a signal with the Fourier basis defined here
+#     is equivalent to computing the discrete Fourier transform on a sliding window of the same size
+#     as that of the basis.
+
+n_samples = 1000
+n_freqs = 20
+
+# define a signal
+signal = np.random.normal(size=n_samples)
+
+# evaluate the basis
+_, eval_basis = nsl.basis.FourierBasis(n_freqs=n_freqs).evaluate_on_grid(n_samples)
+
+# compute the cross-corr with the signal and the basis
+# Note that we are inverting the time axis of the basis because we are aiming
+# for a cross-correlation, while np.convolve compute a convolution which would flip the time axis.
+xcorr = np.array(
+    [
+        np.convolve(eval_basis[::-1, k], signal, mode="valid")[0]
+        for k in range(2 * n_freqs - 1)
+    ]
+)
+
+# compute the power (add back sin(0 * t) = 0)
+fft_complex = np.fft.fft(signal)
+fft_amplitude = np.abs(fft_complex[:n_freqs])
+fft_phase = np.angle(fft_complex[:n_freqs])
+# compute the phase and amplitude from the convolution
+xcorr_phase = np.arctan2(np.hstack([[0], xcorr[n_freqs:]]), xcorr[:n_freqs])
+xcorr_aplitude = np.sqrt(xcorr[:n_freqs] ** 2 + np.hstack([[0], xcorr[n_freqs:]]) ** 2)
+
+fig, ax = plt.subplots(1, 2)
+ax[0].set_aspect("equal")
+ax[0].set_title("Signal amplitude")
+ax[0].scatter(fft_amplitude, xcorr_aplitude)
+ax[0].set_xlabel("FFT")
+ax[0].set_ylabel("cross-correlation")
+
+ax[1].set_aspect("equal")
+ax[1].set_title("Signal phase")
+ax[1].scatter(fft_phase, xcorr_phase)
+ax[1].set_xlabel("FFT")
+ax[1].set_ylabel("cross-correlation")
+plt.tight_layout()

From 497b31f075feb3390d263985123d869415ebb741 Mon Sep 17 00:00:00 2001
From: BalzaniEdoardo <edoardo.balzani87@gmail.com>
Date: Tue, 19 Dec 2023 12:41:46 -0500
Subject: [PATCH 05/12] fixed tests

---
 src/nemos/__init__.py |  10 +--
 src/nemos/basis.py    |   2 +-
 tests/test_basis.py   | 157 +++++++++++++++++++++++-------------------
 3 files changed, 87 insertions(+), 82 deletions(-)

diff --git a/src/nemos/__init__.py b/src/nemos/__init__.py
index 28836560..8bcee80f 100644
--- a/src/nemos/__init__.py
+++ b/src/nemos/__init__.py
@@ -1,11 +1,3 @@
 #!/usr/bin/env python3
 
-from . import (
-    basis,
-    exceptions,
-    glm,
-    observation_models,
-    regularizer,
-    simulation,
-    utils,
-)
+from . import basis, exceptions, glm, observation_models, regularizer, simulation, utils
diff --git a/src/nemos/basis.py b/src/nemos/basis.py
index 0ab22c6b..f0b8b536 100644
--- a/src/nemos/basis.py
+++ b/src/nemos/basis.py
@@ -104,7 +104,7 @@ def _check_evaluate_input(self, *xi: ArrayLike) -> Tuple[NDArray]:
             # make sure array is at least 1d (so that we succeed when only
             # passed a scalar)
             xi = tuple(np.atleast_1d(np.asarray(x, dtype=float)) for x in xi)
-        except TypeError:
+        except (TypeError, ValueError):
             raise TypeError("Input samples must be array-like of floats!")
 
         # check for non-empty samples
diff --git a/tests/test_basis.py b/tests/test_basis.py
index 3ee4cfa5..826c57c9 100644
--- a/tests/test_basis.py
+++ b/tests/test_basis.py
@@ -730,31 +730,31 @@ def test_evaluate_returns_expected_number_of_basis(
         eval_basis = basis_obj.evaluate(np.linspace(0, 1, sample_size))
         assert(eval_basis.shape[1] == 2*n_freqs-1)
 
-    @pytest.mark.parametrize("samples", [[], np.zeros(11)])
-    def test_non_empty_samples(self, samples):
-        if len(samples) == 0:
-            with pytest.raises(
-                ValueError, match="All sample provided must be non empty"
-            ):
-                self.cls(5).evaluate(samples)
-        else:
+    @pytest.mark.parametrize("samples, expectation",
+        [
+            ([], pytest.raises(ValueError, match="All sample provided must be non empty")),
+            (np.zeros(11), does_not_raise())
+        ]
+    )
+    def test_non_empty_samples(self, samples, expectation):
+        with expectation:
             self.cls(5).evaluate(samples)
 
     @pytest.mark.parametrize(
-        "arraylike", [0, [0]*11, (0,)*11, np.array([0]*11), jax.numpy.array([0]*11)]
+        "arraylike, expectation", [
+            (["x"], pytest.raises(TypeError, match="Input samples must be array-like of floats")),
+            ([0]*11, does_not_raise()),
+            ((0,)*11, does_not_raise()),
+            (np.array([0]*11), does_not_raise()),
+            (jax.numpy.array([0]*11), does_not_raise())
+        ]
     )
-    def test_input_to_evaluate_is_arraylike(self, arraylike):
+    def test_input_to_evaluate_is_arraylike(self, arraylike, expectation):
         """
         Checks that the sample size of the output from the evaluate() method matches the input sample size.
         """
         basis_obj = self.cls(n_freqs=5)
-        raise_exception = not isinstance(
-            arraylike, (tuple, list, np.ndarray, jax.numpy.ndarray)
-        )
-        if raise_exception:
-            with pytest.raises(TypeError, match="Input samples must be array-like"):
-                basis_obj.evaluate(arraylike)
-        else:
+        with expectation:
             basis_obj.evaluate(arraylike)
 
     @pytest.mark.parametrize("sample_size", [100, 1000])
@@ -774,43 +774,40 @@ def test_sample_size_of_evaluate_matches_that_of_input(
                 f"The second dimension of the evaluated basis is {eval_basis.shape[0]}",
             )
 
-    @pytest.mark.parametrize("n_freqs", [-1, 0, 1, 3])
-    def test_minimum_number_of_basis_required_is_matched(self, n_freqs):
+    @pytest.mark.parametrize("n_freqs, expectation",
+        [
+            (-1, pytest.raises(ValueError,match=r"Object class FourierBasis requires >= 1 basis elements")),
+            (0, pytest.raises(ValueError,match=r"Object class FourierBasis requires >= 1 basis elements")),
+            (1, does_not_raise()),
+            (3, does_not_raise())
+        ]
+    )
+    def test_minimum_number_of_basis_required_is_matched(self, n_freqs, expectation):
         """
         Verifies that the minimum number of basis functions and order required (i.e., at least 1) and
         order < #basis are enforced.
         """
         n_samples = 10
-        raise_exception = n_freqs < 1
-        if raise_exception:
-            with pytest.raises(
-                    ValueError,
-                    match=rf"Object class {self.cls.__name__} requires >= 1 basis elements",
-            ):
-                basis_obj = self.cls(n_freqs=n_freqs)
-                basis_obj.evaluate(np.linspace(0, 1, n_samples))
-        else:
+        with expectation:
             basis_obj = self.cls(n_freqs=n_freqs)
             basis_obj.evaluate(np.linspace(0, 1, n_samples))
 
-    @pytest.mark.parametrize("n_freqs", [3, 6, 7, 10])
-    def test_minimum_aliasing_detection(self, n_freqs):
+    @pytest.mark.parametrize("n_freqs, expectation",
+                             [
+                                 (3, does_not_raise()),
+                                 (6, does_not_raise()),
+                                 (7, pytest.raises(ValueError,match=rf"Not enough samples, aliasing likely to occur")),
+                                 (10, pytest.raises(ValueError,match=rf"Not enough samples, aliasing likely to occur"))
+                             ]
+                             )
+    def test_minimum_aliasing_detection(self, n_freqs, expectation):
         """
         Verifies that the minimum number of basis functions and order required (i.e., at least 1) and
         order < #basis are enforced.
         """
         n_samples = 10
-        # highest freq * 2 <= n_samples to avoid aliasing.
-        raise_exception = (n_freqs - 1) * 2 > n_samples
-        if raise_exception:
-            with pytest.raises(
-                    ValueError,
-                    match=rf"Not enough samples, aliasing likely to occur",
-            ):
-                basis_obj = self.cls(n_freqs=n_freqs)
-                basis_obj.evaluate(np.linspace(0, 1, n_samples))
-        else:
-            basis_obj = self.cls(n_freqs=n_freqs)
+        basis_obj = self.cls(n_freqs=n_freqs)
+        with expectation:
             basis_obj.evaluate(np.linspace(0, 1, n_samples))
 
     @pytest.mark.parametrize(
@@ -823,56 +820,72 @@ def test_samples_range_matches_evaluate_requirements(self, sample_range: tuple):
         basis_obj = self.cls(n_freqs=5)
         basis_obj.evaluate(np.linspace(*sample_range, 100))
 
-    @pytest.mark.parametrize("n_input", [0, 1, 2, 3])
-    def test_number_of_required_inputs_evaluate(self, n_input):
+    @pytest.mark.parametrize("n_input, expectation", [
+         (0, pytest.raises(TypeError, match=r"FourierBasis\.evaluate\(\) missing 1 required positional argument")),
+         (1, does_not_raise()),
+         (2, pytest.raises(TypeError, match=r"FourierBasis\.evaluate\(\) takes 2 positional arguments but")),
+         (3, pytest.raises(TypeError, match=r"FourierBasis\.evaluate\(\) takes 2 positional arguments but"))
+        ]
+    )
+    def test_number_of_required_inputs_evaluate(self, n_input, expectation):
         """
         Confirms that the evaluate() method correctly handles the number of input samples that are provided.
         """
         basis_obj = self.cls(n_freqs=5)
-        raise_exception = n_input != basis_obj._n_input_dimensionality
         inputs = [np.linspace(0, 1, 20)] * n_input
-        if raise_exception:
-            with pytest.raises(
-                    ValueError,
-                    match="Input dimensionality mismatch. This basis evaluation requires [0-9]+ inputs,",
-            ):
-                basis_obj.evaluate(*inputs)
-        else:
+        with expectation:
             basis_obj.evaluate(*inputs)
 
-    @pytest.mark.parametrize("sample_size", [-1, 0, 10, 11, 100])
-    def test_evaluate_on_grid_meshgrid_size(self, sample_size):
+    @pytest.mark.parametrize(
+        "sample_size, expectation",
+        [
+            (-1, pytest.raises(ValueError, match=r"Invalid input data|All sample counts provided must be greater")),
+            (0, pytest.raises(ValueError, match=r"Invalid input data|All sample counts provided must be greater")),
+            (10, does_not_raise()),
+            (11, does_not_raise()),
+            (100, does_not_raise())
+        ]
+    )
+    def test_evaluate_on_grid_meshgrid_valid_size(self, sample_size, expectation):
         """
         Checks that the evaluate_on_grid() method returns a grid of the expected size.
         """
         basis_obj = self.cls(n_freqs=5)
-        raise_exception = sample_size <= 0
-        if raise_exception:
-            with pytest.raises(
-                    ValueError,
-                    match=r"Invalid input data|"
-                          rf"All sample counts provided must be greater",
-            ):
-                basis_obj.evaluate_on_grid(sample_size)
-        else:
+        with expectation:
+            basis_obj.evaluate_on_grid(sample_size)
+
+    @pytest.mark.parametrize(
+        "sample_size, expectation",
+        [
+            (10, does_not_raise()),
+            (11, does_not_raise()),
+            (100, does_not_raise())
+        ]
+    )
+    def test_evaluate_on_grid_meshgrid_match_size(self, sample_size, expectation):
+        """
+        Checks that the evaluate_on_grid() method returns a grid of the expected size.
+        """
+        basis_obj = self.cls(n_freqs=5)
+        with expectation:
             grid, _ = basis_obj.evaluate_on_grid(sample_size)
             assert grid.shape[0] == sample_size
 
-    @pytest.mark.parametrize("n_input", [0, 1, 2])
-    def test_evaluate_on_grid_input_number(self, n_input):
+    @pytest.mark.parametrize(
+        "n_input, expectation",
+        [
+            (0, pytest.raises(TypeError, match="Input dimensionality mismatch. This basis evaluation")),
+             (1, does_not_raise()),
+            (2, pytest.raises(TypeError, match="Input dimensionality mismatch. This basis evaluation"))
+        ]
+    )
+    def test_evaluate_on_grid_input_number(self, n_input, expectation):
         """
         Validates that the evaluate_on_grid() method correctly handles the number of input samples that are provided.
         """
         basis_obj = self.cls(n_freqs=5)
         inputs = [10] * n_input
-        raise_exception = n_input != basis_obj._n_input_dimensionality
-        if raise_exception:
-            with pytest.raises(
-                    ValueError,
-                    match=r"Input dimensionality mismatch\. This basis evaluation requires [0-9]+ inputs",
-            ):
-                basis_obj.evaluate_on_grid(*inputs)
-        else:
+        with expectation:
             basis_obj.evaluate_on_grid(*inputs)
 
 

From 332d55510dca655aafdcef71f80e4ee003274f43 Mon Sep 17 00:00:00 2001
From: Edoardo Balzani <edoardo.balzani87@gmail.com>
Date: Tue, 19 Dec 2023 17:17:12 -0500
Subject: [PATCH 06/12] Update docs/examples/plot_1D_basis_function.py

Co-authored-by: William F. Broderick <billbrod@gmail.com>
---
 docs/examples/plot_1D_basis_function.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/docs/examples/plot_1D_basis_function.py b/docs/examples/plot_1D_basis_function.py
index 46cdf41a..5e96772e 100644
--- a/docs/examples/plot_1D_basis_function.py
+++ b/docs/examples/plot_1D_basis_function.py
@@ -120,7 +120,7 @@
 #     When computing the cross-correlation of a signal with the Fourier basis functions,
 #     we essentially measure how well the signal correlates with sinusoids of different frequencies,
 #     within a specified temporal window. This process mirrors the operation performed by the Fourier transform.
-#     Therefore, it becomes clear that computing the cross-correlation of a signal with the Fourier basis defined here
+#     Therefore, computing the cross-correlation of a signal with the Fourier basis defined here
 #     is equivalent to computing the discrete Fourier transform on a sliding window of the same size
 #     as that of the basis.
 

From 17a05a97a1a443906005cc3d7d47412bda37a188 Mon Sep 17 00:00:00 2001
From: BalzaniEdoardo <edoardo.balzani87@gmail.com>
Date: Tue, 19 Dec 2023 17:48:47 -0500
Subject: [PATCH 07/12] changed basis

---
 docs/developers_notes/01-basis_module.md |   4 +-
 docs/examples/plot_1D_basis_function.py  |  30 +-
 src/nemos/basis.py                       | 333 ++++++++++++++---------
 3 files changed, 215 insertions(+), 152 deletions(-)

diff --git a/docs/developers_notes/01-basis_module.md b/docs/developers_notes/01-basis_module.md
index 57dd3c86..b8d8a925 100644
--- a/docs/developers_notes/01-basis_module.md
+++ b/docs/developers_notes/01-basis_module.md
@@ -19,9 +19,7 @@ Abstract Class Basis
 │   │
 │   └─ Concrete Subclass CyclicBSplineBasis
 │
-├─ Abstract Subclass RaisedCosineBasis
-│   │
-│   ├─ Concrete Subclass RaisedCosineBasisLinear
+├─ Concrete Subclass RaisedCosineBasisLinear 
 │   │
 │   └─ Concrete Subclass RaisedCosineBasisLog
 │
diff --git a/docs/examples/plot_1D_basis_function.py b/docs/examples/plot_1D_basis_function.py
index 46cdf41a..7d307916 100644
--- a/docs/examples/plot_1D_basis_function.py
+++ b/docs/examples/plot_1D_basis_function.py
@@ -90,13 +90,13 @@
 # Additionally, the Fourier basis has the advantage of being orthogonal, which simplifies the estimation and
 # interpretation of the model parameters, each of which will represent the relative contribution of a specific
 # oscillation frequency to the overall signal.
-
-
+#
 # A Fourier basis can be instantiated with the usual syntax.
 # The user can pass the desired frequencies for the basis or
-# the frequencies will be set to np.arange(n_basis_funcs//2).
+# the frequencies will be set to `np.arange(n_basis_funcs//2)`.
 # The number of basis function is required to be even.
-fourier_basis = nmo.basis.FourierBasis(n_freqs=4)
+
+fourier_basis = nmo.basis.FourierBasis(max_freq=4)
 
 # evaluate on equi-spaced samples
 samples, eval_basis = fourier_basis.evaluate_on_grid(1000)
@@ -112,17 +112,17 @@
 plt.tight_layout()
 
 # %%
-# !!! note "Fourier basis convolution and Fourier transform"
-#     The Fourier transform of a signal $ s(t) $ restricted to a temporal window $ [t_0,\;t_1] $ is
-#     $$ \\hat{x}(\\omega) = \\int_{t_0}^{t_1} s(\\tau) e^{-j\\omega \\tau} d\\tau. $$
-#     where $ e^{-j\\omega \\tau} = \\cos(\\omega \\tau) - j \\sin (\\omega \\tau) $.
+# ## Fourier Basis Convolution and Fourier Transform
+# The Fourier transform of a signal $ s(t) $ restricted to a temporal window $ [t_0,\;t_1] $ is
+# $$ \\hat{x}(\\omega) = \\int_{t_0}^{t_1} s(\\tau) e^{-j\\omega \\tau} d\\tau. $$
+# where $ e^{-j\\omega \\tau} = \\cos(\\omega \\tau) - j \\sin (\\omega \\tau) $.
 #
-#     When computing the cross-correlation of a signal with the Fourier basis functions,
-#     we essentially measure how well the signal correlates with sinusoids of different frequencies,
-#     within a specified temporal window. This process mirrors the operation performed by the Fourier transform.
-#     Therefore, it becomes clear that computing the cross-correlation of a signal with the Fourier basis defined here
-#     is equivalent to computing the discrete Fourier transform on a sliding window of the same size
-#     as that of the basis.
+# When computing the cross-correlation of a signal with the Fourier basis functions,
+# we essentially measure how well the signal correlates with sinusoids of different frequencies,
+# within a specified temporal window. This process mirrors the operation performed by the Fourier transform.
+# Therefore, it becomes clear that computing the cross-correlation of a signal with the Fourier basis defined here
+# is equivalent to computing the discrete Fourier transform on a sliding window of the same size
+# as that of the basis.
 
 n_samples = 1000
 n_freqs = 20
@@ -131,7 +131,7 @@
 signal = np.random.normal(size=n_samples)
 
 # evaluate the basis
-_, eval_basis = nmo.basis.FourierBasis(n_freqs=n_freqs).evaluate_on_grid(n_samples)
+_, eval_basis = nmo.basis.FourierBasis(max_freq=n_freqs).evaluate_on_grid(n_samples)
 
 # compute the cross-corr with the signal and the basis
 # Note that we are inverting the time axis of the basis because we are aiming
diff --git a/src/nemos/basis.py b/src/nemos/basis.py
index f0b8b536..a380bfbd 100644
--- a/src/nemos/basis.py
+++ b/src/nemos/basis.py
@@ -52,7 +52,7 @@ def __init__(self, n_basis_funcs: int) -> None:
         self._check_n_basis_min()
 
     @abc.abstractmethod
-    def evaluate(self, *xi: NDArray) -> NDArray:
+    def evaluate(self, *xi: ArrayLike) -> NDArray:
         """
         Evaluate the basis set at the given samples x1,...,xn using the subclass-specific "evaluate" method.
 
@@ -328,7 +328,7 @@ def __init__(self, basis1: Basis, basis2: Basis) -> None:
     def _check_n_basis_min(self) -> None:
         pass
 
-    def evaluate(self, *xi: NDArray) -> NDArray:
+    def evaluate(self, *xi: ArrayLike) -> NDArray:
         """
         Evaluate the basis at the input samples.
 
@@ -384,7 +384,7 @@ def __init__(self, basis1: Basis, basis2: Basis) -> None:
     def _check_n_basis_min(self) -> None:
         pass
 
-    def evaluate(self, *xi: NDArray) -> NDArray:
+    def evaluate(self, *xi: ArrayLike) -> NDArray:
         """
         Evaluate the basis at the input samples.
 
@@ -506,7 +506,7 @@ def _check_n_basis_min(self) -> None:
 
 
 class MSplineBasis(SplineBasis):
-    """M-spline 1-dimensional basis functions.
+    """M-spline[$^1$](#references) 1-dimensional basis functions.
 
     Parameters
     ----------
@@ -520,15 +520,14 @@ class MSplineBasis(SplineBasis):
 
     References
     ----------
-    [^1]:
-        Ramsay, J. O. (1988). Monotone regression splines in action.
+    1. Ramsay, J. O. (1988). Monotone regression splines in action.
         Statistical science, 3(4), 425-441.
     """
 
     def __init__(self, n_basis_funcs: int, order: int = 2) -> None:
         super().__init__(n_basis_funcs, order)
 
-    def evaluate(self, sample_pts: NDArray) -> NDArray:
+    def evaluate(self, sample_pts: ArrayLike) -> NDArray:
         """Generate basis functions with given spacing.
 
         Parameters
@@ -578,7 +577,7 @@ def evaluate_on_grid(self, n_samples: int) -> Tuple[NDArray, NDArray]:
 
 class BSplineBasis(SplineBasis):
     """
-    B-spline 1-dimensional basis functions.
+    B-spline[$^1$](#references) 1-dimensional basis functions.
 
     Parameters
     ----------
@@ -597,8 +596,7 @@ class BSplineBasis(SplineBasis):
 
     References
     ----------
-    [^2]:
-        Prautzsch, H., Boehm, W., Paluszny, M. (2002). B-spline representation. In: Bézier and B-Spline Techniques.
+    1. Prautzsch, H., Boehm, W., Paluszny, M. (2002). B-spline representation. In: Bézier and B-Spline Techniques.
         Mathematics and Visualization. Springer, Berlin, Heidelberg. https://doi.org/10.1007/978-3-662-04919-8_5
 
     """
@@ -606,7 +604,7 @@ class BSplineBasis(SplineBasis):
     def __init__(self, n_basis_funcs: int, order: int = 2):
         super().__init__(n_basis_funcs, order=order)
 
-    def evaluate(self, sample_pts: NDArray) -> NDArray:
+    def evaluate(self, sample_pts: ArrayLike) -> NDArray:
         """
         Evaluate the B-spline basis functions with given sample points.
 
@@ -694,7 +692,7 @@ def __init__(self, n_basis_funcs: int, order: int = 2):
                 f"order {self.order} specified instead!"
             )
 
-    def evaluate(self, sample_pts: NDArray) -> NDArray:
+    def evaluate(self, sample_pts: ArrayLike) -> NDArray:
         """Evaluate the Cyclic B-spline basis functions with given sample points.
 
         Parameters
@@ -772,34 +770,68 @@ def evaluate_on_grid(self, n_samples: int) -> Tuple[NDArray, NDArray]:
         return super().evaluate_on_grid(n_samples)
 
 
-class RaisedCosineBasis(Basis, abc.ABC):
-    def __init__(self, n_basis_funcs: int) -> None:
+class RaisedCosineBasisLinear(Basis):
+    """Represent linearly-spaced raised cosine basis functions.
+
+    This implementation is based on the cosine bumps used by Pillow et al.[$^1$](#references)
+    to uniformly tile the internal points of the domain.
+
+    Parameters
+    ----------
+    n_basis_funcs :
+        The number of basis functions.
+    width :
+        Width of the raised cosine. By default, it's set to 2.0.
+
+    References
+    ----------
+    1. Pillow, J. W., Paninski, L., Uzzel, V. J., Simoncelli, E. P., & J.,
+        C. E. (2005). Prediction and decoding of retinal ganglion cell responses
+        with a probabilistic spiking model. Journal of Neuroscience, 25(47),
+        11003–11013. http://dx.doi.org/10.1523/jneurosci.3305-05.2005
+    """
+
+    def __init__(self, n_basis_funcs: int, width: float = 2.0) -> None:
         super().__init__(n_basis_funcs)
         self._n_input_dimensionality = 1
+        self._check_width(width)
+        self._width = width
 
-    @abc.abstractmethod
-    def _transform_samples(self, sample_pts: NDArray) -> NDArray:
-        """
-        Abstract method for transforming sample points.
+    @property
+    def width(self):
+        """Return width of the raised cosine."""
+        return self._width
+
+    @staticmethod
+    def _check_width(width: float):
+        """Validate the width value.
 
         Parameters
         ----------
-        sample_pts :
-           The sample points to be transformed, shape (n_samples,).
+        width :
+            The width value to validate.
+
+        Raises
+        ------
+        ValueError
+            If width <= 1 or 2*width is not a positive integer. Values that do not match
+            this constraint will result in:
+            - No overlap between bumps (width < 1).
+            - Oscillatory behavior when summing the basis elements (2*width not integer).
         """
-        pass
+        if width <= 1 or (not np.isclose(width * 2, round(2 * width))):
+            raise ValueError(
+                f"Invalid raised cosine width. "
+                f"2*width must be a positive integer, 2*width = {2 * width} instead!"
+            )
 
-    def evaluate(self, sample_pts: NDArray) -> NDArray:
+    def evaluate(self, sample_pts: ArrayLike) -> NDArray:
         """Generate basis functions with given samples.
 
         Parameters
         ----------
-        sample_pts : (number of samples,)
-            Spacing for basis functions, holding elements on interval [0, 1). A
-            good default is ``nmo.sample_points.raised_cosine_log`` for log
-            spacing (as used in [2]_) or
-            ``nmo.sample_points.raised_cosine_linear`` for linear spacing.
-            Shape (n_samples,).
+        sample_pts :
+            Spacing for basis functions, holding elements on interval [0, 1], Shape (number of samples, ).
 
         Returns
         -------
@@ -816,17 +848,34 @@ def evaluate(self, sample_pts: NDArray) -> NDArray:
         if any(sample_pts < 0) or any(sample_pts > 1):
             raise ValueError("Sample points for RaisedCosine basis must lie in [0,1]!")
 
-        # transform to the proper domain
-        transform_sample_pts = self._transform_samples(sample_pts)
-
-        shifted_sample_pts = (
-            transform_sample_pts[:, None]
-            - (np.pi * np.arange(self.n_basis_funcs))[None, :]
+        peaks = self._compute_peaks()
+        delta = peaks[1] - peaks[0]
+        # generate a set of shifted cosines, and constrain them to be non-zero
+        # over a single period, then enforce the codomain to be [0,1], by adding 1
+        # and then multiply by 0.5
+        basis_funcs = 0.5 * (
+            np.cos(
+                np.clip(
+                    np.pi * (sample_pts[:, None] - peaks[None]) / (delta * self.width),
+                    -np.pi,
+                    np.pi,
+                )
+            )
+            + 1
         )
-        basis_funcs = 0.5 * (np.cos(np.clip(shifted_sample_pts, -np.pi, np.pi)) + 1)
 
         return basis_funcs
 
+    def _compute_peaks(self):
+        """
+        Compute the location of raised cosine peaks.
+
+        Returns
+        -------
+            Peak locations of each basis element.
+        """
+        return np.linspace(0, 1, self.n_basis_funcs)
+
     def evaluate_on_grid(self, n_samples: int) -> Tuple[NDArray, NDArray]:
         """Evaluate the basis set on a grid of equi-spaced sample points.
 
@@ -846,127 +895,140 @@ def evaluate_on_grid(self, n_samples: int) -> Tuple[NDArray, NDArray]:
         """
         return super().evaluate_on_grid(n_samples)
 
-
-class RaisedCosineBasisLinear(RaisedCosineBasis):
-    """Linearly-spaced raised cosine basis functions used by Pillow et al.
-
-    These are "cosine bumps" that uniformly tile the space.
-
-
-    Parameters
-    ----------
-    n_basis_funcs
-        Number of basis functions.
-
-    References
-    ----------
-    [^3]:
-        Pillow, J. W., Paninski, L., Uzzel, V. J., Simoncelli, E. P., & J.,
-        C. E. (2005). Prediction and decoding of retinal ganglion cell responses
-        with a probabilistic spiking model. Journal of Neuroscience, 25(47),
-        11003–11013. http://dx.doi.org/10.1523/jneurosci.3305-05.2005
-
-    """
-
-    def __init__(self, n_basis_funcs: int) -> None:
-        super().__init__(n_basis_funcs)
-
-    def _transform_samples(self, sample_pts: NDArray) -> NDArray:
-        """
-        Linearly map the samples from [0,1] to the the [0, (n_basis_funcs - 1) * pi].
-
-        Parameters
-        ----------
-        sample_pts :
-            The sample points used for evaluating the splines, shape (n_samples,)
-
-        Returns
-        -------
-        :
-            A transformed version of the sample points that matches the Raised Cosine basis domain,
-            shape (number of samples, ).
-        """
-        return sample_pts * np.pi * (self.n_basis_funcs - 1)
-
     def _check_n_basis_min(self) -> None:
         """Check that the user required enough basis elements.
 
-        Check that the number of basis is at least 1.
+        Check that the number of basis is at least 2.
 
         Raises
         ------
         ValueError
-            If an insufficient number of basis element is requested for the basis type
+            If n_basis_funcs < 2.
         """
-        if self.n_basis_funcs < 1:
+        if self.n_basis_funcs < 2:
             raise ValueError(
-                f"Object class {self.__class__.__name__} requires >= 1 basis elements. "
+                f"Object class {self.__class__.__name__} requires >= 2 basis elements. "
                 f"{self.n_basis_funcs} basis elements specified instead"
             )
 
 
-class RaisedCosineBasisLog(RaisedCosineBasis):
-    """Log-spaced raised cosine basis functions used by Pillow et al. [2]_.
+class RaisedCosineBasisLog(RaisedCosineBasisLinear):
+    """Represent log-spaced raised cosine basis functions.
 
-    These are "cosine bumps" that uniformly tile the space.
+    Similar to `RaisedCosineBasisLinear` but the basis functions are log-spaced.
+    This implementation is based on the cosine bumps used by Pillow et al.[$^1$](#references)
+    to uniformly tile the internal points of the domain.
 
     Parameters
     ----------
-    n_basis_funcs
-        Number of basis functions.
+    n_basis_funcs :
+        The number of basis functions.
+    width :
+        Width of the raised cosine. By default, it's set to 2.0.
+    enforce_decay_to_zero:
+        If set to True, the algorithm first constructs a basis with `n_basis_funcs + ceil(width)` elements
+        and subsequently trims off the extra basis elements. This ensures that the final basis element
+        decays to 0.
 
     References
     ----------
-    .. [2] Pillow, J. W., Paninski, L., Uzzel, V. J., Simoncelli, E. P., & J.,
+    1. Pillow, J. W., Paninski, L., Uzzel, V. J., Simoncelli, E. P., & J.,
        C. E. (2005). Prediction and decoding of retinal ganglion cell responses
        with a probabilistic spiking model. Journal of Neuroscience, 25(47),
        11003–11013. http://dx.doi.org/10.1523/jneurosci.3305-05.2005
-
     """
 
-    def __init__(self, n_basis_funcs: int) -> None:
-        super().__init__(n_basis_funcs)
+    def __init__(
+        self,
+        n_basis_funcs: int,
+        width: float = 2.0,
+        time_scaling: float = None,
+        enforce_decay_to_zero: bool = True,
+    ) -> None:
+        super().__init__(n_basis_funcs, width=width)
+        self.enforce_decay_to_zero = enforce_decay_to_zero
+        if time_scaling is None:
+            self._time_scaling = 50.0
+        else:
+            self._check_time_scaling(time_scaling)
+            self._time_scaling = time_scaling
+
+    @property
+    def time_scaling(self):
+        return self._time_scaling
 
-    def _transform_samples(self, sample_pts: NDArray) -> NDArray:
-        """Map the sample domain to log-space.
+    @staticmethod
+    def _check_time_scaling(time_scaling):
+        if time_scaling <= 0:
+            raise ValueError(
+                f"Only strictly positive time_scaling are allowed, {time_scaling} provided instead."
+            )
 
-        Map the equi-spaced samples from [0,1] to log equi-spaced samples [0, (n_basis_funcs - 1) * pi].
+    def _transform_samples(self, sample_pts: NDArray) -> NDArray:
+        """
+        Map the sample domain to log-space.
 
         Parameters
         ----------
         sample_pts :
-            The sample points used for evaluating the splines, shape (n_samples,).
+            Sample points used for evaluating the splines,
+            shape (n_samples, ).
 
         Returns
         -------
-        :
-            A transformed version of the sample points that matches the Raised Cosine basis domain,
-            shape (n_sample_points, ).
+            Transformed version of the sample points that matches the Raised Cosine basis domain,
+            shape (n_samples, ).
         """
-        return (
-            np.power(
-                10,
-                -(np.log10((self.n_basis_funcs - 1) * np.pi) + 1) * sample_pts
-                + np.log10((self.n_basis_funcs - 1) * np.pi),
-            )
-            - 0.1
+        # This log-stretching of the sample axis has the following effect:
+        # - as the time_scaling tends to 0, the points will be linearly spaced across the whole domain.
+        # - as the time_scaling tends to inf, basis will be small and dense around 0 and
+        # progressively larger and less dense towards 1.
+        log_spaced_pts = np.log(self.time_scaling * sample_pts + 1) / np.log(
+            self.time_scaling + 1
         )
+        return log_spaced_pts
 
-    def _check_n_basis_min(self) -> None:
-        """Check that the user required enough basis elements.
+    def _compute_peaks(self):
+        """
+        Peak location of each log-spaced cosine basis element
+
+        Compute the peak location for the log-spaced raised cosine basis.
+        Enforcing that the last basis decays to zero is equivalent to
+        setting the last peak to a value smaller than 1.
+
+        Returns
+        -------
+            Peak locations of each basis element.
+
+        """
+        if self.enforce_decay_to_zero:
+            # compute the last peak location such that the last
+            # basis element decays to zero at the last sample.
+            last_peak = 1 - self.width / (self.n_basis_funcs + self.width - 1)
+        else:
+            last_peak = 1
+        return np.linspace(0, last_peak, self.n_basis_funcs)
+
+    def evaluate(self, sample_pts: ArrayLike) -> NDArray:
+        """Generate log-spaced raised cosine basis with given samples.
+
+        Parameters
+        ----------
+        sample_pts :
+            Spacing for basis functions, holding elements on interval [0, 1].
 
-        Checks that the number of basis is at least 2.
+        Returns
+        -------
+        basis_funcs :
+            Log-raised cosine basis functions, shape (n_samples, n_basis_funcs).
 
         Raises
         ------
         ValueError
-            If an insufficient number of basis element is requested for the basis type
+            If the sample provided do not lie in [0,1].
         """
-        if self.n_basis_funcs < 2:
-            raise ValueError(
-                f"Object class {self.__class__.__name__} requires >= 2 basis elements. "
-                f"{self.n_basis_funcs} basis elements specified instead"
-            )
+        (sample_pts,) = self._check_evaluate_input(sample_pts)
+        return super().evaluate(self._transform_samples(sample_pts))
 
 
 class OrthExponentialBasis(Basis):
@@ -1120,16 +1182,19 @@ def evaluate_on_grid(self, n_samples: int) -> Tuple[NDArray, NDArray]:
 class FourierBasis(Basis):
     """Set of 1D Fourier basis.
 
+    This class defines a cosine and negative sine basis (quadrature pair)
+    with frequencies ranging 0 to max_freq - 1.
+
     Parameters
     ----------
-    n_freqs
-            Number of frequencies. The number of basis function will be 2*n_freqs - 1.
+    max_freq
+            Number of frequencies. The number of basis function will be 2*max_freq - 1.
     """
 
-    def __init__(self, n_freqs: int):
-        super().__init__(n_basis_funcs=2 * n_freqs - 1)
+    def __init__(self, max_freq: int):
+        super().__init__(n_basis_funcs=2 * max_freq - 1)
 
-        self._frequencies = np.arange(n_freqs, dtype=np.float32)
+        self._frequencies = np.arange(max_freq, dtype=float)
         self._n_input_dimensionality = 1
 
     def _check_n_basis_min(self) -> None:
@@ -1140,7 +1205,7 @@ def _check_n_basis_min(self) -> None:
         Raises
         ------
         ValueError
-            If an insufficient number of basis element is requested for the basis type
+            If an insufficient number of basis element is requested for the basis type.
         """
         if self.n_basis_funcs < 1:
             raise ValueError(
@@ -1163,25 +1228,25 @@ def evaluate(self, sample_pts: NDArray) -> NDArray:
 
         Notes
         -----
-            If the frequencies provided are np.arange(n_freq), convolving a signal
-            of length n_samples with this basis is equivalent, but slower,
-            then computing the FFT truncated to the first n_freq components.
+        If the frequencies provided are np.arange(max_freq), convolving a signal
+        of length n_samples with this basis is equivalent, but slower,
+        then computing the FFT truncated to the first max_freq components.
 
-            Therefore, convolving a signal with this basis is equivalent
-            to compute the FFT over sliding window.
+        Therefore, convolving a signal with this basis is equivalent
+        to compute the FFT over a sliding window.
 
         Examples
         --------
-            >>> import nemos as nmo
-            >>> import numpy as np
-            >>> n_samples, n_freqs = 1000, 10
-            >>> basis = nmo.basis.FourierBasis(n_freqs*2)
-            >>> eval_basis = basis.evaluate(np.linspace(0, 1, n_samples))
-            >>> sinusoid = np.cos(3 * np.arange(0, 1000) * np.pi * 2 / 1000.)
-            >>> conv = [np.convolve(eval_basis[::-1, k], sinusoid, mode='valid')[0] for k in range(2*n_freqs-1)]
-            >>> fft = np.fft.fft(sinusoid)
-            >>> print('FFT power:   ', np.round(np.real(fft[:10]), 4))
-            >>> print('Convolution: ', np.round(conv[:10], 4))
+        >>> import nemos as nmo
+        >>> import numpy as np
+        >>> n_samples, max_freq = 1000, 10
+        >>> basis = nmo.basis.FourierBasis(max_freq*2)
+        >>> eval_basis = basis.evaluate(np.linspace(0, 1, n_samples))
+        >>> sinusoid = np.cos(3 * np.arange(0, 1000) * np.pi * 2 / 1000.)
+        >>> conv = [np.convolve(eval_basis[::-1, k], sinusoid, mode='valid')[0] for k in range(2*max_freq-1)]
+        >>> fft = np.fft.fft(sinusoid)
+        >>> print('FFT power:   ', np.round(np.real(fft[:10]), 4))
+        >>> print('Convolution: ', np.round(conv[:10], 4))
         """
         (sample_pts,) = self._check_evaluate_input(sample_pts)
         # assumes equi-spaced samples.

From aa67d9c3044e1e5fd88be3fd67c6267ea44a9ace Mon Sep 17 00:00:00 2001
From: BalzaniEdoardo <edoardo.balzani87@gmail.com>
Date: Tue, 19 Dec 2023 18:19:29 -0500
Subject: [PATCH 08/12] improved dovstrings

---
 docs/examples/plot_1D_basis_function.py | 25 +++++++++++++------------
 src/nemos/basis.py                      | 21 +++++++++++----------
 2 files changed, 24 insertions(+), 22 deletions(-)

diff --git a/docs/examples/plot_1D_basis_function.py b/docs/examples/plot_1D_basis_function.py
index 44e75cfa..6caf5df1 100644
--- a/docs/examples/plot_1D_basis_function.py
+++ b/docs/examples/plot_1D_basis_function.py
@@ -91,12 +91,13 @@
 # interpretation of the model parameters, each of which will represent the relative contribution of a specific
 # oscillation frequency to the overall signal.
 #
-# A Fourier basis can be instantiated with the usual syntax.
-# The user can pass the desired frequencies for the basis or
-# the frequencies will be set to `np.arange(n_basis_funcs//2)`.
-# The number of basis function is required to be even.
+# A Fourier basis can be instantiated with the following syntax:
+# the user can provide the maximum frequency of the cosine and negative
+# sine pairs by setting the `max_freq` parameter.
+# The sinusoidal basis elements will have frequencies from 0 to `max_freq`.
 
-fourier_basis = nmo.basis.FourierBasis(max_freq=4)
+
+fourier_basis = nmo.basis.FourierBasis(max_freq=3)
 
 # evaluate on equi-spaced samples
 samples, eval_basis = fourier_basis.evaluate_on_grid(1000)
@@ -126,13 +127,13 @@
 
 
 n_samples = 1000
-n_freqs = 20
+max_freq = 20
 
 # define a signal
 signal = np.random.normal(size=n_samples)
 
 # evaluate the basis
-_, eval_basis = nmo.basis.FourierBasis(max_freq=n_freqs).evaluate_on_grid(n_samples)
+_, eval_basis = nmo.basis.FourierBasis(max_freq=max_freq).evaluate_on_grid(n_samples)
 
 # compute the cross-corr with the signal and the basis
 # Note that we are inverting the time axis of the basis because we are aiming
@@ -140,17 +141,17 @@
 xcorr = np.array(
     [
         np.convolve(eval_basis[::-1, k], signal, mode="valid")[0]
-        for k in range(2 * n_freqs - 1)
+        for k in range(2 * max_freq + 1)
     ]
 )
 
 # compute the power (add back sin(0 * t) = 0)
 fft_complex = np.fft.fft(signal)
-fft_amplitude = np.abs(fft_complex[:n_freqs])
-fft_phase = np.angle(fft_complex[:n_freqs])
+fft_amplitude = np.abs(fft_complex[:max_freq + 1])
+fft_phase = np.angle(fft_complex[:max_freq + 1])
 # compute the phase and amplitude from the convolution
-xcorr_phase = np.arctan2(np.hstack([[0], xcorr[n_freqs:]]), xcorr[:n_freqs])
-xcorr_aplitude = np.sqrt(xcorr[:n_freqs] ** 2 + np.hstack([[0], xcorr[n_freqs:]]) ** 2)
+xcorr_phase = np.arctan2(np.hstack([[0], xcorr[max_freq+1:]]), xcorr[:max_freq+1])
+xcorr_aplitude = np.sqrt(xcorr[:max_freq+1] ** 2 + np.hstack([[0], xcorr[max_freq+1:]]) ** 2)
 
 fig, ax = plt.subplots(1, 2)
 ax[0].set_aspect("equal")
diff --git a/src/nemos/basis.py b/src/nemos/basis.py
index a380bfbd..363f6211 100644
--- a/src/nemos/basis.py
+++ b/src/nemos/basis.py
@@ -1183,18 +1183,19 @@ class FourierBasis(Basis):
     """Set of 1D Fourier basis.
 
     This class defines a cosine and negative sine basis (quadrature pair)
-    with frequencies ranging 0 to max_freq - 1.
+    with frequencies ranging 0 to max_freq.
 
     Parameters
     ----------
     max_freq
-            Number of frequencies. The number of basis function will be 2*max_freq - 1.
+            Highest frequency of the cosine, negative sine pairs.
+            The number of basis function will be 2*max_freq + 1.
     """
 
     def __init__(self, max_freq: int):
-        super().__init__(n_basis_funcs=2 * max_freq - 1)
+        super().__init__(n_basis_funcs=2 * max_freq + 1)
 
-        self._frequencies = np.arange(max_freq, dtype=float)
+        self._frequencies = np.arange(max_freq + 1, dtype=float)
         self._n_input_dimensionality = 1
 
     def _check_n_basis_min(self) -> None:
@@ -1213,7 +1214,7 @@ def _check_n_basis_min(self) -> None:
                 f"{self.n_basis_funcs} basis elements specified instead"
             )
 
-    def evaluate(self, sample_pts: NDArray) -> NDArray:
+    def evaluate(self, sample_pts: ArrayLike) -> NDArray:
         """Generate basis functions with given spacing.
 
         Parameters
@@ -1228,7 +1229,7 @@ def evaluate(self, sample_pts: NDArray) -> NDArray:
 
         Notes
         -----
-        If the frequencies provided are np.arange(max_freq), convolving a signal
+        The frequencies are set to np.arange(max_freq+1), convolving a signal
         of length n_samples with this basis is equivalent, but slower,
         then computing the FFT truncated to the first max_freq components.
 
@@ -1240,13 +1241,13 @@ def evaluate(self, sample_pts: NDArray) -> NDArray:
         >>> import nemos as nmo
         >>> import numpy as np
         >>> n_samples, max_freq = 1000, 10
-        >>> basis = nmo.basis.FourierBasis(max_freq*2)
+        >>> basis = nmo.basis.FourierBasis(max_freq)
         >>> eval_basis = basis.evaluate(np.linspace(0, 1, n_samples))
         >>> sinusoid = np.cos(3 * np.arange(0, 1000) * np.pi * 2 / 1000.)
-        >>> conv = [np.convolve(eval_basis[::-1, k], sinusoid, mode='valid')[0] for k in range(2*max_freq-1)]
+        >>> conv = [np.convolve(eval_basis[::-1, k], sinusoid, mode='valid')[0] for k in range(2*max_freq+1)]
         >>> fft = np.fft.fft(sinusoid)
-        >>> print('FFT power:   ', np.round(np.real(fft[:10]), 4))
-        >>> print('Convolution: ', np.round(conv[:10], 4))
+        >>> print('FFT power:   ', np.round(np.real(fft[:max_freq]), 4))
+        >>> print('Convolution: ', np.round(conv[:max_freq], 4))
         """
         (sample_pts,) = self._check_evaluate_input(sample_pts)
         # assumes equi-spaced samples.

From adb5efbecc361f39ac5abc13b4df9184e2f31603 Mon Sep 17 00:00:00 2001
From: BalzaniEdoardo <edoardo.balzani87@gmail.com>
Date: Tue, 19 Dec 2023 18:26:05 -0500
Subject: [PATCH 09/12] fixed tests for new parametrization

---
 src/nemos/basis.py  |  2 +-
 tests/test_basis.py | 50 ++++++++++++++++++++++-----------------------
 2 files changed, 26 insertions(+), 26 deletions(-)

diff --git a/src/nemos/basis.py b/src/nemos/basis.py
index 363f6211..c9501e8e 100644
--- a/src/nemos/basis.py
+++ b/src/nemos/basis.py
@@ -1210,7 +1210,7 @@ def _check_n_basis_min(self) -> None:
         """
         if self.n_basis_funcs < 1:
             raise ValueError(
-                f"Object class {self.__class__.__name__} requires >= 1 basis elements. "
+                f"Object class {self.__class__.__name__} requires >= 0 basis elements. "
                 f"{self.n_basis_funcs} basis elements specified instead"
             )
 
diff --git a/tests/test_basis.py b/tests/test_basis.py
index be9df818..9ad6dd9a 100644
--- a/tests/test_basis.py
+++ b/tests/test_basis.py
@@ -849,15 +849,15 @@ def test_decay_rate_size_match_n_basis_func(self, decay_rates, n_basis_func):
 class TestFourierBasis(BasisFuncsTesting):
     cls = basis.FourierBasis
 
-    @pytest.mark.parametrize("n_freqs", [2, 4, 8])
+    @pytest.mark.parametrize("max_freq", [2, 4, 8])
     @pytest.mark.parametrize("sample_size", [20, 1000])
     def test_evaluate_returns_expected_number_of_basis(
-            self, n_freqs, sample_size
+            self, max_freq, sample_size
     ):
         """Tests whether the evaluate method returns the expected number of basis functions."""
-        basis_obj = self.cls(n_freqs=n_freqs)
+        basis_obj = self.cls(max_freq=max_freq)
         eval_basis = basis_obj.evaluate(np.linspace(0, 1, sample_size))
-        assert(eval_basis.shape[1] == 2*n_freqs-1)
+        assert(eval_basis.shape[1] == 2*max_freq+1)
 
     @pytest.mark.parametrize("samples, expectation",
         [
@@ -882,19 +882,19 @@ def test_input_to_evaluate_is_arraylike(self, arraylike, expectation):
         """
         Checks that the sample size of the output from the evaluate() method matches the input sample size.
         """
-        basis_obj = self.cls(n_freqs=5)
+        basis_obj = self.cls(max_freq=5)
         with expectation:
             basis_obj.evaluate(arraylike)
 
     @pytest.mark.parametrize("sample_size", [100, 1000])
-    @pytest.mark.parametrize("n_freqs", [4, 10])
+    @pytest.mark.parametrize("max_freq", [4, 10])
     def test_sample_size_of_evaluate_matches_that_of_input(
-            self, n_freqs, sample_size
+            self, max_freq, sample_size
     ):
         """
         Checks that the sample size of the output from the evaluate() method matches the input sample size.
         """
-        basis_obj = self.cls(n_freqs)
+        basis_obj = self.cls(max_freq)
         eval_basis = basis_obj.evaluate(np.linspace(0, 1, sample_size))
         if eval_basis.shape[0] != sample_size:
             raise ValueError(
@@ -903,39 +903,39 @@ def test_sample_size_of_evaluate_matches_that_of_input(
                 f"The second dimension of the evaluated basis is {eval_basis.shape[0]}",
             )
 
-    @pytest.mark.parametrize("n_freqs, expectation",
+    @pytest.mark.parametrize("max_freq, expectation",
         [
-            (-1, pytest.raises(ValueError,match=r"Object class FourierBasis requires >= 1 basis elements")),
-            (0, pytest.raises(ValueError,match=r"Object class FourierBasis requires >= 1 basis elements")),
+            (-1, pytest.raises(ValueError,match=r"Object class FourierBasis requires >= 0 basis elements")),
+            (0, does_not_raise()),
             (1, does_not_raise()),
             (3, does_not_raise())
         ]
     )
-    def test_minimum_number_of_basis_required_is_matched(self, n_freqs, expectation):
+    def test_minimum_number_of_basis_required_is_matched(self, max_freq, expectation):
         """
         Verifies that the minimum number of basis functions and order required (i.e., at least 1) and
         order < #basis are enforced.
         """
         n_samples = 10
         with expectation:
-            basis_obj = self.cls(n_freqs=n_freqs)
+            basis_obj = self.cls(max_freq=max_freq)
             basis_obj.evaluate(np.linspace(0, 1, n_samples))
 
-    @pytest.mark.parametrize("n_freqs, expectation",
+    @pytest.mark.parametrize("max_freq, expectation",
                              [
                                  (3, does_not_raise()),
-                                 (6, does_not_raise()),
+                                 (6, pytest.raises(ValueError,match=rf"Not enough samples, aliasing likely to occur")),
                                  (7, pytest.raises(ValueError,match=rf"Not enough samples, aliasing likely to occur")),
                                  (10, pytest.raises(ValueError,match=rf"Not enough samples, aliasing likely to occur"))
                              ]
                              )
-    def test_minimum_aliasing_detection(self, n_freqs, expectation):
+    def test_minimum_aliasing_detection(self, max_freq, expectation):
         """
         Verifies that the minimum number of basis functions and order required (i.e., at least 1) and
         order < #basis are enforced.
         """
         n_samples = 10
-        basis_obj = self.cls(n_freqs=n_freqs)
+        basis_obj = self.cls(max_freq=max_freq)
         with expectation:
             basis_obj.evaluate(np.linspace(0, 1, n_samples))
 
@@ -946,7 +946,7 @@ def test_samples_range_matches_evaluate_requirements(self, sample_range: tuple):
         """
         Verifies that the evaluate() method can handle input range.
         """
-        basis_obj = self.cls(n_freqs=5)
+        basis_obj = self.cls(max_freq=5)
         basis_obj.evaluate(np.linspace(*sample_range, 100))
 
     @pytest.mark.parametrize("n_input, expectation", [
@@ -960,7 +960,7 @@ def test_number_of_required_inputs_evaluate(self, n_input, expectation):
         """
         Confirms that the evaluate() method correctly handles the number of input samples that are provided.
         """
-        basis_obj = self.cls(n_freqs=5)
+        basis_obj = self.cls(max_freq=5)
         inputs = [np.linspace(0, 1, 20)] * n_input
         with expectation:
             basis_obj.evaluate(*inputs)
@@ -979,7 +979,7 @@ def test_evaluate_on_grid_meshgrid_valid_size(self, sample_size, expectation):
         """
         Checks that the evaluate_on_grid() method returns a grid of the expected size.
         """
-        basis_obj = self.cls(n_freqs=5)
+        basis_obj = self.cls(max_freq=5)
         with expectation:
             basis_obj.evaluate_on_grid(sample_size)
 
@@ -995,7 +995,7 @@ def test_evaluate_on_grid_meshgrid_match_size(self, sample_size, expectation):
         """
         Checks that the evaluate_on_grid() method returns a grid of the expected size.
         """
-        basis_obj = self.cls(n_freqs=5)
+        basis_obj = self.cls(max_freq=5)
         with expectation:
             grid, _ = basis_obj.evaluate_on_grid(sample_size)
             assert grid.shape[0] == sample_size
@@ -1012,7 +1012,7 @@ def test_evaluate_on_grid_input_number(self, n_input, expectation):
         """
         Validates that the evaluate_on_grid() method correctly handles the number of input samples that are provided.
         """
-        basis_obj = self.cls(n_freqs=5)
+        basis_obj = self.cls(max_freq=5)
         inputs = [10] * n_input
         with expectation:
             basis_obj.evaluate_on_grid(*inputs)
@@ -1426,7 +1426,7 @@ def instantiate_basis(n_basis, basis_class):
                 n_basis_funcs=n_basis, decay_rates=np.arange(1, 1 + n_basis)
             )
         elif basis_class == basis.FourierBasis:
-            basis_obj = basis_class(n_freqs=n_basis)
+            basis_obj = basis_class(max_freq=n_basis)
         elif basis_class == basis.BSplineBasis:
             basis_obj = basis_class(n_basis_funcs=n_basis, order=3)
         elif basis_class == basis.CyclicBSplineBasis:
@@ -1481,7 +1481,7 @@ def test_evaluate_input(self, eval_input):
 
     @pytest.mark.parametrize("n_basis_a", [5, 6])
     @pytest.mark.parametrize("n_basis_b", [5, 6])
-    @pytest.mark.parametrize("sample_size", [10, 1000])
+    @pytest.mark.parametrize("sample_size", [15, 1000])
     @pytest.mark.parametrize(
         "basis_a",
         [class_obj for _, class_obj in utils_testing.get_non_abstract_classes(basis)],
@@ -1697,7 +1697,7 @@ def test_evaluate_input(self, eval_input):
 
     @pytest.mark.parametrize("n_basis_a", [5, 6])
     @pytest.mark.parametrize("n_basis_b", [5, 6])
-    @pytest.mark.parametrize("sample_size", [10, 1000])
+    @pytest.mark.parametrize("sample_size", [15, 1000])
     @pytest.mark.parametrize(
         "basis_a",
         [class_obj for _, class_obj in utils_testing.get_non_abstract_classes(basis)],

From a62712e0052afbbe6d0cd01d4f31cad5d6f61285 Mon Sep 17 00:00:00 2001
From: BalzaniEdoardo <edoardo.balzani87@gmail.com>
Date: Wed, 20 Dec 2023 13:26:17 -0500
Subject: [PATCH 10/12] allows dc term only

---
 src/nemos/basis.py  | 4 ++--
 tests/test_basis.py | 2 +-
 2 files changed, 3 insertions(+), 3 deletions(-)

diff --git a/src/nemos/basis.py b/src/nemos/basis.py
index c9501e8e..8e85f4b1 100644
--- a/src/nemos/basis.py
+++ b/src/nemos/basis.py
@@ -1208,9 +1208,9 @@ def _check_n_basis_min(self) -> None:
         ValueError
             If an insufficient number of basis element is requested for the basis type.
         """
-        if self.n_basis_funcs < 1:
+        if self.n_basis_funcs < 0:
             raise ValueError(
-                f"Object class {self.__class__.__name__} requires >= 0 basis elements. "
+                f"Object class {self.__class__.__name__} requires >= 1 basis elements. "
                 f"{self.n_basis_funcs} basis elements specified instead"
             )
 
diff --git a/tests/test_basis.py b/tests/test_basis.py
index 9ad6dd9a..ee059b53 100644
--- a/tests/test_basis.py
+++ b/tests/test_basis.py
@@ -905,7 +905,7 @@ def test_sample_size_of_evaluate_matches_that_of_input(
 
     @pytest.mark.parametrize("max_freq, expectation",
         [
-            (-1, pytest.raises(ValueError,match=r"Object class FourierBasis requires >= 0 basis elements")),
+            (-1, pytest.raises(ValueError,match=r"Object class FourierBasis requires >= 1 basis elements")),
             (0, does_not_raise()),
             (1, does_not_raise()),
             (3, does_not_raise())

From 07bd9f6e2a7d4271b7cec7fe6c9d5c41ef204830 Mon Sep 17 00:00:00 2001
From: Edoardo Balzani <edoardo.balzani87@gmail.com>
Date: Fri, 23 Feb 2024 08:13:51 -0500
Subject: [PATCH 11/12] Update docs/examples/plot_1D_basis_function.py

Co-authored-by: William F. Broderick <billbrod@gmail.com>
---
 docs/examples/plot_1D_basis_function.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/docs/examples/plot_1D_basis_function.py b/docs/examples/plot_1D_basis_function.py
index 6caf5df1..deddf1da 100644
--- a/docs/examples/plot_1D_basis_function.py
+++ b/docs/examples/plot_1D_basis_function.py
@@ -119,9 +119,9 @@
 # where $ e^{-j\\omega \\tau} = \\cos(\\omega \\tau) - j \\sin (\\omega \\tau) $.
 #
 # When computing the cross-correlation of a signal with the Fourier basis functions,
-# we essentially measure how well the signal correlates with sinusoids of different frequencies,
+# we are measuring how well the signal correlates with sinusoids of different frequencies,
 # within a specified temporal window. This process mirrors the operation performed by the Fourier transform.
-# Therefore, it becomes clear that computing the cross-correlation of a signal with the Fourier basis defined here
+# Therefore, computing the cross-correlation of a signal with the Fourier basis defined here
 # is equivalent to computing the discrete Fourier transform on a sliding window of the same size
 # as that of the basis.
 

From b704c986202f3d8aca6cc14754f13c08bec5b10b Mon Sep 17 00:00:00 2001
From: Edoardo Balzani <edoardo.balzani87@gmail.com>
Date: Fri, 23 Feb 2024 08:14:06 -0500
Subject: [PATCH 12/12] Update docs/examples/plot_1D_basis_function.py

Co-authored-by: William F. Broderick <billbrod@gmail.com>
---
 docs/examples/plot_1D_basis_function.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/docs/examples/plot_1D_basis_function.py b/docs/examples/plot_1D_basis_function.py
index deddf1da..97fceb9f 100644
--- a/docs/examples/plot_1D_basis_function.py
+++ b/docs/examples/plot_1D_basis_function.py
@@ -113,7 +113,7 @@
 plt.tight_layout()
 
 # %%
-# ## Fourier Basis Convolution and Fourier Transform
+# #### Fourier Basis Convolution and Fourier Transform
 # The Fourier transform of a signal $ s(t) $ restricted to a temporal window $ [t_0,\;t_1] $ is
 # $$ \\hat{x}(\\omega) = \\int_{t_0}^{t_1} s(\\tau) e^{-j\\omega \\tau} d\\tau. $$
 # where $ e^{-j\\omega \\tau} = \\cos(\\omega \\tau) - j \\sin (\\omega \\tau) $.