load_file("nrngui.hoc")
//use_mcell_ran4(1) // MUST be used on CINECA
cvode_active(0)

a1=0.1

b1=0.7

a1max=a1

b1max=b1

rap_a1=5
rap_b1=(b1/0.71)*5


numaxon=1
numsoma=1
numbasal=58
numapical=78
numtrunk=57
cont=26  // ******  number of obliques

last=12 //********  n-ple size
ns=26  //********  number of trunk comps with noise

sim4conf=102
double pin[cont]
double cond[cont]
double trunko[cont]
double dendo[cont]
double u2[last+3]
double u[last]

weight=80

xopen("geoc91662.hoc")             // geometry file
xopen("fixnseg.hoc")           

Rm = 28000
RmDend = Rm/1
RmSoma = Rm
RmAx = Rm

Cm    = 1
CmSoma= Cm
CmAx  = Cm
CmDend = Cm*1

RaAll= 150
RaSoma=150  
RaAx = 50

Vrest = -65
dt = 0.1
gna =  .025
AXONM = 5
gkdr = 0.01
celsius = 35.0  
KMULT =  0.03
KMULTP = 0.03
ghd=0.00005
nash=0

objref g, b,c, stim, vbox,vecchio,pc, args, args2, r,ennupla
objref  p, s[cont], rsyn[cont], nc[last], sref, blist[numtrunk],aplist,f[2*cont]
strdef dend, trunk


for i=0, numtrunk-1 {blist[i] = new Vector()}



args = new Vector()
args2 = new Vector()


aplist = new Vector(numapical)

forsec "axon" {insert pas e_pas=Vrest g_pas = 1/RmAx Ra=RaAx cm=CmAx}
forsec "soma" {insert pas e_pas=Vrest g_pas = 1/RmSoma Ra=RaSoma cm=CmSoma}
forsec "dendrite"{insert pas e_pas=Vrest g_pas = 1/RmDend Ra=RaAll cm=CmDend}
forsec "user5" {insert pas e_pas=Vrest g_pas = 1/RmDend Ra=RaAll cm=CmDend}

access soma
freq=50
geom_nseg()
tot=0
forall {tot=tot+nseg}

distance()
maxdist=0
forsec "user5" for(x) {if (distance(x)>maxdist) {maxdist=distance(x)}}
print "total # of segments (50Hz): ",tot, "  max path distance: ", maxdist


//*********mapping bifurcations******************


for i=0, numapical-1 apical_dendrite[i] {
	while (!issection("user5.*")) {
	//print "before  ", i, secname()
		sref = new SectionRef()
		access sref.parent
		sprint(dend, secname())
	}
	//print "apical ",i," ",dend
	for k=0, numtrunk-1 user5[k] {
		sprint(trunk,secname())
		x=strcmp(dend, trunk)
		if (x==0) {blist[k].append(i)aplist.x[i]=k}
	}
}

/*****************lettura da file********************/

vecchio = new File()
vecchio.ropen("obliqui91662.txt")

for i=0,cont-1 {
	trunko[i]=vecchio.scanvar()
	dendo[i]=vecchio.scanvar()
	pin[i]=vecchio.scanvar()
	cond[i]=vecchio.scanvar()
}

vecchio.close()



/********************  fine  **********************/


tstop=30

for z=0, cont-1 {
	s[z] = new NetStims(.5)
	s[z].interval=0.2
	s[z].number = 1
	s[z].start=10
	s[z].noise=0
	s[z].seed(987651119)
	rsyn[z] = new Exp2Syn(0.5)
	rsyn[z].tau1 = 0.4
	rsyn[z].tau2 = 1
	rsyn.e=0
	
}

user5[13] {		// *************cell dependent
	stim= new IClamp(0)
	stim.amp=0
	stim.dur=tstop
	stim.del=0
}

forsec "axon" {   
                insert nax gbar_nax=gna * AXONM	sh_nax=nash
                insert kdr gkdrbar_kdr=gkdr
                insert kap gkabar_kap = KMULTP*1
}

forsec "soma" {   
		insert hd ghdbar_hd=ghd	vhalfl_hd=-73
                insert na3 ar_na3=1 sh_na3=nash gbar_na3=gna
                insert kdr gkdrbar_kdr=gkdr
                insert kap gkabar_kap = KMULTP
		insert ds
}

for i=0, numbasal-1 dendrite[i] {
		insert hd ghdbar_hd=ghd	vhalfl_hd=-73
                insert na3 ar_na3=1 gbar_na3=gna sh_na3=nash
                insert kdr gkdrbar_kdr=gkdr
                insert kap gkabar_kap = KMULTP
}
                

forsec "user5" {
	insert hd ghdbar_hd=ghd
        insert na3 ar_na3=1 gbar_na3=gna sh_na3=nash
        insert kdr gkdrbar_kdr=gkdr
	insert kap gkabar_kap=0
	insert kad gkabar_kad=0
	for (x,0) { 
		xdist = distance(x)
               	ghdbar_hd(x) = ghd*(1+3*xdist/100)
               		if (xdist > 100){
						vhalfl_hd=-81
                    	gkabar_kad(x) = KMULT*(1+xdist/100)
               		} else {
						vhalfl_hd=-73
                    	gkabar_kap(x) = KMULTP*(1+xdist/100)
        			}
	}
}

for i=0, numapical-1 apical_dendrite[i] {
	insert hd 
        insert na3 ar_na3=1 gbar_na3=gna sh_na3=nash
        insert kdr gkdrbar_kdr=gkdr
	insert kap
	insert kad
		gkabar_kad = 1*user5[aplist.x[i]].gkabar_kad(1)
		gkabar_kap = 1*user5[aplist.x[i]].gkabar_kap(1)
		vhalfl_hd = user5[aplist.x[i]].vhalfl_hd
       	ghdbar_hd = user5[aplist.x[i]].ghdbar_hd(1)
}

proc init() {

	t=0
    forall {
    	v=Vrest
    	if (ismembrane("nax") || ismembrane("na3")) {ena=55}
       	if (ismembrane("kdr") || ismembrane("kap") || ismembrane("kad")) {ek=-90}
       	if (ismembrane("hd") ) {ehd_hd=-30}
	}
	finitialize(Vrest)
    fcurrent()

    forall {
		for (x) {
			if (ismembrane("na3")||ismembrane("nax")){e_pas(x)=v(x)+(ina(x)+ik(x))/g_pas(x)}
			if (ismembrane("hd")) {e_pas(x)=e_pas(x)+i_hd(x)/g_pas(x)}
		}
	}
	cvode.re_init()
	cvode.event(tstop)
	access soma
	//g.begin()
	flag=0
	
}


proc advance() {
	fadvance()
	if (vmax_ds(.5)>0) {t=tstop flag=1}
	//g.plot(t)
	//g.flush()
	//p.flush()
	doNotify()
}


func superrun() {local id, numprove, c, k, seme1, seme2, seme3

	a1=$o2.x[0]
	b1=$o2.x[1]
	stim.amp=$o2.x[2]
    
	k=$3
	c=$4
	seme1=(k*sim4conf+c)*31530
	seme2=(k*sim4conf+c)*22210
	seme3=(k*sim4conf+c)*18230

	print "seme1", seme1
	print "seme2", seme2
	print "seme3", seme3

	for kk=0, ns user5[kk] {  // cell-dependent
	f[kk] = new Gfluct2(0.5)
	f[kk].g_e0 = a1
	f[kk].g_i0 = b1
	f[kk].std_e = a1/rap_a1
	f[kk].std_i = b1/rap_b1
	f[kk].new_seed(seme1)
	print " noise #",kk, " at ",secname()
	}
	
	

	
	print k,c
	objref ennupla
	numprove=0
	id = hoc_ac_
	ennupla= new Vector()
	flag=0
	// definiamo ennupla 
	for c=0, cont-1 {
	ennupla.append(c)
	}
	aux=cont
	// print "seme", seme3
	r= new Random()
	r.MCellRan4(seme3)
	for c=0, last-1 {
	ind=r.discunif(0, aux-1)
	cond[ennupla.x[ind]]=80
	aux=aux-1
	print ennupla.x(ind)
	ennupla.remove(ind)
	}
	print "\n"
	for c=0, cont-last-1 {
	cond[ennupla.x(c)]=0
	print ennupla.x(c)

	}
	ennupla.resize(0)

	r= new Random(seme2)
        r.MCellRan4(seme2)

	while(0==0){ 
	numprove=numprove+1
	for c=0, cont-1 {
	ennupla.append(c)
	}
	aux=cont
	// print " \n"
	for c=0, last-1 {
	ind=r.discunif(0, aux-1)
	//print "indice ",ind
	u[c]=ennupla.x[ind]
	u2[c]=dendo[u[c]]
	nc[c] = new NetCon(s[u[c]],rsyn[u[c]],0,0,cond[u[c]]*1.e-3)
	apical_dendrite[u2[c]] rsyn[u[c]].loc(0.5)
	ennupla.remove(ind)
	aux=aux-1
	 
	//print u[c]
	}
	// print " \n"
	ennupla.resize(0)	
			
	// print f[2].g_e0
	run()
	//print numprove*0.025-0.5 ," s "
	if (flag==1) { print "ha sparato dopo",numprove,"prove" spikec=1}
	if(numprove>3600){print "piu di 3600 prove" flag=1 spikec=0}
	
 	}
	pc.post(id, $o2,numprove,seme1,seme2,seme3)
	
	return id
	
}


pc = new ParallelContext()
pc.runworker()


for yy=0, last-1 {u[yy]=0}

order=0
conto=0

print "Begin Simulation"
print "Contesto = ",stim.amp

objref aa1, bb1, aamp 
aa1 = new Vector()
bb1 = new Vector()
aamp = new Vector()
//ennupla = new Vector()

aa1.append(0.02)
//aa1.append(0.25,0.38,0.5,1,2.5,5,7.5)
bb1.append(0.02)

//bb1.append(0.5,0.75,1,2,5,10,15)
//aamp.append(0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5)
contex = 0
proc loop() { local k, first, ind, c
	k=$1
	c=$2
	args.resize(0)      
	args.append(aa1.x[nsim]*a1max)
	args.append(bb1.x[nsim]*b1max)
	args.append(contex)
	pc.submit("superrun", order, args,k,c)
	order=order+1
		
}

for nsim=0, aa1.size()-1 {
	for z=1,sim4conf{
	loop(nsim,z)	    
	}
}

objref outfile, outfile2
outfile = new File()
outfile2 = new File() 

strdef name, name2


total=0

while ((id = pc.working) != 0) {
	pc.take(id)
	args2 = pc.upkvec
	m  = pc.upkscalar
	seme1=pc.upkscalar
	seme2=pc.upkscalar
	seme3=pc.upkscalar
	
	//	args2.printf()
	total = total +1

//	printf("%d configurazioni ", m)
		//for y=0, last-1 {printf(" %g ", args2.x[y])}
//	printf("noise ex %g inh %g context %g semerum %d semeennupla %d \n",args2.x[0]/a1max,args2.x[1]/b1max, args2.x[2],seme1,seme2)
		//write n-uples
	
	
	
	
sprint(name,"le%dple-%dobl-%gecc-%gini-b%g-%d.txt",last,cont,args2.x[0]/a1max,args2.x[1]/b1max,b1max,sim4conf)
	     
sprint(name2,"le%dple-%dobl-%gecc-%gini-b%g-solotempi-%d.txt",last,cont,args2.x[0]/a1max,args2.x[1]/b1max,b1max,sim4conf)	

	outfile.aopen(name)
	outfile.printf(" %d configurazioni ",m)
	args2.printf(outfile, " %g ")	
	outfile.printf(" seme r seme e semel %d %d %d",seme1,seme2,seme3)
	outfile.close()
	
        outfile2.aopen(name2)	
        outfile2.printf(" %d   	%d \n",m,seme1)
        outfile2.close()
}

pc.done