// spatial configuration see also 3d.hoc

// RANDOM
// ======
// procedures to potentiate a random percent of synapses on a branch
objref prnd, pfile
strdef pfilename
objref psynrnd
// indicate lower and upper bounds for pLTP and size of potentiation
proc random_potentiate() { local n, pLTPl, pLTPh, pot, seed, i, j, pLTP, psynLTP
	if (RDP_check()) { 
	n = numarg()
	pLTPl = $1	// lower bound for LTP probability (0,1)
	pLTPh = $2	// higher bound for LTP probability (0,1)
	pot = $3		// size of potentiation ( >1)
	seed = 83659
	if (n == 4) { seed = $4 }
	sprint(pfilename, "extras/spatial/data/random/potentiations/%s-RIP%d-RDP%d.txt", CELL, RIPNUM, RDPNUM)
	pfile = new File()
	pfile.wopen(pfilename)
	pfile.printf("pLTPl=%3.2f\tpLTPh=%3.2f\tpot=%3.2f\n", pLTPl, pLTPh, pot)
	pfile.printf("logsyn\tbrnch\tnsecsyns\tpLTP\tLTP\n")
	prnd = new Random(seed)
	prnd.uniform(pLTPl, pLTPh)	// probability of potentiation in each section
	psynrnd = new Random(seed + 238)
	psynrnd.uniform(0,1)		// probability of a synapse to potentiate
	for i = 0, NLEAFSECS-1 {
		branch = tree.branchlist.object(secLsorti.x(i))
		pLTP = prnd.repick()	// percent of potentiated synapses in section
		nsecsyns = branch.last_node.numlogsyn	// num of syns on the section
		for j = 0, nsecsyns-1 {
			logsyn = branch.last_node.logsynlist.object(j)
			psyn = 1
			psynLTP = psynrnd.repick()
			if (psynLTP < pLTP) {
				psyn = pot
			}
			logsyn.initial_gmax(logsyn.syn.Egmax * psyn)
			pfile.printf("%d\t%d\t%d\t\t%3.2f\t%3.2f\n", logsyn.id, secLsorti.x(i), nsecsyns, pLTP, psyn)
			criterion_mark(logsyn.id, 6)
		}
	}
	write_GMAX0_records() // records.hoc
	pfile.close()
	} else {
		print "RDP locked for writing!"
	}
}

objref pLTP_vec, pLTPsorti_vec, dist_vec, distsorti_vec
objref tmp_vec
objref tmpi_vec

// BETTER GO TO subgroup_potentiate.hoc !!!!!!!!!!!!!!!!!!
// -----------------------------------------------------------------------------

// potentiate the np most distal (proximal) synapses out of all leafsec synapses
// indicate size of potentiation and number of potentiated synapses
proc disperse_potentiate() { local n, pot, type, np, i, j, psyn
	if (RDP_check()) { 
	n = numarg()
	pot = $1		// size of potentiation ( >1)
	np = $2		// number of potentiated synapses
	type = 1		// type: 1=distal ; 2=proximal
	if (n == 3) { type = $3 }
	sprint(pfilename, "extras/spatial/data/random/potentiations/%s-RIP%d-RDP%d.txt", CELL, RIPNUM, RDPNUM)
	pfile = new File()
	pfile.wopen(pfilename)
	pfile.printf("np=%3.2f\tpot=%3.2f\n", np, pot)
	pfile.printf("logsyn\tbrnch\tnsecsyns\tdist\tLTP\n")
	dist_vec = new Vector()
	tmpi_vec = new Vector()
	distsorti_vec = new Vector()
	pLTP_vec = new Vector()
	for i = 0, NLEAFSECS-1 {
		branch = tree.branchlist.object(secLsorti.x(i))
		nsecsyns = branch.last_node.numlogsyn	// num of syns on the section
		for j = 0, nsecsyns-1 {
			logsyn = branch.last_node.logsynlist.object(j)
			dist_vec.append(abs(logsyn.dist))
		}
	}
	tmpi_vec = dist_vec.sortindex()
	if (type == 1) { tmpi_vec.reverse() }
	distsorti_vec = tmpi_vec.sortindex()
	k = 0
	for i = 0, NLEAFSECS-1 {
		branch = tree.branchlist.object(secLsorti.x(i))
		nsecsyns = branch.last_node.numlogsyn	// num of syns on the section
		for j = 0, nsecsyns-1 {
			logsyn = branch.last_node.logsynlist.object(j)
			psyn = 1
			if (distsorti_vec.x(k) < np) {
				psyn = pot
			}
			pLTP_vec.append(psyn)
			logsyn.initial_gmax(logsyn.syn.Egmax * psyn)
			pfile.printf("%d\t%d\t%d\t\t%3.2f\t%3.2f\n", logsyn.id, secLsorti.x(i), nsecsyns, logsyn.dist, psyn)
			criterion_mark(logsyn.id, 6)
			k += 1
		}
	}
	write_GMAX0_records() // records.hoc
	pfile.close()
	} else {
		print "RDP locked for writing!"
	}
}

np=0	// number of potentiated synapses
proc show_random_potentiate() { local check, i, j, k, id, psyn
	sprint(pfilename, "extras/spatial/data/random/potentiations/%s-RIP%d-RDP%d.txt", CELL, RIPNUM, RDPNUM)
	pfile = new File()
	check = pfile.ropen(pfilename)
	if (check) {
		for i = 0, NLEAFSECS-1 {
			branch = tree.branchlist.object(secLsorti.x(i))
			nsecsyns = branch.last_node.numlogsyn	// num of syns on the section
			for j = 0, nsecsyns-1 {
				id = pfile.scanvar()
				for k = 1, 3 {
					pfile.scanvar()
				}
				psyn = pfile.scanvar()
				if (psyn > 1) {
					np+=1
					criterion_mark(id, 6)
				}
			}
		}
		pfile.close()
	} else {
		print "Could not find ", pfilename
	}
}

// =============================================

// create a new sugroup with all synapses belonging to leafsecs so that they can
// be controlled as a group
proc leafsec_subgroup() { local ok, i
	if(RDP_check()) {
		if (NSUBGS && SUBG.nlogsyns==0) {
			ok=dbox.dialog("This procedure will create a new subgroup!")
			if (ok) {
				for i = 0, NLEAFSECS-1 {
					branch = tree.branchlist.object(secLsorti.x(i))
					nsecsyns = branch.last_node.numlogsyn	// num of syns on the section
					for j = 0, nsecsyns-1 {
						logsyn = branch.last_node.logsynlist.object(j)
						add_logsyn_to_subg(1,logsyn)
					}
				}
			}
		} else { dbox.dialog("First create a new subgroup for this RDP!") }
	} else { dbox.dialog("Current RDP is locked for writing!") }
}

// =============================================
// plot for each potentiated synapse its residual potentiation vs. spatial segregation

objref pvec, p0vec, ssvec, ss0vec
objref ss0pvec, ss0dvec, ppvec, pdvec, ss0bvec, pbvec
objref meanbvec
strdef cnfgs

// indicate RIP and RDP for potentiated and for control simulations
// as well as the synaptic spatial interval dx, and alpha (averaging space constant)
proc residual_random_potentiation_vs_SS() { local i,j,k,l, id, lambda, lseg, l, epsilon, n, norm, x0, x1, y0, y1, dx, alpha
	// step 1: compute initial and residual potentiation
	plotter = plotterlist.object(0)
	p0vec = new Vector()
	pvec = new Vector()
	sprint(cnfgs, "potentiated: RIP%d RDP%d; contrl: RIP%d RDP%d", $1, $2, $3, $4)
	// change to potentiated configuration
	change_configuration($1,$2)  // configurations.hoc
	// change to gmax0
	plotter.gtype = 0
	plotter.ordinate_selection(0)
	for i = 0, NLOGSYNS-1 {
		p0vec.append(logsynlist.object(i).syn.Egmax0)
	}
	// change to gmax
	plotter.gtype = 0
	plotter.ordinate_selection(1)
	for i = 0, NLOGSYNS-1 {
		pvec.append(logsynlist.object(i).syn.Egmax)
	}
	// change to control configuration
	change_configuration($3,$4)  // configurations.hoc
	// change to gmax0
	plotter.gtype = 0
	plotter.ordinate_selection(0)
	for i = 0, NLOGSYNS-1 {
		p0vec.x(i) /= logsynlist.object(i).syn.Egmax0
	}
	// change to gmax
	plotter.gtype = 0
	plotter.ordinate_selection(1)
	for i = 0, NLOGSYNS-1 {
		pvec.x(i) /= logsynlist.object(i).syn.Egmax
	}
	// step 2: compute spatial segregation of p0vec to plot pvec against
	print cnfgs
	ss0vec = new Vector(NLOGSYNS, 0)
	ssvec = new Vector(NLOGSYNS, 0)
	meanbvec = new Vector(NLEAFSECS, 0)
	epsilon = 0.01	// smallest distance weight
	g = plotter.graph_export
	g.erase_all
	dx = $5	// spatial interval between every two synapses
	alpha = $6	// averaging space constant
	n = int(-alpha/dx*log(epsilon))	// space window half-size
//	n = 40
	norm = (1 + exp(-dx/alpha) - 2*exp(-(n+1)*dx/alpha)) / (1 - exp(-dx/alpha)) // normalization
	for i = 0, NLEAFSECS-1 {
		branch = tree.branchlist.object(secLsorti.x(i))
		nsecsyns = branch.last_node.numlogsyn	// num of syns on the section
/*
		// it turns out that a constant space constant (alpha) produces much
		// better result than a passive lambda-dependent one
		branch.last_node.sec.sec {
			lseg = L / nseg
			l = 0
			for(x,0) {
				l += lseg / sqrt(1e+4*diam(x) / (4*g_pas(x)*Ra))
			}
			lambda = L / l
		}
		norm = (1 + exp(-dx/(alpha*lambda)) - 2*exp(-(n+1)*dx/(alpha*lambda))) / (1 - exp(-dx/(alpha*lambda))) // normalization
*/
		for j = 0, nsecsyns-1 {
			logsyn = branch.last_node.logsynlist.object(j)
			for k = j-n, j+n {
				f = exp(-abs(k-j)*dx / alpha)
//				f = exp(-abs(k-j)*dx / (alpha*lambda))
				l = k
				if (l > nsecsyns-1) { // beyond edge of branch - mirror
					l = 2*nsecsyns - k-1
				}
				if (l < 0) { // assume no potentiation outside the section
					ss0vec.x(logsyn.id) += p0vec.x(logsyn.id) / 1 * f
					ssvec.x(logsyn.id) += pvec.x(logsyn.id) / 1 * f
				} else {
					ss0vec.x(logsyn.id) += p0vec.x(logsyn.id) / p0vec.x(branch.last_node.logsynlist.object(l).id) * f
					ssvec.x(logsyn.id) += pvec.x(logsyn.id) / pvec.x(branch.last_node.logsynlist.object(l).id) * f
				}
			}
			meanbvec.x(i) += abs(ss0vec.x(logsyn.id)/norm-1)
		}
		meanbvec.x(i) /= nsecsyns
	}
	ss0vec.div(norm)
	ssvec.div(norm)
	ss0pvec = new Vector()
	ss0dvec = new Vector()
	ss0bvec = new Vector()
	ppvec = new Vector()
	pdvec = new Vector()
	pbvec = new Vector()
	for i = 0, NLOGSYNS-1 {
		logsynlist.object(i).range = ss0vec.x(i)
		if (ss0vec.x(i)) {	// only special branches
			pbvec.append(pvec.x(i))
			ss0bvec.append(ss0vec.x(i))
			if (p0vec.x(i) > 1) {
				ppvec.append(pvec.x(i))
				ss0pvec.append(ss0vec.x(i))
			} else {
				pdvec.append(pvec.x(i))
				ss0dvec.append(ss0vec.x(i))
			}
		}
	}
	x0 = 0.5
	y0 = 0.5
	x1 = 1.5
	y1 = 1.5
	g.size(x0, x1, y0, y1)
	ppvec.mark(g, ss0pvec, "O", 3, 2, 1)
	g.color(2)
	linear_fit(ss0pvec, ppvec)
	sprint(lbl, "r_2 = %3.2f", r_2)
	g.align(0.5, 1)
	g.label(0.3,0.9, lbl)
	g.beginline
	g.line(x0, x0*a+b)
	g.line(x1, x1*a+b)
	g.flush
	print "alpha = ", alpha, "dx = ", dx
	print "potentiated: slope = ", a, "intrcpt = ", b, "r^2 = ", r_2
	ap = a
	bp = b
	r_2p = r_2
	pdvec.mark(g, ss0dvec, "O", 3, 3, 1)
	g.color(3)
	linear_fit(ss0dvec, pdvec)
	sprint(lbl, "r_2 = %3.2f", r_2)
	g.label(0.7,0.2, lbl)
	g.beginline
	g.line(x0, x0*a+b)
	g.line(x1, x1*a+b)
	g.flush
	print "non-potentiated: slope = ", a, "intrcpt = ", b, "r^2 = ", r_2
}

objref electrotonic_vec

proc space_constant() { local i, l, lseg, avg
	electrotonic_vec = new Vector()
	g = plotterlist.object(0).graph_export
	g.erase_all
	avg = 0
	for i = 0, NLEAFSECS-1 {
		branch = tree.branchlist.object(secLsorti.x(i))
		nsecsyns = branch.last_node.numlogsyn	// num of syns on the section
		branch.last_node.sec.sec {
			lseg = L / nseg
			l = 0
			for(x,0) {
				l += lseg / sqrt(1e+4*diam(x) / (4*g_pas(x)*Ra))
			}
			lambda = L / l
			avg += lambda
			electrotonic_vec.append(L/lambda)
//			g.mark(L, lambda, "O", 3, 1, 1)
			print nsecsyns, L, lambda, L/lambda, nsecsyns/lambda
		}
	}
	avg /= NLEAFSECS
	print avg
}

objref xv, yv
proc linear_fit() { local c
	xv = $o1
	yv = $o2
	SS_xx = xv.dot(xv) - xv.size*xv.mean^2
	SS_yy = yv.dot(yv) - yv.size*yv.mean^2
	SS_xy = xv.dot(yv) - xv.size*xv.mean*yv.mean
	a = SS_xy / SS_xx
	b = yv.mean - a*xv.mean
	r_2 = SS_xy^2 / (SS_xx * SS_yy)
}

// loop over alpha values and write a file with a, b, r^2
// indicate RIP and RDP for potentiated and for control simulations
// and dx, alpha min max and step size
proc alpha_loop() { local min, max, step, alpha, dx
	dx = $5
	min = $6
	max = $7
	step = $8
	sprint(pfilename, "extras/spatial/data/random/alphas/%s-RIP%d-RDP%d-RIP%d-RDP%d.txt", CELL, $1, $2, $3, $4)
	pfile = new File()
	pfile.wopen(pfilename)
	pfile.printf("\tpotentiated\t\tothers\n")
	pfile.printf("alpha\tslope\tintrcpt\tr^2\tslope\tintrcpt\tr^2\n")
	for (alpha = min; alpha <= max; alpha += step) {
		residual_random_potentiation_vs_SS($1, $2, $3, $4, dx, alpha)
		pfile.printf("%3.2g\t%3.4g\t%3.4g\t%3.4g\t%3.4g\t%3.4g\t%3.4g\n", alpha, ap, bp, r_2p, a, b, r_2)
	}
	pfile.close()
}

objref alphvec
objref apvec, bpvec, rpvec
objref advec, bdvec, rdvec

// indicate RIP and RDP for potentiated and for control simulations
proc plot_alphas() { local check
	sprint(pfilename, "extras/spatial/data/random/alphas/%s-RIP%d-RDP%d-RIP%d-RDP%d.txt", CELL, $1, $2, $3, $4)
	pfile = new File()
	check = pfile.ropen(pfilename)
	if (check) {
		plotter = plotterlist.object(0)
		g = plotter.graph_export
		g.erase_all
		alphvec = new Vector()
		apvec = new Vector()
		bpvec = new Vector()
		rpvec = new Vector()
		advec = new Vector()
		bdvec = new Vector()
		rdvec = new Vector()
		while(1) {
			if (pfile.eof()) { break }
			alphvec.append(pfile.scanvar())
			apvec.append(pfile.scanvar())
			bpvec.append(pfile.scanvar())
			rpvec.append(pfile.scanvar())
			advec.append(pfile.scanvar())
			bdvec.append(pfile.scanvar())
			rdvec.append(pfile.scanvar())
		}
		rpvec.plot(g, alphvec, 2, 2)
//		apvec.plot(g, alphvec, 2, 1)
		rdvec.plot(g, alphvec, 3, 2)
//		advec.plot(g, alphvec, 3, 1)
	} else { 
		print "file not found"
	}
}
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

// vector saying for each synapse whether it is not on one of the special branches (0)
// whether it is and was potentiated (1) or not potentiated (-1)
// and logsyn id list of potentiated (psynid) and depressed (dsynid) synapses
objref syntype, psynid, dsynid
// select relevant potentiation configuration
proc prepare_syntype() { local check, i, j, k, id, psyn
	syntype = new Vector(NLOGSYNS,0)
	psynid = new Vector()
	dsynid = new Vector()
	sprint(pfilename, "extras/spatial/data/random/potentiations/%s-RIP%d-RDP%d.txt", CELL, RIPNUM, RDPNUM)
	pfile = new File()
	check = pfile.ropen(pfilename)
	if (check) {
		for i = 0, NLEAFSECS-1 {
			branch = tree.branchlist.object(secLsorti.x(i))
			nsecsyns = branch.last_node.numlogsyn	// num of syns on the section
			for j = 0, nsecsyns-1 {
				id = pfile.scanvar()
				for k = 1, 3 {
					pfile.scanvar()
				}
				psyn = pfile.scanvar()
				if (psyn > 1) {
//				if (psyn < 1) {
					syntype.x(id) = 1
					psynid.append(id)
				} else {
					syntype.x(id) = -1
					dsynid.append(id)
				}
			}
		}
		pfile.close()
	} else {
		print "Could not find ", prfilename
	}
}


/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

// =============================================
// compute correlation between gmax profile at t0 and t1
// for each representative branch

objref t0_vec, t1_vec, bcc_vec, bt0_vec, bt1_vec

proc corr_random_potentiate() { local i, j, std0, std1, mu0, mu1, n
	t0_vec = new Vector()
	t1_vec = new Vector()
	bcc_vec = new Vector()
	// step 1: compute initial and residual potentiation
	plotter = plotterlist.object(0)
	p0vec = new Vector()
	pvec = new Vector()
	sprint(cnfgs, "potentiated: RIP%d RDP%d; contrl: RIP%d RDP%d", $1, $2, $3, $4)
	// change to potentiated configuration
	change_configuration($1,$2)  // configurations.hoc
	// change to gmax0
	plotter.gtype = 0
	plotter.ordinate_selection(0)
	for i = 0, NLOGSYNS-1 {
		p0vec.append(logsynlist.object(i).syn.Egmax0)
	}
	// change to gmax
	plotter.gtype = 0
	plotter.ordinate_selection(1)
	for i = 0, NLOGSYNS-1 {
		pvec.append(logsynlist.object(i).syn.Egmax)
	}
	// change to control configuration
	change_configuration($3,$4)  // configurations.hoc
	// change to gmax0
	plotter.gtype = 0
	plotter.ordinate_selection(0)
	for i = 0, NLOGSYNS-1 {
		p0vec.x(i) /= logsynlist.object(i).syn.Egmax0
	}
	// change to gmax
	plotter.gtype = 0
	plotter.ordinate_selection(1)
	for i = 0, NLOGSYNS-1 {
		pvec.x(i) /= logsynlist.object(i).syn.Egmax
	}
	for i = 0, NLEAFSECS-1 {
		bt0_vec = new Vector()
		bt1_vec = new Vector()
		branch = tree.branchlist.object(secLsorti.x(i))
		nsecsyns = branch.last_node.numlogsyn	// num of syns on the section
		for j = 0, nsecsyns-1 {
			logsyn = branch.last_node.logsynlist.object(j)
			bt0_vec.append(p0vec.x(logsyn.id))
			bt1_vec.append(pvec.x(logsyn.id))
		}
		std0 = bt0_vec.stdev
		std1 = bt1_vec.stdev
		mu0 = bt0_vec.mean
		mu1 = bt1_vec.mean
		n = bt0_vec.size
		bt0_vec.sub(mu0)
		bt1_vec.sub(mu1)
		bcc_vec.append(bt0_vec.dot(bt1_vec)/(n*std0*std1))
	}
}


objref nsecsyns_vec, pLTP_vec
// load data about each branch
proc branch_data() { local check, i, j, k, id, psyn, pLTP, max
	sprint(pfilename, "extras/spatial/data/random/potentiations/%s-RIP%d-RDP%d.txt", CELL, RIPNUM, RDPNUM)
	pfile = new File()
	check = pfile.ropen(pfilename)
	if (check) {
		nsecsyns_vec = new Vector()
		pLTP_vec = new Vector()
		dist_vec = new Vector(NLEAFSECS, 0)
		max = 1
		for i = 0, NLEAFSECS-1 {
			branch = tree.branchlist.object(secLsorti.x(i))
			nsecsyns = branch.last_node.numlogsyn	// num of syns on the section
			nsecsyns_vec.append(nsecsyns)
			for j = 0, nsecsyns-1 {
				id = pfile.scanvar()
				logsyn = branch.last_node.logsynlist.object(j)
				for k = 1, 2 {
					pfile.scanvar()
				}
				pLTP = pfile.scanvar()
				psyn = pfile.scanvar()
				if (psyn > max) { max = psyn }
				dist_vec.x(i) += abs(logsyn.dist)
			}
			dist_vec.x(i) /= nsecsyns
			pLTP_vec.append(pLTP)
		}
		pfile.close()
		rp0 = max
//		rp0 = 0.5
		print "spatial dispersion index = ", pLTP_vec.dot(dist_vec)/dist_vec.sum
	} else {
		print "Could not find ", pfilename
	}
}

objref avgp_vec, avgd_vec, err_vec
// compute rp measures for each branch
// specify RIP and RDP for potentiated and control
proc branch_rp() { local check, i, j, id, psyn, avgp, ctp, avgd, ctd
	avgp_vec = new Vector()
	avgd_vec = new Vector()
	pvec = new Vector()
	err_vec = new Vector()
	// change to potentiated
	change_configuration($1,$2)  // configurations.hoc
	prepare_syntype($1,$2)
	branch_data()
	// change to gmax
	plotter.gtype = 0
	plotter.ordinate_selection(1)
	for i = 0, NLOGSYNS-1 {
		logsyn = logsynlist.object(i)
		pvec.append(logsyn.syn.Egmax)
	}
	// change to control
	change_configuration($3,$4)  // configurations.hoc
	// change to gmax
	plotter.gtype = 0
	plotter.ordinate_selection(1)
	for i = 0, NLOGSYNS-1 {
		logsyn = logsynlist.object(i)
		pvec.x(i) /= logsyn.syn.Egmax
	}
	for i = 0, NLEAFSECS-1 {
		branch = tree.branchlist.object(secLsorti.x(i))
		nsecsyns = branch.last_node.numlogsyn	// num of syns on the section
		avgp = 0
		ctp = 0
		avgd = 0
		ctd = 0
		for j = 0, nsecsyns-1 {
			id = branch.last_node.logsynlist.object(j).id
			psyn = pvec.x(id)
			if (syntype.x(id) == 1) { 
				avgp += psyn
				ctp += 1
			} else if (syntype.x(id) == -1) {
				avgd += psyn
				ctd += 1
			}
		}
		if (ctp) { avgp_vec.append(avgp/ctp) }
		if (ctd) { avgd_vec.append(avgd/ctd) }
	}
	g = plotter.graph_export
	g.erase_all
	err_vec = pLTP_vec.c
	err_vec.apply("rp_vs_pLTP")
	err_vec.sub(avgp_vec)
	r_2 = err_vec.var
	r_2 = 1 - r_2 / avgp_vec.var
	avgp_vec.mark(g, pLTP_vec, "O", 3, 2, 1)
	g.color(2)
	sprint(lbl, "r_2 = %3.2f", r_2)
	g.align(0.5, 1)
	g.label(0.3,0.9, lbl)
	theoretical_rp_vs_pLTP(rp0, 2)
//	plot_fit(pLTP_vec, avgp_vec, 2)
//	plot_fit(pLTP_vec, avgd_vec, 3)
}

// theoretical curves
// ===========
// rp0 must be a global variable
func rp_vs_pLTP() { local pLTP, rp
	pLTP = $1
	rp = rp0 / ( pLTP*(rp0 - 1) + 1)
	return rp
}

// rp0 must be a global variable
func rp_vs_N() { local N, rp
	N = $1
	rp = rp0 / ( 1/N*(rp0 - 1) + 1)
	return rp
}

objref rpvec, pLTPvec
// specify rp0 and color
proc theoretical_rp_vs_pLTP() { local c
	rp0 = $1
	c = $2
	g = plotterlist.object(0).graph_export
	rpvec = new Vector()
	pLTPvec = new Vector()
	pLTPvec = new Vector()
	pLTPvec.indgen(0,1,0.1)
	rpvec = pLTPvec.c
	rpvec.apply("rp_vs_pLTP")
	g.color(c)
	rpvec.line(g, pLTPvec)
}

// rp0 must be a global variable
func single_synapse() { local kappa, Vtrg, Esyn, rp, T
	rp = $1
	kappa = 1
	Vtrg = 5
	Esyn = 65
//	T = kappa/(Vtrg - Esyn)^2 *(Vtrg*log(rp0/rp) - Esyn*log((rp0-1)/(rp-1)))
	T = kappa/(Vtrg*(Vtrg-Esyn)) * ((Vtrg-Esyn)*log(rp/rp0) + Esyn*log((1-rp)/(1-rp0)))
	return T
}

objref T_vec, rp_vec
// specify rp0 and color
proc theoretical_single_synapse() { local c
	rp0 = $1	// global variable
	c = $2
	g = plotterlist.object(0).graph_export
	rp_vec = new Vector()
	T_vec = new Vector()
//	rp_vec.indgen(1.001,rp0,0.001)
	rp_vec.indgen(rp0,0.999,0.001)
	T_vec = rp_vec.c
	T_vec.apply("single_synapse")
	g.color(c)
	rp_vec.line(g, T_vec)
}

/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////


// plot residual potentiation at t2 vs t1
objref pt0, pt1, pt2, dt0, dt1, dt2

// indicate RIP and RDP for potentiated and for control simulations for t1 and for t2
proc potentiation_t2_vs_t1() { local i, kp, kd, x0, x1, y0, y1
	plotter = plotterlist.object(0)
	pt0 = new Vector()
	pt1 = new Vector()
	pt2 = new Vector()
	dt0 = new Vector()
	dt1 = new Vector()
	dt2 = new Vector()
	// change to t1 potentiated
	change_configuration($1,$2)  // configurations.hoc
	prepare_syntype($1,$2)
	// change to gmax
	plotter.gtype = 0
	plotter.ordinate_selection(1)
	for i = 0, NLOGSYNS-1 {
		if (syntype.x(i)) {
			logsyn = logsynlist.object(i)
			if (syntype.x(i) == 1) {
				pt1.append(logsyn.syn.Egmax)
				pt0.append(logsyn.syn.Egmax0)
			} else {
				dt1.append(logsyn.syn.Egmax)
				dt0.append(logsyn.syn.Egmax0)
			}
		}
	}
	// change to t1 control
	change_configuration($3,$4)  // configurations.hoc
	// change to gmax
	plotter.gtype = 0
	plotter.ordinate_selection(1)
	kp = 0
	kd = 0
	for i = 0, NLOGSYNS-1 {
		if (syntype.x(i)) {
			logsyn = logsynlist.object(i)
			if (syntype.x(i) == 1) {
				pt0.x(kp) /= logsyn.syn.Egmax0
				pt1.x(kp) /= logsyn.syn.Egmax
				kp += 1
			} else {
				dt0.x(kd) /= logsyn.syn.Egmax0
				dt1.x(kd) /= logsyn.syn.Egmax
				kd += 1
			}
		}
	}
	// change to t2 potentiated
	change_configuration($5,$6)  // configurations.hoc
	// change to gmax
	plotter.gtype = 0
	plotter.ordinate_selection(1)
	for i = 0, NLOGSYNS-1 {
		if (syntype.x(i)) {
			logsyn = logsynlist.object(i)
			if (syntype.x(i) == 1) {
				pt2.append(logsyn.syn.Egmax)
			} else {
				dt2.append(logsyn.syn.Egmax)
			}
		}
	}
	// change to t2 control
	change_configuration($7,$8)  // configurations.hoc
	// change to gmax
	plotter.gtype = 0
	plotter.ordinate_selection(1)
	kp = 0
	kd = 0
	for i = 0, NLOGSYNS-1 {
		if (syntype.x(i)) {
			logsyn = logsynlist.object(i)
			if (syntype.x(i) == 1) {
				pt2.x(kp) /= logsyn.syn.Egmax
				kp += 1
			} else {
				dt2.x(kd) /= logsyn.syn.Egmax
				kd += 1
			}
		}
	}
	// plot
	g.erase_all
	plot_fit(pt1, pt2, 2)
	plot_fit(dt1, dt2, 3)
}


// indicate RIP and RDP for potentiated and for control simulations 
proc potentiation_t1_vs_t0() { local i, kp, kd, x0, x1, y0, y1
	plotter = plotterlist.object(0)
	pt0 = new Vector()
	pt1 = new Vector()
	dt0 = new Vector()
	dt1 = new Vector()
	// change to potentiated
	change_configuration($1,$2)  // configurations.hoc
	prepare_syntype($1,$2)
	// change to gmax
	plotter.gtype = 0
	plotter.ordinate_selection(1)
	for i = 0, NLOGSYNS-1 {
		if (syntype.x(i)) {
			logsyn = logsynlist.object(i)
			if (syntype.x(i) == 1) {
				pt1.append(logsyn.syn.Egmax)
				pt0.append(logsyn.syn.Egmax0)
			} else {
				dt1.append(logsyn.syn.Egmax)
				dt0.append(logsyn.syn.Egmax0)
			}
		}
	}
	// change to control
	change_configuration($3,$4)  // configurations.hoc
	// change to gmax
	plotter.gtype = 0
	plotter.ordinate_selection(1)
	kp = 0
	kd = 0
	for i = 0, NLOGSYNS-1 {
		if (syntype.x(i)) {
			logsyn = logsynlist.object(i)
			if (syntype.x(i) == 1) {
				pt0.x(kp) /= logsyn.syn.Egmax0
				pt1.x(kp) /= logsyn.syn.Egmax
				kp += 1
			} else {
				dt0.x(kd) /= logsyn.syn.Egmax0
				dt1.x(kd) /= logsyn.syn.Egmax
				kd += 1
			}
		}
	}
	// plot
	g.erase_all
	plot_fit(pt0, pt1, 2)
	plot_fit(dt0, dt1, 3)
}

objref X,Y

proc plot_fit() { local x0, x1, y0, y1, color
	n = numarg()
	X = $o1
	Y = $o2
	color = $3
	plotter = plotterlist.object(0)
	g = plotter.graph_export
	x0 = X.min
	y0 = Y.min
	x1 = X.max
	y1 = Y.max
	g.size(x0, x1, y0, y1)
	Y.mark(g, X, "O", 3, color, 1)
	linear_fit(X, Y)
	g.color(color)
	sprint(lbl, "r_2 = %3.2f", r_2)
	g.align(0.5, 1)
	g.label(0.3,0.9, lbl)
	g.beginline
	g.line(x0, x0*a+b)
	g.line(x1, x1*a+b)
	g.flush
}