[Enh] add example and handle one-dimensional arrays in CSD.generate l…

…fp (NeuralEnsemble#594)

* add example to documentation
* add handling for 1-dimensional arrays
* add regression test, testing for one- and two-dimensional arrays

---------

Co-authored-by: Michael Denker <[email protected]>
Co-authored-by: Alexander Kleinjohann <[email protected]>

elephant/current_source_density.py

@@ -71,16 +71,16 @@ def estimate_csd(lfp, coordinates='coordinates', method=None,
 
     Parameters
     ----------
-    lfp : neo.AnalogSignal
+    lfp : :class:`neo.core.AnalogSignal`
         Positions of electrodes can be added as an array annotation
     coordinates : array-like Quantity or string
         Specifies the corresponding spatial coordinates of the electrodes.
         Coordinates can be directly supplied by a NxM array-like Quantity
         with dimension of space, where M is the number of signals in 'lfp',
         and N is equal to the dimensionality of the method.
         Alternatively, if coordinates is a string, the function will fetch the
-        coordinates, supplied in the same format, as annotation of 'lfp' by that
-        name.
+        coordinates, supplied in the same format, as annotation of 'lfp' by
+        that name.
         Default: 'coordinates'
     method : string
         Pick a method corresponding to the setup, in this implementation
@@ -182,14 +182,14 @@ def estimate_csd(lfp, coordinates='coordinates', method=None,
                 # All iCSD methods explicitly assume a source
                 # diameter in contrast to the stdCSD  that
                 # implicitly assume infinite source radius
-                raise ValueError("Parameter diam must be specified for iCSD \
-                                  methods: {}".format(", ".join(icsd_methods)))
+                raise ValueError(f"Parameter diam must be specified for iCSD "
+                                 f"methods: {', '.join(icsd_methods)}")
 
         if 'f_type' in kwargs:
             if (kwargs['f_type'] != 'identity') and  \
                (kwargs['f_order'] is None):
-                raise ValueError("The order of {} filter must be \
-                                  specified".format(kwargs['f_type']))
+                raise ValueError(f"The order of {kwargs['f_type']} filter must"
+                                 f" be specified")
 
         csd_method = getattr(icsd, method)  # fetch class from icsd.py file
         csd_estimator = csd_method(lfp=lfp.T.magnitude * lfp.units,
@@ -228,12 +228,15 @@ def generate_lfp(csd_profile, x_positions, y_positions=None, z_positions=None,
         2D : large_source_2D and small_source_2D
         3D : gauss_3d_dipole
     x_positions : np.ndarray
-        Positions of the x coordinates of the electrodes
+        A 2D column vector (N x 1 array) containing the positions of the x
+        coordinates of the electrodes
     y_positions : np.ndarray, optional
-        Positions of the y coordinates of the electrodes
+        A 2D column vector (N x 1 array) containing the positions of the y
+        coordinates of the electrodes
         Defaults to None, use in 2D or 3D cases only
     z_positions : np.ndarray, optional
-        Positions of the z coordinates of the electrodes
+        A 2D column vector (N x 1 array) containing the positions of the z
+        coordinates of the electrodes
         Defaults to None, use in 3D case only
     x_limits : list, optional
         A list of [start, end].
@@ -253,10 +256,31 @@ def generate_lfp(csd_profile, x_positions, y_positions=None, z_positions=None,
 
     Returns
     -------
-    LFP : neo.AnalogSignal
+    LFP : :class:`neo.core.AnalogSignal`
        The potentials created by the csd profile at the electrode positions.
        The electrode positions are attached as an annotation named
        'coordinates'.
+
+    Examples
+    --------
+    >>> import numpy as np
+    >>> from elephant.current_source_density import generate_lfp, estimate_csd
+    >>> from elephant.current_source_density_src.utility_functions import gauss_1d_dipole  # noqa
+    >>> # 1. Define an array xs to x coordinate values with a length of 2304
+    >>> xs=np.linspace(0, 10, 2304)
+
+    >>> # 2. Run generate_lfp(gauss_1d_dipole, xs)
+    >>> lfp = generate_lfp(gauss_1d_dipole, xs)
+
+    >>> # 3. Run estimate_csd(lfp, method="StandardCSD")
+    >>> csd = estimate_csd(lfp, method="StandardCSD")  #doctest: +ELLIPSIS
+    discrete ...
+    >>> # 4. Print the results
+    >>> print(f"LFPs: {lfp}")
+    LFPs: [[-0.01483716 -0.01483396 -0.01483075 ...  0.01219233  0.0121911
+       0.01218986]] mV
+    >>> print(f"CSD estimate: {csd}")  #doctest: +ELLIPSIS
+    CSD estimate: [[-1.00025592e-04 -6.06684588e-05  ...
     """
 
     def integrate_1D(x0, csd_x, csd, h):
@@ -297,6 +321,10 @@ def integrate_3D(x, y, z, csd, xlin, ylin, zlin, X, Y, Z):
     sigma = 1.0
     h = 50.
     if dim == 1:
+        # Handle one dimensional case,
+        # see https://github.com/NeuralEnsemble/elephant/issues/546
+        if len(x_positions.shape) == 1:
+            x_positions = np.expand_dims(x_positions, axis=1)
         chrg_x = x
         csd = csd_profile(chrg_x)
         pots = integrate_1D(x_positions, chrg_x, csd, h)

elephant/test/test_csd.py → elephant/test/test_current_source_density.py

@@ -15,6 +15,7 @@
 import quantities as pq
 from elephant import current_source_density as csd
 import elephant.current_source_density_src.utility_functions as utils
+from elephant.current_source_density import generate_lfp
 
 available_1d = ['StandardCSD', 'DeltaiCSD', 'StepiCSD', 'SplineiCSD', 'KCSD1D']
 available_2d = ['KCSD2D', 'MoIKCSD']
@@ -153,5 +154,29 @@ def test_kcsd2d_init(self):
         self.assertEqual(len(result.times), 1)
 
 
+class GenerateLfpTestCase(unittest.TestCase):
+    @classmethod
+    def setUpClass(cls) -> None:
+        cls.one_dimensional = np.linspace(0, 10, 2304)
+        cls.two_dimensional = np.linspace(0, 10, 2304
+                                          ).reshape(2304, 1)
+
+    def test_generate_lfp_one_dimensional_array(self):
+        """
+        Regression test for Issue #546,
+        see: https://github.com/NeuralEnsemble/elephant/issues/546
+        """
+        # this should raise NOT an error
+        generate_lfp(utils.gauss_1d_dipole, self.one_dimensional)
+
+    def test_generate_lfp_two_dimensional_array(self):
+        """
+        Regression test for Issue #546,
+        see: https://github.com/NeuralEnsemble/elephant/issues/546
+        """
+        # this should NOT raise an error
+        generate_lfp(utils.gauss_1d_dipole, self.two_dimensional)
+
+
 if __name__ == '__main__':
     unittest.main()