# -*- coding: utf-8 -*-
"""
Created on Wed Jul 27 16:28:09 2016
@author: tomasz, Bartosz, Romain
"""
from __future__ import print_function
import numpy as np
import numba
from scipy import stats
import math
import pylab as plt
#from PoissonPointProcess import SNCP_corr, corr_poisson, CCVF
def generate_poisson_old(freq, tmax):
"""Draw spike times from Poisson distribution"""
t = np.cumsum(stats.expon.rvs(scale=1/freq, size=round(freq * tmax)))
t = t[:t.searchsorted(tmax)]
return t
# this function might be slower than the previous implementation - please check
def generate_poisson(freq, tmax):
return np.sort(np.random.uniform(0,tmax,np.random.poisson(lam = freq*tmax)))
def generate_correlated(mother, r):
"""Generate correlated spike train from a mother spike train
process"""
rn = np.random.rand(len(mother))
return mother[rn<r]
def jitter_spikes(spt, tau):
"""add normaly distributed jitter to spike trains"""
delta_t = stats.expon.rvs(scale=tau, size=len(spt))
# to have symmetric distribution around 0 uncomment the following lines
# sign = (np.random.rand(len(spt)) > 0.5) * 2 - 1
# delta_t = sign * delta_t
return np.sort(spt + delta_t)
def spike_trains_hierarch(nb_comp, nb_syn, rate, tmax, corr_coef_L, corr_coef_G, jtr):
N = rate * tmax
if corr_coef_L == 0. and corr_coef_G == 0.:
spike_trains = [[np.sort(generate_poisson(rate, tmax)) for i in range(nb_syn)] for j in range(nb_comp)]
return spike_trains
if corr_coef_G == 0.:
r_LC = corr_coef_L
M = round(N / r_LC)
freq_mother = float(M) / tmax
mother_loc = [np.sort(generate_poisson(freq_mother, tmax)) for j in range(nb_comp)]
spike_trains = [[jitter_spikes(generate_correlated(mother_loc[j], r_LC), jtr) for i in range(nb_syn)] for j
in range(nb_comp)]
return spike_trains
else:
if corr_coef_L < 0. or corr_coef_L > 1.:
raise ValueError('\n\n *** Local correlation outside [0,1]! ***')
if corr_coef_G < 0. or corr_coef_G > 1.:
raise ValueError('\n\n *** Global correlation outside [0,1]! ***')
if corr_coef_L < corr_coef_G:
raise ValueError('\n\n *** Global correlation cannot be lower than local correlation! ***')
if corr_coef_G == 0.:
"print generation"
r_LC = corr_coef_L
M = round(N / r_LC)
freq_mother = float(M) / tmax
mother_loc = [np.sort(generate_poisson(freq_mother, tmax)) for j in range(nb_comp)]
spike_trains = [[jitter_spikes(generate_correlated(mother_loc[j], r_LC), jtr) for i in range(nb_syn)] for j
in range(nb_comp)]
return spike_trains
r_LC = corr_coef_L
r_GL = corr_coef_G / corr_coef_L
M = round(N / (r_GL * r_LC))
freq_mother = float(M) / tmax
mother_glob = generate_poisson(freq_mother, tmax)
mother_loc = [generate_correlated(mother_glob, r_GL) for j in range(nb_comp)]
spike_trains = [[jitter_spikes(generate_correlated(mother_loc[j], r_LC), jtr) for i in range(nb_syn)] for j in range(nb_comp)]
return spike_trains
# def spike_trains_hierarch_ind_global_bad(nc, ns, rate, tmax, rL, rG, jtr):
#
# N = rate * tmax
#
# M = round(N * nc * ns / (1 + rG * (nc - 1) + rL * (ns - 1)))
#
# fM = M / tmax
#
# ts = generate_poisson(fM, tmax)
#
# C = np.random.randint(nc, size=len(ts))
#
# S = np.random.randint(ns, size=len(ts))
#
# ts_loc = [[[ts[i]], [[C[i], S[i]]]] for i in range(len(ts))]
#
#
# for i in range(len(ts_loc)):
#
# for s in range(ns):
#
# if s != ts_loc[i][1][0][1]:
#
# r = np.random.rand()
#
# if r < rL:
#
# ts_loc[i][1] = ts_loc[i][1] + [[ts_loc[i][1][0][0],s]]
#
# for c in range(nc):
#
# if c != ts_loc[i][1][0][0]:
#
# r = np.random.rand()
#
# if r < rG:
#
# ts_loc[i][1] = ts_loc[i][1] + [[c, np.random.randint(ns)]]
#
#
# ts_loc_a = [[[] for i in range(ns)] for j in range(nc)]
#
# total_number = 0
#
# for i in range(len(ts_loc)):
#
# for j in range(len(ts_loc[i][1])):
# ts_loc_a[ts_loc[i][1][j][0]][ts_loc[i][1][j][1]] = ts_loc_a[ts_loc[i][1][j][0]][ts_loc[i][1][j][1]] + ts_loc[i][0]
#
# total_number = total_number + len(ts_loc[i][0])
#
# print('Number of spikes: ' + str(total_number))
# print('Number of spikes per synapse: ' + str(total_number/(nc * ns)))
# print(N)
#
# return ts_loc_a
def connect_matrix(t_length, nc, ns, rL, rG):
from random import shuffle
C = np.random.randint(nc, size=t_length)
mat0 = np.zeros((t_length, nc * ns))
glob_array = [rG] + [0] * (ns - 1)
loc_array = [1.] + [rL] * (ns - 1)
for i in range(t_length):
shuffle(glob_array)
shuffle(loc_array)
mat0[i, :] = glob_array * nc
mat0[i, (C[i])*ns : (C[i]+1) * ns] = loc_array
return mat0
def generate_correlated_new(mother, vec):
rn = np.random.rand(len(mother))
return mother[rn<vec]
def spike_trains_hierarch_ind_global(nc, ns, rate, tmax, rL, rG, jtr):
N = rate * tmax
M = round(N * nc * ns / (1 + rG * (nc - 1) + rL * (ns - 1)))
fM = M / tmax
ts = generate_poisson(fM, tmax)
cm = connect_matrix(len(ts), nc, ns, rL, rG)
spike_train = [[jitter_spikes(generate_correlated_new(ts, cm[:,j*ns + i]), jtr) for i in range(ns)] for j in range(nc)]
return spike_train
# import time
#
# nc = 3
# ns = 2
#
# tmax = 1
#
# rate = 1
#
# rL = 1.
#
# rG = 0.5
#
# jtr = 0.1
#
#
# t0 = time.time()
#
# sp1 = spike_trains_hierarch_ind_global(nc, ns, rate, tmax, rL, rG, jtr)
#
# print('new new method: ' + str(time.time() - t0))
#
# print(len(sp1[0][0]))
#
#
#
#
# import matplotlib.pyplot as plt
#
# plt.figure()
#
# for i in range(len(sp1)):
#
# for j in range(len(sp1[0])):
#
# plt.plot([i, sp1[i][j][0]],'+')
#
# plt.show()
#
#
#
#
#
# t0 = time.time()
#
# sp2 = spike_trains_hierarch_ind_global_bad(nc, ns, rate, tmax, rL, rG, jtr)
#
# print('new method: ' + str(time.time() - t0))
#
# print(len(sp2[0][0]))
#
#
#
# t0 = time.time()
#
# sp3 = spike_trains_hierarch(nc, ns, rate, tmax, rL, rG, jtr)
#
# print('old method: ' + str(time.time() - t0))
#
# print(len(sp3[0][0]))