# -*- coding: utf-8 -*-

# Not yet functionning with brian2 !!

from brian2 import *

tau = 0.2 * msecond
tmax = 800 * ms

# Tonic spiking
a_count = 0.02
b_count = 0.2
c_count = -50 * mV
d_count = 2 * mV

# Phasic bursting
a_onset = 0.02
b_onset = 0.25
c_onset = -55 * mV
d_onset = 0.05 * mV
I0 = 0.6 * mV

eqs= '''
dv/dt = (0.04*10**3 * volt**-1 * v**2 + 5 * v + 140*mvolt - u + I)/tau : volt
du/dt = a * (b * v - u)/tau : volt
I : volt
a : 1
b : 1
c : volt
d : volt
'''

eqs_count = '''
dv/dt = (0.04*10**3 * volt**-1 * v**2 + 5 * v + 140*mvolt - u + ge + gi)/tau : volt
du/dt = a * (b * v - u)/tau : volt
dge/dt = -ge/(15*ms) : volt
dgi/dt = -gi/(15*ms) : volt
a : 1
b : 1
c : volt
d : volt
'''

reset= '''
v=c
u=u+d
'''

G_onset = NeuronGroup(N=1, model=eqs,
              threshold='v >= 30*mV', reset=reset)

G_count = NeuronGroup(N = 2*10 , model = eqs_count, threshold ='v >= 30*mV', reset=reset)

# Définition de la séquence d'entrée
I_values = [0 * mV]
I_times = [0*ms]
for i in range(12):
    I_values.append(I0)
    I_values.append(0)
    t0 = 20*ms + i * 60 * ms + (2 * rand() - 1.0) * 20 * ms 
    I_times.append(t0)
    I_times.append(t0+20*ms)
print(I_times)
G_onset.I = TimedArray(I_values, times=array(I_times)) 

# On initialise le potentiel au potentiel de reset
G_onset.v = -65*mV
G_onset.u = b_onset * G_onset.v
G_onset.a = a_onset
G_onset.b = b_onset
G_onset.c = c_onset
G_onset.d = d_onset

# Définition du groupe pour compter
G_count.v = -65*mV
G_count.u = b_count * G_count.v
G_count.a = a_count
G_count.b = b_count
G_count.c = c_count
G_count.d = d_count

# On connecte le G_onset au premier neurone de G_count
S = Synapses(G_onset, G_count, 'ge', pre='V_post += w',connect=True)
for i in range(1,10):
    S.ge[0,2*i] = 0.25 * mV
S.ge[0,0] = 0.75 * mV

C_count_exci = Connection(G_count, G_count, 'ge', delay=2*ms)
C_count_inhi = Connection(G_count, G_count, 'gi', delay=2*ms)
for i in range(10):
    C_count_exci[2*i,2*i+1] = 1.5 * mV
    C_count_exci[2*i+1,2*i] = 1.5 * mV

for i in range(9):
    C_count_exci[2*i+1, 2*(i+1)] = 0.5 * mV
    C_count_inhi[2*(i+1)+1, 2*i] = -3.0 * mV

for i in range(1,9):
    C_count_inhi[2*i+1, 0] = -3.0 * mV

# On ferme la boucle
C_count_exci[19, 0] = 0.5* mV
C_count_inhi[1, 18] = -3.0* mV

# On monitore ses spikes et son potentiel
vmon_onset = StateMonitor(G_onset, 'v', record=0)
inputmon_onset = StateMonitor(G_onset, 'I', record=0)

spikemon_count = SpikeMonitor(G_count)
vmon_count = StateMonitor(G_count, 'v', record=True)


run(tmax)

# Détecteur d'onset
figure()
subplot(411)
title('Input potential')
inputmon_onset.plot()
ylim([-I0, 2*I0])
ylabel('Potential (V)')
xlabel('Time(s)')
legend(['I'])

subplot(412)
vmon_onset.plot()
ylabel('Potential (V)')
xlabel('Time(s)')
legend(['V'])

subplot(413)
raster_plot(spikemon_count)
xlim(0*ms, tmax * 10**3)

subplot(414)
vmon_count.plot()

savefig('../spike_counting.png')

show()