neuroinflab / kcsd-python Goto Github PK

kCSD/kcsd/utility_functions.py:189: FutureWarning: arrays to stack must be passed as a "sequence" type such as list or tuple. Support for non-sequence iterables such as generators is deprecated as of NumPy 1.16 and will raise an error in the future.
  unique_elec_pos = np.vstack({tuple(row) for row in elec_pos})

numpy: 1.16.2
python: 3.7.2

As we do not test package for Python <3.8, I think backward compatibility code like:

Across all the files, the headers must be made consistent.

Information on the authorship should be removed, and a simpler maximum 3 line description of the contents must be included.

All authorship contributions must be made to https://github.com/Neuroinflab/kCSD-python/blob/master/docs/AUTHORS.rst

increase figure size

See branch

This branch contains forward models pre-computed and saved as pickles on the fly in the home folder, for future use. This would work in python3, and not in python2, where I save the interpolated function as a pickle and resurrect it on each subsequent call. This would fail at when changing parameters of MoI kcsd method, as it has not been implemented in include the exception of sigma_tissue, sigma_saline and number of moi iterations.

If this enhancement is warranted is an open question. I notice an improvement of 6-8 seconds in 30 second simulation runs when using crossvalidation. Perhaps with L-curve we can do better. Systematic quantification of this improvement is necessary before merger of this enhancement (if it is).

A possible fix for python2 is to store the numpy arrays instead and then perform the interpolation on the fly. I noticed that this is not as efficient. So perhaps, skip this optimization for python2.

I think

    super(sKCSD, self).__init__(ele_pos, pots, **kwargs)

Right now it is skipping call to KCSD1D.__init__() which I doubt to be the intended behaviour.
Although it is mostly harmless due to the current implementation of KCSD1D class, it may result in surprising errors in case the implementation changes.

Start from master now (tagged as version 1.2), create a branch that fixes all these issues and PR back, please. please comment below if you find any of these suggested changes redundant for the sake of documentation.

1. Rename https://github.com/Neuroinflab/kCSD-python/blob/master/figures/sKCSD_paper/functions.py
2. Is this (https://github.com/Neuroinflab/kCSD-python/tree/master/figures/sKCSD_paper/morphology) a duplicate of https://github.com/Neuroinflab/kCSD-python/tree/master/data/morphology
3. https://github.com/Neuroinflab/kCSD-python/blob/master/figures/sKCSD_paper/ElcoordsDomi14.txt file belongs elsewhere?

4. kCSD-python/kcsd/__init__.py

   Line 6 in 7dd7cdf

   from .sKCSD import sKCSD, sKCSDcell

   is sKCSDCell used somewhere other than within the class itself - what is the difference between the sKCSD and sKCSDCell classes?
5. obj = Data("Data/gang_7x7_200") would fail, please use kcsd.sample_data_path instead.
6. please change this function name

   kCSD-python/kcsd/validation/plotting_functions.py

   Line 215 in 7dd7cdf

   def plot(ax_i, what, **kwargs):

7. Need atleast 1 line of docstring in these functions https://github.com/Neuroinflab/kCSD-python/blob/master/kcsd/validation/plotting_functions.py and https://github.com/Neuroinflab/kCSD-python/blob/master/kcsd/utility_functions.py

8. kCSD-python/kcsd/utility_functions.py

   Line 41 in 7dd7cdf

   sys.path.insert(0, os.path.abspath(os.path.join(os.path.dirname(__file__),

   why do this?

Cheers!

https://github.com/Neuroinflab/kCSD-python/blob/826be34c081dd6af4ffdb8439589c2e3f163b614/kcsd/KCSD.py#L425

utility_functions.py imports kcsd
kcsd imports utility_functions

This creates an error when importing uninstalled code, and is a bad practice as far as I am aware.

Inconsistent point_errors formulations in tutorial_basic.ipynb, ValidateKCSD.py and the kCSD paper.

Fix the text - new section for tutorial, skcsd paper citation is wrong.

As pointed by Asia, the figures do not have the units

Could not solve for environment specs
The following packages are incompatible
├─ libgcc-ng ==11.2.0 h1234567_1 is requested and can be installed;
└─ libsqlite is not installable because it requires
└─ libgcc-ng >=12 , which conflicts with any installable versions previously reported.
time: 47.690
Removing intermediate container 0a2f091d0a1f
The command '/bin/sh -c TIMEFORMAT='time: %3R' bash -c 'time ${MAMBA_EXE} env update -p ${NB_PYTHON_PREFIX} --file "binder/environment.yml" && time ${MAMBA_EXE} clean --all -f -y && ${MAMBA_EXE} list -p ${NB_PYTHON_PREFIX} '' returned a non-zero code: 1

Matrix inversion is considered numerically harmful, but KCSD.values() uses it.

I have cooked a ground truth CSD which demonstrates the issue (code below).

The pictures has been generated with the following code:

import numpy as np
import matplotlib.pyplot as plt
from kcsd import KCSD2D, utility_functions, csd_profile
from scipy import integrate, linalg
import matplotlib.cm as cm

def make_plot(xx, yy, zz, title, cmap=cm.bwr):
    fig = plt.figure(figsize=(7, 7))
    ax = plt.subplot(111)
    ax.set_aspect('equal')
    t_max = np.max(np.abs(zz))
    levels = np.linspace(-1 * t_max, t_max, 32)
    im = ax.contourf(xx, yy, zz, levels=levels, cmap=cmap)
    ax.set_xlabel('X (mm)')
    ax.set_ylabel('Y (mm)')
    ax.set_title(title)
    ticks = np.linspace(-1 * t_max, t_max, 3, endpoint=True)
    plt.colorbar(im, orientation='horizontal', format='%.2g', ticks=ticks)
    return ax

def integrate_2d(csd_at, true_csd, ele_pos, h, csd_lims):
    csd_x, csd_y = csd_at
    xlin = csd_lims[0]
    ylin = csd_lims[1]
    Ny = ylin.shape[0]
    m = np.sqrt((ele_pos[0] - csd_x)**2 + (ele_pos[1] - csd_y)**2)
    m[m < 0.0000001] = 0.0000001
    y = np.arcsinh(2 * h / m) * true_csd
    integral_1D = np.zeros(Ny)
    for i in range(Ny):
        integral_1D[i] = integrate.simps(y[:, i], ylin)

    integral = integrate.simps(integral_1D, xlin)
    return integral

def forward_method(ele_pos, csd_at, true_csd):
    pots = np.zeros(ele_pos.shape[0])
    xlin = csd_at[0, :, 0]
    ylin = csd_at[1, 0, :]
    h = 50. # distance between the electrode plane and the CSD plane
    conductivity = 1.0 # S/m
    for ii in range(ele_pos.shape[0]):
        pots[ii] = integrate_2d(csd_at, true_csd,
        [ele_pos[ii][0], ele_pos[ii][1]], h,
        [xlin, ylin])
    return pots / (2 * np.pi * conductivity)


xmin = 0.0
xmax = 1.0
ymin = 0.0
ymax = 1.0
n_src_init = 1000
R_init = 1.
ext_x = 0.0
ext_y = 0.0
h = 50. # distance between the electrode plane and the CSD plane
conductivity = 1.0 # S/m

def do_kcsd(ele_pos, pots):
    pots = pots.reshape((len(ele_pos), 1)) # first time point
    return KCSD2D(ele_pos, pots, h=h, sigma=conductivity,
                  xmin=xmin, xmax=xmax,
                  ymin=ymin, ymax=ymax,
                  n_src_init=n_src_init,
                  src_type='gauss',
                  R_init=R_init)

csd_at = np.mgrid[0.:1.:100j,
                  0.:1.:100j]
csd_x, csd_y = csd_at
csd_pos = np.vstack((csd_x.flatten(), csd_y.flatten())).T

est_to_ele = slice(35, 65, 3)
ele_x, ele_y = csd_at[:, est_to_ele, est_to_ele]
ele_pos = np.vstack((ele_x.flatten(), ele_y.flatten())).T

k = do_kcsd(ele_pos, np.zeros(ele_pos.shape[0]))
kernel = k.k_pot
ck = k.k_interp_cross
true_csd = k.process_estimate(np.dot(ck, np.dot(kernel, linalg.svd(np.dot(np.linalg.inv(kernel), kernel))[2][0].reshape((-1, 1)))))[:, :, 0]

pots = forward_method(ele_pos, csd_at, true_csd)
true_pots = pots.reshape(ele_x.shape)

k = do_kcsd(ele_pos, pots)
est_csd = k.values('CSD')
kcsd_pot = k.values('POT')

make_plot(csd_x, csd_y, true_csd, 'True CSD')
make_plot(ele_x, ele_y, true_pots, 'True POT', cmap=cm.PRGn)

make_plot(csd_x, csd_y, est_csd[:,:,0], 'kCSD CSD')
make_plot(ele_x, ele_y, kcsd_pot[est_to_ele, est_to_ele, 0], 'kCSD POT', cmap=cm.PRGn)

make_plot(csd_x, csd_y, true_csd - est_csd[:,:,0], 'True - kCSD CSD')
make_plot(ele_x, ele_y, true_pots - kcsd_pot[est_to_ele, est_to_ele, 0], 'True - kCSD POT', cmap=cm.PRGn)
plt.show()

Tutorials on binders don't seem to be working at the moment. verify and tag.

from kcsd import generate as utils

This is old API, need to be updated with an equivalent one to the one in the advanced tutorial

val = config.ValidateKCSD(dim=int(dim_select.value[0]))

When calling KCSD.py , for example python kCSD-python/figures/kcsd_properties/tutorial1.py, why are there multiple processes being called open. Where is this coming from?

try htop when running this file. Why are there multiple nodes opening up? This should not happen.

• figures/kCSD_properties/L_curve_simulation.py - throws error: NameError: name 'noise_lvl' is not defined. I think function make_plot_perf (line 23 from figure_LCandCVperformance.py) should be rewritten and take 'noise_lvl' as an argument to solve the issue;
  line 121: polish name is used

• figures/kCSD_properties/figure_LC.py - it throws error: NameError: name 'df' is not defined. Uncomment lines 103 and 106 and change name from 'data_fig4_and_fig13_lc_noise25.0.npz' to 'data_fig4_and_fig13_LC_noise25.0.npz' in line 106;
  line 108: polish name is used

check_for_duplicated_electrodes function from kcsd/utility_functions.py gives a FutureWarning:
arrays to stack must be passed as a "sequence" type such as list or tuple. Support for non-sequence iterables such as generators is deprecated as of NumPy 1.16 and will raise an error in the future.
unique_elec_pos = np.vstack({tuple(row) for row in elec_pos})

1. estimate='CSD'is a parameter of L_curve, it is not used anywhere.
2. L_curve does not return anything despite the statement in the docstrings.
3. decimal=3 in prints, skcsd uses source widths of the order of tens of microns.
4. Rs is converted to an array, lambdas isn't, it would be quite useful, if it was.

#35 provides useful scripts

As reported by a Christoph H via email

When running pip install . in windows (untested in other OS) this error shows up.

File 'setup.py', line 12, in readme
return f.read()
'charmap codec can't decode byte 0x81 in position 154 ...'

This is already at the setup.py level

Suggested fix, to add a line at the very beginning of setup.py

# -- coding: utf-8 --

Any chance that the 2.0 release could be published to PyPI? I have a package foo that depends on kcsd, but I cannot upload foo to PyPI, because PyPI will not accept git URL dependencies. Of course I could fork and publish myself, but I maybe a project author would prefer to? Thanks -Graham

we need to improve the code coverage of the tests.

Files from ipynb_test branch are not compatible with Python 3:

• widget_helpers.py (dict.keys() return views instead of lists in Python 3)
• tutorial_helpers.py (replace xrange() with range())

Pertains to ipynb_tests branch

1. Whats the difference between kCSD-python/data/ball_and_stick_8 and kCSD/tutorials/Data/ball_and_stick_8 ? All data now moved to kCSD-python/data/* so use this accordingly. [Same for gang_7x7_200] - urgent.

2. Why kCSD-python/kcsd/sKCSD.py and kCSD-python/kcsd/sKCSDCell.py - I am merging these into one single file.

3. kcsd/plotting_functions.py . See if this can me moved / added into the kcsd/validation/ part of the code. Please avoid plotting functions here. @m-kowalska perhaps you can look into this?

4. All figures now moved to kCSD-python/figures/sKCSD_paper, tutorials folder reserved for ipython notebooks. If you'd like to write a tutorial that'd be great. Like wise all the figures for the current paper would in this folder. - Also would be useful to generate these said figures and include them here as png - low priority. Add a short readme pointing to the doi of the skcsd paper

5. kCSD-python/figures/sKCSD/functions.py - is a poor choice for a file name.

6. pep8 convention is missing. Please fix.

7. In the class definition of sKCSD.py inherits KCSD1D but the documentation reads differently. Which is it? There also seems to be some functions from the original class re-written here - is there some reason for this (for example def forward_model).

8. Remove all sys.path.append for importing kcsd. It is now a package. Use, import kcsd or from kcsd import KCSD1D etc.

More later.
@asiaszmek @danek8317

I think that the magnitudes should be shown as here, as an inset plot - as most of this is otherwise an empty, which does not contain much information. Figure B, previously is now Figure A, inset instead. Open to criticism.

Can you look into the alignment of the y axis labels, plus the titles compatible with 3.5 python version f'string' doesn't work in 3.5.

Cheers!

https://numpy.org/devdocs/release/1.20.0-notes.html#deprecations

Either numpy requirement should be pinned <1.20.0, or code should be updated to use builtin types (e.g. complex rather than np.complex).

1. Remove splines on top and right for all the figures
2. Remove box around the legend - it legend should be below the figures, without a bounding box.
3. Follow the https://github.com/Neuroinflab/kCSD-python/blob/master/figures/kCSD_properties/figure_properties.py and use gridspec as in https://github.com/Neuroinflab/kCSD-python/blob/master/figures/kCSD_properties/figure_vectors.py - consistency in font size and font styles across all figures in the paper.
4. Recommend rearranging so that Topleft is potentials, bottom left is truecsd/ kCSD, top right is L curve bottom right is CV.
5. Increase the gap between the rows of the figure
6. Align the y-axis labels
7. Increase the resolution to 300 dpi

0. Use figure properties and gridspec as above
1. Remove splines top and right for both subplots. Top subplot has a shared x-axis with the bottom - consider something like figure_vectors.py where only the last row has the necessary x-axis units and labels.
2. Legend must be below the figures, without a bounding box.
3. Align the y axis labels.

1. Why does the cbar here have both positive and negative values when the countour has only positive values? Move the cbar below, like in the error map figures.
2. For this subplot 1 and 2 - combine these into one figure - use gridspect and figure properties for the font size and font styles.

1. Move cbar below the figure.
2. The red-ellipse is not as apparent.

try:
from . import sKCSD, sKCSDcell

Do we want to include this case? It requires to include additional package in the requirements and installation file.

possible solution:
try:
from skmonaco import mcmiser
skmonaco_available = True
import multiprocessing
num_cores = multiprocessing.cpu_count()
except ImportError:
skmonaco_available = False

According to:

In make_column_plot function:
1st, 2nd and 4th subplots show all of the currents
3rd - potentials, show potentials in moment defined by the variable time_pt_interest - that is 3000, in the middle of 6000 long signal.

ax1 = plt.subplot(241, aspect='equal')
plot_all_currents(ax1, xmin, xmax, ymin, ymax, all_x, all_y, all_z, all_val, letter='A')

ax2 = plt.subplot(242, aspect='equal')
plot_csd_slice(ax2, xmin, xmax, ymin, ymax, all_x, all_y, all_z, all_val, letter='B')

ax3 = plt.subplot(243, aspect='equal')
plot_dense_potentials(ax3, h, pop_names, time_pts, time_pt_interest,
letter='C', filt=False)

ax4 = plt.subplot(244, aspect='equal')
xx, yy = np.mgrid[xmin:xmax:np.complex(0, true_csd.shape[0]),
ymin:ymax:np.complex(0, true_csd.shape[1])]
plot_csd_smooth(ax4, xmin, xmax, ymin, ymax, true_csd[:, :],
xx, yy, letter='D')

Later, in the same function:

plot_csd_smooth(ax_list1[i], xmin, xmax, ymin, ymax, kcsd[:, :, 750], x, y,
letter=letter_list1[i], ele_pos=ele_pos)

The kcsd is plotted in time point of 750.

1. Inconsistent documentation of the class 'nr_basis' for instance. 'mask', 'kcsd_xlims', is not mentioned. 'ele_xlims' is in the documentation but not in the kwargs pops / also only self.ele_lims seems to be used.
2. Why are the remaining parameters passed in the make_reconstruction function as opposed to the original function - what is the justification this separation?
3. Whats with the csd_seed being passed in make_reconstruction and also in the class init for ValidateKCSD1D, 2D, 3D etc? Please consider writing them as separate functions if you want to update them. Like def update_csd_seed(): , def update_ele_lims()., csd_profile is another. It is confusing to pass these values twice.
4. consider changing the name of function 'recon' - its similar to make_reconstruction - and both call the KCSD classes in the end.

This class needs more code review.

I think the following imports are unnecessary:

According to https://docs.python.org/3/library/builtins.html, the builtins module:

can be useful in modules that provide objects with the same name as a built-in value, but in which the built-in of that name is also needed

which seems not to be the case.

https://github.com/Neuroinflab/kCSD-python/blob/master/__init__.py

Needs to be removed.

dim_selection: after pressing 1D or 2D button and then 3D button 'UnboundLocalError: local variable 'big_wid' referenced before assignment' occurs

Hello all,

Please install the latest version and raise issues regarding any code quality/bugs and anything you find about this package.

see where the documentation could be better etc.

Any comments on the tutorials would also help. Perhaps get some naive users to try?

I'll leave this for Marta to coordinate.

All through use 300 dpi resolution

1. Please upload the scripts to generate these figures
2. Use gridspec and move the cbar to occupy the whole two rows.
3. Use figure_properties and for using consistent font size and font styles.
4. The 'x labels' 'depth' can be omitted for the first row
5. The subplot labels 'A', 'B' etc must be figure labels and not inside the title itself.

Side question - are these same as the error calculation plots in the tutorial? If so we should use consistent cmap.

1 )The numbers here on the contours are not visible at all - either remove them or use a different color - like blue.
2) The cbar needs units
3) Can you try transparency instead of overlay of two different contours?

1. Use figure properties for the font size and font consistency
2. D you can reduce the cbar range so that the darks are more visible. You can use colorbar's extend parameter (https://matplotlib.org/gallery/images_contours_and_fields/image_masked.html#sphx-glr-gallery-images-contours-and-fields-image-masked-py)

Can we get the documentation to be autogenerated like in numpy?

see numpy-docs package that we include.

If we add numba as a requirement instead of skmonaco, we can then then deploy with just-in-time compiler
https://numba.pydata.org/numba-doc/dev/user/5minguide.html
Adding this line above all functions in basis_functions.py and before forward_method() and int_pot1d() / 2d() / 3d() in KCSD.py

@jit(nopython=True)

already improves the performance drastically, I am able to knock off ~7 secs from a 30 sec tests.

is IMHO mathematically equivalent with:

estimation[:, t] = np.dot(estimation_table, beta)

which should be faster.

However, the obtained results may differ due to numerical issues.

Hello
I am trying to use the KCSD package on some toy data. I used the code below and for a LFP vector size of (9,1), an estimated CSD size of (100,1) was given as an output. Why is that? Should nt they be of the same size?

from kcsd import KCSD1D
import numpy as np
from matplotlib import pyplot as plt

elec_pos = np.array([[-0.5], [0], [1], [1.5], [3.5], [4.1], [5.0], [7.0], [8.0]])

pots = np.array([[-0.1], [0.3], [-0.4], [0.2], [0.8], [0.5], [0.2], [0.5], [0.6]])

def do_kcsd(elec_pos, pots):
h = 1.
sigma = 1.0
pots = pots.reshape((len(elec_pos), 1)) # first time point
k = KCSD1D(elec_pos, pots,h=h, sigma=sigma,
xmin=0.0, xmax=1.0, src_type='gauss', R_init=1.)
return k

k = do_kcsd(elec_pos, pots)
est_csd = k.values('CSD')
est_pot = k.values('POT')

I think the docstrings should be inspected and updated.

For example docstring for KCSD.create_lookup() suggests it has a parameter dist_table_density (it has not; it is provided as attribute):