{load_file("nrngui.hoc")}

initialization_tstart = startsw()

strdef modelname
strdef cmd
strdef fileName
objref fileObj

fileObj = new File()

rdSeed = 1
Nmc = 100
con_rad = 5	
beta = 0.1
isolat =0

frstate	= 4 // 0, low firing rate; 1, medium firing rate; 2, high firing rate

if (frstate == 0) {
	fr_pf = 1
	fr_in = 1
} else if (frstate == 1) {
	fr_pf = 6
	fr_in = 2
} else if (frstate == 2) {
	fr_pf = 8
	fr_in = 2
} else if (frstate == 3) {
	fr_pf = 10
	fr_in = 2
} else {
	fr_pf = 12
	fr_in = 2
}

condition = 0

connectivity = 1 

celsius = 34

tstop = 30000
dt = 0.02
step_per_ms= 1/dt

v_init = -70


NsynI = 10
Econ = 0.0005
Icon = 0.001/NsynI 

NcontE = 2
NcontI = 5 
Ncontbs = 4*NsynI
Ncontst = 6*NsynI

NGrc = 4000
Nst = 30*0.6
Nbs = 4

NsynE = NGrc*NcontE
NsynIst = Nst*Ncontst
NsynIbs = Nbs*Ncontbs

rateE = 0.2700*fr_pf
rateIst = 14.4000*fr_in
rateIbs = 14.4000*fr_in

xopen("pkj.hoc")

objref tempvec
tempvec = new Vector()
{tempvec.append(Nmc)}

{Ncells2save = tempvec.min()}

//=================== creating neuron ================================

objref pkj_TTC
objref sl //synaptic locations list

//=================== random variables ================================

objref rds1,rds2,rds3,rds4
{rds1 = new Random(1000*rdSeed)}
{rds2 = new Random(2000*rdSeed)}
{rds1.uniform(0,1)}
{rds2.uniform(0,1)}

objref conMat,cons,cont,vecN,switchEdge,target,conInd,switchInd,switchInd2,si
conMat = new Matrix(Nmc,Nmc) 
cons = new Matrix(Nmc,con_rad)
cont = new Matrix(Nmc,con_rad)
vecN = new Vector()
vecN.indgen(0,Nmc-1,1)
for (i=0;i<con_rad;i+=1) {
	cons.setcol(i,vecN)
	cont.setcol(i,vecN.c.add(i+1))
}
for(i=0;i<Nmc;i+=1){
	for(j=0;j<con_rad;j+=1){
		cont.x[i][j] = (cont.x[i][j])%Nmc
		conMat.x[cons.x[i][1]][cont.x[i][j]] = 1
	}
}
conMat.add(conMat.c.transpose())


switchEdge = new Vector(2*con_rad)
switchInd = new Vector()
switchInd2= new Vector()
target = new Vector(Nmc)
conInd = new Vector()

for (i = 0;i<Nmc;i+=1) {
	switchEdge.setrand(rds1)
	target.setrand(rds2)
	target.x(i) = 0
	switchInd = switchEdge.c.indvwhere("<",beta)
	switchInd2 = switchEdge.c.indvwhere(">=",beta)
	conInd = conMat.getcol(i).c.indvwhere(">",0)
	if (switchInd.size()*conInd.size()) {
		for (j = 0;j<switchInd2.size();j+=1) {
			target.x(conInd.x(switchInd2.x(j))) = 0
		}
		si = target.sortindex.reverse()
		for (k = 0;k<switchInd.size();k+=1) {
			conMat.x[conInd.x(switchInd.x[k])][i]=0
		}
		for (k = 0;k<switchInd.size();k+=1) {
			conMat.x[si.x[k]][i]=1
		}
	}
}
if (isolat) {
conMat.zero()
con_rad = 0
}

//==================== presynaptic spike trains ====================

objref preTrainList

preTrainList = new List()

for(i=0;i<Nmc;i+=1){
	{rds2 = new Random(1000*rdSeed+i)}
	{rds3 = new Random(1000*rdSeed+i)}
	{rds4 = new Random(1000*rdSeed+2*Nmc+i)}
	{rds2.negexp(1/rateE)}
	{rds3.negexp(1/rateIst)}
	{rds4.negexp(1/rateIbs)}

	{preTrainList.append(new List())}
	
	for(i2=0;i2<(NGrc+Nst+Nbs);i2+=1){
		pst=0 
		if (i2<NGrc) {
			for j =0,NcontE-1 {
				preTrainList.o[i].append(new Vector())
			}
			while (pst<tstop) {
				pst+= 1000*rds2.repick()
				for j = 0,NcontE-1 {preTrainList.o[i].o[preTrainList.o[i].count()-1-j].append(pst)}
			}
		} else if  (i2<(NGrc+Nst)){
			for j =0,Ncontst-1 {
				preTrainList.o[i].append(new Vector())
			}
			while (pst<tstop) {
				pst+= 1000*rds3.repick()
				for j = 0,Ncontst-1 {preTrainList.o[i].o[preTrainList.o[i].count()-1-j].append(pst)}
			}
		} else {
			for j =0,Ncontbs-1 {
				preTrainList.o[i].append(new Vector())
			}
			while (pst<tstop) {
				pst+= 1000*rds4.repick()
				for j = 0,Ncontbs-1 {preTrainList.o[i].o[preTrainList.o[i].count()-1-j].append(pst)}
			}		
		}
	}
		
}
	

{load_file("netparmpi.hoc")}
objref epnm

epnm = new ParallelNetManager(Nmc)
{epnm.round_robin()}

for(i=0;i<Nmc;i+=1){
	if (epnm.gid_exists(i)) {
		pkj_TTC = new pkj_neuron() 
		epnm.register_cell(i,pkj_TTC)
		epnm.pc.gid2cell(i).initRand(1000*rdSeed+i)
		epnm.pc.gid2cell(i).setnetworkparameters(Econ,Icon,NsynE,NsynIst,NsynIbs,NcontI)
		epnm.pc.gid2cell(i).distributeSyn()
		epnm.pc.gid2cell(i).setpretrains(preTrainList.o[i])
	}
}

//========================== microcircuit connections =============================================

objref syninds

for(i=0;i<Nmc;i+=1){
	if (epnm.gid_exists(i)) {
		epnm.pc.gid2cell(i).insertMCcons(conMat.getcol(i))  

	}
}

syninds = new Vector()
for(i=0;i<Nmc;i+=1){
	syninds.append(NsynE+NsynIst+NsynIbs)
}


for(i=0;i<Nmc;i+=1){
	for(j=0;j<Nmc;j+=1){
		if (conMat.x[j][i] != 0){  
			for(jj=0;jj<NcontI;jj+=1){
				epnm.nc_append(j,i,syninds.x[i],1e-3/NcontI,1.5)  
				syninds.x[i] +=1   
			}
		}
	}
}


for(i=0;i<Nmc;i+=1){
	if (epnm.gid_exists(i)) {
		{epnm.pc.gid2cell(i).queuePreTrains()} 
	}1
}


strdef outDir

if (frstate==0) {
	modelname = "Purkinje_lowfr"
} else if (frstate==1) {
	modelname = "Purkinje_midfr"
} else {
	modelname = "Purkinje_highfr"
}

sprint(outDir,"simdata/%s/N%d_c%d_cond%d_seed%d",modelname,Nmc,connectivity,condition,rdSeed)

sprint(cmd,"system(\"mkdir -p %s\")",outDir)
execute(cmd)

if(epnm.gid_exists(0)){
	sprint(fileName,"%s/conmat.txt",outDir)
	fileObj.wopen(fileName)
	conMat.fprint(0,fileObj)
	fileObj.close()
}


//============================================================================
//==================== recording settings ==========================


objref vSomaList
objref apcvecList, apcList

apcvecList = new List()
apcList = new List()
vSomaList = new List()

// record Vsoma
for(i=0;i<Ncells2save;i+=1){
	if (epnm.gid_exists(i)) {
		vSomaList.append(new Vector())
		access epnm.pc.gid2cell(i).somaA
		vSomaList.o[vSomaList.count()-1].record(&v(0.5),dt)
	}
}

//record APCount
for(i=0;i<Nmc;i+=1){
	if (epnm.gid_exists(i)) {
		access epnm.pc.gid2cell(i).somaA
		apcList.append(new APCount(0.5))
		apcvecList.append(new Vector())
		apcList.o[apcList.count()-1].thresh= -30
		apcList.o[apcList.count()-1].record(apcvecList.o[apcList.count()-1])
	}
}

{epnm.set_maxstep(100)}
finitialize(-65)            //not stdinit() unless input dt and step_per_ms= 1/dt again!

if (epnm.gid_exists(0)) {
	print "\n"
	sim_tstart = startsw()
	initializationtime = (sim_tstart-initialization_tstart)/3600
	print "Initialization completed. Initialization took ", initializationtime, " hours\n"
	print "Starting simulation\n"
	print "\n"
}

{epnm.psolve(tstop)}

if (epnm.gid_exists(0)) {
	simruntime = (startsw() - sim_tstart)/3600
	print "Simulation took ", simruntime, " hours\n"
}

i2=0
for(i=0;i<Ncells2save;i+=1){
	if (epnm.gid_exists(i)) {
		tempvec = new Vector()
        tempvec = vSomaList.o[i2].c
		sprint(fileName,"%s/vTraceN%d.txt",outDir,i)
		fileObj.wopen(fileName)
		tempvec.printf(fileObj,"%2.6e\n")
		fileObj.close()
		i2+=1
	}
}

i2=0
for(i=0;i<Nmc;i+=1){
	if (epnm.gid_exists(i)) {
		tempvec = new Vector()
        tempvec = apcvecList.o[i2].c
		sprint(fileName,"%s/spikeTimesN%d.txt",outDir,i)
		fileObj.wopen(fileName)
		tempvec.printf(fileObj,"%2.6e\n")
		fileObj.close()
        
        tempvec = new Vector()   
        Nspike =  apcList.o[i2].n
        tempvec = apcvecList.o[i2].c(1,Nspike-1).sub(apcvecList.o[i2].c(0,Nspike-2))
        cv = tempvec.stdev()/tempvec.mean()
        printf("the mean ISI of cell%d is %f\n",i,tempvec.mean)
        printf("the ISI CV of cell%d is %f\n",i,cv)
        
		i2+=1
	}
}

{epnm.pc.runworker()}
{epnm.pc.done()}