import sys
sys.path.insert(0, '..')
import time
import matplotlib.pyplot as plt
%matplotlib notebook  
import numpy as np
import scipy.stats

from Configuration import Configuration
from NeuralTract import NeuralTract

conf = Configuration('confNeuralTractSpikes.rmto')

t = np.arange(0.0, conf.simDuration_ms, conf.timeStep_ms)

pools = dict()
pools[0] = NeuralTract(conf, 'CMExt')

tic = time.time()
for i in xrange(0,len(t)-1):
    pools[0].atualizePool(t[i], 1000/12.0, 10)
toc = time.time()
print str(toc - tic) + ' seconds'

pools[0].listSpikes()

"""
Explanation: Simulation of 10000 ms of 400 independent descending commands following a gamma distribution with mean of 142 ms and order 10
End of explanation
"""


plt.figure()
plt.plot(pools[0].poolTerminalSpikes[:, 0],
         pools[0].poolTerminalSpikes[:, 1]+1, '.')
plt.xlabel('t (ms)')
plt.ylabel('Descending Command index')

"""
Explanation: The spike times of all descending commands along the 10000 ms of simulation is shown in Fig. \ref{fig:spikesDesc}.
End of explanation
"""


plt.figure()
plt.plot(pools[0].poolTerminalSpikes[pools[0].poolTerminalSpikes[:, 0]>9000, 0],
         pools[0].poolTerminalSpikes[pools[0].poolTerminalSpikes[:, 0]>9000, 1]+1, '.')
plt.xlabel('t (ms)')
plt.ylabel('Descending Command index')

"""
Explanation: The spike times of all descending commands during the last 1000 ms of the simulation is shown in Fig. \ref{fig:spikesDescLast}.
End of explanation
"""


ISI = np.array([])
for i in xrange(0,len(pools[0].unit)):
    ISI = np.append(ISI, np.diff(np.reshape(np.array(pools[0].unit[i].terminalSpikeTrain), (-1,2))[:,0]))
    

plt.figure()
plt.hist(ISI)
plt.xlabel('ISI (ms)')
plt.ylabel('Counts')

"""
Explanation: The histogram of the interspike intevals of all the descending commands is shown in Fig. \ref{fig:hist}. Note that the peak is in the specified ISI at the beginning of the simulation.
End of explanation
"""


SD = np.std(ISI)
M = np.mean(ISI)
SK = scipy.stats.skew(ISI)
CV = SD / M

print 'ISI Mean = ' + str(M) + ' ms'
print 'ISI Standard deviation = ' + str(SD) + ' ms'
print 'ISI CV = ' + str(CV)

M_FR = 1000.0 / M
SD_FR = np.sqrt((SD**2) * 1000 / (M**3) + 1/6.0 + (SD**4) / (2*M**4) - SK/(3*M**3))


print 'Firing rate mean = ' + str(M_FR) + ' Hz'
print 'Firing rate standard deviation = ' + str(SD_FR) + ' Hz'

CV_FR = SD_FR / M_FR
print 'CV of Firing rate = ' + str(CV_FR)

"""
Explanation: Below different statistics of the interspike intervals and firing rate are obtained.
End of explanation
"""


