Skip to content

Commit

Permalink
Merge pull request #166 from espenhgn/cythonfix
Browse files Browse the repository at this point in the history 
fix Cython FutureWarning
  • Loading branch information
@espenhgn
espenhgn committed Mar 25, 2020
2 parents cfeeb04 + 2bd53b1 commit b4af4e0
Show file tree
Hide file tree
Showing 2 changed files with 62 additions and 61 deletions.
31 changes: 16 additions & 15 deletions LFPy/alias_method.pyx
View file
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
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

0 comments on commit b4af4e0

Please sign in to comment.