# -*- coding: utf-8 -*-
"""
Created on Tue Nov 8 19:12:05 2011
@author: chris
"""
from __future__ import division
#import sys
#sys.path.insert(0, "/home/chris/lib/python") # for new matplotlib!!!
from pylab import *
from numpy import round, random, any
from units import *
import time
from NeuroTools import stgen
import h5py
import os
import shutil
def noclip(ax):
"Turn off all clipping in axes ax; call immediately before drawing"
ax.set_clip_on(False)
artists = []
artists.extend(ax.collections)
artists.extend(ax.patches)
artists.extend(ax.lines)
artists.extend(ax.texts)
artists.extend(ax.artists)
for a in artists:
a.set_clip_on(False)
def make_colormap(seq):
import matplotlib.colors as mcolors
"""Return a LinearSegmentedColormap
seq: a sequence of floats and RGB-tuples. The floats should be increasing
and in the interval (0,1).
"""
seq = [(None,) * 3, 0.0] + list(seq) + [1.0, (None,) * 3]
cdict = {'red': [], 'green': [], 'blue': []}
for i, item in enumerate(seq):
if isinstance(item, float):
r1, g1, b1 = seq[i - 1]
r2, g2, b2 = seq[i + 1]
cdict['red'].append([item, r1, r2])
cdict['green'].append([item, g1, g2])
cdict['blue'].append([item, b1, b2])
return mcolors.LinearSegmentedColormap('CustomMap', cdict)
def create_colnoise(t, sexp, cutf = None, seed = None, onf = None):
'''x = create_colnoise(t, sexp, cutf)
Make coloured noise signal
t = vector of times
sexp = spectral exponent - Power ~ 1 / f^sexp
cutf = frequency cutoff - Power flat (white) for f <~ cutf
cutf == None OR sexp == 0: white noise
output: mean = 0, std of signal = 1/2 => max 95% is 1
'''
from scipy import signal
nt = len(t)
dt = (t[-1] - t[0]) / (nt - 1)
random.seed(seed)
#if sexp == 0: # no smooth cutoff
# N = 10001 # number of filter taps
# xp = int((N-1)/2)
# x = random.standard_normal(size = len(t)+xp)
#else:
# x = random.standard_normal(size = shape(t))
x = random.standard_normal(size = shape(t))/2
#print std(x)
if cutf == None:
pass
else:
if sexp == 0: # no cutoff
pass
# OLD: sharp cutoff
#x = fft(x)
#f = fftfreq(nt, dt)
#x[nonzero(f == 0)] = 0 # remove zero frequency contribution
##i = nonzero(f != 0) # find indices i of non-zero frequencies
##x[i]=x[i] / (cutf ** 2 + f[i] ** 2) ** (sexp / 4) # using i allows cutf = 0
#x = real(ifft(x)) # ignore imaginary part (numerical error)
#x = x / std(x) # std of signal = 1/2, mean = 0
#fs = 1/dt # sampling rate
#fc = cutf/(0.5*fs) # cutoff frequency
#a = 1 # filter denominator
#b = signal.firwin(N, cutoff=fc, window='hamming') #'blackmanharris') # filter numerator
#x = signal.lfilter(b, a, x) # filtered output
#x = x[xp:]
##w,h = signal.freqz(b,a)
##h_dB = 20 * log10 (abs(h))
##figure(77)
##semilogx(w/max(w),h_dB)
##show()
elif sexp == -1: # sharp cutoff
x, freq, freq_wp, freq_used = create_multisines(t, freq_used=None, cutf=cutf, onf=onf)
#print std(x)
x = x / std(x) / 2 # std of signal = 1, mean = 0
else: # smooth cutoff
x = fft(x)
f = fftfreq(nt, dt)
x[nonzero(f == 0)] = 0 # remove zero frequency contribution
i = nonzero(f != 0) # find indices i of non-zero frequencies
x[i]=x[i] / (cutf ** 2 + f[i] ** 2) ** (sexp / 4) # using i allows cutf = 0
x = real(ifft(x)) # ignore imaginary part (numerical error)
x = x / std(x) / 2 # std of signal = 1/2, mean = 0
return x
def create_multisines(t, freq_used=array([1]), cutf = None, onf = None):
"""
This function will produce a colored noise signal using the time points
as defined in array t. The signal will be constructed using sinosoids with
frequencies as defined in the array freq_used and randomized phases.
The output consists of signal x
"""
tstop = t[-1]
df = 1 / tstop # frequency stepsize
dt = t[2] - t[1]
data_points = len(t) # number of time or frequency steps
vector_freq = zeros(data_points)
vector_phase = zeros(data_points)
if cutf != None:
f = arange(0,data_points)*df
if onf != None:
freq_used = f[nonzero((f <= cutf) & (f >= onf))]
#print freq_used
else:
freq_used = f[nonzero((f <= cutf) & (f > 0))]
index_f_used = round(freq_used / df).astype('int') # get indices of used frequencies in frequency vector
index_fneg_used = (data_points - index_f_used).astype('int') # indices of negative frequencies
index_fall_used = concatenate((index_f_used, index_fneg_used)) # indices of pos+neg frequencies
vector_freq[index_fall_used] = data_points / 2 # each frequency used ???
phase = 2*pi*(random.rand(len(freq_used),1)-0.5) # pick phases randomly shifted by +-pi (sould there by another 2* to shift +-2pi???)
vector_phase[index_f_used] = phase # assign positive phases to full vector
vector_phase[index_fneg_used] = -phase # assign negative phases to full vector
freqtemp = vector_freq * exp(1j * vector_phase) # generate frequency domain response
x = real(ifft(freqtemp)) # convert into time domain
#print "- Number of msine frequencies: " + str(2 * std(x) ** 2)
noise_data = x/max(abs(x)) # scale so that signal amplitude is 1
freq = fftfreq(data_points, dt)[ 0:round(data_points / 2) ] # positive frequency vector
noise_power = abs(fft(noise_data))[ 0:round(data_points / 2) ] # compute noise power
freq_wp = find(noise_power > 2 * std(noise_power)) # threshold to discriminate the indexes of peak frequencies
freq_used = freq[freq_wp] # vector of used frequencies [Hz]
return noise_data, freq, freq_wp, freq_used
def create_singlesine(fu = 5, amp = 0.5, ihold = 1, dt = 0.025*ms, periods = 10, minlength = 1*s, t_prestim = 2*s, l_zeros = 2):
"""
This function will produce a single sine signal of frequency fu with holding current ihold
Signal has at least the length periods*T (s) or minlength (s).
Use stimulate with pre stimulus of length t_prestim (s)
"""
fs = 1 / dt # sampling rate
tnext = 0
# delay for no noise input
start_zeros = zeros(l_zeros * fs)
t_zeros = tnext + arange(0, l_zeros, dt)
tnext = t_zeros[-1] + dt
l_pre_signal = ceil(t_prestim / (1. / fu)) * 1. / fu # length of pre stimulus should to be at least t_prestim seconds but with length of full periods
t_pre_signal = arange(0, l_pre_signal, dt) # create pre time vector
pre_signal = amp * sin(2 * pi * t_pre_signal * fu) # create pre signal vector
t_pre_signal = t_pre_signal + tnext
tnext = t_pre_signal[-1] + dt
l_t = max(minlength, periods * 1 / fu) # length of input_signal: stimulate for at least periods*T or minlength
t_input_signal = arange(0, l_t, dt) # create stimulus time vector
#window = sin(2 * pi * t_input_signal * 1/l_t/2) # not really good if nonlinear membrane function!!!!
input_signal = amp * sin(2 * pi * t_input_signal * fu)
t_input_signal = t_input_signal + tnext
i_start = len(start_zeros) + len(pre_signal)
i_stop = len(start_zeros) + len(pre_signal) + len(input_signal)
tnext = t_input_signal[-1] + dt
l_post_signal = 1 # length of post stimulus should only be 1 s, equivalent to 1 Hz lower bound for spike cutoff
t_post_signal = arange(0, l_post_signal, dt) # create pre time vector
post_signal = amp * sin(2 * pi * t_post_signal * fu) # create pre signal vector
t_post_signal = t_post_signal + tnext
all_data = concatenate((start_zeros, pre_signal, input_signal, post_signal)) # combine all
t = concatenate((t_zeros, t_pre_signal, t_input_signal, t_post_signal)) # combine all
t1 = arange(0, size(all_data) * dt, dt) # time vector of stimulus [s]
i_startstop = array([i_start, i_stop])
t_startstop = array([t[i_start], t[i_stop]])
iholdvec = concatenate((zeros(1 * fs), ones(len(all_data) - 1 * fs) * ihold)) # create holding current vector
#iholdvec = concatenate((zeros(1 * fs), ones(len(all_data) - 2 * fs) * ihold, zeros(1 * fs))) # create holding current vector
stimulus = all_data + iholdvec # create full stimulus vector
return t, stimulus, i_startstop, t_startstop
def create_ZAP(dt=0.025*ms, ihold=1*nA, amp=0.1*nA, fmax=100, t_stim=30*s, ex=2):
t = arange(0,t_stim,dt)
zap = amp*sin(2*pi*((fmax*t**ex)/((ex+1)*t_stim**ex))*t) + ihold
f=(fmax*t**ex)/(t_stim**ex)
return t, zap, f
def aiftransfer(freq = arange(0,1000,1), tau = 20*ms, f0 = 100*Hz, i0 = 0*nA, rm = 10*MOhm, Vreset = 0, Vth = 1, Vrest = 0, delta_t = 0):
"""
Create theoretical transfer function of a simplified integrate and fire neuron
from Knight 1972
"""
gamma = 1/tau
# theoretical transfer function H(2*pi*i*freq)
z = (gamma + 2 * pi * freq * 1j) / f0
if i0 == 0:
s0 = gamma / (1 - exp(-gamma / f0)) # theoretical current for f0
H = f0/s0 * exp(gamma / f0) * (1 - exp(-z)) / z * exp(- 2*pi*delta_t*1j)
else:
s0 = (i0*rm + (Vrest - Vreset)) / (Vth - Vreset)
H = rm/(Vth - Vreset) * f0/s0 * exp(gamma / f0) * (1 - exp(-z)) / z * exp(- 2*pi*delta_t*1j)
# useful to return TF for zero gamma
z = (2 * pi * freq * 1j) / f0
H0 = f0/s0 * (1 - exp(-z)) / z
return H, H0
def syn_kernel(t, tau1 = 5, tau2 = 10):
if tau1 == 0:
G = exp(-t/tau2)
else:
if (tau1/tau2 > .9999):
tau1 = .9999*tau2
tp = (tau1*tau2)/(tau2 - tau1) * log(tau2/tau1)
factor = -exp(-tp/tau1) + exp(-tp/tau2)
factor = 1/factor
G = factor * (exp(-t/tau2) - exp(-t/tau1))
return G
#def pop_giftransfer(freq = arange(0,1000,1), C = 0.5*1e-9, g = 0.025*1e-6, gr = 0.025*1e-6, tau_r = 100*1e-3, Vreset = -56*1e-3, Vth = -50*1e-3):
# """
# Create theoretical population transfer function of a general integrate and fire neuron with resonance
# from Magnus et al. 2003
# """
#
# omega = 2*pi*freq
#
# gamma = gr/g
# tau_m = C/g
#
# def Y(omega):
# (1 - g * tau_m * (Vth - Vreset) * A_LIF(omega)) / (1 + 1j * omega * tau_m)
#
# def A_LIF(omega):
#
# def A(omega):
# return (1 - gamma * Y(omega) / (1 + gamma * Y(omega) + 1j * omega * tau_r)) * A_LIF(omega)
#
# return A(omega)
def fit_aiftransfer(freq, H0, f0, i0):
"""
Fit theoretical transfer function of a simplified integrate and fire neuron to actual transfer function
freq = used frequencies
H = transfer function
"""
from scipy.optimize import fmin, leastsq
def peval(freq, f0, i0, p):
H, _ = aiftransfer(freq, tau = p[0], f0 = f0, i0 = i0, delta_t = p[2])
return p[1] * H
def residuals(p, H0, freq, f0, i0):
H = peval(freq, f0, i0, p)
#err = sum( abs(abs(H0) - abs(H)) + abs(unwrap(angle(H0)) - unwrap(angle(H0))) * 180/pi )**2 #
err = sum( ( concatenate((real(H0), imag(H0))) - concatenate((real(H), imag(H))) )**2 )
#err = sum( (abs(H0) - abs(H))**2 )
return err
p0 = array([20*ms, 1, 0]) # tau, rm
plsq = fmin(residuals, p0, args=(H0, freq, f0, i0))
#print plsq
tau = plsq[0] #[0] #plsq[0][0]
scale = plsq[1] #[1]
delta_t = plsq[2]
H, _ = aiftransfer(freq, tau = tau, f0 = f0, i0 = i0, delta_t = delta_t)
H = scale*H
return tau, scale, H, delta_t
def fit_sinusoid(f, t0, u0, t = None):
"""
Fit sinusoid of known frequency to signal
f = known frequency
t0 = timebase
u0 = signal
t = optional interpolation timebase
"""
e = ones(size(t0))
c = sin(2 * pi * f * t0)
s = cos(2 * pi * f * t0)
a = array((e, c, s)).T
b = linalg.lstsq(a, u0.T)[0]
umean = b[0]
uamp = sqrt(b[1] ** 2 + b[2] ** 2)
uphase = math.atan2(b[2], b[1]) # convention: cos(2*pi*f*t+uphase)
if t is not None:
u = umean + uamp * sin(2 * pi * f * t + uphase)
else:
u = umean + uamp * sin(2 * pi * f * t0 + uphase)
t = t0
# from scipy import io
#if f == 1:
# figure(9)
# plot(t0, u0, 'b.', t, u, 'r')
# #ylim(25, 40)
# from datetime import datetime
# idate = datetime.now().strftime('%Y%m%d_%H%M%S') # %S
# savefig("./figs/" + idate + "-sin_" + str(f) + ".png", dpi = 300) # save it
# txt=array((t0,u0,u)).T
# savetxt("./figs/" + idate + "-sin_" + str(f) + ".txt", txt)
# show()
# clf()
return umean, uamp, uphase, u
def fit_sinusoid2(f, t0, u0, t = None):
"""
Fit sinusoid of known frequency to signal
f = known frequency
t0 = timebase
u0 = signal
t = optional interpolation timebase
"""
from scipy.optimize import leastsq
def peval(t, f, p):
return p[0] + p[1] * sin(2 * pi * f * t + p[2])
def residuals(p, u, t, f):
err = u-peval(t, f, p)
return err
p0 = array([0, 0, 0])
plsq = leastsq(residuals, p0, args=(u0, t0, f))
umean = plsq[0][0]
uamp = plsq[0][1]
uphase = plsq[0][2]
if t is not None:
u = umean + uamp * sin(2 * pi * f * t + uphase)
else:
u = umean + uamp * sin(2 * pi * f * t0 + uphase)
t = t0
return umean, uamp, uphase, u
def fit_sinusoid_fft(f, t0, u0, t = None, i = None):
"""
Fit sinusoid of known frequency to signal using fft
f = known frequency
t0 = timebase
u0 = signal
t = optional interpolation timebase
i = optional stimulus current (will return magnitude)
"""
e = ones(size(t0))
c = sin(2 * pi * f * t0)
s = cos(2 * pi * f * t0)
a = array((e, c, s)).T
b = linalg.lstsq(a, u0.T)[0]
umean = b[0]
umean = mean(u0)
if i is not None:
b = linalg.lstsq(a, i.T)[0]
imean = b[0]
expb = exp(-2 * pi * 1j * f * t0) # exponetial basis for Fourier series
G = sum((u0 - umean) * expb) / sum((i - imean) * expb)
uamp = abs(G)
uphase = angle(G)
u = []
else:
expb = exp(-2 * pi * 1j * f * t0 / len(t0)) # exponetial basis for Fourier series
G = sum((u0 - umean) * expb)
uamp = abs(G)
uphase = angle(G)-pi/2
if t is not None:
u = umean + uamp * sin(2 * pi * f * t + uphase)
else:
u = umean + uamp * sin(2 * pi * f * t0 + uphase)
t = t0
return umean, uamp, uphase, u
def fit_exp(t0, u0, p0 = array([0, 0]), delay = 0):
"""
Fit exponential to signal
t0 = timebase
u0 = signal
"""
from scipy.optimize import leastsq
def peval(t, p, delay = 0):
dt = t[1]-t[0]
u = p[0] * (1 - exp(-(t-delay)/p[1]))
u[0:delay/dt]=0
return u
def residuals(p, u, t, delay = 0):
err = u-peval(t, p, delay)
return err
plsq = leastsq(residuals, p0, args=(u0, t0, delay))
udecay = plsq[0][0]
utau = plsq[0][1]
u = peval(t0, array([udecay,utau]), delay)
return udecay, utau, u
def fit_dualexp(t0, u0, p0 = array([0, 0, 0]), delay = 0):
"""
Fit exponential to signal
t0 = timebase
u0 = signal
"""
from scipy.optimize import leastsq
def peval(t, p, delay = 0):
dt = t[1]-t[0]
tp = (p[1]*p[2])/(p[2] - p[1]) * log(p[2]/p[1])
factor = 1 / (-exp(-tp/p[1]) + exp(-tp/p[2]))
u = factor * p[0] * ( -exp(-(t-delay)/p[1]) + exp(-(t-delay)/p[2]) )
u[0:delay/dt]=0
return u
def residuals(p, u, t, delay = 0):
err = u-peval(t, p, delay)
return err
plsq = leastsq(residuals, p0, args=(u0, t0, delay))
umax = plsq[0][0]
utau1 = plsq[0][1]
utau2 = plsq[0][2]
u = peval(t0, array([umax,utau1,utau2]), delay)
return umax, utau1, utau2, u
def fit_twodualexp(t0, u0, p0 = array([0, 0, 0, 0, 0, 0]), delay = 0):
"""
Fit exponential to signal
t0 = timebase
u0 = signal
"""
from scipy.optimize import leastsq
def peval(t, p, delay = 0):
dt = t[1]-t[0]
u = p[0] * ( -exp(-(t-delay)/p[1]) + exp(-(t-delay)/p[2]) )
v = p[3] * ( -exp(-(t-delay)/p[4]) + exp(-(t-delay)/p[5]) )
tp = (p[1]*p[2])/(p[2] - p[1]) * log(p[2]/p[1])
factor = 1 / (-exp(-tp/p[1]) + exp(-tp/p[2]))
u = factor * p[0] * ( -exp(-(t-delay)/p[1]) + exp(-(t-delay)/p[2]) )
tp = (p[4]*p[5])/(p[5] - p[4]) * log(p[5]/p[4])
factor = 1 / (-exp(-tp/p[4]) + exp(-tp/p[5]))
v = p[3] * ( -exp(-(t-delay)/p[4]) + exp(-(t-delay)/p[5]) )
u = u+v
u[0:delay/dt]=0
return u
def residuals(p, u, t, delay = 0):
err = u-peval(t, p, delay)
return err
plsq = leastsq(residuals, p0, args=(u0, t0, delay))
umax = plsq[0][0]
utau1 = plsq[0][1]
utau2 = plsq[0][2]
vmax = plsq[0][3]
vtau1 = plsq[0][4]
vtau2 = plsq[0][5]
u = peval(t0, array([umax,utau1,utau2,vmax,vtau1,vtau2]), delay)
return umax, utau1, utau2, u, vmax,vtau1,vtau2
def fit_threedualexp(t0, u0, p0 = array([0, 0, 0, 0, 0, 0, 0, 0, 0]), delay = 0):
"""
Fit exponential to signal
t0 = timebase
u0 = signal
"""
from scipy.optimize import leastsq
def peval(t, p, delay = 0):
dt = t[1]-t[0]
tp = (p[1]*p[2])/(p[2] - p[1]) * log(p[2]/p[1])
factor = 1 / (-exp(-tp/p[1]) + exp(-tp/p[2]))
u = factor * p[0] * ( -exp(-(t-delay)/p[1]) + exp(-(t-delay)/p[2]) )
tp = (p[4]*p[5])/(p[5] - p[4]) * log(p[5]/p[4])
factor = 1 / (-exp(-tp/p[4]) + exp(-tp/p[5]))
v = p[3] * ( -exp(-(t-delay)/p[4]) + exp(-(t-delay)/p[5]) )
tp = (p[7]*p[8])/(p[8] - p[7]) * log(p[8]/p[7])
factor = 1 / (-exp(-tp/p[7]) + exp(-tp/p[8]))
z = p[6] * ( -exp(-(t-delay)/p[7]) + exp(-(t-delay)/p[8]) )
u = u+v+z
u[0:delay/dt]=0
return u
def residuals(p, u, t, delay = 0):
err = u-peval(t, p, delay)
return err
plsq = leastsq(residuals, p0, args=(u0, t0, delay))
umax = plsq[0][0]
utau1 = plsq[0][1]
utau2 = plsq[0][2]
vmax = plsq[0][3]
vtau1 = plsq[0][4]
vtau2 = plsq[0][5]
zmax = plsq[0][6]
ztau1 = plsq[0][7]
ztau2 = plsq[0][8]
u = peval(t0, array([umax,utau1,utau2,vmax,vtau1,vtau2,zmax,ztau1,ztau2]), delay)
return umax, utau1, utau2, u, vmax, vtau1, vtau2, zmax, ztau1, ztau2
def fit_tripleexp(t0, u0, p0 = array([0, 0, 0, 0]), delay = 0):
"""
Fit exponential to signal
t0 = timebase
u0 = signal
"""
from scipy.optimize import leastsq
def peval(t, p, delay = 0):
dt = t[1]-t[0]
u = p[0] * ( -exp(-(t-delay)/p[1]) + exp(-(t-delay)/p[2]) + exp(-(t-delay)/p[3]) )
u[0:delay/dt]=0
return u
def residuals(p, u, t, delay = 0):
err = u-peval(t, p, delay)
return err
plsq = leastsq(residuals, p0, args=(u0, t0, delay))
umax = plsq[0][0]
utau1 = plsq[0][1]
utau2 = plsq[0][2]
utau3 = plsq[0][3]
u = peval(t0, array([umax,utau1,utau2,utau3]), delay)
return umax, utau1, utau2, utau3, u
def shannon_interp(t0, y0, t, tu=None, chunk_size=0):
"""
function y = shannon_interp(t0, y0, t, tu)
band-passed interpolation of non-uniform sampled data to
done by brute-force n x n matrix inversion, so only viable for
small sample numbers n.
t0 = non-uniform sample times
with t0 = [] uses linspaced data with same range as tu
y0 = sample data
tu = uniform sampling times (same length as t0)
t = timebase for interpolation
y = values at interpolated times t
"""
#print size(y0)
y = []
tn = []
it_start = 0
if chunk_size == 0:
chunk_size = len(t0)
for i in xrange(0, len(t0), chunk_size): # split with given chunk_size
t0_chunk = t0[i:i + chunk_size]
y0_chunk = y0[i:i + chunk_size]
if i == int((ceil((len(t0) / double(chunk_size)) - 1) * chunk_size)): # if this the last run
it_end = len(t) # use full rest of t
else:
it_end = t.searchsorted(t0_chunk[-1]) + 1
## test searchsorted()
# t = array([1,2,3,3.9,4.1,5,6])
# t.searchsorted(4)
# print str(it_start) + " " + str(it_end)
t_offset = t[it_start]
t0_chunk = t0_chunk - t_offset
t_chunk = t[it_start:it_end+1] - t_offset
it_start = it_end + 1
nsamp = len(t0_chunk)
if tu == None:
tu = linspace(t0_chunk[0], t0_chunk[-1], nsamp)
T = tu[2]-tu[1]
G = zeros((nsamp, nsamp))
for l in range(nsamp):
for k in range(nsamp):
G[l, k] = sinc((t0_chunk[l] - tu[k]) / T)
yu = linalg.solve(G,y0_chunk)
y_chunk = 0
for l in range(nsamp):
y_chunk = y_chunk + yu[l] * sinc((t_chunk - tu[l]) / T)
y.append(y_chunk)
tn.append(t_chunk)
y = concatenate(y)
tn = concatenate(tn)
#plot(t0, y0, "*b", t, y, "-r")
#axis(ymin=min(y0), ymax=max(y0))
#show()
## check missed values
#print set(t).difference(set(tn))
return y
def get_spikes(v, vthres, t):
'''
s, ts = get_spikes(v, vthres, t)
Spike detection implemented as voltage threshold crossings
v = membrane potential
vthres = threshold
t = time vector (in s)
-
s = spike 'bit' vector s
ts = (interpolated) spike times
'''
s = zeros(len(v))
s[1:] = ((v[:-1] <= vthres) & (v[1:] > vthres)).astype(float)
i = nonzero(s)[0]
ts = t[i-1]+(t[i]-t[i-1])*(vthres-v[i-1])/(v[i]-v[i-1])
return s, ts
def get_spikefreq(spike_times, stimlength = 0, compute_mean = 0, change_factor = 0.3):
'''
Compute spike frequency, set to second spike for causality,
also returns mean and onset frequency
Note: time in s now!
'''
if len(np.array(spike_times)) >= 2: # compute only if at least 2 spikes
spike_diff = diff(spike_times) # differentiate for spike intervals!
spike_freq = 1 / spike_diff # spike_times is in s
freq_times = spike_times[1:] # Set frequency to second spike for causality!
freq_onset = spike_freq[0] # Onset frequency only
a = [0]
if compute_mean > 1:
a = where(freq_times > (stimlength - 1))
print a
elif change_factor == 0:
freq_mean = mean(spike_freq[int(np.ceil(len(spike_freq)*9/10.)):])
elif any(spike_times > (stimlength / 2)) & compute_mean: # compute only if there is sustained spiking and a request is actually made
freq_change = abs(diff(spike_freq)) # check how the firing rate is changing
a = where(freq_change < change_factor) # check if rate is adapting
if any(a): # only compute mean firing rate if spiking adapts
freq_stable = a[0][0] # get first tuple and first element. here the firing rate has adapted
freq_mean = mean(spike_freq[freq_stable:]) # compute mean firing rate
else:
freq_mean = float('NaN')
else:
freq_times = np.array([])
spike_freq = np.array([])
freq_mean = float('NaN')
freq_onset = float('NaN')
return freq_times, spike_freq, freq_mean, freq_onset
def get_magphase(input_signal, t, output_signal, t_out, method = "peak", f = None, deb = 0):
"""
Compute the mean amplitude, mean magnitude and mean phase from an input and
output signal with a single frequency.
Default method *peak": peak detection
Magnitude is computed by: full oscillation out / full oscillation in
Minimum length of signal and response has to be 1.5 * 1 / f
Other method "fit": least-square fit, f has to be given!
"""
if method == "peak": # t and t_out have to be identical!
if size(t) != size(t_out):
raise ValueError('Input and output need same timebase! Use method "fit" instead!')
dt = t[2] - t[1]
input_signal_temp = input_signal - mean(input_signal) # remove holding current for amplitude detection
pos_halfwaves = concatenate((where(input_signal_temp > 0.75 * max(input_signal_temp))[0], array([len(input_signal_temp)]))) # detect the positive half waves
iend = pos_halfwaves[where(diff(pos_halfwaves) > 2)[0]] # detect where the half waves end, limit of 2 is arbitrary, 1 would also work
f = 1 / (mean(diff(iend)) * dt) # compute frequency
peak_vec_in = zeros(len(iend)-1) # construct vectors for peak values
peak_vec_out = zeros(len(iend)-1) # construct vectors for peak values
peak_vec_in_index = zeros(len(iend)-1, int) # construct vectors for peak indexes
peak_vec_out_index = zeros(len(iend)-1, int) # construct vectors for peak indexes
amptwice_vec_in = zeros(len(iend)-1) # construct vectors for size of full oscillation
amptwice_vec_out = zeros(len(iend)-1) # construct vectors for size of full oscillation
mean_vec_out = zeros(len(iend)-1) # construct vectors for size of full oscillation
for i in range(len(iend)-1): # iterate through all sections
wave_cut_in = input_signal[iend[i]:iend[i+1]] # cut out input section
peak_vec_in[i], i_in, min_in = wave_cut_in.max(0), wave_cut_in.argmax(0), wave_cut_in.min(0) # get peak positions, indexes and minimas
amptwice_vec_in[i] = peak_vec_in[i] - min_in
peak_vec_in_index[i] = int(iend[i] + i_in) # convert to global indexes
wave_cut_out = output_signal[iend[i]:iend[i+1]] # cut out output section
peak_vec_out[i], i_out, min_out = wave_cut_out.max(0), wave_cut_out.argmax(0), wave_cut_out.min(0) # get peak positions, indexes and minimas
amptwice_vec_out[i] = peak_vec_out[i] - min_out
peak_vec_out_index[i] = int(iend[i] + i_out) # convert to global indexes
mean_vec_out[i] = amptwice_vec_out[i]*0.5 + min_out
phase_mean = mean((peak_vec_in_index-peak_vec_out_index) * dt / (1 / f) * 360) # compute phase
mag_mean = mean(amptwice_vec_out / amptwice_vec_in) # compute magnitude
amp_mean = mean(peak_vec_out) # just compute mean of amplitude
umean = mean(mean_vec_out[i])
if method == "fit":
print "- fitting"
imean, iamp, iphase, i = fit_sinusoid(f, t, input_signal)
umean, uamp, uphase, u = fit_sinusoid(f, t_out, output_signal, t)
phase_mean = (uphase - iphase) * (180 / pi)
mag_mean = uamp / iamp
amp_mean = uamp + umean # for voltage amplitude DAngelo 2001
if method == "fft":
if any(abs(t_out - t) > 0.00000001): # bloody floating point precision
raise ValueError('Time base of input (t) and output (t_out) have to be identical')
imean, iamp, iphase, i = fit_sinusoid(f, t, input_signal)
umean, ramp, rphase, u = fit_sinusoid_fft(f, t_out, output_signal, i = input_signal)
phase_mean = rphase * (180 / pi)
mag_mean = ramp
amp_mean = ramp * iamp
print f, amp_mean
#plt.figure('FIT')
#ax99 = plt.subplot(1,1,1)
#ax99.plot(t_out, output_signal,'b')
#ax99.plot(t_out, u,'r')
#plt.savefig("./figs/fit_" + str(f) + ".pdf", dpi = 300, transparent=True) # save it
#plt.clf()
return amp_mean, mag_mean, phase_mean, umean, u
def construct_Stimulus(noise_data, fs, amp = 1, ihold = 0, tail_points = 2, delay_baseline = 4):
"""
Construct Stimulus from cnoise/msine input and other parameters.
"""
#inin = 8 # stimulate before with 10 s of signal
inin = np.array((len(noise_data)/fs)*0.1).clip(max=8)
stim_data = concatenate((noise_data[-inin*fs:], noise_data)) # increase length of stimulus # no normalization here: / max(abs(noise_data))
stimulus = concatenate((concatenate((zeros(round(delay_baseline*fs)), amp * stim_data)), zeros(round(tail_points*fs)))) # construct stimulus
iholdvec = concatenate((zeros(round(fs)), ones(round(len(stimulus) - 1 * fs)) * ihold))
stimulus = stimulus + iholdvec
dt = 1 / fs
t = arange(0, len(stimulus) * dt,dt) # time vector of stimulus [s]
t_startstop = np.array([inin+delay_baseline, inin+delay_baseline+len(noise_data)/fs])
return stimulus, t, t_startstop
def construct_Pulsestim(dt = 0.025e-3, pulses = 1, latency = 10e-3, stim_start = 0.02, stim_end = 0.02, len_pulse = 0.5e-3, amp_init = 1, amp_next = None):
"""
Construct a pulse stimulus in the form of |---stim_start---|-len_pulse-|--(latency-len_pulse)--|...|-len_pulse-|--(latency-len_pulse)--|--stim_end--|
For stim_end shorter than pulse: stim_end = stim_end + len_pulse
"""
#print dt
fs = 1 / dt
if len(np.shape(amp_next)) > 0:
pulses = len(amp_next)
amp_vec = amp_next
else:
if amp_next == None:
amp_vec = np.ones(pulses)*amp_init
else:
amp_vec = np.ones(pulses)*amp_next
amp_vec[0] = amp_init
if len(np.shape(latency)) > 0:
pass
else:
latency = np.ones(pulses)*latency
if len(amp_vec) != len(latency):
raise ValueError('amp_vec and latency vectors do not have the same size!!!')
if stim_end < len_pulse:
print "From construct_Pulsestim: stim_end shorter than pulse, setting stim_end = stim_end + len_pulse"
stim_end = stim_end + len_pulse
ivec = zeros(round((stim_start + sum(latency) + stim_end)*fs)) # construct zero vector to begin with
t = arange(0, len(ivec))*dt # construct time vector
ivec[round(stim_start*fs):round((stim_start+len_pulse)*fs)] = amp_vec[0]
for i in range(1, pulses):
ivec[round((stim_start+sum(latency[0:i]))*fs):round((stim_start+sum(latency[0:i])+len_pulse)*fs)] = amp_vec[i]
#ivec_new = concatenate((zeros(round((stim_start + (i - 1) * latency[i-1]) * fs)), ones(round(len_pulse * fs)) * amp)) # stimulus delay + spike
#
#ivec_new = concatenate((ivec_new, zeros(round(((pulses - i + 1) * latency[i-1] - len_pulse + stim_end) * fs)))) # rest of stimulus
#
#if len(ivec) > len(ivec_new):
# ivec_new = concatenate(( ivec_new, zeros(len(ivec)-len(ivec_new)) ))
#
#if len(ivec) < len(ivec_new):
# rem = -1*(len(ivec_new)-len(ivec))
# ivec_new = ivec_new[:rem]
#ivec = ivec + ivec_new
#if len(t) != len(ivec):
#raise ValueError('Both vectors do not have the same size!!!')
return t, ivec
def compute_Impedance(voltage, current, t1, stimulus, t, noise_data_points, freq_wp = None, do_csd = 0, w_length = 5):
"""
Compute the impedance using input and output and cross spectral density method,
cut out relevant parts with knowledge of noise_data_points
w_length: desired length of csd window
"""
from matplotlib.mlab import csd
response = voltage
#stimulus_out = current
# no interpolation is best !!
#stimulus_out_ds = interp(t,t1[:-1],stimulus_out[:-1]) # interpolate (downsample) to be eqivalent with input
response_ds = interp(t,t1[:-1],response[:-1]) # interpolate (downsample) to be eqivalent with input
#stimulus_out_ds = stimulus_out[:end-1] # no interpolation
#response_ds = response[:end-1]
#noiseinput = stimulus_out_ds[:end-1] # interplation introduces error!!!
noiseinput = stimulus # interplation introduces error: use original input
dt = t[2] - t[1]
fdataend = max(find(noiseinput - noiseinput[round(1.5 / dt)])) # get rid of constant input for detection
noiseinput_cut = noiseinput[(fdataend - noise_data_points):fdataend] # cut out relevant part
response_ds_cut = response_ds[(fdataend - noise_data_points):fdataend] # cut out relevant part
if do_csd:
# power spectrum & transfer function using cross spectral density
fs = 1 / dt
nfft = 2 ** int(ceil(log2(w_length * fs))) # nfft is the number of data points used
# in each block for the FFT. Must be even;
# a power 2 is most efficient.
#w_true_length = nfft * dt # true length window
c_in_out, freq = csd(noiseinput_cut, response_ds_cut, Fs = fs, NFFT = nfft, noverlap = nfft/4)
c_in_in, freq = csd(noiseinput_cut, noiseinput_cut, Fs = fs, NFFT = nfft, noverlap = nfft/4)
#c_in_out = abs(c_in_out) # phas einformation is lost by doing this
#c_in_in = abs(c_in_in)
zdata_pos = c_in_out / c_in_in
pwr = abs(zdata_pos) # Impedance Power
pha = angle(zdata_pos, deg = True) # Impedance Phase in degree
if freq_wp == None: # check if specific frequencies are used
fstart = find(freq >= 1. / w_length)[0] # lowest frequency depends on window length
fend = max(find(freq < 1000)) # only use frequencies smaller than 1kHz
magz = pwr[fstart:fend] # don't use freq = 0
phaz = pha[fstart:fend] # don't use freq = 0
freq_used = freq[fstart:fend] # don't use freq = 0
#print "- number of cnoise frequencies: " + str(len(freq_used)) + " Hz"
else: # multi sine stimulation
freq_was = fftfreq(noise_data_points, dt)[ 0:round(noise_data_points / 2) ] # full frequency vector without windowing
freq_used = freq_was[freq_wp] # the frequencies that were used
#freq_wp_new = freq.searchsorted(freq_used)
freq_wp_new = [abs(freq-i).argmin() for i in freq_used] # find corresponding indexes to sued frequencies
error = sum(abs(freq_used - freq[freq_wp_new]))
print "- frequency conversion error due to csd: " + str(error) + " Hz"
magz = pwr[freq_wp_new]
phaz = pha[freq_wp_new]
freq_used = freq[freq_wp_new] # freq_used is changed here to new values !
else:
noiseinput_f = fft(noiseinput_cut) # transform into frequency domain
response_f = fft(response_ds_cut) # transform into frequency domain
zdata = response_f / noiseinput_f # mV/nA = MOhms
zdata_pos = zdata[0:round(len(noiseinput_f) / 2)] # Take only the first half of the power to avoid redundancy
pwr = abs(zdata_pos) # Impedance Power
pha = angle(zdata_pos, deg=True) # Impedance Phase in degree
freq = fftfreq(noise_data_points, dt)[0:round(noise_data_points / 2)] # positive frequency vector
if freq_wp == None: # check if specific frequencies are used
fstart = find(freq >= 1. / (noise_data_points * dt))[0] # lowest frequency depends on window length
fend = max(find(freq < 1000)) # only use frequencies smaller than 1kHz
magz = pwr[fstart:fend] # don't use freq = 0
phaz = pha[fstart:fend] # don't use freq = 0
freq_used = freq[fstart:fend] # don't use freq = 0
#print "- number of cnoise frequencies: " + str(len(freq_used)) + " Hz"
else: # multi sine stimulation
magz = pwr[freq_wp]
phaz = pha[freq_wp]
freq_used = freq[freq_wp]
ca = 1 / (magz * exp(pi / 180 * 1j * phaz))
#figure; plot(t1, current); show()
#figure; semilogx(freq_used, magz); show()
return magz, phaz, ca, freq, freq_used
def compute_Transfer(spike_freq, freq_times, stimulus, t, noise_data_points, gsyn, gsyn_in = None, freq_wp = None, method_interpol = array(["linear", "quadratic", "shannon", "syn"]), do_csd = 0, nc_delay = 0, w_length = 2*s, t_kernel = 0, t_qual = 0, K_mat_old = np.array([]), t_startstop = [], give_psd = False):
"""
Compute the transfer function using cross spectral density method,
cut out relevant parts with knowledge of noise_data_points
w_length: desired length of csd window
t_kernel: length of final kernel, default = full length
t_qual: length of quality comparison, default: 0*s comparison
"""
from matplotlib.mlab import csd
from scipy.interpolate import InterpolatedUnivariateSpline
import neuronpy.util.spiketrain
import nitime.algorithms as tsa
#from csdmt import csdmt
dt = t[2] - t[1]
shift = 0 # only important for syn if nc_delay was used
fmean = None
for i, ii in enumerate(method_interpol):
# Here we have to interpolate!!
if "linear" in ii:
# linear interpolation
response_ds = interp(t, freq_times, spike_freq, left = 0, right = 0) # interpolate (upsample) to be eqivalent with input, set zero beginning and end
if "bin" in ii:
dt0 = dt
t0 = t
response_ds = spike_freq
t = freq_times
dt = t[2] - t[1]
if dt != dt0: #adjust stimulus
stimulus = interp(t,t0,stimulus)
print "Adjust Stimulus", len(response_ds), len(stimulus)
if "dslin" in ii:
# linear interpolation
response_ds = interp(t, freq_times, spike_freq, left = 0, right = 0) # interpolate (upsample) to be eqivalent with input, set zero beginning and end
# test downsampling
dt2 = 5e-3
dt_fac = dt/dt2
t2 = arange(0, t[-1], dt2)
response_ds = interp(t2,t,response_ds)
# downsample to bin size dt!
stimulus = interp(t2,t,stimulus)
t = t2
dt = dt2
noise_data_points = noise_data_points*dt_fac
#figure(44)
#plot(freq_times, ones(len(freq_times)), '*')
#plot(t, response_ds)
if "binary" in ii:
# binary spike code
dt2 = 1e-3
dt_fac = dt/dt2
t2 = arange(0, t[-1], dt2) # define edges, otherwise get_histogram does not define edges consistently!!!
[response_ds, _] = neuronpy.util.spiketrain.get_histogram(freq_times, bins = t2)
response_ds = concatenate((zeros(1),response_ds))
# downsample to bin size dt!
stimulus = interp(t2,t,stimulus)
t = t2
dt = dt2
noise_data_points = noise_data_points*dt_fac
if "dt" in ii:
# binary spike code using dt as bin size!
[response_ds, _] = neuronpy.util.spiketrain.get_histogram(freq_times, bins = t)
response_ds = concatenate((zeros(1),response_ds))
if "quadratic" in ii:
# quadratic interpolation
sfun = InterpolatedUnivariateSpline(freq_times, spike_freq, k=2)
response_ds = sfun(t)
if "shannon" in ii:
# shannon interpolation
response_ds = shannon_interp(freq_times, spike_freq, t)
if "syn" in ii:
# synaptic integration
response_ds = gsyn
shift = int(nc_delay/dt) # shift response by the nc delay to remove offset
if "gsyn_in" in ii:
# synaptic integration
response_ds = gsyn_in
shift = (2*ms)/dt # shift response by the delay of 1*ms spike play + 1*ms synapse nc to remove offset
noiseinput = stimulus # interplation introduces error: use original input
if len(t_startstop) < 2:
fdataend = max(find(noiseinput-noiseinput[round(1.5 / dt)]))+1 # get rid of constant input for detection
noiseinput_cut = noiseinput[(fdataend-noise_data_points):fdataend] # cut out relevant part
response_ds_cut = response_ds[(fdataend-noise_data_points)+shift:fdataend+shift] # cut out relevant part
else:
print "- t_startstop given!"
#i_start = mlab.find(t >= t_startstop[0])
#i_stop = mlab.find(t <= t_startstop[1])
#noiseinput_cut = noiseinput[i_start:i_stop] # cut out relevant part
#response_ds_cut = response_ds[i_start+shift:i_stop+shift] # cut out relevant part
noiseinput_cut = noiseinput[int(t_startstop[0]/dt):int(t_startstop[1]/dt)] # cut out relevant part
response_ds_cut = response_ds[int(t_startstop[0]/dt)+shift:int(t_startstop[1]/dt)+shift] # cut out relevant part
if "linear" in ii:
fmean = mean(response_ds_cut) # take fmean from linear interpolation
rmean = mean(response_ds_cut)
smean = mean(noiseinput_cut)
if do_csd: # Use CSD for colored/white noise AND multi sine estimation
print "- using csd for transfer function estimation! \r"
# power spectrum & transfer function using cross spectral density
fs = 1 / dt
nfft = 2 ** int(ceil(log2(w_length * fs))) # nfft is the number of data points used
# in each block for the FFT. Must be even;
# a power 2 is most efficient.
w_true_length = nfft * dt # true length window
r = response_ds_cut-rmean # response, no mean
s = noiseinput_cut-smean # signal, no mean
if "NOT_USED" in ii:
sr = np.array([s,r])
BW = None #0.05
#NW = BW / (2 * fs) * len(s)
#print NW
freq, P = tsa.multi_taper_csd(sr, Fs=fs, BW=BW, low_bias=True, adaptive=True, sides='twosided')
P_ss = P[0,0,:]
P_rr = P[1,1,:]
P_sr = P[0,1,:]
P_rs = P[1,0,:]
# results are returned like a = array([0, 1, 2, -3, -2, -1])
# already ordered for ifft!
# but freq is wrong!!
else:
noverlap = nfft / 4
#print np.shape(r)
P_rr, freq = csd(r, r, Fs = fs, NFFT = nfft, noverlap = noverlap, sides = 'twosided') # , window = mlab.window_none
P_sr, freq = csd(s, r, Fs = fs, NFFT = nfft, noverlap = noverlap, sides = 'twosided')
P_rs, freq = csd(r, s, Fs = fs, NFFT = nfft, noverlap = noverlap, sides = 'twosided')
#P_rs = concatenate([array([0]),P_sr[1:][::-1]]) # P_rs(f) = P_sr(-f)
#P_rs, freq = csd(r, s, Fs = fs, NFFT = nfft, noverlap = nfft / 4, sides = 'twosided')
P_ss, freq = csd(s, s, Fs = fs, NFFT = nfft, noverlap = noverlap, sides = 'twosided')
# results are returned like a = array([-3, -2, -1, 0, 1, 2])
# for ifft it should be: ifftshift(a) = array([ 0, 1, 2, -3, -2, -1])
# reorder here!!
P_rr = ifftshift(P_rr)
P_sr = ifftshift(P_sr)
P_rs = ifftshift(P_rs)
P_ss = ifftshift(P_ss)
freq = ifftshift(freq)
if sum(real(P_ss))<=0: # there is no power in the input, just check the output!
give_psd = True
P_ss = np.ones(len(P_ss))
P_sr = np.ones(len(P_sr))
P_rs = np.ones(len(P_rs))
print "- return PSD of output only!"
P_nn = (P_ss - abs(P_sr)**2/P_rr)
SNR = (P_ss / P_nn)[0:len(freq)/2] # only positive frequencies
VAFf = (abs(P_sr)**2/(P_rr*P_ss))[0:len(freq)/2]
#fend = max(find(freq[0:len(freq)/2] < 100)) # only frequencies up to 100 Hz
#print "Total VAF: " + str(sum(( abs(P_sr)**2 / P_ss )[0:len(freq)/2][0:fend]) / sum(P_rr[0:len(freq)/2][0:fend]))
if give_psd:
fdata_pos = (P_rr)[0:len(freq)/2] # only positive frequencies
else:
fdata_pos = (P_sr/P_ss)[0:len(freq)/2] # only positive frequencies
Kf0 = (P_rs/P_rr)
P_ss = P_ss[0:len(freq)/2]
freq = freq[0:len(freq)/2] # only positive frequencies
#figure(33)
#plot(freq, P_ss[0:len(P_ss)/2])
#figure(34)
#plot(freq, P_rr[0:len(P_rr)/2])
#figure(35)
#plot(freq, P_sr[0:len(P_sr)/2])
#figure(36)
#plot(freq, P_rs[0:len(P_rs)/2])
#figure(37)
#semilogx(freq, 10*log10(SNR))
#show()
pwr = abs(fdata_pos) # Transfer Power
pha = angle(fdata_pos, deg = True) # Transfer Phase in degree
if t_kernel == 0:
n_kernel = len(Kf0)
else:
n_kernel = int(t_kernel/dt)
tk = arange(n_kernel)*dt
if freq_wp == None: # colored/white noise: check if specific frequencies are used
#fstart = find(freq >= 1. / w_length)[0] # lowest frequency depends on window length
fstart = 1
fend = max(find(freq < 1000)) # only use frequencies smaller than 1kHz
#freq_used = freq[fstart:fend] # don't use freq = 0
freq_used = freq[fstart:fend]
if i == 0:
# create matrix to hold all different interpolation methods:
ca_mat = zeros((len(method_interpol),len(fdata_pos)), dtype=complex)
mag_mat = zeros((len(method_interpol),len(freq_used))) # impedance magnitude vector
pha_mat = zeros((len(method_interpol),len(freq_used))) # imedance phase vector
K_mat = zeros((len(method_interpol),n_kernel), dtype=complex)
SNR_mat = zeros((len(method_interpol),len(freq_used)))
VAFf_mat = zeros((len(method_interpol),len(freq_used)))
Qual_mat = zeros((len(method_interpol),3,3))
mag_mat[i,:] = pwr[fstart:fend].T # don't use freq = 0
pha_mat[i,:] = pha[fstart:fend].T # don't use freq = 0
# USE FULL FREQ FOR THIS
ca_mat[i,:] = fdata_pos.T #[fstart:fend].T
SNRi = 10*log10(real(SNR[fstart:fend].T))
SNR_mat[i,:] = SNRi
VAFi = real(VAFf[fstart:fend].T)
VAFf_mat[i,:] = VAFi
P_ss = P_ss[fstart:fend]
#K_mat[i,:] = ifftshift(ifft(Kf0,n_kernel))
#K_mat[i,:] = ifftshift(ifft(concatenate((Kf0,Kf0[::-1][:-1])),n_kernel))
K_mat[i,:] = ifftshift(ifft(Kf0))[len(Kf0)//2-n_kernel//2:len(Kf0)//2-n_kernel//2+n_kernel] # cutting is better than "downsampling"
#print "- number of cnoise frequencies: " + str(len(freq_used))
else: # multi sine stimulation
freq_was = fftfreq(noise_data_points, dt)[ 0:round(noise_data_points / 2) ] # full frequency vector without windowing
freq_used = freq_was[freq_wp] # the frequencies that were used
#freq_wp_new = freq.searchsorted(freq_used)
freq_wp_new = [abs(freq-i).argmin() for i in freq_used] # find corresponding indexes to sued frequencies
error = sum(abs(freq_used - freq[freq_wp_new]))
print "- frequency conversion error due to csd: " + str(error) + " Hz"
print freq_used
print freq[freq_wp_new]
freq_used = freq[freq_wp_new] # freq_used is changed here to new values !!
if i == 0:
# create matrix to hold all different interpolation methods:
ca_mat = zeros((len(method_interpol),len(freq_used_new)), dtype=complex)
mag_mat = zeros((len(method_interpol),len(freq_used_new))) # impedance magnitude vector
pha_mat = zeros((len(method_interpol),len(freq_used_new))) # imedance phase vector
K_mat = zeros((len(method_interpol),n_kernel), dtype=complex)
SNR_mat = zeros((len(method_interpol),len(freq_used_new)))
VAFf_mat = zeros((len(method_interpol),len(freq_used_new)))
Qual_mat = zeros((len(method_interpol),3,3))
mag_mat[i,:] = pwr[freq_wp_new]
pha_mat[i,:] = pha[freq_wp_new]
ca_mat[i,:] = fdata_pos[freq_wp_new]
K_mat[i,:] = ifftshift(ifft(Kf0,n_kernel))[freq_wp_new]
SNRi = 10*log10(real(SNR[freq_wp_new]))
SNR_mat[i,:] = SNRi
VAFi = real(VAFf[freq_wp_new])
VAFf_mat[i,:] = VAFi
P_ss = []
# get SNR and VAF quality
SNR0 = mean(SNRi[0:3])
iSNRc = find(SNRi <= (SNR0-3))
if len(iSNRc) == 0:
SNR_cutff = float('NaN')
else:
SNR_cutff = freq_used[iSNRc[0]];
SNR_mean = mean(SNRi)
Qual_mat[i,0,0] = SNR0
Qual_mat[i,0,1] = SNR_cutff
Qual_mat[i,0,2] = SNR_mean
#print 'SNR0:', SNR0
VAF0 = mean(VAFi[0:3])
iVAFc = find(VAFi <= 0.9*VAF0)
if len(iVAFc) == 0:
VAF_cutff = float('NaN')
else:
VAF_cutff = freq_used[iVAFc[0]];
VAF_mean = mean(VAFi)
Qual_mat[i,1,0] = VAF0
Qual_mat[i,1,1] = VAF_cutff
Qual_mat[i,1,2] = VAF_mean
if t_qual > 0: # compare signal to estimated signal
resp = r[0:t_qual/dt] # used response
stim = s[0:t_qual/dt] # used signal
t2 = t[0:t_qual/dt]
if len(K_mat_old) > 0:
K_mat = K_mat_old
tc, resp_cc, stim_cc, stim_re_cc, noise_cc, CF, VAF = reconstruct_Stimulus(K_mat[i,:], resp, stim, t2)
if i == 0:
# create matrix to hold all different interpolation methods:
stim_re_mat = zeros((len(method_interpol),len(stim_re_cc)))
resp_mat = zeros((len(method_interpol),len(resp_cc)))
noise_mat = zeros((len(method_interpol),len(noise_cc)))
CF_mat = zeros(len(method_interpol))
VAF_mat = zeros(len(method_interpol))
stim_re_mat[i,:] = stim_re_cc
stim = stim_cc
resp_mat[i,:] = resp_cc
noise_mat[i,:] = noise_cc
CF_mat[i] = CF
VAF_mat[i] = VAF
Qual_mat[i,2,0] = CF
Qual_mat[i,2,1] = VAF
else: # no comparison
stim_re_mat = []
stim = []
resp_mat = []
noise_mat = []
CF_mat = []
VAF_mat = []
tc = []
else: # NO CSD
noiseinput_f = fft(noiseinput_cut) # transform into frequency domain
response_f = fft(response_ds_cut) # transform into frequency domain
fdata = response_f / noiseinput_f # Hz / nA
fdata_pos = fdata[0:round(len(noiseinput_f) / 2)] # Take only the first half of the power to avoid redundancy
pwr = abs(fdata_pos) # Transfer Power
pha = angle(fdata_pos, deg = True) # Transfer Phase in degree
freq = fftfreq(noise_data_points, dt)[0:round(noise_data_points / 2)] # positive frequency vector
if freq_wp == None: # colored/white noise: check if specific frequencies are used
fstart = find(freq >= 1. / (noise_data_points*dt))[0] # lowest frequency depends on window length
fend = max(find(freq < 1000)) # only use frequencies smaller than 1kHz
freq_used = freq[fstart:fend] # don't use freq = 0
if i == 0:
# create matrix to hold all different interpolation methods:
ca_mat = zeros((len(method_interpol),len(freq_used)), dtype=complex)
mag_mat = zeros((len(method_interpol),len(freq_used))) # impedance magnitude vector
pha_mat = zeros((len(method_interpol),len(freq_used))) # imedance phase vector
K_mat = zeros((len(method_interpol)), dtype=complex)
tk = []
mag_mat[i,:] = pwr[fstart:fend] # don't use freq = 0
pha_mat[i,:] = pha[fstart:fend] # don't use freq = 0
ca_mat[i,:] = 1 / (mag_mat[i,:] * exp(pi / 180 * 1j * pha_mat[i,:]))
#print "- number of cnoise frequencies: " + str(len(freq_used)) + " Hz"
else: # multi sine stimulation
freq_used = freq[freq_wp]
if i == 0:
# create matrix to hold all different interpolation methods:
ca_mat = zeros((len(method_interpol),len(freq_used)), dtype=complex)
mag_mat = zeros((len(method_interpol),len(freq_used))) # impedance magnitude vector
pha_mat = zeros((len(method_interpol),len(freq_used))) # imedance phase vector
K_mat = zeros((len(method_interpol)), dtype=complex)
tk = []
mag_mat[i,:] = pwr[freq_wp]
pha_mat[i,:] = pha[freq_wp]
ca_mat[i,:] = 1 / (mag_mat[i,:] * exp(pi / 180 * 1j * pha_mat[i,:]))
stim = []
stim_re_mat = []
resp_mat = []
noise_mat = []
CF_mat = []
VAF_mat = []
VAFf_mat = []
SNR_mat = []
tc = []
P_ss = []
Qual_mat = []
SNR_mat = (freq_used, SNR_mat)
VAFf_mat = (freq_used, VAFf_mat)
return {'mag_mat':mag_mat,'pha_mat':pha_mat,'ca_mat':ca_mat,'freq':freq,'freq_used':freq_used,
'fmean':fmean,'tk':tk,'K_mat':K_mat,'SNR_mat':SNR_mat,'Qual_mat':Qual_mat,
'stim':stim,'stim_re_mat':stim_re_mat,'resp_mat':resp_mat,
'noise_mat':noise_mat,'CF_mat':CF_mat,'VAF_mat':VAF_mat, 'VAFf_mat':VAFf_mat, 'tc':tc, 'P_ss':P_ss }
def reconstruct_Stimulus(Ker, resp, stim, t):
dt = t[2] - t[1]
# if len(Ker)%2==0:
# Ker_conv = Ker[0:-1-1]
# else:
# Ker_conv = Ker
Ker_conv = Ker
n_kernel = len(Ker_conv) # NO SHIFT NECESSARY!!!
start_time = time.time()
stim_re_ = real(convolve(Ker_conv, resp)) # reconstructed signal
process_time = time.time() - start_time
print "Convolution took ", process_time
#n_shift = floor(n_kernel/2)
#n_shift = 1 # JUST SHIFT ONE TIMESTEP, IDKY???
#stim_re = concatenate([stim_re_[n_shift:], zeros(n_shift)]) # acausal filter was used: shift to reconstrcuted signal
#stim_re = concatenate([zeros(n_shift), stim_re_[0:(len(stim_re_)-n_shift)]]) # acausal filter was used: shift to reconstrcuted signal
istart = int(len(Ker_conv)//2) #+ 0.01/dt
stim_re = stim_re_
resp_cc = resp[int(n_kernel):int(len(resp))-int(n_kernel)] # cut at beginning and end to eliminate conv artefacts
stim_cc = stim[int(n_kernel):int(len(resp))-int(n_kernel)] # cut at beginning and end to eliminate conv artefacts
stim_re_cc = stim_re[istart:istart+int(len(resp))][int(n_kernel):int(len(resp))-int(n_kernel)]
#figure(99)
#plot(stim_cc, 'b')
#plot(stim_re_cc, 'g')
#show()
# old
#stim_re_ = real(convolve(Ker, resp, 'same')) # reconstructed signal
#stim_re = stim_re_
#resp_cc = resp[2*n_kernel:len(resp)-2*n_kernel]
#stim_cc = stim[2*n_kernel:len(stim)-2*n_kernel]
#stim_re_cc = stim_re[2*n_kernel:len(stim_re)-2*n_kernel]
#resp_cc = resp_cc[1/dt:len(resp_cc)-1/dt]
#stim_cc = stim_cc[1/dt:len(stim_cc)-1/dt]
#stim_re_cc = stim_re_cc[1/dt:len(stim_re_cc)-1/dt]
# if n_kernel%2==0:
# Ker_conv = np.r_[np.zeros(len(Ker)-1),Ker[0:-1-1]]
# else:
# Ker_conv = np.r_[np.zeros(len(Ker)),Ker[0:-1]]
#
# window_len = len(Ker_conv)
#
# print "Kernel length:", window_len, n_kernel
#
# resp_conv = np.r_[resp[window_len-1:0:-1],resp,resp[-1:-window_len:-1]]
# resp_conv = resp
#
# start_time = time.time()
# stim_re_ = np.convolve(Ker_conv,resp_conv,mode='valid')
# process_time = time.time() - start_time
#
# print "Convolution took ", process_time
#
# #stim_re = stim_re_[window_len//2:len(stim_re_)-window_len//2]
# stim_re = stim_re_#[window_len:len(stim_re_)-window_len]
#
# resp_cc = resp[window_len//2:len(resp)-window_len//2]
# stim_cc = stim[window_len//2:len(stim)-window_len//2]
# stim_re_cc = stim_re
noise_cc = stim_re_cc - stim_cc # noise
#P_nn_, freq_ = csd(noise_cc, noise_cc, Fs = fs, NFFT = nfft, noverlap = nfft / 4, sides = 'twosided')
#P_ss_, freq_ = csd(stim_cc, stim_cc, Fs = fs, NFFT = nfft, noverlap = nfft / 4, sides = 'twosided')
#SNR_ = P_ss_/P_nn_
err = mean(noise_cc**2)
CF = 1 - (sqrt(err)/std(stim_cc))
VAF = 1 - (err/var(stim_cc))
tc = arange(len(stim_cc))*dt
return tc, resp_cc, stim_cc, stim_re_cc, noise_cc, CF, VAF
def est_quality(t, f, R, amp, pha, fmean):
dt = t[2]-t[1]
n = len(R)
FT = fft(R)[0:round(n / 2)] # transform into frequency domain
freq = fftfreq(n, dt)[0:round(n / 2)] # positive frequency vector
# figure(66)
# subplot(2,1,1)
# plot(freq[1:100], abs(FT[1:100])**2,'*')
# subplot(2,1,2)
# plot(t, R)
# show()
i_f = abs(freq-f).argmin() # find corresponding indexes to frequencies
i_3f = abs(freq-3*f).argmin() # find corresponding indexes to frequencies
NI = log10(abs(FT[i_3f])**2 / abs(FT[i_f])**2)
# Variance Accounted For
R_est = amp*sin(2*pi*f*t+pha)+fmean
N = (R-R_est)
VAF = 1 - (mean(N**2)/var(R))
#figure(99)
#plot(t, R, t, R_est)
#show()
return {'NI':NI, 'VAF':VAF, 'R_est':R_est, 'N':N}
def if_extrap(x, xp, yp):
"""interp function with linear extrapolation"""
#print type(x), type(xp), type(yp)
#print shape(x), shape(xp), shape(yp)
if np.sum((np.diff(xp) < 0)) > (len(xp)/2): # reverse interpolation
xp = xp[::-1]
yp = yp[::-1]
y = interp(x, xp, yp, left = 0, right = 0)
m = mean(diff(yp)) / mean(diff(xp))
y = where(x > xp[-1], yp[-1]+(x-xp[-1])*m, y)
return y
#def adjust_spines(ax,spines, color='k'):
# for loc, spine in ax.spines.iteritems():
# if loc in spines:
# spine.set_position(('outward',10)) # outward by 10 points
# #spine.set_smart_bounds(True)
# else:
# spine.set_color('none') # don't draw spine
#
# # turn off ticks where there is no spine
# if 'left' in spines:
# ax.yaxis.set_ticks_position('left')
# elif 'right' not in spines:
# # no yaxis ticks
# ax.yaxis.set_ticks([])
#
# if 'right' in spines:
# ax.yaxis.set_ticks_position('right')
# s = ax.transAxes.transform((1, 1)) # transform into display coordinates
# e = ax.transAxes.transform((1, 0)) # transform into display coordinates
# line = Line2D((s[0]+10,e[0]+10), (s[1],e[1]), color=color, linewidth=rcParams['axes.linewidth'], transform=None) # ax.transAxes
# line.set_clip_on(False) # show line!
# ax.add_line(line)
#
# #second_right = matplotlib.spines.Spine(ax, 'right', ax.spines['left']._path)
# #second_right.set_position(('outward', 10))
# #ax.spines['second_right'] = second_right
# #ax.spines['second_right'].set_color('k')
# #ax.spines['right'].set_color('k')
#
#
# if 'bottom' in spines:
# ax.xaxis.set_ticks_position('bottom')
# else:
# # no xaxis ticks
# ax.xaxis.set_ticks([])
# ax.axes.get_xaxis().set_visible(False)
def if_spike_train(frequency=100, duration=1000, start_time=0, dt=0.1, seed=None, noise=0, refrac=0, tau=0):
"""
generate constant spike train with memoryless integrate and fire, add noise
"""
ni = noise/2. # adjusted so that noise = +-2*std
t_refrac = -1
rangen = np.random.mtrand.RandomState()
if seed is not None:
rangen.seed(seed)
if np.size(frequency) > 1:
duration = frequency[0][-1]
start_time = frequency[0][0]
time_vec = frequency[0]
freq_vec = frequency[1]
dt = time_vec[2]-time_vec[1]
else:
time_vec = np.arange(start_time,duration,dt)
freq_vec = frequency * np.ones(len(time_vec))
train = np.array([])
u = np.zeros(len(time_vec))
u[0] = rangen.uniform() # start potential
n = np.zeros(len(time_vec))
n[0] = 0
thresh_v = np.zeros(len(time_vec))
if u[0] == 1: # neuron at threshold
train = np.append(train, time_vec[0]) # record spike times
u[0] = 0 # reset
#if tau < 0:
# thresh = 1 + noise * rangen.normal()
#else:
# thresh = 1
thresh = 1
for i, t in enumerate(time_vec):
freq = freq_vec[i]
if tau == 0:
n[i] = freq/1000. + (ni/1000.) * rangen.normal()
elif tau < 0:
n[i] = freq/1000.
else:
n[i] = n[i-1] + (1. - exp(-dt/tau)) * (freq/1000. - n[i-1]) + sqrt(1.-exp(-2.*dt/tau)) * (ni/1000.) * rangen.normal()
u[i] = u[i-1] + dt*n[i]
if t < t_refrac:
if tau < 0:
u[i] = noise * rangen.normal()
else:
u[i] = 0
#if tau < 0:
# thresh = 1 + noise * rangen.normal()
# #if thresh < 0: thresh = 1e-6
elif (u[i] >= thresh): # treshold crossing (u[i-1] < thresh) &
#dts = (u[i]-1)/(u[i]-u[i-1])*dt # interpolate spike time
t_spike = t#+dts
train = np.append(train, t_spike) # spike times
#u[i] = 0
t_refrac = t + refrac
if tau < 0:
u[i] = noise * rangen.normal()
#if u[i] > 1: u[i] = 1-1e-6
else:
u[i] = 0
#if tau < 0:
# thresh = 1 + noise * rangen.normal()
# #if thresh < 0: thresh = 1e-12
#thresh_v[i] = thresh
#print std(n)*1000
#plt.figure(255)
#plt.plot(time_vec, n*1000)
#import matplotlib.pyplot as plt
#freq_times, spike_freq, freq_mean, freq_interp = get_spikefreq(train, time_vec)
#plt.figure(1)
#plt.plot(time_vec,freq_interp, time_vec, freq_vec)
#plt.figure(2)
#plt.plot(time_vec, u)
#plt.show()
#print 'freq_mean: ',freq_mean
#n = thresh_v
return train, u, n # *1000 # freq_vec
def if_spike_train2(frequency=100, duration=1000, start_time=0, dt=0.1, seed=None, noise=0, refrac=0, tau=0):
"""
generate constant spike train with memoryless integrate and fire, add noise
"""
ni = noise/2 # adjusted so that noise = +-2*std
t_refrac = -1
rangen = np.random.mtrand.RandomState()
if seed is not None:
rangen.seed(seed)
if np.size(frequency) > 1:
duration = frequency[0][-1]
start_time = frequency[0][0]
time_vec = frequency[0]
freq_vec = frequency[1]
dt = time_vec[2]-time_vec[1]
else:
time_vec = np.arange(start_time,duration,dt)
freq_vec = frequency * np.ones(len(time_vec))
#print time_vec
#print freq_vec
train = np.array([])
u = np.zeros(len(time_vec))
cur_u = rangen.uniform() # start potential
n = np.zeros(len(time_vec))
cur_n = 0
if cur_u == 1: # neuron at threshold
train = np.append(train, time_vec[0]) # record spike times
cur_u = 0 # reset
u[0] = cur_u
for i, cur_time in enumerate(time_vec):
#freq = freq_vec[mlab.find(time_vec >= cur_time)[0]]/1e3
freq = freq_vec[i]
if tau == 0:
next_n = freq/1000. + (ni/1000.) * rangen.normal()
else:
next_n = cur_n + (1. - exp(-dt/tau)) * (freq/1000. - cur_n) + sqrt(1.-exp(-2.*dt/tau)) * (ni/1000.) * rangen.normal()
stim = cur_n
#stim = freq/1000. + (ni/1000.) * rangen.normal()
next_u = cur_u + dt*stim
if cur_time < t_refrac:
next_u = 0
elif (cur_u < 1) & (next_u >= 1): # treshold crossing
dts = (next_u-1)/(next_u-cur_u)*dt # interpolate spike time
train = np.append(train, cur_time+dts) # spike times
next_u = 0
t_refrac = cur_time + refrac
cur_u = next_u
u[i] = next_u
#print next_n
cur_n = next_n
n[i] = next_n
#print std(n)*1000
#plt.figure(255)
#plt.plot(time_vec, n*1000)
#import matplotlib.pyplot as plt
#freq_times, spike_freq, freq_mean, freq_interp = get_spikefreq(train, time_vec)
#plt.figure(1)
#plt.plot(time_vec,freq_interp, time_vec, freq_vec)
#plt.figure(2)
#plt.plot(time_vec, u)
#plt.show()
#print 'freq_mean: ',freq_mean
return train, u, n*1000 # freq_vec
def spike_train(frequency=100, duration=1000, start_time=0, seed=None, noise=0, jitter=0):
"""
generate constant spike train and add noise and/or jitter
"""
rangen = np.random.mtrand.RandomState()
if seed is not None:
rangen.seed(seed)
if np.size(frequency) > 1:
duration = frequency[0][-1]
start_time = frequency[0][0]
time_vec = frequency[0]
freq_vec = frequency[1]
freq_vec_max = max(freq_vec)/1000.
else:
freq_vec_max = freq_vec/1000.
train = np.array([])
cur_time = duration-10*rangen.uniform() # start with end
while cur_time >= start_time:
train = np.append(train,np.array([cur_time]))
if np.size(frequency) > 1:
freq = freq_vec[mlab.find(time_vec >= cur_time)[0]]/1000.
else:
freq = frequency/1000.
isi = 1/freq
cur_time -= (1. - noise)*isi + noise*isi*rangen.exponential()
train = train[::-1] # reverse
if jitter > 0: # jitter in ms!!!
x = rangen.normal(0, jitter, len(train))
train = train + x
train = np.sort(train)
train = np.delete(train, [i for i,x in enumerate(train) if x < start_time or x > duration])
#print '- method: spike_train'
return train
def mod_spike_train(modulation, noise = 0, seed=None, noise_tau = 0, noise_a=1e9):
"""
generate sinusoidal modulated spike train with and without noise
modulation = (time, rate) with time in seconds, rate in Hz!!
"""
modulation = (modulation[0]/ms, modulation[1].clip(min=0)) # has to be in ms! No negative values!
u = []
if noise_a >= 1e9:
#print "--- train generate from integrate and fire"
train, u, n = if_spike_train(frequency=modulation, seed=seed, noise=noise, tau=noise_tau/ms, refrac=1)
elif noise_a == 1:
#print "--- train generate poisson"
gen = stgen.StGen(rng = np.random.mtrand.RandomState(), seed=seed)
train = gen.inh_poisson_generator(rate=modulation[1], t=modulation[0], t_stop=modulation[0][-1], array=True)
else:
#print "--- train generate gamma"
a = np.ones_like(modulation[0])*noise_a
b = 1/a/modulation[1]
t = modulation[0]
#print shape(a),shape(b),shape(t)
gen = stgen.StGen(rng = np.random.mtrand.RandomState(), seed=seed)
train = gen.inh_gamma_generator(a, b, t=t, t_stop=t[-1], array=True)
return train, u # output in ms
def oscill_spike_train(sor = 4, spike_prob = 0.25, noise_fraction = 4, end_time = 200.e3, seed = 1):
print sor, spike_prob, noise_fraction, end_time
# Subthreshold oscill rate
sor_isi = 1000/sor
spt = np.arange(sor_isi,end_time,sor_isi)
np.random.seed(seed)
rnd = np.random.uniform(0.,1.,spt.size)
spt = np.array([s for (s,r) in zip(spt,rnd) if r < spike_prob])
# Add uniformely distributed noise to spt times
# noise mean = 0 added on the peak of oscillation
# noise amplitude = 1/4 of cycle
spt = np.array([s+np.random.uniform(-1.,1.)*sor_isi/noise_fraction/2 for s in spt])
return spt # output in ms
def slugify(value):
"""
Normalizes string, converts to lowercase, removes non-alpha characters,
and converts spaces to hyphens.
From Django's "django/template/defaultfilters.py".
"""
#import string
#valid_chars = "-_.() %s%s" % (string.ascii_letters, string.digits)
#return ''.join(c for c in value if c in valid_chars)
for c in r'[]/\;,><&*:%=+@!#^\'()$| ?^':
value = value.replace(c,'')
return value
#import re
#_slugify_strip_re = re.compile(r'[^\w\s-]')
#_slugify_hyphenate_re = re.compile(r'[-\s]+')
#import unicodedata
#if not isinstance(value, unicode):
# value = unicode(value)
#value = unicodedata.normalize('NFKD', value).encode('ascii', 'ignore')
#value = unicode(_slugify_strip_re.sub('', value).strip().lower())
#return _slugify_hyphenate_re.sub('_', value)
def gauss_func(x, a, x0, sigma):
return a*np.exp(-(x-x0)**2/(2*sigma**2))
def nans(shape, dtype=float):
a = np.empty(shape, dtype)
a.fill(np.nan)
return a
def nanmean(a, dim=0):
return np.mean(np.ma.masked_array(a,np.isnan(a)),dim).filled(np.nan)
def r_squared(actual, ideal):
actual_mean = np.mean(actual)
ideal_dev = np.sum([(val-actual_mean)**2])
actual_dev = np.sum([(val-actual_mean)**2])
return ideal_dev / actual_dev
#slope, intercept, r_value, p_value; std_err = scipy.stats.linregress(actual, ideal)
#return r_value**2
def rw_hdf5(filepath = "data_dict.hdf5", data_dict = None, export = False):
if data_dict == None: # load
data_dict = {}
print "rw_hdf5 load:", filepath
f = h5py.File(filepath, 'r')
for i in f.items():
data_dict[i[0]] = np.array(i[1])
f.close()
if export:
if not os.path.exists(export):
os.makedirs(export)
shutil.copy(filepath, export)
else:
f = h5py.File(filepath, 'w')
for name in data_dict:
data0 = data_dict[name]
if (type(data0) == np.ndarray):
pass
else:
data0 = np.array(data0)
#print data0
if len(np.shape(data0)) < 1:
f.create_dataset(name, data=data0)
else:
f.create_dataset(name, data=data0, compression='lzf')
f.close()
return data_dict
# test code
if __name__ == '__main__':
test = array(["get_magphase", "fit_sinusoid", "create_multisines", "create_colnoise"])
test = array(["fit_exp"])
test = array(["fit_dualexp"])
test = array(["fit_tripleexp"])
#test = array(["shannon_interp"])
test = array(["construct_Pulsestim"])
test = array(["construct_Test"])
#test = array(["fit_aiftransfer"])
test = array(["create_colnoise"])
#test = array(["test_amp"])
#test = array(["create_ZAP"])
#test = array(["est_quality"])
#test = array(["create_multisines"])
test = array(["syn_kernel"])
#test = array(["aiftransfer"])
#test = array(["create_singlesine"])
#test = array(["if_spike_train_compare_dttest"])
#test = array(["if_spike_train_compare"])
#test = array(["get_spikefreq"])
#test = array(["test_syn"])
#test = array(["hdf5"])
#test = array(["fit_expnew"])
if "test_syn" in test:
duration=1000
start_time=0
dt=0.1
time_vec = np.arange(start_time,duration,dt)
g = np.zeros(len(time_vec))
g[0] = 1
cur_g = 1
tau = 2.86
for i, cur_time in enumerate(time_vec):
s = 0
if (cur_time >= 500) and (cur_time <= 600):
s = 10
next_g = cur_g + dt*(s*(1-cur_g)-cur_g*tau)
cur_g = next_g
g[i] = next_g
figure(1)
plot(time_vec, g)
show()
if "get_spikefreq" in test:
spike_times = np.array([0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1,1.1,1.2,1.3,1.4,1.5,1.6,1.7,1.8,1.9,2,2.2,2.4,2.6,2.8,3,3.2])
freq_times, spike_freq, freq_mean, freq_onset = get_spikefreq(spike_times, stimlength = 3, compute_mean = 2, change_factor = 0.3)
print freq_times, spike_freq, freq_mean, freq_onset
if "get_magphase" in test:
tstart = 0; dt = 0.025*ms; tstop = 7.5 # s
fs = 1/dt # sampling rate
fmax = fs/2 # maximum frequency (nyquist)
t = arange(tstart,tstop,dt) # create stimulus time vector
fu = 0.2 # [Hz]
input_signal = sin(2 * pi * t * fu)
t_out = sort(random.uniform(tstart, tstop, size(t)/10000)) # different timebase! only fit possible!
output_signal = (10 * sin(2 * pi * t_out * fu + pi / 2)) - 70
amp_mean, mag_mean, phase_mean, umean = get_magphase(input_signal, t, output_signal, t_out, method = "fit", f = fu)
print "mean amplitude output signal " + str(amp_mean)
print "mean magnitude " + str(mag_mean) + " ,mean phase = " + str(phase_mean)
plot(t, input_signal ,'--', t_out, output_signal, '*-')
show()
if "fit_sinusoid" in test:
a0 = 3; b0 = 2; c0 = 0; f = 3
t0 = linspace(0, 1, 1000)
noise = rand(1000).T * 0
u0 = a0 + b0 * sin(2 * pi * f * t0 + c0) + noise
t = linspace(0, 1, 1000)
a, b, c, u = fit_sinusoid_fft(f, t0, u0, t)
print a
print b
print c
figure(1)
plot(t0, u0, '.', t, u)
show()
if "fit_exp" in test:
t0 = arange(0,1,0.001)
decay = -10
u0 = decay * (1 - exp(-t0/20e-3))
un0 = u0 + 0.2*np.random.normal(size=len(t0))
udecay, utau, u = fit_exp(t0, un0, p0 = array([-1, 1e-3]))
print udecay, utau
figure(1)
plot(t0, un0, '.', t0, u0, '--', t0, u)
show()
if "fit_dualexp" in test:
dt = 0.025*ms
t0 = arange(0,1,dt)
tau1 = 1*ms
tau2 = 100*ms
maxim = 1
delay = 100*ms
tau_rise=(tau1*tau2)/(tau1-tau2)
tau_fall=tau1
B = 1/( ((tau2/tau1)**(tau_rise/tau1))-((tau2/tau1)**(tau_rise/tau2)) )
print B
u0 = ( maxim*( -exp(-(t0-delay)/tau1)+exp(-(t0-delay)/tau2) ) )
u0[0:delay/dt]=0;
#F = F + ydata(floor(delay/dt))
un0 = u0 + 0.1*np.random.normal(size=len(t0))
umax, utau1, utau2, u = fit_dualexp(t0, un0, p0 = array([1, 1*ms, 1*ms]),delay=delay)
print umax, utau1, utau2
figure(1)
plot(t0, un0, '.', t0, u0, '--', t0, u)
show()
if "fit_tripleexp" in test:
dt = 0.025*ms
t0 = arange(0,0.2,dt)
tau1 = 1*ms
tau2 = 10*ms
tau3 = 1000*ms
maxim = 1
delay = 100*ms
u0 = ( maxim*( -exp(-(t0-delay)/tau1)+exp(-(t0-delay)/tau2)+exp(-(t0-delay)/tau3) ) )
u0[0:delay/dt]=0;
#F = F + ydata(floor(delay/dt))
un0 = u0 + 0.1*np.random.normal(size=len(t0))
umax, utau1, utau2, utau3, u = fit_tripleexp(t0, un0, p0 = array([1, 1*ms, 1*ms, 1*ms]),delay=delay)
print umax, utau1, utau2, utau3
figure(1)
plot(t0, un0, '.', t0, u0, '--', t0, u)
show()
if "fit_expnew" in test:
tau = 500e-3
cutf = 20
sexp = -1
dt = 0.001
t_noise = arange(0, 100, dt)
noise_data = create_colnoise(t_noise, sexp, cutf)
stimulus, t, t_startstop = construct_Stimulus(noise_data, 1/dt, amp=1, ihold = 0, tail_points = 0, delay_baseline = 10)
t_kernel = arange(0, tau*10, dt)
kernel = syn_kernel(t_kernel, 0, tau)
kernel = np.concatenate((zeros(len(kernel)-1),kernel))
filtered = np.convolve(stimulus, kernel, mode='same')
filtered1 = filtered[int(t_startstop[0]/dt):int(t_startstop[1]/dt)]
stimulus1 = stimulus[int(t_startstop[0]/dt):int(t_startstop[1]/dt)]
t1 = t[int(t_startstop[0]/dt):int(t_startstop[1]/dt)]
w_length = 2
fs = 1 / dt
nfft = 2 ** int(ceil(log2(w_length * fs)))
w_true_length = nfft * dt # true length window
noverlap = nfft / 4
r = stimulus1
P_rr, freq = psd(r, Fs = 1/dt, NFFT = nfft, noverlap = noverlap, sides = 'twosided') # , window = mlab.window_none
P_rr = ifftshift(P_rr)
fdata_pos = (P_rr)[1:len(freq)/2] # only positive frequencies
pwr_s = abs(fdata_pos) # Transfer Power
pha_s = angle(fdata_pos, deg = True) # Transfer Phase in degree
r = filtered1
P_rr, freq = psd(r, Fs = 1/dt, NFFT = nfft, noverlap = noverlap, sides = 'twosided') # , window = mlab.window_none
P_rr = ifftshift(P_rr)
fdata_pos = (P_rr)[1:len(freq)/2] # only positive frequencies
pwr_f = abs(fdata_pos) # Transfer Power
pha_f = angle(fdata_pos, deg = True) # Transfer Phase in degree
freq = ifftshift(freq)
freq = freq[1:len(freq)/2]
figure(1)
subplot(2,1,1)
plot(t1, stimulus1/max(stimulus1), 'k-', t1, filtered1/max(filtered1), 'r-')
subplot(2,1,2)
semilogx(freq, pwr_s/pwr_s[0], 'k-', freq, pwr_f/pwr_f[0], 'r-')
show()
if "shannon_interp" in test:
a0 = 3; b0 = 2; c0 = 1.234; freq = 10;
t0 = array([0.00459605074693903,0.00726375896667564,0.00784973441063397,0.0127224903712386,0.0171270655461151,0.0393663919020232,0.0609900368303785,0.0635371421307821,0.0654594032207050,0.0999845349934256,0.119188624525719,0.131158588842870,0.134423875905300,0.138989575327831,0.139516447247645,0.182021984615913,0.189528321172140,0.201667933001159,0.202604832168432,0.218847864233656,0.219479711282702,0.226652337702447,0.228111424603092,0.239969946662225,0.246105341192989,0.250848554299275,0.254989299897020,0.275066518258535,0.296418212747334,0.301828987940621,0.332499026325882,0.348105182140276,0.351266689715571,0.354722452499086,0.359590326763347,0.374299226239647,0.383307652959454,0.383364191948955,0.389899994370318,0.404459485959467,0.421426947556085,0.442676684844202,0.449368223144463,0.456522496631530,0.469306274142778,0.490669133318145,0.492682694168636,0.504719391132661,0.507205621367988,0.518527193799826,0.521561514445141,0.534135096228315,0.568392189742432,0.568484472347038,0.573023005985662,0.577029656618231,0.582603276061774,0.584055978860854,0.587796282178598,0.588501277932966,0.608464044116267,0.625893464986574,0.626008271991061,0.627958771815969,0.630538464601416,0.648639994155523,0.651808804576273,0.660791197220462,0.680866786193684,0.690596392893720,0.699105430999923,0.704644674745085,0.707497758866246,0.724810615824970,0.726322828528511,0.752762163986078,0.758617951982347,0.774559534114423,0.780486968166248,0.791867080004023,0.799655776201573,0.802171594064305,0.848515015994186,0.856540805722978,0.864847431324938,0.876360704800671,0.877835683926097,0.883152342556065,0.884554440010395,0.885558670718519,0.929895247208519,0.941655353501460,0.945441426048917,0.946800026376379,0.948273065241326,0.957117920800973,0.978062788376808,0.989045885911051,0.994412760780681,0.995238107661316])
t = linspace(0, 1, 1000)
u0 = a0 + b0 * sin(2 * pi * freq * t0 + c0)
u = shannon_interp(t0, u0, t, chunk_size=10)
print len(u)
plot(t0, u0, "*b", t, u, "-r")
axis(ymin=0, ymax=6)
show()
if "create_multisines" in test:
fstart = 1; fstop = 100; fsteps = 1 # Hz
freq_used = arange(fstart, fstop, fsteps) # frequencies we want to examine
tstart = 0; t_stim = 10; dt = 0.025*ms # s
fs = 1/dt # sampling rate
t_noise = arange(tstart, t_stim, dt) # create stimulus time vector, make sure stimulus is even!!!
noise_data, f1, freq_wp, f_used_check = create_multisines(t_noise, freq_used) # create multi sine signal
pxx, freq = psd(noise_data, Fs = 1 / dt) # compute power spectrum
data_points = len(t_noise)
noise_power = abs(fft(noise_data))[ 0:round(data_points / 2) ] # compute noise power
figure(1); clf()
subplot(3, 1, 1)
plot(t_noise, noise_data)
subplot(3, 1, 2)
loglog(freq, pxx / pxx[1])
subplot(3, 1, 3)
semilogx(f1, noise_power)
show()
if "create_colnoise" in test:
dt = 0.025e-3
fs = 1/dt
t = arange(0, 100, dt)
# sexp = 0 # reasonably smooth
# cutf= 10 # 0.01kHz = 10Hz frequency cutoff
# #cutf= None
# x = create_colnoise(t, sexp, cutf)
#
# cutf= 10
# dt = t[2] - t[1]
#
# # display power spectrum
# figure(2)
# pxx, f = psd(x, Fs=1/dt)
# figure(1); #clf()
# subplot(2, 1, 1)
# plot(t, x)
# subplot(2, 1, 2)
# semilogx(f, pxx / pxx[1],
# [f[1], cutf], [1, 1], 'r:', # white
# [cutf, f[-1]], [1, (cutf / f[-1]) ** sexp], 'r:') # coloured
w_length = 1
nfft = 2 ** int(ceil(log2(w_length * fs)))
sexp = 0 #-1
cutf= 10 # 0.01kHz = 10Hz frequency cutoff
#cutf= None
sexp = 4 #-1
cutf= 20 # 0.01kHz = 10Hz frequency cutoff
x = create_colnoise(t, sexp, cutf)
figure(3)
# the histogram of the data
n, bins, patches = hist(x, 100, normed=True, facecolor='green', alpha=0.75)
# add a 'best fit' line
y = normpdf(bins, 0, 1)
plot(bins, y, 'r--', linewidth=1)
# display power spectrum
figure(2)
P, f = psd(x, Fs = fs, NFFT = nfft, noverlap = nfft / 4, sides = 'twosided')
P = ifftshift(P)[0:len(f)/2]
f = ifftshift(f)[0:len(f)/2] # only positive frequencies
figure(1)
subplot(2, 1, 1)
print shape(t), shape(x)
plot(t, x)
subplot(2, 1, 2)
semilogx(f,P)
sexp = -1
cutf= 20 # 0.01kHz = 10Hz frequency cutoff
#cutf= None
onf = 10
x = create_colnoise(t, sexp, cutf, onf=onf)
figure(3)
n, bins, patches = hist(x, 100, normed=True, facecolor='blue', alpha=0.75)
# display power spectrum
figure(2)
P, f = psd(x, Fs = fs, NFFT = nfft, noverlap = nfft / 4, sides = 'twosided')
P = ifftshift(P)[0:len(f)/2]
f = ifftshift(f)[0:len(f)/2] # only positive frequencies
figure(1)
subplot(2, 1, 1)
plot(t, x)
subplot(2, 1, 2)
semilogx(f,P)
show()
#from scipy import signal
#FS = 1000.0 # sampling rate
#FC = 0.05/(0.5*FS) # cutoff frequency at 0.05 Hz
#N = 1001 # number of filter taps
#a = 1 # filter denominator
#b = signal.firwin(N, cutoff=FC, window='hamming') # filter numerator
#
#M = FS*60 # number of samples (60 seconds)
#n = arange(M) # time index
#x1 = cos(2*pi*n*0.025/FS) # signal at 0.025 Hz
#x = x1 + 2*rand(M) # signal + noise
#y = signal.lfilter(b, a, x) # filtered output
#
#plot(n/FS, x); plot(n/FS, y, 'r') # output in red
#grid()
if "test_amp" in test:
dt = 0.025e-3
fs = 1/dt
t1 = arange(0, 5, dt)
sexp = -1
cutf= 100 # 0.01kHz = 10Hz frequency cutoff
#sexp = 0
#cutf= 0 # 0.01kHz = 10Hz frequency cutoff
x = create_colnoise(t1, sexp, cutf)
print std(x)
print var(x)
print std(x)**2
a = 0.1
x = x*a
t, stimulus, i_startstop, t_startstop = create_singlesine(fu = 5, amp = a, ihold = 0, dt = dt, periods = 1, minlength = 1, t_prestim = 1)
figure(1)
plot(t1, x, t, stimulus)
show()
if "create_singlesine" in test:
t, stimulus, i_startstop, t_startstop = create_singlesine(fu = 100, amp = 0.5, ihold = 1, dt = 0.1*ms, periods = 1, minlength = 1, t_prestim = 1)
plot(t, stimulus)
show()
if "aiftransfer" in test:
freq = arange(0,1000,1)
H, H0 = aiftransfer(freq = freq, tau = 10*ms, f0 = 100, i0 = 1*nA, rm = 10*MOhm, Vreset = -70, Vth = -60, Vrest = -70)
H1, H0 = aiftransfer(freq = freq, tau = 10*ms, f0 = 100, i0 = 1*nA, rm = 10*MOhm, Vreset = -70, Vth = -60, Vrest = -70, delta_t = 10e-3)
magfA = abs(H)
phafA = unwrap(angle(H)) * (180 / pi)
magfA1 = abs(H1)
phafA1 = unwrap(angle(H1)) * (180 / pi)
subplot(2,1,1)
semilogx(freq, magfA, '--', label = "analytical model")
semilogx(freq, magfA1, 'r--', label = "analytical model")
subplot(2,1,2)
semilogx(freq, phafA, '--', label = "analytical model")
semilogx(freq, phafA1, 'r--', label = "analytical model")
figure(2)
H = H/H0
mag = abs(H)
loglog(freq, mag)
show()
if "syn_kernel" in test:
t = arange(-2,2,1e-3)
G = syn_kernel(t, 0.2, 0.2)
G[0:2/1e-3] = 0
plot(t,G)
show()
if "fit_aiftransfer" in test:
freq = arange(0,1000,1)
tau0 = 10*ms
f0 = 80*Hz
i0 = 0.010*nA
H0 = aiftransfer(freq = freq, tau = tau0, f0 = f0, i0 = i0, rm = 10*MOhm, Vreset = -70, Vth = -60, Vrest = -70)
tau, scale, H = fit_aiftransfer(freq, H0, f0, i0)
print tau0
print tau
print scale
magfA0 = abs(H0)
phafA0 = unwrap(angle(H0)) * (180 / pi)
magfA = abs(H)
phafA = unwrap(angle(H)) * (180 / pi)
subplot(2,1,1)
semilogx(freq, magfA0, freq, magfA, '--', label = "analytical model")
subplot(2,1,2)
semilogx(freq, phafA0, freq, phafA, '--', label = "analytical model")
show()
if "construct_Pulsestim" in test:
#t, ivec = construct_Pulsestim(pulses = 2, dt = 2.5e-05, latency = 0.0091, stim_start = 0.01, stim_end = 0.01, len_pulse = 0.0002, amp_init = 0.64, amp_next = 0.24)
t, ivec = construct_Pulsestim(dt = 1e-3, latency = [1,2,3,4], stim_start = 1.5, stim_end = 0.5, len_pulse = 0.1, amp_next = [1,-1,0.5,0.3])
print len(t)
print len (ivec)
plot(t, ivec)
show()
# t = arange(0,1,1e-3)
# x = zeros(len(t))
# x[0:100] = 1
# x[900:990] = -1
# plot(t,x)
#
# X = fft(x)
# plot(fftshift(X))
#
# x2 = ifft(X)
# plot(x2)
#
# show()
if "construct_Test" in test:
import scipy
from scipy import io
dt = 1*ms
pulses = 50
istep = np.random.rand(pulses)-0.5
print istep
tstep = np.arange(0,pulses*0.5,0.5)
tstep = tstep + np.random.rand(pulses)*(0.3)
t, ivec = construct_Pulsestim(dt = dt, latency = np.append(diff(tstep), 0), stim_start = 0*s, stim_end = 1*s, len_pulse = 0.01*s, amp_next = istep)
ivec = ivec + 0.5
istep = -1*istep
t, ivec2 = construct_Pulsestim(dt = dt, latency = np.append(diff(tstep), 0), stim_start = 0*s, stim_end = 1*s, len_pulse = 0.01*s, amp_next = istep)
ivec2 = ivec2 + 0.5
#cutf = 20
#sexp = -1
#t = arange(0, 12, dt)
#seed = 3
#ivec = create_colnoise(t, sexp, cutf, seed)*0.2 + 0.5
tau = 0.06*s
t_kernel = arange(0, 0.5*s, dt)
kernel = syn_kernel(t_kernel, tau, tau)
kernel = np.concatenate((zeros(len(kernel)-1),kernel))
ovec = np.convolve(ivec, kernel, mode='same')
ivec = ivec[int(1/dt):int(21/dt)]
ivec2 = ivec2[int(1/dt):int(21/dt)]
ovec = ovec[int(1/dt):int(21/dt)]
t = t[int(1/dt):int(21/dt)]
ovec = ovec/max(ovec)
print len(ovec), len(ivec)
print "Saving .mat"
data = {}
data['inputSequence'] = np.array([ones(len(ivec)), ivec]) #
data['outputSequence'] = ovec
scipy.io.savemat('./input.mat',data)
plot(t, ivec, 'k', t, ovec, 'r',)
show()
if "create_ZAP" in test:
ihold = 1
t, zap, f = create_ZAP(ihold=ihold, ex=2)
# Find all indices right before a rising-edge zero crossing
zap_ = zap - ihold
indices = find((zap_[1:] >= 0) & (zap_[:-1] < 0))
# More accurate, using linear interpolation to find intersample
# zero-crossings (Measures 1000.000129 Hz for 1000 Hz, for instance)
crossings = [i - zap_[i] / (zap_[i+1] - zap_[i]) for i in indices]
f_est = 1/diff(t[crossings])
t_est = t[crossings[:-1]]
plot(t, zap)
twinx()
plot(t, f, 'b', t_est, f_est, 'r')
show()
if "if_spike_train" in test:
train, u = if_spike_train(frequency=100, duration=10/ms, start_time=0, dt=0.025, seed=30, noise=10, tau=0.1)
rate = 1/diff(train*ms)
figure(1)
plot(train[:-1], rate, 'b')
figure(2)
plot(train, ones(len(train)), 'b|')
#train, u = if_spike_train(frequency=100, duration=10/ms, start_time=0, dt=0.0025, seed=30, noise=10, tau=100)
#rate = 1/diff(train*ms)
#figure(1)
#plot(train[:-1], rate, 'r')
#figure(2)
#plot(train, ones(len(train)), 'r|')
show()
if "if_spike_train_compare_dttest" in test:
import numpy
from NeuroTools import stgen, signals
from pylab import zeros_like, plot
import sys
import matplotlib.pyplot as plt
import time
tsim = 100000
fmean = 40
noise = 0.01*40
tau = 5
#noise = 0.36*40
#tau = 5
#noise = 1*40
#tau = 5
#noise = 40
#tau = 0
noise = 0.1
tau = -1
dt = 0.025 #10.0
t = numpy.arange(0,tsim,dt)
rate = numpy.ones_like(t)*fmean
modulation = (t, rate) # has to be in ms!
t0 = time.time()
train, u, n = if_spike_train(frequency=modulation, seed=1, noise=noise, tau=tau)
tsolve = time.time() - t0
values, xaxis = numpy.histogram(train, bins=50)
fano = numpy.var(values)/numpy.mean(values)
print numpy.mean(1/diff(train)), numpy.std(diff(train))/numpy.mean(diff(train)), tsolve, fano
plt.figure(1)
plt.plot(t, n, 'b')
plt.figure(2)
plt.plot(train[1:], 1e3/diff(train),"*-b")
figure(4)
values, xaxis = numpy.histogram(diff(train), bins=50)
values = values/float(values.sum())
xaxis = xaxis[:-1]
plt.plot(xaxis,values, color="b")
dt = 0.1 #10.0
t = numpy.arange(0,tsim,dt)
rate = numpy.ones_like(t)*fmean
modulation = (t, rate) # has to be in ms!
t0 = time.time()
train, u, n = if_spike_train(frequency=modulation, seed=1, noise=noise, tau=tau)
tsolve = time.time() - t0
values, xaxis = numpy.histogram(train, bins=50)
fano = numpy.var(values)/numpy.mean(values)
print numpy.mean(1/diff(train)), numpy.std(diff(train))/numpy.mean(diff(train)), tsolve, fano
plt.figure(1)
plt.plot(t, n, 'r')
plt.figure(2)
plt.plot(train[1:], 1e3/diff(train),"*-r")
figure(4)
values, xaxis = numpy.histogram(diff(train), bins=50)
values = values/float(values.sum())
xaxis = xaxis[:-1]
plt.plot(xaxis,values, color="r")
fmean = 60
dt = 0.01 #10.0
t = numpy.arange(0,tsim,dt)
rate = numpy.ones_like(t)*fmean
modulation = (t, rate) # has to be in ms!
t0 = time.time()
train, u, n = if_spike_train(frequency=modulation, seed=1, noise=noise, tau=tau)
tsolve = time.time() - t0
values, xaxis = numpy.histogram(train, bins=50)
fano = numpy.var(values)/numpy.mean(values)
print numpy.mean(1/diff(train)), numpy.std(diff(train))/numpy.mean(diff(train)), tsolve, fano
plt.figure(1)
plt.plot(t, n, 'm')
plt.figure(2)
plt.plot(train[1:], 1e3/diff(train),"*-m")
figure(4)
values, xaxis = numpy.histogram(diff(train), bins=50)
values = values/float(values.sum())
xaxis = xaxis[:-1]
plt.plot(xaxis,values, color="m")
# NEURON FASTER:
# import numpy as np
# from neuron import h
# from cells.IfCell import *
#
# dt = 0.1
#
# t0 = time.time()
#
# cell = IfCell(C = 1, R = 1e12, e = 0, thresh = 1, vrefrac = 0)
#
# gid = 0
# noiseRandObj = h.Random() # provides NOISE with random stream
# fluct = h.Ifluct2(cell.soma(0.5))
# fluct.m = 0
# fluct.s = noise/2
# fluct.tau = tau
# fluct.noiseFromRandom(noiseRandObj) # connect random generator!
# noiseRandObj.MCellRan4(1, gid+1) # set lowindex to gid+1, set highindex to > 0
# noiseRandObj.normal(0,1)
#
# stim = h.IClamp(cell.soma(0.5))
# stim.amp = fmean
# stim.delay = 0
# stim.dur = tsim
#
# tstop = tsim
#
# rec_t = h.Vector()
# rec_t.record(h._ref_t)
#
# rec_v = h.Vector()
# rec_v.record(cell.soma(0.5)._ref_v)
#
# rec_s = h.Vector()
# nc = cell.connect_target(None) # threshold is set in neuron definition, or here!
# nc.record(rec_s) # record indexes of the positive zero crossings
#
#
# h.load_file("stdrun.hoc")
# h.celsius = 0
# h.init()
# h.tstop = tstop
# h.dt = dt
# h.steps_per_ms = 1 / dt
#
# h.v_init = 0
#
# h.run()
#
# t1 = np.array(rec_t)
# voltage = np.array(rec_v)
# train = np.array(rec_s)
#
# tsolve = time.time() - t0
#
# values, xaxis = numpy.histogram(train, bins=50)
# fano = numpy.var(values)/numpy.mean(values)
#
# print numpy.mean(1/diff(train)), numpy.std(diff(train))/numpy.mean(diff(train)), tsolve, fano
#
# #plt.figure(1)
# #plt.plot(t, n, 'g')
#
# plt.figure(2)
# plt.plot(train[1:], 1e3/diff(train),"*-g")
#
# plt.figure(3)
# plt.plot(t1, voltage,"g")
#
# figure(4)
# values, xaxis = numpy.histogram(diff(train), bins=50)
# values = values/float(values.sum())
# xaxis = xaxis[:-1]
# plt.plot(xaxis,values, color="g")
#
plt.show()
if "if_spike_train_compare" in test:
# Generate the PSTH for an inhomogeneous gamma renewal process
# with a step change in the rate (b changes, a stays fixed)
# This script generates Figure 5 in:
# Eilif Muller, Lars Buesing, Johannes Schemmel, and Karlheinz Meier
# Spike-Frequency Adapting Neural Ensembles: Beyond Mean Adaptation and Renewal Theories
# Neural Comput. 2007 19: 2958-3010.
# i.e. the PSTH for a 2D adapting markov process undergoing
# a step change in statistics due to a step stimulus.
import numpy
from NeuroTools import stgen, signals
from pylab import zeros_like, plot
import sys
import matplotlib.pyplot as plt
import time
tsim = 10000.0
dt = 0.025 #10.0
t = numpy.arange(0,tsim,dt)
bins = numpy.arange(0,tsim,10)
trials = 1 #5000
fmean = 40
# stepup
i_start = t.searchsorted(4000.0,'right')-1
i_end = t.searchsorted(6000.0,'right')-1
rate = numpy.ones_like(t)*fmean #11.346 #*20.0
#rate[i_start:i_end] = 2*fmean #92.25 #40.0
psth = zeros_like(bins)
stg = stgen.StGen()
# POISSON
a = numpy.ones_like(t)*1
b = 1/a/rate
print "Running %d trials of %.2f milliseconds" % (trials, tsim)
for i in xrange(trials):
if i%100==0:
print "%d" % i,
sys.stdout.flush()
stg.seed(i)
t0 = time.time()
#st = stg.inh_gamma_generator(a,b,t,tsim,array=True)
st = stg.inh_poisson_generator(rate, t, tsim, array=True)
tsolve = time.time() - t0
s1 = signals.SpikeTrain(st)
isi = s1.isi()
fano_factor_isi = numpy.var(isi)/numpy.mean(isi)
cv_isi = numpy.std(isi)/numpy.mean(isi)
figure(3)
values, xaxis = numpy.histogram(s1.isi(), bins=50)
values = values/float(values.sum())
xaxis = xaxis[:-1]
plt.plot(xaxis,values, color="k")
values, xaxis = numpy.histogram(s1.spike_times, bins=50)
fano = numpy.var(values)/numpy.mean(values)
print "POISSON, CV: ", cv_isi, "MEAN: ", s1.mean_rate(), "FANO:", fano, "tsolve:", tsolve, "s"
psth[1:]+=numpy.histogram(st,bins)[0]
figure(2)
plt.plot(st, 0*numpy.ones_like(st),"|k")
print "\n"
# normalize
psth = psth.astype(float)
psth/= dt*float(trials)/1000.0
# this is for correct 'steps' plotting only
psth[0] = psth[1]
figure(1)
plt.plot(bins,psth,linestyle='steps', color="k")
psth = zeros_like(bins)
stg = stgen.StGen()
# LOW NOISE IF
noise = 0.01
noise_tau = -1
print "Running %d trials of %.2f milliseconds" % (trials, tsim)
for i in xrange(trials):
if i%100==0:
print "%d" % i,
sys.stdout.flush()
#stg.seed(i)
#st = stg.inh_gamma_generator(a,b,t,1000.0,array=True)
t0 = time.time()
st, u = mod_spike_train(modulation=[t*ms, rate], noise = noise, seed=i, noise_tau = noise_tau)
tsolve = time.time() - t0
s1 = signals.SpikeTrain(st)
figure(3)
values, xaxis = numpy.histogram(s1.isi(), bins=50)
values = values/float(values.sum())
xaxis = xaxis[:-1]
plt.plot(xaxis,values, color="b")
values, xaxis = numpy.histogram(s1.spike_times, bins=50)
fano = numpy.var(values)/numpy.mean(values)
print "Low Noise IF, NO TAU, CV: ", s1.cv_isi(), "MEAN: ", s1.mean_rate(), "FANO:", fano, "tsolve:", tsolve, "s"
psth[1:]+=numpy.histogram(st,bins)[0]
figure(2)
plt.plot(st, 1*numpy.ones_like(st),"|b")
figure(4)
plt.plot(u, "b")
#plt.plot(n, "b--")
print "\n"
# normalize
psth = psth.astype(float)
psth/= dt*float(trials)/1000.0
# this is for correct 'steps' plotting only
psth[0] = psth[1]
figure(1)
plt.plot(bins,psth,linestyle='steps', color="b")
psth = zeros_like(bins)
stg = stgen.StGen()
# HIGH NOISE IF, NO TAU
noise = 10
noise_tau = -1
print "Running %d trials of %.2f milliseconds" % (trials, tsim)
for i in xrange(trials):
if i%100==0:
print "%d" % i,
sys.stdout.flush()
t0 = time.time()
st, u = mod_spike_train(modulation=[t*ms, rate], noise = noise, seed=i, noise_tau = noise_tau)
tsolve = time.time() - t0
s1 = signals.SpikeTrain(st)
figure(3)
values, xaxis = numpy.histogram(s1.isi(), bins=50)
values = values/float(values.sum())
xaxis = xaxis[:-1]
plt.plot(xaxis,values, color="r")
values, xaxis = numpy.histogram(s1.spike_times, bins=50)
fano = numpy.var(values)/numpy.mean(values)
print "High Noise IF, NO TAU, CV: ", s1.cv_isi(), "MEAN: ", s1.mean_rate(), "FANO:", fano, "tsolve:", tsolve, "s"
psth[1:]+=numpy.histogram(st,bins)[0]
figure(2)
plt.plot(st, 2*numpy.ones_like(st),"|r")
figure(4)
plt.plot(u, "r")
#plt.plot(n, "r--")
print "\n"
# normalize
psth = psth.astype(float)
psth/= dt*float(trials)/1000.0
# this is for correct 'steps' plotting only
psth[0] = psth[1]
figure(1)
plt.plot(bins,psth,linestyle='steps', color="r")
psth = zeros_like(bins)
stg = stgen.StGen()
# MEDIUM NOISE IF, NO TAU
noise = 0.5
noise_tau = -1
print "Running %d trials of %.2f milliseconds" % (trials, tsim)
for i in xrange(trials):
if i%100==0:
print "%d" % i,
sys.stdout.flush()
t0 = time.time()
st, u = mod_spike_train(modulation=[t*ms, rate], noise = noise, seed=i, noise_tau = noise_tau)
tsolve = time.time() - t0
s1 = signals.SpikeTrain(st)
figure(3)
values, xaxis = numpy.histogram(s1.isi(), bins=50)
values = values/float(values.sum())
xaxis = xaxis[:-1]
plt.plot(xaxis,values, color="m")
values, xaxis = numpy.histogram(s1.spike_times, bins=50)
fano = numpy.var(values)/numpy.mean(values)
print "Medium Noise IF, NO TAU, CV: ", s1.cv_isi(), "MEAN: ", s1.mean_rate(), "FANO:", fano, "tsolve:", tsolve, "s"
psth[1:]+=numpy.histogram(st,bins)[0]
figure(2)
plt.plot(st, 3*numpy.ones_like(st),"|m")
figure(4)
plt.plot(u, "g")
#plt.plot(n, "r--")
print "\n"
# normalize
psth = psth.astype(float)
psth/= dt*float(trials)/1000.0
# this is for correct 'steps' plotting only
psth[0] = psth[1]
figure(1)
plt.plot(bins,psth,linestyle='steps', color="m")
psth = zeros_like(bins)
stg = stgen.StGen()
# # GAMMA
# a = numpy.ones_like(t)*100
# b = 1/a/rate
#
# print "Running %d trials of %.2f milliseconds" % (trials, tsim)
# for i in xrange(trials):
# if i%100==0:
# print "%d" % i,
# sys.stdout.flush()
# stg.seed(i)
#
# t0 = time.time()
#
# st = stg.inh_gamma_generator(a,b,t,tsim,array=True)
#
# tsolve = time.time() - t0
#
# s1 = signals.SpikeTrain(st)
#
# isi = s1.isi()
# fano_factor_isi = numpy.var(isi)/numpy.mean(isi)
# cv_isi = numpy.std(isi)/numpy.mean(isi)
#
# figure(3)
# values, xaxis = numpy.histogram(s1.isi(), bins=50)
# values = values/float(values.sum())
# xaxis = xaxis[:-1]
# plt.plot(xaxis,values, color="g")
#
# values, xaxis = numpy.histogram(s1.spike_times, bins=50)
# fano = numpy.var(values)/numpy.mean(values)
#
# print "GAMMA, CV: ", cv_isi, "MEAN: ", s1.mean_rate(), "FANO:", fano, "tsolve:", tsolve, "s"
#
#
# psth[1:]+=numpy.histogram(st,bins)[0]
# figure(2)
# plt.plot(st, 4*numpy.ones_like(st),"|g")
#
# print "\n"
#
# # normalize
# psth = psth.astype(float)
# psth/= dt*float(trials)/1000.0
# # this is for correct 'steps' plotting only
# psth[0] = psth[1]
# figure(1)
# plt.plot(bins,psth,linestyle='steps', color="g")
# psth = zeros_like(t)
#
# a = numpy.ones_like(t)*40 # 11.346
# bq = numpy.ones_like(t)*0.1231*14.48
# a[i_start:i_end] = 40# 92.25
# bq[i_start:i_end] = 0.09793*14.48
#
# tau_s = 110.0
# tau_r = 1.97
# qrqs = 221.96
# t_stop = 1000.0
#
# #a is in units of Hz. Typical values are available
# #in Fig. 1 of Muller et al 2007, a~5-80Hz (low to high stimulus)
#
# #bq here is taken to be the quantity b*q_s in Muller et al 2007, is thus
# #dimensionless, and has typical values bq~3.0-1.0 (low to high stimulus)
#
# #qrqs is the quantity q_r/q_s in Muller et al 2007,
# #where a value of qrqs = 3124.0nS/14.48nS = 221.96 was used.
#
# #tau_s has typical values on the order of 100 ms
# #tau_r has typical values on the order of 2 ms
#
#
# print "Running %d trials of %.2f milliseconds" % (trials, tsim)
# for i in xrange(trials):
# if i%100==0:
# print "%d" % i,
# sys.stdout.flush()
# #stg.seed(i)
# st = stg.inh_2Dadaptingmarkov_generator(a,bq,tau_s,tau_r,qrqs,t,t_stop,array=True)
#
# psth[1:]+=numpy.histogram(st,t)[0]
# figure(3)
# plt.plot(st, i*numpy.ones_like(st),"|r")
#
# print "\n"
#
# # normalize
#
# psth = psth.astype(float)
# psth/= dt*float(trials)/1000.0 #psth/= dt*10000.0/1000.0
#
# # this is for correct 'steps' plotting only
# psth[0] = psth[1]
#
# figure(1)
# plt.plot(t,psth,linestyle='steps', color="r")
plt.show()
if "plot_train" in test:
t = arange(0, 1, 0.001)
noise = 40*create_colnoise(t, -1, 20, seed=4)+40
noise = ones(len(t))*6
modulation = (t, noise)
train, _ = mod_spike_train(modulation)
print train
figure()
subplot(1,2,1)
plot(t, noise, 'b')
subplot(1,2,2)
plot(train*1e-3, ones(len(train)), 'b|')
if "plot_train2" in test:
t = arange(0, 1, 0.001)
noise = create_colnoise(t, -1, 20, seed=4)
addnoise = create_colnoise(t, -1, 100, seed=4)
noise1 = 10*40*noise+40
noise2 = -10*40*noise+40
noise1 = noise1.clip(min=0)
noise2 = noise2.clip(min=0)
#noise = ones(len(t))*6
noise_a = 1; noise_syn = 0; noise_syn_tau = 0 # poisson
noise_a = 1e9; noise_syn = 0.5; noise_syn_tau = -1 # medium noise
modulation1 = (t, noise1)
train1, _ = mod_spike_train(modulation1, noise = noise_syn, seed = 1, noise_tau = noise_syn_tau, noise_a = noise_a)
modulation2 = (t, noise2)
train2, _ = mod_spike_train(modulation2, noise = noise_syn, seed = 1, noise_tau = noise_syn_tau, noise_a = noise_a)
tau = 0.1
t_kernel = np.arange(0, tau*6, 0.001)
kernel1 = syn_kernel(t_kernel, 0, tau)
kernel = np.concatenate((np.zeros(len(kernel1)-1),kernel1))
print "- Basis convolution"
noise0 = concatenate([noise,noise,noise,noise,noise])
noise_conv0 = np.convolve(noise0, kernel, mode='same')
noise_conv = noise_conv0[3000:4000]
figure()
subplot(4,2,1)
plot(t, noise, 'b')
plt.axis([0.25, 0.8, -1.5, 1.5])
subplot(4,2,3)
plot(t, noise1, 'b')
plt.axis([0.25, 0.8, -10, 650])
subplot(4,2,5)
plot(train1*1e-3, ones(len(train1)), 'b|')
plt.axis([0.25, 0.8, 0.94, 1.06])
subplot(4,2,4)
plot(t, noise2, 'b')
plt.axis([0.25, 0.8, -10, 650])
subplot(4,2,6)
plot(train2*1e-3, ones(len(train2)), 'b|')
plt.axis([0.25, 0.8, 0.94, 1.06])
subplot(4,2,7)
noisy=concatenate([noise_conv[0:600]+6*addnoise[0:600],noise_conv[600:1000]+20*addnoise[600:1000]])
plot(t, noise_conv, 'b', t, noisy, 'r')
plt.axis([0.25, 0.8, -60, 60])
subplot(4,2,8)
plot(t_kernel, kernel1, 'b')
plt.axis([0, 0.55, 0, 1.1])
if "est_quality" in test: # Shows that SNR depends on noise!!!
tstart = 0; dt = 1*ms; tstop = 100 # s
fs = 1/dt # sampling rate
fmax = fs/2 # maximum frequency (nyquist)
t = arange(tstart,tstop,dt) # create stimulus time vector
fu = 10 # [Hz]
amp = 0.2
ihold = 0
n_in = 0.2*create_colnoise(t, sexp = 4, cutf = None)
addmean = 0 #40
ampmult = 10
# SINGLE SINE, NOISE N from fit
input_signal = amp*sin(2 * pi * t * fu) + ihold
output_signal = n_in+(ampmult * amp * sin(2 * pi * t * fu + 0*(pi / 2)) + ihold + addmean)
amp_mean, mag_mean, phase_mean, fmean = get_magphase(input_signal, t, output_signal, t, method = "fft", f = fu)
print "mean amplitude output signal " + str(amp_mean)
print "mean magnitude " + str(mag_mean) + " ,mean phase = " + str(phase_mean)
results = est_quality(t, fu, output_signal, amp*mag_mean, phase_mean/ (180 / pi), fmean)
NI, VAF, R_est, N = results.get('NI'), results.get('VAF'), results.get('R_est'), results.get('N')
#amp_mean, mag_mean, phase_mean, fmean = get_magphase(input_signal, t, N, t, method = "fft", f = fu)
SNR1 = mean(input_signal**2) / std(N)
SNR2 = mean(input_signal**2) / mean(N**2)
print "- NI: " + str(NI) + ", VAF: " + str(VAF) + ", SNR2: " + str(SNR2) + ", SNR1: " + str(SNR1)
figure(11)
Ps, f = psd(input_signal, Fs=1/dt)
Pn, f = psd(N, Fs=1/dt)
SNRf = Ps/Pn
nfft = 2 ** int(ceil(log2(2 * fs)))
P_nn, freq = csd(N, N, Fs = fs, NFFT = nfft, noverlap = nfft / 4, sides = 'twosided')
P_ss, freq = csd(input_signal, input_signal, Fs = fs, NFFT = nfft, noverlap = nfft / 4, sides = 'twosided')
SNRx = ifftshift(P_ss / P_nn)[0:len(freq)/2] # 10*log10()
freq = ifftshift(freq)[0:len(freq)/2] # only positive frequencies
figure(12)
semilogx(f, SNRf, 'r', freq, SNRx, 'g')
print mean(SNRf)
print mean(real(SNRx))
figure(1)
plot(t, input_signal ,'g--', t, output_signal, 'k-', t, R_est, 'b-', t, N, 'r-')
# NOISE, NOISE N from reconstruction
method_interpol = array(["linear"])
input_signal = amp*create_colnoise(t, sexp = 4, cutf = None)
output_signal = n_in + (ampmult * input_signal + addmean)
noise_data_points = len(input_signal)
w_length = 5*s
results = compute_Transfer(output_signal, t, input_signal, t, noise_data_points, gsyn = output_signal, freq_wp = None, method_interpol = method_interpol, do_csd = 1, nc_delay = 0, w_length = w_length, t_kernel = 0, t_qual = 990)
freq, tc, stim_re_mat, stim, resp_mat, noise_mat, CF_mat, SNR_mat = results.get('freq_used'), results.get('tc'), results.get('stim_re_mat'), results.get('stim'), results.get('resp_mat'), results.get('noise_mat'), results.get('CF_mat'), results.get('SNR_mat')
tk, K_mat = results.get('tk'), results.get('K_mat')
SNR = SNR_mat[1][0,:]
figure(2)
semilogx(freq, SNR, 'b')
print "- SNR(f=" + str(freq[(freq>=9.8) & (freq<=10)]) + "): " + str(real(SNR[(freq>=9.8) & (freq<=10)]))
print mean(SNR)
w_length = 0.5*s
results = compute_Transfer(output_signal, t, input_signal, t, noise_data_points, gsyn = output_signal, freq_wp = None, method_interpol = method_interpol, do_csd = 1, nc_delay = 0, w_length = w_length, t_kernel = 0, t_qual = 990)
freq, tc, stim_re_mat, stim, resp_mat, noise_mat, CF_mat, SNR_mat = results.get('freq_used'), results.get('tc'), results.get('stim_re_mat'), results.get('stim'), results.get('resp_mat'), results.get('noise_mat'), results.get('CF_mat'), results.get('SNR_mat')
tk, K_mat = results.get('tk'), results.get('K_mat')
figure(4)
plot(tk, K_mat[0,:])
SNR = SNR_mat[1][0,:]
figure(2)
semilogx(freq, SNR, 'g')
print "- SNR(f=" + str(freq[(freq>=9.5) & (freq<=10.5)]) + "): " + str(real(SNR[(freq>=9.5) & (freq<=10)]))
print mean(SNR)
figure(3)
tc, resp_cc, stim_cc, stim_re_cc, noise_cc, CF, VAF = reconstruct_Stimulus(K_mat[0,:], resp_mat[0,:], stim, t)
plot(tc, stim_cc, 'k')
plot(tc, stim_re_cc, 'b')
plot(tc, noise_cc, 'r')
print "- CF=" + str(CF) + " VAF=" + str(VAF)
#show()
#SNR CHECK!
figure(31)
nfft = 2 ** int(ceil(log2(w_length * fs)))
P_nn, freq = csd(noise_cc, noise_cc, Fs = fs, NFFT = nfft, noverlap = nfft / 4, sides = 'twosided')
P_ss, freq = csd(stim_cc, stim_cc, Fs = fs, NFFT = nfft, noverlap = nfft / 4, sides = 'twosided')
SNR = ifftshift(P_ss / P_nn)[0:len(freq)/2] # 10*log10()
freq = ifftshift(freq)[0:len(freq)/2] # only positive frequencies
figure(2)
semilogx(freq, SNR, 'r')
figure(32)
semilogx(freq, SNR, 'r')
figure(33)
pxx, f = psd(input_signal, Fs=1/dt)
clf()
semilogx(f, pxx / pxx[1])
show()
if "hdf5" in test:
filename = 'text.hdf5'
y = []
x = np.arange(0, 1, 0.1)
y.append(x)
y.append(x)
d = {}
d['y'] = y
d['b'] = 'array'
print d
rw_hdf5(filename, d)
d = rw_hdf5(filename)
print d