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.

Sign up for GitHub

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

Jump to bottom

fix Cython FutureWarning #166

Merged
merged 1 commit into from
Mar 25, 2020
Merged

fix Cython FutureWarning #166

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
31 changes: 16 additions & 15 deletions LFPy/alias_method.pyx
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,5 +1,6 @@
#!/usr/bin/env python
# -*- coding: utf-8 -*-
# cython: language_level=2

from __future__ import division
import numpy as np
Expand All @@ -22,7 +23,7 @@ cpdef np.ndarray[long, ndim=1, negative_indices=False] alias_method(np.ndarray[l
"""
Alias method for drawing random numbers from a discrete probability
distribution. See http://www.keithschwarz.com/darts-dice-coins/

Parameters
----------
idx : np.ndarray
Expand All @@ -42,31 +43,31 @@ cpdef np.ndarray[long, ndim=1, negative_indices=False] alias_method(np.ndarray[l
assert idx.size == probs.size
except AssertionError as ae:
raise ae('length of idx and probs arrays must be equal')

#C-declare variables
cdef np.ndarray[long, ndim=1, negative_indices=False] J, spc
cdef np.ndarray[DTYPE_t, ndim=1, negative_indices=False] q
cdef np.ndarray[DTYPE_t, ndim=2, negative_indices=False] rands
cdef int nn, j, ad, K, kk

# Construct the table.
J, q = alias_setup(probs)

#output array
spc = np.zeros(nsyn, dtype=int)

#prefetch random numbers, alias_draw needs nsyn x 2 numbers
rands = np.random.rand(nsyn, 2)
K = J.size

K = J.size
# Generate variates using alias draw method
for nn in range(nsyn):
kk = floor(rands[nn, 0]*K)
if rands[nn, 1] < q[kk]:
spc[nn] = idx[kk]
else:
spc[nn] = idx[J[kk]]

return spc


Expand All @@ -75,7 +76,7 @@ cpdef np.ndarray[long, ndim=1, negative_indices=False] alias_method(np.ndarray[l
cpdef alias_setup(np.ndarray[DTYPE_t, ndim=1, negative_indices=False] probs):
"""Set up function for alias method.
See http://www.keithschwarz.com/darts-dice-coins/

Parameters
----------
probs : np.ndarray
Expand All @@ -95,7 +96,7 @@ cpdef alias_setup(np.ndarray[DTYPE_t, ndim=1, negative_indices=False] probs):
cdef long K
cdef int small, large, kk, s_i, l_i
cdef DTYPE_t prob

K = probs.size
q = probs*K
J = np.zeros(K, dtype=int)
Expand All @@ -113,25 +114,25 @@ cpdef alias_setup(np.ndarray[DTYPE_t, ndim=1, negative_indices=False] probs):
else:
larger[l_i] = kk
l_i += 1

s_i -= 1
l_i -= 1

# Loop though and create little binary mixtures that
# appropriately allocate the larger outcomes over the
# overall uniform mixture.
while s_i >= 0 and l_i >= 0:
small = smaller[s_i]
large = larger[l_i]

J[small] = large
q[large] = q[large] + q[small] - 1

s_i -= 1

if q[large] < 1:
s_i += 1
l_i -= 1
smaller[s_i] = large

return J, q
92 changes: 46 additions & 46 deletions LFPy/run_simulation.pyx
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,5 +1,6 @@
#!/usr/bin/env python
# -*- coding: utf-8 -*-
# cython: language_level=2
"""Copyright (C) 2012 Computational Neuroscience Group, NMBU.

This program is free software: you can redistribute it and/or modify
Expand Down Expand Up @@ -38,22 +39,22 @@ def _run_simulation(cell, cvode, variable_dt=False, atol=0.001):
cvode.atol(atol)
else:
cvode.active(0)

#re-initialize state
neuron.h.finitialize(cell.v_init)

#initialize current- and record
if cvode.active():
cvode.re_init()
else:
neuron.h.fcurrent()
neuron.h.frecord_init()

#Starting simulation at t != 0
neuron.h.t = cell.tstart

cell._loadspikes()

#print sim.time at intervals
cdef int counter = 0
cdef double interval
Expand All @@ -65,7 +66,7 @@ def _run_simulation(cell, cvode, variable_dt=False, atol=0.001):
interval = 1000. / cell.dt
else:
interval = 100. / cell.dt

while neuron.h.t < tstop:
neuron.h.fadvance()
counter += 1
Expand All @@ -88,7 +89,7 @@ def _run_simulation_with_electrode(cell, cvode, electrode=None,
electrode argument used to determine coefficient
matrix, and calculate the LFP on every time step.
"""

#c-declare some variables
cdef int i, j, tstep#, ncoeffs
#cdef int totnsegs = cell.totnsegs
Expand All @@ -103,7 +104,7 @@ def _run_simulation_with_electrode(cell, cvode, electrode=None,
cdef np.ndarray[DTYPE_t, ndim=2, negative_indices=False] coeffs
cdef np.ndarray[DTYPE_t, ndim=2, negative_indices=False] current_dipole_moment
cdef np.ndarray[DTYPE_t, ndim=2, negative_indices=False] midpoints

#check if h5py exist and saving is possible
try:
import h5py
Expand All @@ -112,12 +113,12 @@ def _run_simulation_with_electrode(cell, cvode, electrode=None,
to_file = False
file_name = None


# Use electrode object(s) to calculate coefficient matrices for LFP
# calculations. If electrode is a list, then
if cell.verbose:
print('precalculating geometry - LFP mapping')

#put electrodecoeff in a list, if it isn't already
if dotprodcoeffs is not None:
if type(dotprodcoeffs) != list:
Expand All @@ -126,44 +127,44 @@ def _run_simulation_with_electrode(cell, cvode, electrode=None,
else:
#create empty list if no dotprodcoeffs are supplied
dotprodcoeffs = []

#just for safekeeping
lendotprodcoeffs0 = len(dotprodcoeffs)
#access electrode object and append mapping

#access electrode object and append mapping
if electrode is not None:
#put electrode argument in list if needed
if type(electrode) == list:
electrodes = electrode
else:
electrodes = [electrode]

for el in electrodes:
el.calc_mapping(cell)
dotprodcoeffs.append(el.mapping)
elif electrode is None:
electrodes = None

# Initialize NEURON simulations of cell object

# Initialize NEURON simulations of cell object
neuron.h.dt = dt

#don't know if this is the way to do, but needed for variable dt method
if cell.dt <= 1E-8:
cvode.active(1)
cvode.atol(atol)

#re-initialize state
neuron.h.finitialize(cell.v_init)
neuron.h.frecord_init() # wrong voltages t=0 for tstart < 0 otherwise
neuron.h.fcurrent()

#Starting simulation at t != 0
neuron.h.t = cell.tstart

#load spike times from NetCon
cell._loadspikes()

#print sim.time at intervals
counter = 0
tstep = 0
Expand All @@ -173,7 +174,7 @@ def _run_simulation_with_electrode(cell, cvode, electrode=None,
interval = 1000. / dt
else:
interval = 100. / dt

#temp vector to store membrane currents at each timestep
imem = np.zeros(cell.totnsegs)
#LFPs for each electrode will be put here during simulation
Expand All @@ -200,7 +201,7 @@ def _run_simulation_with_electrode(cell, cvode, electrode=None,
current_dipole_moment = cell.current_dipole_moment.copy()
cell.current_dipole_moment = np.array([[]])
midpoints = np.c_[cell.xmid, cell.ymid, cell.zmid]

#run fadvance until time limit, and calculate LFPs for each timestep
while neuron.h.t < tstop:
if neuron.h.t >= 0:
Expand All @@ -216,7 +217,7 @@ def _run_simulation_with_electrode(cell, cvode, electrode=None,
if to_memory:
for j, coeffs in enumerate(dotprodcoeffs):
electrodesLFP[j][:, tstep] = np.dot(coeffs, imem)

if to_file:
for j, coeffs in enumerate(dotprodcoeffs):
el_LFP_file['electrode{:03d}'.format(j)
Expand All @@ -232,7 +233,7 @@ def _run_simulation_with_electrode(cell, cvode, electrode=None,
rtfactor))
t0 = time()
ti = neuron.h.t

try:
#calculate LFP after final fadvance()
i = 0
Expand All @@ -254,11 +255,11 @@ def _run_simulation_with_electrode(cell, cvode, electrode=None,

except:
pass

# update current dipole moment values
if rec_current_dipole_moment:
cell.current_dipole_moment = current_dipole_moment

# Final step, put LFPs in the electrode object, superimpose if necessary
# If electrode.perCellLFP, store individual LFPs
if to_memory:
Expand All @@ -278,34 +279,34 @@ def _run_simulation_with_electrode(cell, cvode, electrode=None,
electrodes[j-lendotprodcoeffs0].CellLFP = []
electrodes[j-lendotprodcoeffs0].CellLFP.append(LFP)
electrodes[j-lendotprodcoeffs0].electrodecoeff = dotprodcoeffs[j]

if to_file:
el_LFP_file.close()


cpdef _collect_geometry_neuron(cell):
"""Loop over allseclist to determine area, diam, xyz-start- and
endpoints, embed geometry to cell object"""


cdef np.ndarray[DTYPE_t, ndim=1, negative_indices=False] areavec = np.zeros(cell.totnsegs)
cdef np.ndarray[DTYPE_t, ndim=1, negative_indices=False] diamvec = np.zeros(cell.totnsegs)
cdef np.ndarray[DTYPE_t, ndim=1, negative_indices=False] lengthvec = np.zeros(cell.totnsegs)

cdef np.ndarray[DTYPE_t, ndim=1, negative_indices=False] xstartvec = np.zeros(cell.totnsegs)
cdef np.ndarray[DTYPE_t, ndim=1, negative_indices=False] xendvec = np.zeros(cell.totnsegs)
cdef np.ndarray[DTYPE_t, ndim=1, negative_indices=False] ystartvec = np.zeros(cell.totnsegs)
cdef np.ndarray[DTYPE_t, ndim=1, negative_indices=False] yendvec = np.zeros(cell.totnsegs)
cdef np.ndarray[DTYPE_t, ndim=1, negative_indices=False] zstartvec = np.zeros(cell.totnsegs)
cdef np.ndarray[DTYPE_t, ndim=1, negative_indices=False] zendvec = np.zeros(cell.totnsegs)

cdef DTYPE_t gsen2, secL
cdef LTYPE_t counter, nseg, n3d, i


cdef np.ndarray[DTYPE_t, ndim=1, negative_indices=False] L, x, y, z, segx, segx0, segx1


counter = 0

#loop over all segments
Expand All @@ -325,17 +326,17 @@ cpdef _collect_geometry_neuron(cell):
x[i] = neuron.h.x3d(i)
y[i] = neuron.h.y3d(i)
z[i] = neuron.h.z3d(i)

#normalize as seg.x [0, 1]
L /= secL

#temporary store position of segment midpoints
segx = np.zeros(nseg)
i = 0
for seg in sec:
segx[i] = seg.x
i += 1

#can't be >0 which may happen due to NEURON->Python float transfer:
#segx0 = (segx - gsen2).round(decimals=6)
#segx1 = (segx + gsen2).round(decimals=6)
Expand All @@ -345,10 +346,10 @@ cpdef _collect_geometry_neuron(cell):
#fill vectors with interpolated coordinates of start and end points
xstartvec[counter:counter+nseg] = np.interp(segx0, L, x)
xendvec[counter:counter+nseg] = np.interp(segx1, L, x)

ystartvec[counter:counter+nseg] = np.interp(segx0, L, y)
yendvec[counter:counter+nseg] = np.interp(segx1, L, y)

zstartvec[counter:counter+nseg] = np.interp(segx0, L, z)
zendvec[counter:counter+nseg] = np.interp(segx1, L, z)

Expand All @@ -359,18 +360,17 @@ cpdef _collect_geometry_neuron(cell):
lengthvec[counter] = secL/nseg

counter += 1


#set cell attributes
cell.xstart = xstartvec
cell.ystart = ystartvec
cell.zstart = zstartvec

cell.xend = xendvec
cell.yend = yendvec
cell.zend = zendvec

cell.area = areavec
cell.diam = diamvec
cell.length = lengthvec