# -*- 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 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
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 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 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]
return t, ivec
def slugify(value):
"""
Normalizes string, converts to lowercase, removes non-alpha characters,
and converts spaces to hyphens.
From Django's "django/template/defaultfilters.py".
"""
for c in r'[]/\;,><&*:%=+@!#^\'()$| ?^':
value = value.replace(c,'')
return value
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 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__':
pass