# Plot an alpha beta table generated by VoltageDepTabChannel etc.

# Copyright 2007 John L Baker. All rights reserved.
# This software is provided AS IS under the terms of the Open Source
# MIT License. See http://www.opensource.org/licenses/mit-license.php.

# Example invocation:
# plot_alpha_beta("test-rallpack-n.txt",vlim=c(-90,40))

plot_alpha_beta<-function(
	file,				# File name with alpha-beta values
	vlim=c(-Inf,Inf),		# Range of voltages to plot (in mV)
	plotFit=TRUE,		# Plot fitted energy barrier results
	drive="C:",			# Drive (or file system) containing data	
	dirPath=file.path(drive,"BNSF_1_0","Data","Default")) # Data directory
{
	# Read alpha-bet values in standard units (volts, seconds, etc.)
	AB<-read.table(file.path(dirPath,file),header=FALSE,sep=',');

	if (ncol(AB)>3) {
		AB<-AB[,1:3];  # discard any unused columns
	}
	if (ncol(AB)>=3) {
		colnames(AB)<-c("v","alpha","beta");
		AB$tau<-1/(AB$alpha+AB$beta);
		AB$xinf<-AB$alpha/(AB$alpha+AB$beta);
	}
	else {
		stop("Unexpected number of columns read");
	}

	# Rescale alpha and beta to 1/ms and also	
	# rescale voltage to mV and tau to ms.
	AB$alpha<-AB$alpha/1000;
	AB$beta<-AB$beta/1000;
	AB$v<-AB$v*1000;
	AB$tau<-AB$tau*1000;

	# Get maximum and minimum tau values for later use
	tauMaxIdx<-which.max(AB$tau);
	tauMaxV<-AB$v[tauMaxIdx];
	tauMin<-min(AB$tau);
	tauMax<-AB$tau[tauMaxIdx];

	# Locate half activation point and estimate
	# parameters for a single energy barrier model.
	# Note that this fitting method can be of limited
	# accuracy since it only looks at values for a few
	# different (but key) voltage values.

	# Note also that there is an assumption that no gating
	# variable exponent was already included in the alpha-beta
	# values (see VoltageDepTabChannel::xinfExponent).
	
	nAB<-nrow(AB);
	xinf.monotonic<-TRUE;

	if (	(AB$xinf[1]<AB$xinf[nAB] && any(AB$xinf[1:(nAB-1)]>AB$xinf[2:nAB])) ||
		(AB$xinf[1]>AB$xinf[nAB] && any(AB$xinf[1:(nAB-1)]<AB$xinf[2:nAB])) ) {
		xinf.monotonic<-FALSE;
		warning("X.infinity values are not monotonic");
	}
	else if (plotFit) {


		# The energy barrier model fit here is of the form:
		#
		# xinf = 1/(1+exp(-(v-vhalf)/k))
		# tau = tmin+1/rate*1/(exp(g*(v-vhalf)/k)+(exp(-(1-g)*(v-vhalf)/k)))
		#
		# where rate is value set to ensure a specified maximum for tau.
		#
		# Except for the tmin correction factor (ala Borg-Graham)
		# an equivalent formulation is:
		#
		# alpha = k0*exp(zeta*gamma*(v-vhalf)*F/(R*T))
		# beta =  k0*exp(-zeta*(1-gamma)*(v-vhalf)*F/(R*T))
		#
		# where k0 (ie rate) is also set to ensure a maximum tau value.

		dv<-(AB$v[nAB]-AB$v[1])/(nAB-1);
		vhalf<-approx(AB$xinf,AB$v,0.5,rule=2)$y;
		k<-(2*dv)/diff(approx(AB$v,AB$xinf,c(vhalf-dv,vhalf+dv))$y)/4;

		g<-1/(1+exp((tauMaxV-vhalf)/k));
		if (tauMax-tauMin>1e-6) {
			rate<-1/((exp(g*(tauMaxV-vhalf)/k)+exp(-(1-g)*(tauMaxV-vhalf)/k))*(tauMax-tauMin));
		}
		else {
			rate<-Inf;
		}

		cat("\nEstimated energy barrier model fit");
		cat("\n half-activation voltage (vhalf) =",vhalf,"(mV)");
		cat("\n response voltage sensitivity (k) =",k,"(mV)");
		cat("\n voltage sensor relative location (gamma) =",g);
		cat("\n minimum tau (tauMin) =",tauMin,"(ms)");
		cat("\n maximum tau (tauMax) =",tauMax,"(ms)");
		cat("\n rate =",rate,"(1/ms)\n");
	}

	# Plot each value in turn as a function of voltage
	# after selecting the range requested.
	ABmap<-AB$v>=vlim[1] & AB$v<=vlim[2];
	AB<-AB[ABmap,];

	# Get energy barrier estimates for xinf and tau if needed
	if (plotFit & xinf.monotonic) {
		eb.xinf<-1/(1+exp(-(AB$v-vhalf)/k));
		if (is.finite(rate)) {
			eb.tau<-tauMin+1/(rate*(exp(g*(AB$v-vhalf)/k)+exp(-(1-g)*(AB$v-vhalf)/k)));
		}
		else {
			eb.tau<-rep(tauMax,length(AB$v));
		}
	}

	# Plot alpha, beta, xinf, and tau all on one screen
	layout(matrix(c(1,2,3,4),2,2,byrow=TRUE))
	
	plot(AB$v,AB$alpha,type='l',col='black',
		xlab="Voltage (mV)",
		ylab="Alpha (1/ms)");

	plot(AB$v,AB$beta,type='l',col='black',
		xlab="Voltage (mV)",
		ylab="Beta (1/ms)");

	plot(AB$v,AB$xinf,type='l',col='black',
		ylim=c(0,1.5),
		xlab="Voltage (mV)",
		ylab="X.infinity");
	if (plotFit & xinf.monotonic) {
		lines(AB$v,eb.xinf,col='blue',lty=2);
		legend("topright",c("Xinf","EB fit"),
			lty=c(1,2),col=c("black","blue"))
	}
	
	plot(AB$v,AB$tau,type='l',col='black',
		ylim=c(0,1.5*tauMax),
		xlab="Voltage (mV)",
		ylab="Tau (ms)");
	if (plotFit & xinf.monotonic) {
		lines(AB$v,eb.tau,col='blue',lty=2);
		legend("topright",c("Tau","EB fit"),
			lty=c(1,2),col=c("black","blue"))
	}

	# Reset layout
	layout(matrix(1,1,1))
}