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

from brian2 import *

defaultclock.dt = 0.01*msecond

tau = 20 * msecond        # constant de temps 
Vt = -55 * mvolt          # seuil pour spiker
Vr = -75 * mvolt          # potentiel de reset
Vs = -65 * mvolt          # potentiel de repos
I0 =  15 * mvolt
tmax = 500 * msecond

# Notre groupe contient un neurone intègre et tire 
eqs = '''
dV/dt = (-(V-Vs) + I)/tau : volt
I : volt
'''
G = NeuronGroup(N=1,
                model=eqs,
                threshold='V >= Vt',
                reset='V = Vr',
                method='linear')

# On monitore ses spikes et son potentiel
spikemon = SpikeMonitor(G)
statemon = StateMonitor(G, 'V', record=0)

# On initialise le potentiel aléatoirement entre le reset et le seuil
G.V = Vr + (Vt - Vr) * rand(len(G))
G.I = I0

# On lance la simulation
run(tmax)

# On affiche combien de spikes ont été émis
print("Nombre de spikes émis : ", spikemon.num_spikes)
print("Liste des temps de spikes :")
print(spikemon.spike_trains())
print(', '.join(["{} ".format(ts) for ts in spikemon.spike_trains()[0]]))
print("Temps inter-spikes (s.):")
print(', '.join(["{} ".format(t1-t0) for t0, t1 in zip(spikemon.spike_trains()[0][:-1], spikemon.spike_trains()[0][1:])]))

predicted_isi = tau * log((Vr - (Vs + I0))/(Vt - (Vs + I0)))
print("Intervalle interspike {}".format(predicted_isi))
print("Fréquence de décharge : {}".format(1.0/predicted_isi))

# On trace nos mesures
figure()
subplot(211)
for i, ts in zip(spikemon.i, spikemon.t):
    plot([ts/msecond, ts/msecond], [i, i+1], 'k')
ylim([-1,2])
xlabel("Temps (ms)")
ylabel(u"Numéro du neurone")

subplot(212)
plot(statemon.t, statemon.V.T)
ylabel('Potentiel (V)')
xlabel('Temps(s.)')

savefig('ex1.png')

show()