/*
WHAT THIS DOES

It allows you to check the phase lag between two similar vectors
(usually voltage spikes, but other vectors may work).
It will run the simulation until steady state is reached.
The phase lag and period (in ms) will be displayed as the run takes place.

You can customize...
a) Interval:	How often it checks the phase lag.
b) Count:	How many equal phase lag values before it is considered steady state.
c) Threshold:	It will only check the values above the threshold.
***Note: Interval should be large enough to allow for 2 spikes (for each signal).
	 A safe rule is Interval > 3*Period.

HOW TO USE

If using this tool individually (not via extra.hoc),
place this file along with other model hoc files.
In the oc> prompt, type...
load_file("phaselag.hoc")

Find the PhaseLag window in 'Tools' -> 'Miscellaneous'


 -- Leo Ng
    June 15, 2004
*/


begintemplate PhaseLag

// recA and recB record values of vectors A and B, respectively.
// thres stores the values of recA or recB above the threshold value.
// slope stores the derivatives of signal A or signal B.
// peakA[2] and peakB[2] store when the last 2 peaks of signal A and B occured, respectively.
// diff[1] store the last phase lag and diff[0] store the second last phase lag.
public Start, interval, count, threshold, recA, recB, period, phase
external tstop, run, continuerun
objref this, recA, recB, thres, slope, scA, scB
strdef sig, cmd, info
double peakA[2], peakB[2], diff[2]

// Sets the default values and creates the vectors.
proc init() {
	period = 0
	phase = -1
	count = 10
	interval = 300
	threshold = 30

	recA = new Vector(interval/dt)
	recB = new Vector(interval/dt)
	thres = new Vector()
	slope = new Vector()
	scA = new SymChooser("Choose Vector A")
	scB = new SymChooser("Choose Vector B")

	CreatePanel()
}

// Creates the GUI.
proc CreatePanel() {
	xpanel("PhaseLag")
	xpvalue("Interval length (ms)", &interval, 1)
	xpvalue("Steady State Count", &count, 1)
	xpvalue("Spike Threshold (units of vector)", &threshold, 1)
	xvalue("Period (ms)", "period", 2)
	xvalue("Phase (ms)", "phase", 2)
	xvarlabel(info)
	xvarlabel("")
	xbutton("Start", "Start()")
	xpanel()
}

// Finds where the last two peaks occurred.
// Returns 1 if successful, 0 if failed.
func FindPeak() { local i, j
	thres.clear()
	slope.clear()

	thres.indvwhere($o1, ">", threshold)
	if (thres.size() == 0) return 0		// Threshold too high.
	slope.deriv($o1, dt)
	$&2[1] = 0
	$&2[0] = 0

	// Gets the last peak.
	// If slope changes sign (+ to -), we have a peak.
	for (i = thres.size() - 1; i >= 0; i -= 1)  {
		if ((slope.x[thres.x[i]] < 0) && (slope.x[thres.x[i - 1]] >= 0)) {
			$&2[1] = thres.x[i] * dt
			break
		}
	}

	// Get the second last peak.
	for (j = i - 1; j >= 0; j -= 1)  {
		if ((slope.x[thres.x[j]] < 0) && (slope.x[thres.x[j - 1]] >= 0)) {
			$&2[0] = thres.x[j] * dt
			break
		}
	}

	if (($&2[1] == 0) || ($&2[0] == 0)) return 0 else return 1	// Can't find peaks.
}

// Finds the phase lag from the two peaks.
// Phase should be +ve.
proc CalcLag() { local tmp1, tmp2
	period = peakA[1] - peakA[0]
	tmp1 = abs(peakA[1] - peakB[1])
	if (tmp1 > period) tmp1 -= period
	tmp2 = abs(peakB[1] - peakA[0])
	if (tmp2 > period) tmp2 -= period
	if (tmp1 == tmp2) tmp2 = period - tmp2
	if (tmp1 < tmp2) diff[1] = tmp1	else diff[1] = tmp2
	phase = diff[1]		// Takes the smaller one.
}

// Keeps the simulation running until steady state is reached.
// Steady state is reached when the same phase values occurred 'count' times.
proc CheckSteady() { local tmpCount, tmpA, tmpB
	tmpCount = count - 1
	while (tmpCount != 0) {
		recA.clear()
		recB.clear()
		continuerun(t + interval)
		tmpA = FindPeak(recA, &peakA)
		tmpB = FindPeak(recB, &peakB)
		if (tmpA == 1 && tmpB == 1) {
			diff[0] = diff[1]
			CalcLag()
			if (diff[0] == diff[1]) tmpCount -= 1 else tmpCount = count - 1
		} else {
			info = "Error: Interval too small or Threshold too high!!!"
			return
		}
	}
	info = "Steady State reached."
}

// Loads the vector choosing menu.
// Starts the run.
proc Start() {
	// Opens the SymChooser menu.
	// Stores the chosen variable in 'sig'.
	// Tells recA and recB to record 'sig'.
	while (!scA.run()) print "Must choose something for Vector A!"
	scA.text(sig)
	sprint(cmd, "%s.recA.record(&%s)", this, sig)
	execute(cmd)
	while (!scB.run()) print "Must choose something Vector B!"
	scB.text(sig)
	sprint(cmd, "%s.recB.record(&%s)", this, sig)
	execute(cmd)

	info = "Waiting for steady state..."
	tstop = interval
	diff[0] = -2
	diff[1] = -1
	recA.clear()
	recB.clear()

	run()	
	FindPeak(recA, &peakA)
	FindPeak(recB, &peakB)
	CalcLag()
	CheckSteady()
}

endtemplate PhaseLag

// Adds window to NEURON's main menu.
objref phaseLag
nrnmainmenu_.miscellaneous_add("PhaseLag", "phaseLag = new PhaseLag()")