Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

[Enh] add example and handle one-dimensional arrays in CSD.generate lfp #594

Merged
merged 13 commits into from
Oct 24, 2023
41 changes: 33 additions & 8 deletions elephant/current_source_density.py
Original file line number Diff line number Diff line change
Expand Up @@ -79,8 +79,8 @@ def estimate_csd(lfp, coordinates='coordinates', method=None,
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
Expand Down Expand Up @@ -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,
Expand Down Expand Up @@ -227,7 +227,7 @@ def generate_lfp(csd_profile, x_positions, y_positions=None, z_positions=None,
1D : gauss_1d_dipole
2D : large_source_2D and small_source_2D
3D : gauss_3d_dipole
x_positions : np.ndarray
x_positions : np.ndarray # TODO: add dimensionality for x-y-z_positions
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I assume this TODO should be addressed in another future PR and not here?

Copy link
Member Author

@Moritz-Alexander-Kern Moritz-Alexander-Kern Oct 20, 2023

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thank you for the reminder, I should have made this TODO more explicit.

One might consider to further discuss the expected input shape in a future PR or Issue (?)

I've now (fe7f9ec) updated the docstring to explicitly specify the expected shape of the inputs as a 2D column vector with dimensions 'N x 1 array.' In the example mentioned in the issue, an error occurred because a flat array (1D) was passed as input.

For now, documenting the expected input is a reasonable solution ?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yup, sounds good!

Positions of the x coordinates of the electrodes
y_positions : np.ndarray, optional
Positions of the y coordinates of the electrodes
Expand All @@ -253,10 +253,31 @@ def generate_lfp(csd_profile, x_positions, y_positions=None, z_positions=None,

Returns
-------
LFP : neo.AnalogSignal
LFP : :class:`neo.core.AnalogSignal`
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This change should also apply in line 74, right?

lfp : neo.AnalogSignal

(Can't directly suggest the change here, because it's too far from any lines changed in this PR.

Copy link
Member Author

@Moritz-Alexander-Kern Moritz-Alexander-Kern Oct 20, 2023

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Indeed, thanks for catching this, changed accordingly.

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):
Expand Down Expand Up @@ -297,6 +318,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)
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -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']
Expand Down Expand Up @@ -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 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 raise an error
generate_lfp(utils.gauss_1d_dipole, self.two_dimensional)


if __name__ == '__main__':
unittest.main()