<jupyter_start><jupyter_text>Code released on GitHub with the lecture "Analysis of e-phys data" during the neuroscience track core course at IST Austria. https://github.com/michnard/SpikeSortingTutorial<jupyter_code>import numpy as np
import matplotlib.pyplot as plt
import scipy.ndimage as nd
# generate spikes on 4 different tetrodes
# each spike will have different amplitude, spike form, etc
# one spike trace
def spike(ampl, width, bef, aft, sigm):
    w = int(width)
    a=np.zeros(64)
    a[30-w:32-w]=bef *ampl
    a[32]=-ampl
    a[32+w:35+w] = aft * ampl
    return nd.gaussian_filter1d(a,sigm)
plt.plot(spike(3,3,0.01,0.02,2))
# ten cells: get random params
N = 10
np.random.seed(25)
ampls = 2*np.random.rand(N) + 0.5 # this then will be multiplied by
befs = np.random.rand(N)/5
afts = np.random.rand(N)/5
wids = np.random.randint(low=1, high = 4, size=N, dtype=int)
sigs = np.random.randn(N) + 3
# channel amplifs
chans = [np.random.rand(4)*2 + 0.5 for i in range(N)]
# visualize the spike shapes
plt.figure(figsize=(6,6))
cols = ['tab:blue', 'tab:orange', 'tab:green', 'tab:red', 'tab:purple', 'tab:brown', 'tab:pink', 'tab:gray', 'tab:olive', 'tab:cyan']
for i in range(N):
    spikes = [spike(ampls[i]*chans[i][j], wids[i], befs[i],afts[i],sigs[i]) for j in range(4)]
    for j in range(4):
        plt.subplot(N, 4, i*4 + j+1)
        plt.plot(spikes[j],c=cols[i])
        plt.axis('off')
        plt.ylim([np.min(spikes)-0.1, np.max(spikes)+0.1])
# generate 10 random firing fields on linear track 0~100
np.random.seed(123)
cents = np.random.rand(N)
x=np.linspace(0,1,100)
ws = np.random.randint(low=10,high=30,size=N)
frs = np.random.randint(low=5,high=50,size=N) + np.random.rand(N)
tun = np.array([frs[ic]*np.exp(-ws[ic]*(x - c)**2) for ic,c in enumerate(cents)])
for i in range(N):
    plt.plot(tun[i])
# average speed on each part of the track
avg_sp = 20*np.exp(-30*(x - 0.5)**2) + 8
# simulate running - let's say 20 times back and forths
locs = []
pos = 0
for tri in range(20):
    while pos<99:
        # jumps forward with a probability prop to avg_sp - let's say this is every 1ms
        pos = pos + int(np.random.rand()<avg_sp[pos]/1000)
        locs.append(pos)
    # random wait
    for i in range(np.random.randint(1000, 10000)):
        locs.append(pos)
    while pos>0:
        # jumps backward with a probability prop to avg_sp - let's say this is every 1ms
        pos = pos - int(np.random.rand()<avg_sp[pos]/1000)
        locs.append(pos)
    # random wait
    for i in range(np.random.randint(1000, 10000)):
        locs.append(pos)
plt.plot(locs)
plt.xlabel('Time')
plt.ylabel('Location')
# save "behavior"
np.save('locs.npy',locs)
# simulate firing - add spikes to "LFP" trace
lfp=np.zeros([4, len(locs)*20],dtype=float)

spikes = np.array([[spike(ampls[i]*chans[i][j], 
                 wids[i], 
                 befs[i],
                 afts[i],
                 sigs[i]) for j in range(4)] for i in range(N)])
# poisson firing of each cell based on tuning
fir = np.random.poisson(tun[:,locs]/1000)
# number of spikes for each cell
np.sum(fir,1)
# add spikes + noise to the lfp
for i in range(N):
    dels = np.random.randint(1,31,int(np.sum(fir[i])))
    for ii,n in enumerate(np.where(fir[i])[0]):
        ind = n*20 + dels[ii]
        lfp[:, ind-31:ind+33] = spikes[i] + np.random.randn(4,64)*0.001
np.save('fir.npy',fir)
# add oscillations and other, spurious things
osc_fr = 20000 / np.array([1,4,10,40,80,120,250])
ws = osc_fr/20000
osc = np.sum([0.2*np.sin(np.arange(len(lfp[0])) / o) for io,o in enumerate(osc_fr)],0)
from scipy import interpolate
# add background stuff
# noise = np.zeros_like(lfp)
for i in range(4):
    b = np.fft.irfft(np.random.randn(int(len(lfp[0])/8)), n=int(len(lfp[0])/4))
    f = interpolate.interp1d(np.arange(int(len(lfp[0])/4))*8,b,kind='cubic',fill_value='extrapolate')
    b_int = f(np.arange(len(lfp[0])))
    lfp[i] += osc/5 + b_int*25
plt.plot(lfp[0][:5000])
# resize from 0 to 65536 -> it can be saved then as a uint16 (16bit unsigned)
for i in range(4):
    lfp[i] = lfp[i] - np.min(lfp[i])
    lfp[i] = lfp[i] / np.max(lfp[i]) * 65536
lfp = lfp.astype(np.uint16)
# save binary file
binwrite=open('testint.in','wb')
lfp.tofile(binwrite)
binwrite.close()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.17290594291062236, 'fraction_numerical': 0.048198408984557795, 'mean_word_length': 3.723756906077348, 'pattern_counts': {'":': 0, '<': 8, '<?xml version=': 0, '>': 7, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '1', 'dir': 'jupyter-structured-clean-dedup', 'id': '334819', 'license': 'MIT', 'n_tokens_mistral': 1890, 'n_tokens_neox': 1736, 'n_words': 407, 'path': 'Generate_Data.ipynb', 'repo_name': 'michnard/SpikeSortingTutorial', 'size': '4274'}