/* --- Initialize Variables --- */		//in NEURON variables are global by default
strdef basepath, esynfullpath, isynfullpath, syndir
strdef esynfile, isynfile, syndir
objref LGMD, Str
objref esyn, isyn
objref syn, sref
objref synvec[1]
objref savdata
objref t_rec, rvec[1]
objref rec_matrix

objref simsecs
strdef filename, savename, datadir
strdef tmpstr, strtmp, sect, filename
strdef cmd, cwd

t0=0
tstart=0
n=0
ts = 0
count=0
val=0
save_model_view=0
e4AP=0.60
eZD=0.98
eXE=0.98


/* --- Functions and Procedures for simulations --- */

func isclass() {
// boolean = isclass(object, class type)
   classname($o1, strtmp)
   return strcmp(strtmp, $s2) == 0
}

proc add4AP() {
	forall {
		if (ismembrane("KD_ca3")) gmax_KD_ca3 = gmax_KD_ca3*(1-e4AP)
	}
	v_init = -62
}
proc wash4AP() {
	forall {
		if (ismembrane("KD_ca3")) gmax_KD_ca3 = gmax_KD_ca3/(1-e4AP)
	}
	v_init = -65
}
proc addZD7288() {
	forall {
		if (ismembrane("hcn")) gmax_hcn = gmax_hcn*(1-eZD)
	}
	v_init = -71
}
proc washZD7288() {
	forall {
		if (ismembrane("hcn")) gmax_hcn = gmax_hcn/(1-eZD)
	}
	v_init = -65
}


proc MakeStringList() { local i,count localobj sl
// MakeStringList( list, string1, string2, ...)
	
	count = numarg()
	sl = $o1
	for i=2,count {
		sl.append(new String($si))
	}
	
}


proc MakeSecList() { local i,count localobj sl, sref
//MakeSecList( sectionlist, char list)
// create a new sectionlist from an input list of section names. uses greedy regexp ('' matches all)
	count = numarg()
	sl = $o1
	for i=2,count {
		forall ifsec $si sl.append()
	}
}

proc MechList() { local i localobj ml, mt
//MechList( string list )
// creates a list of strings with the names of all membrane mechanisms in cas
	ml = $o1
	mt = new MechanismType(0)

	for i=0,mt.count()-1 {
		mt.select(i)
		mt.selected(strtmp)
		ml.append(new String(strtmp))
	}
}


func Rm() { local i,flag localobj gls
//Rm( mechanisms, Rm/Gm flag ) measure membrane resistivity ohm*cm^2 or conductivity S/cm^2
	
	flag = 0
	RmGtot=0

	if (numarg() > 0) {
		if (argtype(1) == 1) {
			gls = $o1
		} else if (argtype(1)==2) {
			gls = new List()
			MakeStringList(gls,$s1)
		} else if (argtype(1)==0) {
			gls = new List()
			getglist(gls)
		}
	} else {
		gls = new List()
		getglist(gls)
	}
	
	if (numarg() > 1) flag = $2
	
	for i=0,gls.count()-1 {
		if (ismembrane(gls.o(i).s)) {
			sprint( strtmp, "RmGtot += g_%s", gls.o(i).s)
			//if (verbosity > 3) printf("%s", strtmp)
			execute(strtmp)
		}
	}
	if (flag ==1) {
		return RmGtot
	} else {
		if (RmGtot == 0) return 1e30 else return 1e-3/RmGtot
	}
}

func meanRm() { local argn,i,flag, Asum,Gsum localobj gls, sl
//Rm = meanRm( sections, mechanisms, Rm/Gm flag)

	flag = 0
	Gsum=0
	Asum = 0
	argn = numarg()
	
	if (verbosity > 3) for i=1,argn print argtype(i)

	if (argn >= 1) {
		if (argtype(1)==2)	{
			sl = new SectionList()
			forall ifsec $s1 sl.append()
			if (verbosity > 2) sl.printnames()
		} else if (argtype(1)==1)	sl = $o1
	}
	

	if (argn >= 2) {
		if (argtype(2)==1) {
			gls = $o2
		} else if (argtype(2)==2) {
			gls = new List()
			MakeStringList(gls,$s2)
		} else if (argtype(2)==0) {
			gls = new List()
			getglist(gls)
		}
	} else {
		gls = new List()
		getglist(gls)
	}
	
	if ((argn >= 3) && (argtype(3)==0) ) flag = $3

	if (object_id(sl) != 0) {
		forsec sl for (x,0) {
			Asum += area(x)
			Gsum += Rm(gls,1)*area(x)
		}
	} else {
		forall for (x,0) {
			Asum += area(x)
			Gsum += Rm(gls,1)*area(x)
		}
	}

	if (flag ==1) {
		return Gsum/Asum
	} else {
		if (Gsum == 0) return 1e30 else return Asum*1e-3/Gsum
	}
}


proc SynFileRead() { local i,flag,Ca,ts,gm,isp localobj inFile, synlist
// SynFileRead( synlist, filename, Ca )
// generates a list of synapse objects from text file "filename".
// If Ca is false creates AlphSynapses; if true creates AlphSynapseCa (see mod file)

	count=0
	syntau=0
	syng=0
	syne=0
	spont_rate=0
	gm=-1
	Ca=0
	
	synlist=$o1
	synlist.remove_all()
	filename = $s2
	if (numarg() > 2) Ca=$3
	inFile = new File(filename)
	flag = inFile.ropen()
	if (flag==0) {
		sprint( filename, "%s%s", syndir, filename)
		inFile = new File(filename)
		flag = inFile.ropen()
		if (flag==0) {
			sprint( filename, "%s/%s", syndir, filename)
			inFile = new File(filename)
			flag = inFile.ropen()
		}
		if (flag==0) {
			printf("Error: File not found. SynFileRead")
			return
		}
	}

	flag = inFile.gets(tmpstr)
	num = sscanf(tmpstr, "Time of Collision: %f", &t0)
	if (verbosity > 1) printf("t0 = %f\n", t0)
	flag = inFile.gets(tmpstr)
	num = sscanf(tmpstr, "Visual Stimulation Start Time: %f", &tstart)
	if (verbosity > 1) printf("start time = %f\n", tstart)
	flag = inFile.gets(tmpstr)
	num = sscanf(tmpstr, "Number of %s Synapses: %i", tmpstr, &count)
	if (verbosity > 1) printf("# of synapses = %g\n", count)
	flag = inFile.gets(tmpstr)
	num = sscanf(tmpstr, "Time Constant: %f", &syntau)
	if (verbosity > 1) printf("synapse tau = %f\n", syntau)
	flag = inFile.gets(tmpstr)
	num = sscanf(tmpstr, "Max Conductance: %f", &syng)
	if (verbosity > 1) printf("synapse gmax = %f\n", syng)
	flag = inFile.gets(tmpstr)
	num = sscanf(tmpstr, "Reversal Potential: %f", &syne)
	if (verbosity > 1) printf("synapse reversal = %f\n", syne)
	flag = inFile.gets(tmpstr)
	num = sscanf(tmpstr, "Kinematic parameters: %s", strtmp)
	if (verbosity > 1) printf("Kinematic Parameters = %s\n", strtmp)
	flag = inFile.gets(tmpstr)
	num = sscanf(tmpstr, "Spontaneous rate: %f", &spont_rate)
	if (verbosity > 1) printf("Spontaneous rate = %f\n", spont_rate)
	
	if (issplit == 1) {
		isp = 1
		stopPar()
	}
	for n=0,count-1 {
		flag = inFile.gets(tmpstr)
		num = sscanf( tmpstr, "%s %f %f", sect, &ts, &gm)
		if ((verbosity > 1) && (n==0)) printf("section: %s, time:%f, gmax: %f\n", sect, ts, gm)
		if (flag != -1) {
			if ( Ca == 1 ) { sprint(strtmp, "%s syn = new AlphaSynapseCa(0.5)", sect)
			} else sprint(strtmp, "%s syn = new AlphaSynapse(0.5)", sect)
			execute(strtmp)
			syn.onset=ts
			syn.tau=syntau
			if (gm==-1) {
				syn.gmax=syng
			} else syn.gmax=gm
			syn.e=syne
			synlist.append(syn)

		} else if (verbosity > 0) printf("failed to read synapse")
	}

	inFile.close()
	objref syn
	if (isp==1) startPar()
}



proc adjparams() {local mult	localobj tmplist
// adjparams(synlist, synapse parameter, new value, multiply)

	mult = 0
	tmplist = $o1
	val = $3
	if (numarg()>3) mult = $4
	
	for i=0,tmplist.count()-1 {
		
		syn = tmplist.o(i)
		if (mult == 0) {
			sprint(tmpstr, "syn.%s = val", $s2)
		} else sprint(tmpstr, "syn.%s = syn.%s*val", $s2, $s2)
		execute(tmpstr)
		
	}

}


proc findobj() {local p0,p1,x,rng	localobj out, in, ms
// findobj(subset, full list, parameter, value)
// findobj(subset, full list, parameter, min value, max value)
//
// finds objects within List that has parameter equal to value or between min and max vales
// returns a List with found objects

	out = $o1
	in = $o2
	tmpstr = $s3
	p0 = $4
	rng = 0
	
	if (numarg()>4) {
		p1 = $5
		rng = 1
	}
	
	classname(in.o(0), strtmp)
	ms = new MechanismStandard(strtmp)
	
	for i=0,in.count()-1 {
		
		ms.in(in.o(i))
		x = ms.get(tmpstr)
		if (rng == 0) {
			if (x == p0)   out.append(in.o(i))
		} else {
			if (( x >= p0) && (x <= p1))  out.append(in.o(i))
		}
	}

}

proc genbursts() {local ns,k,a,scale	localobj sub
// genbursts( val, rescale )
// val is multiple of change within burst periods
// scale is multiple of percent of change to cancel out in non-burst periods (0-1)
// used to produce onset responses to new coarse pixels matching the data. smells like defeat

	scale=1
	a=2
	if (numarg() > 0) a=$1
	if (numarg() > 1) scale=$2
	sub = new List()
	ns = esyn.count()
	
	findobj( sub, esyn, "onset", 2040, 2140)
	findobj( sub, esyn, "onset", 2590, 2640)
	findobj( sub, esyn, "onset", 4160, 4210)
	findobj( sub, esyn, "onset", 4620, 4650)
	findobj( sub, esyn, "onset", 4730, 4755)
	findobj( sub, esyn, "onset", 4880, 4900)
	//findobj( sub, esyn, "onset", 4860, 4880)
	findobj( sub, esyn, "onset", 5040, 5060)
	
	k = sub.count()
	
	if (scale>0)  adjparams(esyn, "gmax", (ns-k*scale*(a-1))/ns, 1)
	adjparams(sub,"gmax", a, 1)
}


func totalarea() { local sum
// area = totalarea( )
// calculate total surface area of all compartments

  finitialize()
  sum = 0
  forall for (x,0) sum += area(x)
  return sum
}


func sectionarea() { local sum localobj sl
// area = sectionarea( sectionlist )
// calculate total surface area of compartments within sectionlist

	finitialize()
	sum = 0

	if (argtype(1)==2)	{
		sl = new SectionList()
		forall ifsec $s1 sl.append()
		if (verbosity > 2) sl.printnames()
	} else if (argtype(1)==1)	sl = $o1

	forsec sl for (x,0) sum += area(x)
	return sum
}


proc FindEnds() { local i,count localobj sl, sref
// FineEnds( sectionlist, string )
// finds the end compartment of every branch within compartments specified by ifsec string.
// if no input argument finds branch ends for entire neuron
// 
	count = numarg()
	sl = $o1
	tmpstr=""
	if (count > 1) tmpstr=$s2
	
	forall {
		ifsec tmpstr {
			sref = new SectionRef()
			if (sref.nchild==0) {
				sl.append()
			}
		}
	}
}

proc Zratio_g() { local zf,grel,up,pwr,max,g,mt localobj zz, chan, sl
// Zratio_g( section list, init section, chan, grel, frequency, increase/decrease, pwr, max )
// 
// sets conductance within sections based on the ratio of impedance between that section 
// and a reference location

zf=100
grel=1e-5
up=1
pwr=1
max=1

finitialize(v_init)

chan = new String()
if (argtype(1)==2)	{
	sl = new SectionList()
	forall ifsec $s1 sl.append()
	if (verbosity > 2) sl.printnames()
} else if (argtype(1)==1)	sl = $o1

if (numarg()>1) {
	if (argtype(2)==2)	{
		sprint(tmpstr, "%s sref = new SectionRef()", $s2)
		if (verbosity > 2) printf("%s\n", tmpstr)
		execute(tmpstr)
	}
	if (argtype(2)==1)	sref = $o2
} else soma[0] sref = new SectionRef()
if (numarg()>2) {
	if (argtype(3)==2)	chan.s = $s3
	if (argtype(3)==1)	chan = $o3
	if (verbosity > 2) printf("%s\n", chan.s)
}
if (numarg()>3) grel = $4
if (numarg()>4) zf = $5
if (numarg()>5) up = $6
if (numarg()>6) pwr = $7
if (numarg()>7) max = $8

if (issplit) {
	mt = 1
	stopPar()
}
	if (verbosity > 2) printf("grel = %g\tzf = %g\tup = %g\tpwr = %g\tmax = %g\n", grel,zf,up,pwr,max)

	zz = new Impedance()
	sref.sec { zz.loc(0.5) }
	zz.compute(zf)
	forsec sl {
		if (ismembrane(chan.s)) {
			
			if (grel == 0) {
				sprint(tmpstr, "val = gmax_%s", chan.s)
				execute(tmpstr)
			} else val = grel

			if (up==0) {
				g = val*zz.ratio(0.5)^pwr
				if (g>max) g=max
				sprint(tmpstr, "gmax_%s = %g", chan.s, g)
			}
			if (up>0) {
				g = val/zz.ratio(0.5)^pwr
				if (g>max) g=max
				sprint(tmpstr, "gmax_%s = %g", chan.s, g)
			}
			execute(tmpstr)
		}
	}
	
	if (mt==1) startPar()
}

proc visualstim() { local i, n, count, newesyn, newisyn, ims,sca, sv
// visualstim( esynfile, isynfile, section list, savename,tstop)
// visualstim( 0, 0, section list, ...)		use already loaded synapse files

	n = numarg()
	if (argtype(1)==2)	{
		newesyn=1
		esynfile = $s1
		esyn = new List()
	} else newesyn=0	// if excitatory synapse files already loaded 
	if (argtype(2)==2)	{
		newisyn=1
		isynfile = $s2
		isyn = new List()
	} else newisyn=0	// if inhibitory synapse files already loaded 
	simsecs = $o3
	if (n > 3) {
		savename = $s4
		sv = 1
	} else sv=0
	
	count=0
	sca=0
	forsec simsecs { count+=1 }

	t_rec = new Vector()
	t_rec.record(&t)

	objref rvec[count]
	i=0
	forsec simsecs {
		rvec[i] = new Vector()
		rvec[i].record(&v(0.5))
		i=i+1
	}

	if (verbosity > 1) printf("loading excitation ... \n")
	if (synCa) sca=1
	if (newesyn) SynFileRead(esyn,esynfile,sca)
	if (verbosity > 1) printf("done.\n")

	if (verbosity > 1) printf("loading inhibition ... \n")
	if (newisyn) SynFileRead(isyn,isynfile)
	if (verbosity > 1) printf("done.\n")

	if (n > 4) {
		tstop = $5
	} else tstop = t0+250

	if (save_model_view) {
		if (verbosity > 0) printf("saving param file ... \n")
		ims = issplit
		if (ims) stopPar()
		LGMD = new ModelView(0)
		if (verbosity > 1) printf("ModelView created ... \n")
		sprint(strtmp, "%s_prm.txt", savename )
		LGMD.text(strtmp)
		LGMD.destroy()
		objref LGMD
		if (ims) startPar()
	}
	
	if (verbosity > 0) printf("running simulation ... \n")
	
	finitialize(v_init)
	run()

	if (sv==1) {
		if (verbosity > 0) printf("saving data ... \n")
	
		sprint(filename, "%s", savename )
		savdata = new File()
		savdata.wopen(filename)
		if (verbosity > 1) printf("opened file %s for writing ... \n", filename)

		rec_matrix = new Matrix()
		rec_matrix.resize(t_rec.size(),count+1)
		rec_matrix.setcol(0,t_rec)
		for (i=1; i<=count; i=i+1) {
			rec_matrix.setcol(i,rvec[i-1])
		}
	
		savdata.printf("Visual stimulus simulation\nt0=%g\ttstart=%g\nExcitation: %s\nInhibition: %s\n",t0,tstart,esynfile,isynfile)
		rec_matrix.fprint(savdata,"%-1.8g ")
		savdata.printf( "time " )
		forsec simsecs {savdata.printf( "%s ", secname() )}
		savdata.printf( "\n" )
		savdata.close()
	}

	if (verbosity > 0) {
		printf("visual stimuli finished\n")
		if (sv==1) printf("Data was saved to %s\n", filename) else printf("Data was not saved\n")
	}
}