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

gLeak = 4.0
thresh = 0.0
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 = []

area = 'PFC'
if len(sys.argv) > 1:
  area = sys.argv[1]

unpicklefile = open('taus_cortical_subjectwise_'+area+'.sav', 'rb')
unpickledlist = pickle.load(unpicklefile,encoding='bytes')
unpicklefile.close()
tau_interpolated,AUCdata,tau_interpolated_pop,isubjs_HC,isubjs_SCZ = unpickledlist[0],unpickledlist[1],unpickledlist[2],unpickledlist[3],unpickledlist[4]
tauNeur = tau_interpolated[:]

for iamp in range(0,len(amps)):
  stimAmp = amps[iamp]
  if exists('fIs/fI_fittedtaus_'+area+gLeak_addition+thresh_addition+'_amp'+str(stimAmp)+'.mat'):
    A = scipy.io.loadmat('fIs/fI_fittedtaus_'+area+gLeak_addition+thresh_addition+'_amp'+str(stimAmp)+'.mat')
  else:
    from brian2 import *

    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(len(tauNeur), 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_fittedtaus_'+area+gLeak_addition+thresh_addition+'_amp'+str(stimAmp)+'.mat',  {'spikes': [PopulationSpikeMonitor.t/msecond, array(PopulationSpikeMonitor.i)]})