#cp sim_mmns_savespikesonly.py sim_fI.py
#cp sim_mmns_savespikesonly.py sim_mmns_sep_savespikesonly.py #Add slowly excited population
from brian2 import *
from pylab import *
import scipy.io
import time
from os.path import exists


tauNeur = 10.0
if len(sys.argv) > 1:
  tauNeur = float(sys.argv[1])
gLeak = 2.0
if len(sys.argv) > 2:
  gLeak = float(sys.argv[2])
thresh = -40
if len(sys.argv) > 3:
  thresh = float(sys.argv[3])
gLeak_addition = '_gLeak'+str(gLeak) if gLeak != 2.0 else ''
thresh_addition = '_thresh'+str(thresh) if thresh != -40 else ''

amps = [100.0*i for i in range(0,11)]

nSps = []
f,axarr = subplots(2,1)

for iamp in range(0,len(amps)):
  stimAmp = amps[iamp]
  if exists('fIs/fI_highamps_tau'+str(tauNeur)+gLeak_addition+thresh_addition+'_amp'+str(stimAmp)+'.mat'):
    A = scipy.io.loadmat('fIs/fI_highamps_tau'+str(tauNeur)+gLeak_addition+thresh_addition+'_amp'+str(stimAmp)+'.mat')
  else:
    print('fIs/fI_highamps_tau'+str(tauNeur)+gLeak_addition+thresh_addition+'_amp'+str(stimAmp)+'.mat not found, simulating...')
    start_scope()

    eqs = '''
dv/dt = (stimulus(t) + g_leak*(e_leak-v))/tau : 1 (unless refractory)
I : 1
tau : second
e_leak : 1
g_leak : 1
'''

    '''
stimulus:responsible for the  signal
'''

    stimList = [0]*10+[stimAmp]*100
    stimulus = TimedArray(array(stimList), dt=100*ms)

    Population = NeuronGroup(1, eqs, threshold='v>'+str(thresh), reset='v = -80', refractory=2*ms, method='rk4')

    # Configuration for Population
    Population.v = -80
    Population.tau = [tauNeur]*ms #[10]*(2*Nperpop)*ms
    Population.g_leak = [gLeak]
    Population.e_leak = [-80]

    tstop = (10+100)*100
    print("tstop = "+str(tstop))

    V_outputPopulation = StateMonitor(Population, 'v', record=True)
    PopulationSpikeMonitor = SpikeMonitor(Population)

    timenow = time.time()
    run(tstop*ms)
    print("Simulation run in "+str(int(time.time()-timenow))+" seconds, amp="+str(stimAmp))

    scipy.io.savemat('fIs/fI_highamps_tau'+str(tauNeur)+gLeak_addition+thresh_addition+'_amp'+str(stimAmp)+'.mat',  {'spikes': [PopulationSpikeMonitor.t/msecond, array(PopulationSpikeMonitor.i)]})