You can not select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
734 lines
25 KiB
734 lines
25 KiB
from __future__ import print_function |
|
|
|
#!/usr/bin/python |
|
# |
|
# utilities for NEURON, in Python |
|
# Module neuron for cnmodel |
|
|
|
import numpy as np |
|
import numpy.ma as ma # masked array |
|
import re, sys, gc, collections |
|
|
|
import neuron |
|
|
|
|
|
_mechtype_cache = None |
|
|
|
|
|
def all_mechanism_types(): |
|
"""Return a dictionary of all available mechanism types. |
|
|
|
Each dictionary key is the name of a mechanism and each value is |
|
another dictionary containing information about the mechanism:: |
|
|
|
mechanism_types = { |
|
'mech_name1': { |
|
'point_process': bool, |
|
'artificial_cell': bool, |
|
'netcon_target': bool, |
|
'has_netevent': bool, |
|
'internal_type': int, |
|
'globals': {name:size, ...}, |
|
'parameters': {name:size, ...}, |
|
'assigned': {name:size, ...}, |
|
'state': {name:size, ...}, |
|
}, |
|
'mech_name2': {...}, |
|
'mech_name3': {...}, |
|
... |
|
} |
|
|
|
* point_process: False for distributed mechanisms, True for point |
|
processes and artificial cells. |
|
* artificial_cell: True for artificial cells, False otherwise |
|
* netcon_target: True if the mechanism can receive NetCon events |
|
* has_netevent: True if the mechanism can emit NetCon events |
|
* internal_type: Integer specifying the NEURON internal type index of |
|
the mechanism |
|
* globals: dict of the name and vector size of the mechanism's global |
|
variables |
|
* parameters: dict of the name and vector size of the mechanism's |
|
parameter variables |
|
* assigned: dict of the name and vector size of the mechanism's |
|
assigned variables |
|
* state: dict of the name and vector size of the mechanism's state |
|
variables |
|
|
|
|
|
Note: The returned data structure is cached; do not modify it. |
|
|
|
For more information on global, parameter, assigned, and state |
|
variables see: |
|
http://www.neuron.yale.edu/neuron/static/docs/help/neuron/nmodl/nmodl.html |
|
""" |
|
global _mechtype_cache |
|
if _mechtype_cache is None: |
|
_mechtype_cache = collections.OrderedDict() |
|
mname = neuron.h.ref("") |
|
# Iterate over two mechanism types (distributed, point/artificial) |
|
for i in [0, 1]: |
|
mt = neuron.h.MechanismType(i) |
|
nmech = int(mt.count()) |
|
# Iterate over all mechanisms of this type |
|
for j in range(nmech): |
|
mt.select(j) |
|
mt.selected(mname) |
|
|
|
# General mechanism properties |
|
name = mname[0] # convert hoc string ptr to python str |
|
|
|
desc = { |
|
"point_process": bool(i), |
|
"netcon_target": bool(mt.is_netcon_target(j)), |
|
"has_netevent": bool(mt.has_net_event(j)), |
|
"artificial_cell": bool(mt.is_artificial(j)), |
|
"internal_type": int(mt.internal_type()), |
|
} |
|
|
|
# Collect information about 4 different types of variables |
|
for k, ptype in [ |
|
(-1, "globals"), |
|
(1, "parameters"), |
|
(2, "assigned"), |
|
(3, "state"), |
|
]: |
|
desc[ptype] = {} # collections.OrderedDict() |
|
ms = neuron.h.MechanismStandard(name, k) |
|
for l in range(int(ms.count())): |
|
psize = ms.name(mname, l) |
|
pname = mname[0] # parameter name |
|
desc[ptype][pname] = int(psize) |
|
|
|
# Assemble everything in one place |
|
_mechtype_cache[name] = desc |
|
|
|
return _mechtype_cache |
|
|
|
|
|
def reset(raiseError=True): |
|
"""Introspect the NEURON kernel to verify that no objects are left over |
|
from previous simulation runs. |
|
""" |
|
# Release objects held by an internal buffer |
|
# See https://www.neuron.yale.edu/phpBB/viewtopic.php?f=2&t=3221 |
|
neuron.h.Vector().size() |
|
|
|
# Make sure nothing is hanging around in an old exception or because of |
|
# reference cycles |
|
|
|
# sys.exc_clear() |
|
gc.collect(2) |
|
neuron.h.Vector().size() |
|
numsec = 0 |
|
|
|
remaining = [] |
|
n = len(list(neuron.h.allsec())) |
|
|
|
if n > 0: |
|
remaining.append((n, "Section")) |
|
|
|
n = len(neuron.h.List("NetCon")) |
|
if n > 0: |
|
remaining.append((n, "NetCon")) |
|
|
|
# No point processes or artificial cells left |
|
for name, typ in all_mechanism_types().items(): |
|
if typ["artificial_cell"] or typ["point_process"]: |
|
n = len(neuron.h.List(name)) |
|
if n > 0: |
|
remaining.append((n, name)) |
|
|
|
if ( |
|
len(remaining) > 0 and raiseError |
|
): # note that not raising the error leads to memory leak |
|
msg = "Cannot reset--old objects have not been cleared: %s" % ", ".join( |
|
["%d %s" % rem for rem in remaining] |
|
) |
|
raise RuntimeError(msg) |
|
|
|
|
|
def custom_init(v_init=-60.0): |
|
""" |
|
Perform a custom initialization of the current model/section. |
|
|
|
This initialization follows the scheme outlined in the |
|
NEURON book, 8.4.2, p 197 for initializing to steady state. |
|
|
|
N.B.: For a complex model with dendrites/axons and different channels, |
|
this initialization will not find the steady-state the whole cell, |
|
leading to initial transient currents. In that case, this initialization |
|
should be followed with a 0.1-5 second run (depends on the rates in the |
|
various channel mechanisms) with no current injection or active point |
|
processes to allow the system to settle to a steady- state. Use either |
|
h.svstate or h.batch_save to save states and restore them. Batch is |
|
preferred |
|
|
|
Parameters |
|
---------- |
|
v_init : float (default: -60 mV) |
|
Voltage to start the initialization process. This should |
|
be close to the expected resting state. |
|
""" |
|
inittime = -1e10 |
|
tdt = neuron.h.dt # save current step size |
|
dtstep = 1e9 |
|
neuron.h.finitialize(v_init) |
|
neuron.h.t = inittime # set time to large negative value (avoid activating |
|
# point processes, we hope) |
|
tmp = neuron.h.cvode.active() # check state of variable step integrator |
|
if tmp != 0: # turn off CVode variable step integrator if it was active |
|
neuron.h.cvode.active(0) # now just use backward Euler with large step |
|
neuron.h.dt = dtstep |
|
n = 0 |
|
while neuron.h.t < -1e9: # Step forward |
|
neuron.h.fadvance() |
|
n += 1 |
|
# print('advances: ', n) |
|
if tmp != 0: |
|
neuron.h.cvode.active(1) # restore integrator |
|
neuron.h.t = 0 |
|
if neuron.h.cvode.active(): |
|
neuron.h.cvode.re_init() # update d(state)/dt and currents |
|
else: |
|
neuron.h.fcurrent() # recalculate currents |
|
neuron.h.frecord_init() # save new state variables |
|
neuron.h.dt = tdt # restore original time step |
|
|
|
|
|
# routine to convert conductances from nS as given elsewhere |
|
# to mho/cm2 as required by NEURON 1/28/99 P. Manis |
|
# units: nano siemens, soma area in um^2 |
|
# |
|
def nstomho(ns, somaarea, refarea=None): |
|
if refarea == None: |
|
return 1e-9 * float(ns) / float(somaarea) |
|
else: |
|
return 1e9 * float(ns) / float(refarea) |
|
|
|
|
|
def mho2ns(mho, somaarea): |
|
return float(mho) * somaarea / 1e-9 |
|
|
|
|
|
def spherearea(dia): |
|
""" |
|
given diameter in microns, return sphere area in cm2 |
|
""" |
|
r = dia * 1e-4 # convert to cm |
|
return 4 * np.pi * r ** 2 |
|
|
|
|
|
def get_sections(h): |
|
""" |
|
go through all the sections and find the names of the sections and all of their |
|
parts (ids). Returns a dict, of sec: [id0, id1...] |
|
|
|
""" |
|
secnames = {} |
|
resec = re.compile("(\w+)\[(\d*)\]") |
|
for sec in h.allsec(): |
|
g = resec.match(sec.name()) |
|
if g.group(1) not in secnames.keys(): |
|
secnames[g.group(1)] = [int(g.group(2))] |
|
else: |
|
secnames[g.group(1)].append(int(g.group(2))) |
|
return secnames |
|
|
|
|
|
def all_objects(): |
|
""" Return a dict of all objects known to NEURON. |
|
|
|
Keys are 'Section', 'Segment', 'Mechanism', 'Vector', 'PointProcess', |
|
'NetCon', ... |
|
""" |
|
objs = {} |
|
objs["Section"] = list(h.all_sec()) |
|
objs["Segment"] = [] |
|
for sec in objs["Section"]: |
|
objs["Segment"].extend(list(sec.allseg())) |
|
objs["PointProcess"] = [] |
|
for seg in objs["Segment"]: |
|
objs["PointProcess"].extend(list(seg.point_processes())) |
|
|
|
return objs |
|
|
|
|
|
def alpha(alpha=0.1, delay=1, amp=1.0, tdur=50.0, dt=0.010): |
|
tvec = np.arange(0, tdur, dt) |
|
aw = np.zeros(tvec.shape) |
|
i = 0 |
|
for t in tvec: |
|
if t > delay: |
|
aw[i] = ( |
|
amp * (t - delay) * (1.0 / alpha) * np.exp(-(t - delay) / alpha) |
|
) # alpha waveform time course |
|
else: |
|
aw[i] = 0.0 |
|
i += 1 |
|
return (aw, tvec) |
|
|
|
|
|
def syns( |
|
alpha=0.1, |
|
rate=10, |
|
delay=0, |
|
dur=50, |
|
amp=1.0, |
|
dt=0.020, |
|
N=1, |
|
mindur=120, |
|
makewave=True, |
|
): |
|
""" Calculate a poisson train of alpha waves |
|
with mean rate rate, with a delay and duration (in mseco) dt in msec. |
|
N specifies the number of such independent waveforms to sum """ |
|
deadtime = 0.7 |
|
if dur + delay < mindur: |
|
tvec = np.arange(0.0, mindur, dt) |
|
else: |
|
tvec = np.arange(0.0, dur + delay, dt) |
|
npts = len(tvec) |
|
ta = np.arange(0.0, 20.0, dt) |
|
aw = ta * alpha * np.exp(-ta / alpha) / alpha # alpha waveform time course |
|
spt = [[]] * N # list of spike times |
|
wave = np.array([]) # waveform |
|
sptime = [] |
|
for j in range(0, N): |
|
done = False |
|
t = 0.0 |
|
nsp = 0 |
|
while not done: |
|
a = np.random.sample(1) |
|
if t < delay: |
|
t = delay |
|
continue |
|
if t >= delay and t <= (delay + dur): |
|
ti = -np.log(a) / ( |
|
rate / 1000.0 |
|
) # convert to exponential distribution with rate |
|
if ti < deadtime: |
|
continue |
|
t = t + ti # running time |
|
if t > delay + dur: |
|
done = True |
|
continue |
|
if nsp is 0: |
|
sptime = t |
|
nsp = nsp + 1 |
|
else: |
|
sptime = np.append(sptime, t) |
|
nsp = nsp + 1 |
|
if j is 0: |
|
wavej = np.zeros(len(tvec)) |
|
for i in range(0, len(sptime)): |
|
st = int(sptime[i] / dt) |
|
wavej[st] = wavej[st] + 1 |
|
spt[j] = sptime |
|
|
|
if makewave: |
|
w = np.convolve(wavej, aw / max(aw)) * amp |
|
if len(w) < npts: |
|
w = np.append(w, np.zeros(npts - len(w))) |
|
if len(w) > npts: |
|
w = w[0:npts] |
|
if j is 0: |
|
wave = w |
|
else: |
|
wave = wave + w |
|
return (spt, wave, tvec, N) |
|
|
|
|
|
def an_syn( |
|
alpha=0.1, |
|
spont=10, |
|
driven=100, |
|
delay=50, |
|
dur=100, |
|
post=20, |
|
amp=0.1, |
|
dt=0.020, |
|
N=1, |
|
makewave=True, |
|
): |
|
# constants for AN: |
|
deadtime = 0.7 # min time between spikes, msec |
|
trise = 0.2 # rise rate, ms |
|
tfall = 0.5 # fall rate, ms |
|
rss = driven / 1000.0 # spikes/millisecond |
|
rr = 3 * rss # transient driven rate |
|
rst = rss |
|
taur = 3 # rapid decay, msec |
|
taust = 10 |
|
ton = delay # msec |
|
stim_end = ton + dur |
|
trace_end = stim_end + post |
|
tvec = np.arange(0.0, trace_end, dt) # dt is in msec, so tvec is in milliseconds |
|
ta = np.arange(0.0, 20.0, dt) |
|
aw = ta * alpha * np.exp(-ta / alpha) / alpha # alpha waveform time course |
|
spt = [[]] * N # list of spike times |
|
wave = [[]] * N # waveform |
|
for j in range(0, N): # for each |
|
done = False |
|
sptime = [] |
|
qe = 0 |
|
nsp = 0 |
|
i = int(0) |
|
if spont <= 0: |
|
q = 1e6 |
|
else: |
|
q = 1000.0 / spont # q is in msec (spont in spikes/second) |
|
t = 0.0 |
|
while not done: |
|
a = np.random.sample(1) |
|
if t < ton: |
|
if spont <= 0: |
|
t = ton |
|
continue |
|
ti = -(np.log(a) / (spont / 1000.0)) # convert to msec |
|
if ti < deadtime: # delete intervals less than deadtime |
|
continue |
|
t = t + ti |
|
if t > ton: # if the interval would step us to the stimulus onset |
|
t = ton # then set to to the stimulus onset |
|
continue |
|
if t >= ton and t < stim_end: |
|
if t > ton: |
|
rise = 1.0 - np.exp(-(t - ton) / trise) |
|
else: |
|
rise = 1.0 |
|
ra = rr * np.exp(-(t - ton) / taur) |
|
rs = rst * np.exp(-(t - ton) / taust) |
|
q = rise * (ra + rs + rss) |
|
ti = -np.log(a) / (q + spont / 1000) # random.negexp(1000/q) |
|
if ti < deadtime: |
|
continue |
|
t = t + ti |
|
if t > stim_end: # only include interval if it falls inside window |
|
t = stim_end |
|
continue |
|
if t >= stim_end and t <= trace_end: |
|
if spont <= 0.0: |
|
t = trace_end |
|
continue |
|
if qe is 0: # have not calculated the new qe at end of stimulus |
|
rise = 1.0 - np.exp(-(stim_end - ton) / trise) |
|
ra = rr * np.exp(-(stim_end - ton) / taur) |
|
rs = rst * np.exp(-(stim_end - ton) / taust) |
|
qe = rise * (ra + rs + rss) # calculate the rate at the end |
|
fall = np.exp(-(t - stim_end) / tfall) |
|
q = qe * fall |
|
ti = -np.log(a) / ( |
|
q + spont / 1000.0 |
|
) # keeps rate from falling below spont rate |
|
if ti < deadtime: |
|
continue |
|
t = t + ti |
|
if t >= trace_end: |
|
done = True |
|
continue |
|
# now add the spike time to the list |
|
if nsp is 0: |
|
sptime = t |
|
nsp = nsp + 1 |
|
else: |
|
sptime = np.append(sptime, t) |
|
nsp = nsp + 1 |
|
# end of for loop on i |
|
if j is 0: |
|
wavej = np.zeros(len(tvec)) |
|
for i in range(0, len(sptime)): |
|
st = int(sptime[i] / dt) |
|
wavej[st] = wavej[st] + 1 |
|
spt[j] = sptime |
|
npts = len(tvec) |
|
if makewave: |
|
w = np.convolve(wavej, aw / max(aw)) * amp |
|
wave[j] = w[0:npts] |
|
return (spt, wave, tvec, N) |
|
|
|
|
|
def findspikes(t, v, thresh): |
|
""" findspikes identifies the times of action potential in the trace v, with the |
|
times in t. An action potential is simply timed at the first point that exceeds |
|
the threshold. |
|
""" |
|
tm = np.array(t) |
|
s0 = ( |
|
np.array(v) > thresh |
|
) # np.where(v > thresh) # np.array(v) > thresh # find points above threshold |
|
|
|
# print ('v: ', v) |
|
dsp = tm[s0] |
|
if dsp.shape[0] == 1: |
|
dsp = np.array(dsp) |
|
sd = np.append(True, np.diff(dsp) > 1.0) # find first points of spikes |
|
if len(dsp) > 0: |
|
sp = dsp[sd] |
|
else: |
|
sp = [] |
|
return sp # list of spike times. |
|
|
|
|
|
def measure(mode, x, y, x0, x1): |
|
""" return the mean and standard deviation of y in the window x0 to x1 |
|
""" |
|
xm = ma.masked_outside(x, x0, x1) |
|
ym = ma.array(y, mask=ma.getmask(xm)) |
|
if mode == "mean": |
|
r1 = ma.mean(ym) |
|
r2 = ma.std(ym) |
|
if mode == "max": |
|
r1 = ma.max(ym) |
|
r2 = 0 |
|
if mode == "min": |
|
r1 = ma.min(ym) |
|
r2 = 0 |
|
if mode == "median": |
|
r1 = ma.median(ym) |
|
r2 = 0 |
|
if mode == "p2p": # peak to peak |
|
r1 = ma.ptp(ym) |
|
r2 = 0 |
|
return (r1, r2) |
|
|
|
|
|
def mask(x, xm, x0, x1): |
|
xmask = ma.masked_outside(xm, x0, x1) |
|
xnew = ma.array(x, mask=ma.getmask(xmask)) |
|
return xnew.compressed() |
|
|
|
|
|
def vector_strength(spikes, freq): |
|
""" |
|
Calculate vector strength and related parameters from a spike train, for the specified frequency |
|
:param spikes: Spike train, in msec. |
|
:param freq: Stimulus frequency in Hz |
|
:return: a dictionary containing: |
|
|
|
r: vector strength |
|
n: number of spikes |
|
R: Rayleigh coefficient |
|
p: p value (is distribution not flat?) |
|
ph: the circularized spike train over period of the stimulus freq, freq, in radians |
|
d: the "dispersion" computed according to Ashida et al., 2010, etc. |
|
""" |
|
|
|
per = 1e3 / freq # convert from Hz to period in msec |
|
ph = 2 * np.pi * np.fmod(spikes, per) / (per) # convert to radians within a cycle |
|
c = np.sum(np.cos(ph)) ** 2 |
|
s = np.sum(np.sin(ph)) ** 2 |
|
vs = (1.0 / len(ph)) * np.sqrt(c + s) # standard vector strength computation |
|
n = len(spikes) |
|
R = n * vs # Raleigh coefficient |
|
Rp = np.exp(-n * vs * vs) # p value for n > 50 (see Ashida et al. 2010). |
|
d = np.sqrt(2.0 * (1 - vs)) / (2 * np.pi * freq) |
|
return {"r": vs, "n": n, "R": R, "p": Rp, "ph": ph, "d": d} |
|
|
|
|
|
def isi_cv2(splist, binwidth=1, t0=0, t1=300, tgrace=25): |
|
""" compute the cv and regularity according to Young et al., J. Neurophys, 60: 1, 1988. |
|
Analysis is limited to isi's starting at or after t0 but before t1, and ending completely |
|
before t1 + tgrace(to avoid end effects). t1 should correspond to the |
|
the end of the stimulus |
|
VERSION using dictionary for cvisi |
|
""" |
|
cvisit = np.arange(0, t1, binwidth) # build time bins |
|
cvisi = {} # isi is dictionary, since each bin may have different length |
|
for i in range(0, len(splist)): # for all the traces |
|
isit = splist[i] # get the spike times for this trial [1:-1] |
|
if len(isit) <= 1: # need 2 spikes to get an interval |
|
continue |
|
isib = np.floor(isit[0:-2] / binwidth) # discreetize |
|
isii = np.diff(splist[i]) # isis. |
|
for j in range(0, len(isib)): # loop over possible start time bins |
|
if ( |
|
isit[j] < t0 or isit[j] > t1 or isit[j + 1] > t1 + tgrace |
|
): # start time and interval in the window |
|
continue |
|
if isib[j] in cvisi: |
|
print("spike in bin: %d" % (isib[j])) |
|
cvisi[isib[j]] = np.append( |
|
cvisi[isib[j]], isii[j] |
|
) # and add the isi in that bin |
|
else: |
|
cvisi[isib[j]] = isii[j] # create it |
|
cvm = np.array([]) # set up numpy arrays for mean, std and time for cv analysis |
|
cvs = np.array([]) |
|
cvt = np.array([]) |
|
for i in cvisi.keys(): # for each entry (possible bin) |
|
c = [cvisi[i]] |
|
s = c.shape |
|
# print c |
|
if len(s) > 1 and s[1] >= 3: # require 3 spikes in a bin for statistics |
|
cvm = np.append(cvm, np.mean(c)) |
|
cvs = np.append(cvs, np.std(c)) |
|
cvt = np.append(cvt, i * binwidth) |
|
return (cvisit, cvisi, cvt, cvm, cvs) |
|
|
|
|
|
def isi_cv(splist, binwidth=1, t0=0, t1=300, tgrace=25): |
|
""" compute the cv and regularity according to Young et al., J. Neurophys, 60: 1, 1988. |
|
Analysis is limited to isi's starting at or after t0 but before t1, and ending completely |
|
before t1 + tgrace(to avoid end effects). t1 should correspond to the |
|
the end of the stimulus |
|
Version using a list of numpy arrays for cvisi |
|
""" |
|
cvisit = np.arange(0, t1, binwidth) # build time bins |
|
cvisi = [[]] * len(cvisit) |
|
for i in range(0, len(splist)): # for all the traces |
|
if len(splist[i]) < 2: # need at least 2 spikes |
|
continue |
|
isib = np.floor( |
|
splist[i][0:-2] / binwidth |
|
) # begining spike times for each interval |
|
isii = np.diff(splist[i]) # associated intervals |
|
for j in range(0, len(isib)): # loop over spikes |
|
if ( |
|
splist[i][j] < t0 or splist[i][j] > t1 or splist[i][j + 1] > t1 + tgrace |
|
): # start time and interval in the window |
|
continue |
|
cvisi[int(isib[j])] = np.append( |
|
cvisi[int(isib[j])], isii[j] |
|
) # and add the isi in that bin |
|
cvm = np.array([]) # set up numpy arrays for mean, std and time for cv analysis |
|
cvs = np.array([]) |
|
cvt = np.array([]) |
|
for i in range(0, len(cvisi)): # for each entry (possible bin) |
|
c = cvisi[i] |
|
if len(c) >= 3: # require 3 spikes in a bin for statistics |
|
cvm = np.append(cvm, np.mean(c)) |
|
cvs = np.append(cvs, np.std(c)) |
|
cvt = np.append(cvt, i * binwidth) |
|
return (cvisit, cvisi, cvt, cvm, cvs) |
|
|
|
|
|
if __name__ == "__main__": |
|
|
|
test = "isicv" |
|
|
|
if test == "isicv": |
|
""" this test is not perfect. Given an ISI, we calculate spike times |
|
by drawing from a normal distribution whose standard deviation varies |
|
with time, from 0 (regular) to 1 (irregular). As coded, the standard |
|
deviation never reaches the target value because spikes fall before or |
|
after previous spikes (thus reducing the stdev). Nonetheless, this shows |
|
that the CV calculation works correctly. """ |
|
nreps = 500 |
|
# cv will be 0 for first 50 msec, 0.5 for next 50 msec, and 1 for next 50 msec |
|
d = [[]] * nreps |
|
isi = 5.0 # mean isi |
|
# we create 100 msec of data where the CV goes from 0 to 1 |
|
maxt = 100.0 |
|
for i in range(nreps): |
|
for j in range(int(maxt / isi) + 1): |
|
t = float(j) * isi |
|
sd = float(j) / isi |
|
if sd == 0.0: |
|
d[i] = np.append(d[i], t) |
|
else: |
|
d[i] = np.append(d[i], np.random.normal(t, sd, 1)) |
|
for j in range(1, 10): # add more intervals at the end |
|
te = t + float(j) * isi |
|
d[i] = np.append(d[i], np.random.normal(te, sd, 1)) |
|
d[i] = np.sort(d[i]) |
|
# print d[i] |
|
# print diff(d[i]) |
|
sh = np.array([]) |
|
for i in range(len(d)): |
|
sh = np.append(sh, np.array(d[i])) |
|
(hist, bins) = np.histogram(sh, bins=250, range=(0, 250), new=True) |
|
if len(bins) > len(hist): |
|
bins = bins[0 : len(hist)] |
|
|
|
pl.figure(1) |
|
pl.subplot(2, 2, 1) |
|
pl.plot(bins, hist) |
|
|
|
(cvisit, cvisi, cvt, cvm, cvs) = isi_cv( |
|
d, binwidth=0.5, t0=0, t1=100, tgrace=25 |
|
) |
|
order = np.argsort(cvt) |
|
cvt = cvt[order] |
|
cvs = cvs[order] |
|
cvm = cvm[order] |
|
pl.subplot(2, 2, 2) |
|
pl.plot(cvt, cvm) |
|
pl.hold(True) |
|
pl.plot(cvt, cvs) |
|
pl.subplot(2, 2, 4) |
|
pl.plot(cvt, cvs / cvm) |
|
pl.show() |
|
|
|
if test == "measure": |
|
x = np.arange(0, 100, 0.1) |
|
s = np.shape(x) |
|
y = np.random.randn(s[0]) |
|
for i in range(0, 4): |
|
print("\ni is : %d" % (i)) |
|
x0 = i * 20 |
|
x1 = x0 + 20 |
|
(r0, r1) = measure("mean", x, y, x0, x1) |
|
print("mean: %f std: %f [0, 20]" % (r0, r1)) |
|
(r0, r1) = measure("max", x, y, x0, x1) |
|
print("max: %f std: %f [0, 20]" % (r0, r1)) |
|
(r0, r1) = measure("min", x, y, x0, x1) |
|
print("min: %f std: %f [0, 20]" % (r0, r1)) |
|
(r0, r1) = measure("median", x, y, x0, x1) |
|
print("median: %f std: %f [0, 20]" % (r0, r1)) |
|
(r0, r1) = measure("p2p", x, y, x0, x1) |
|
print("peak to peak: %f std: %f [0, 20]" % (r0, r1)) |
|
|
|
if test == "an_syn": |
|
(s, w, t, n) = an_syn(N=50, spont=50, driven=150, post=100, makewave=True) |
|
sh = np.array([]) |
|
for i in range(len(s)): |
|
sh = np.append(sh, np.array(s[i])) |
|
(hist, bins) = np.histogram(sh, bins=250, range=(0, 250), new=True) |
|
if len(bins) > len(hist): |
|
bins = bins[0 : len(hist)] |
|
|
|
import pylab as pl |
|
|
|
pl.figure(1) |
|
pl.subplot(2, 2, 1) |
|
pl.plot(bins, hist) |
|
|
|
pl.subplot(2, 2, 3) |
|
for i in range(len(w)): |
|
pl.plot(t, w[i]) |
|
pl.hold = True |
|
(cvisit, cvisi, cvt, cvm, cvs) = isi_cv(s) |
|
order = np.argsort(cvt) |
|
cvt = cvt[order] |
|
cvs = cvs[order] |
|
cvm = cvm[order] |
|
pl.subplot(2, 2, 2) |
|
pl.plot(cvt, cvs / cvm) |
|
pl.show() |
|
|
|
if test == "syns": |
|
(s, w, t, n) = syns(rate=20, delay=0, dur=100.0, N=5, makewave=True) |
|
sh = np.array([]) |
|
for i in range(len(s)): |
|
sh = np.append(sh, np.array(s[i])) |
|
(hist, bins) = np.histogram(sh, bins=250, range=(0, 250), new=True) |
|
if len(bins) > len(hist): |
|
bins = bins[0 : len(hist)] |
|
|
|
import pylab as pl |
|
|
|
pl.figure(1) |
|
pl.subplot(2, 2, 1) |
|
pl.plot(bins, hist) |
|
|
|
pl.subplot(2, 2, 3) |
|
pl.plot(t, w) |
|
pl.hold = True |
|
(cvisit, cvisi, cvt, cvm, cvs) = isi_cv(s) |
|
order = np.argsort(cvt) |
|
cvt = cvt[order] |
|
cvs = cvs[order] |
|
cvm = cvm[order] |
|
pl.subplot(2, 2, 2) |
|
pl.plot(cvt, cvs / cvm) |
|
pl.show()
|
|
|