################################################################
# function codes -- sort of like m-files.  we just have to put
# this in once

# comp will give you the smoothing spline fit
# to get smoothing parameter values, type summary(*)

align<-function(data,alignment='justified'){
	for (w in levels(data$word)) {
		w.ntoken <- max(subset(data,word==w)$token)

		if(alignment=='right')
		{
			for (i in 1:w.ntoken)
			{
				temp <- subset(data,word==w & token==i);
				temp$Time <- temp$Time - temp$Time[dim(temp)[1]];
				data[data$word==w & data$token==i,]$Time <- temp$Time
			}
		} else { if(alignment=='left')
		{
			for (i in 1:w.ntoken)
			{
				temp <- subset(data,word==w & token==i);
				temp$Time <- temp$Time - temp$Time[1];
				data[data$word==w & data$token==i,]$Time <- temp$Time
			}
		} else { if(alignment=='center')
		{
			for (i in 1:w.ntoken)
			{
				temp <- subset(data,word==w & token==i);
				temp$Time <- temp$Time - (temp$Time[dim(temp)[1]] + temp$Time[1])/2;
				data[data$word==w & data$token==i,]$Time <- temp$Time
			}
		} else { if(alignment=='justified')
		{
			for (i in 1:w.ntoken)
			{
				temp <- subset(data,word==w & token==i);
				temp$Time <- (temp$Time - temp$Time[1])/(temp$Time[dim(temp)[1]] - temp$Time[1]);
				data[data$word==w & data$token==i,]$Time <- temp$Time
			}
		}}}}

	}
	data
}

compareformants<-function(w1,w2,data,yrange=c(0,4000),bw=FALSE,weight=3,output='screen',BCI=FALSE,f1=TRUE,f2=TRUE,f3=TRUE,csv=FALSE){
	require(assist)
	options(memory=1000000000)

	w1w2<-rbind(subset(data,word == w1),subset(data,word == w2));
	a<-levels(w1w2$word)
	word3<-a[a!=w1 & a!=w2]
	levels(w1w2$word)<-list(word1=c(w1,word3),word2=w2)

	if(f1)
	{
		f1model<-ssr(F1~word*Time,rk=list(cubic(Time),rk.prod(cubic(Time),shrink1(word))),data=w1w2,scale=T)
		print("The following summary is for F1");
		print(summary(f1model));
	}
	else
	{
		f1model<-NULL;
	}
	if(f2)
	{
		f2model<-ssr(F2~word*Time,rk=list(cubic(Time),rk.prod(cubic(Time),shrink1(word))),data=w1w2,scale=T)
		print("The following summary is for F2");
		print(summary(f2model));
	}
	else
	{
		f2model<-NULL;
	}
	if(f3)
	{
		f3model<-ssr(F3~word*Time,rk=list(cubic(Time),rk.prod(cubic(Time),shrink1(word))),data=w1w2,scale=T)
		print("The following summary is for F3");
		print(summary(f3model));
	}
	else
	{
		f3model<-NULL;
	}

	if(BCI==FALSE)
	{
		if(output=='pdf')
		{
			pdf(paste(w1,' ',w2,'.pdf',sep=''));
			comp.plotformants(f1model,f2model,f3model,w1,w2,yrange=yrange,bw,weight,f1,f2,f3,csv=csv)
			dev.off()
		} else if(output=='emf')
		{
			win.metafile(paste(w1,' ',w2,'.emf',sep=''));
			comp.plotformants(f1model,f2model,f3model,w1,w2,yrange=yrange,bw,weight,f1,f2,f3,csv=csv)
			dev.off()
		}
		else
		{
			comp.plotformants(f1model,f2model,f3model,w1,w2,yrange=yrange,bw,weight,f1,f2,f3,csv=csv)
		}
	}
	else
	{
		if(output=='pdf')
		{
			pdf(paste(w1,' ',w2,'.pdf',sep=''));
			comp.plotbci(f1model,f2model,f3model,w1,w2,yrange=yrange,bw,weight,f1,f2,f3)
			dev.off()
		} else if(output=='emf')
		{
			win.metafile(paste(w1,' ',w2,'.emf',sep=''));
			comp.plotbci(f1model,f2model,f3model,w1,w2,yrange=yrange,bw,weight,f1,f2,f3)
			dev.off()
		}
		else
		{
			comp.plotbci(f1model,f2model,f3model,w1,w2,yrange=yrange,bw,weight,f1,f2,f3)
		}
	}

}

comp.plotformants<-function(f1model,f2model,f3model,word.1,word.2,yrange,bw,weight,f1,f2,f3,csv=FALSE){
	if(bw==TRUE){
		w1.col <- 1;
		w2.col <- 3;
	} else
	{
		w1.col <- 2;
		w2.col <- 3;
	}
	
	par(mfrow=c(1,1))

	isplotted <- FALSE;

	# F1
	if(f1)
	{
		w1<-f1model$data[f1model$data$word=="word1",]
		w2<-f1model$data[f1model$data$word=="word2",]
		n1<-nrow(w1)
		ntot<-nrow(w1)+nrow(w2)
		f1model.pred<-predict(f1model)
		f1model.w1<-data.frame(fit=f1model.pred$fit[1:n1],pstd=f1model.pred$pstd[1:n1])
		f1model.w2<-data.frame(fit=f1model.pred$fit[n1+1:ntot],pstd=f1model.pred$pstd[n1+1:ntot])
		order1<-order(w1$Time)
		order2<-order(w2$Time)
		plot(f1model$data$Time,f1model$data$F1,type="n",xlab="Time (s)",ylab="Frequency (Hz)",ylim=yrange)
		isplotted<-TRUE;
		lines(f1model$data$Time,f1model$data$F1,type="n")
		lines(w1$Time[order1],f1model.w1$fit[order1],col=w1.col,lwd=weight)
		lines(w1$Time[order1],f1model.w1$fit[order1]+3*f1model.w1$pstd[order1],col=w1.col,lwd=1,lty=3)
		lines(w1$Time[order1],f1model.w1$fit[order1]-3*f1model.w1$pstd[order1],col=w1.col,lwd=1,lty=3)
		lines(w2$Time[order2],f1model.w2$fit[order2],col=w2.col,lwd=weight)
		lines(w2$Time[order2],f1model.w2$fit[order2]+3*f1model.w2$pstd[order2],col=w2.col,lwd=1,lty=3)
		lines(w2$Time[order2],f1model.w2$fit[order2]-3*f1model.w2$pstd[order2],col=w2.col,lwd=1,lty=3)

		if(csv)
		{
			table <- unique(data.frame(time=w1$Time[order1],spline=f1model.w1$fit[order1],upper=f1model.w1$fit[order1]+3*f1model.w1$pstd[order1],lower=f1model.w1$fit[order1]-3*f1model.w1$pstd[order1]));
			write.csv(table,file=paste(word.1,'f1.csv',sep='_'),row.names=FALSE,quote=FALSE);

			table <- unique(data.frame(time=w2$Time[order2],spline=f1model.w2$fit[order2],upper=f1model.w2$fit[order2]+3*f1model.w2$pstd[order2],lower=f1model.w2$fit[order2]-3*f1model.w2$pstd[order2]));
			write.csv(table,file=paste(word.2,'f1.csv',sep='_'),row.names=FALSE,quote=FALSE);
		}
	}

	if(f2)
	{
		w1<-f2model$data[f2model$data$word=="word1",]
		w2<-f2model$data[f2model$data$word=="word2",]
		n1<-nrow(w1)
		ntot<-nrow(w1)+nrow(w2)
		f2model.pred<-predict(f2model)
		f2model.w1<-data.frame(fit=f2model.pred$fit[1:n1],pstd=f2model.pred$pstd[1:n1])
		f2model.w2<-data.frame(fit=f2model.pred$fit[n1+1:ntot],pstd=f2model.pred$pstd[n1+1:ntot])
		order1<-order(w1$Time)
		order2<-order(w2$Time)
		if(isplotted)
		{
			lines(f2model$data$Time,f2model$data$F2,type="n")
		}
		else
		{
			plot(f2model$data$Time,f2model$data$F2,type="n",xlab="Time (s)",ylab="Frequency (Hz)",ylim=yrange)
			isplotted<-TRUE;
		}
		lines(w1$Time[order1],f2model.w1$fit[order1],col=w1.col,lwd=weight)
		lines(w1$Time[order1],f2model.w1$fit[order1]+3*f2model.w1$pstd[order1],col=w1.col,lwd=1,lty=3)
		lines(w1$Time[order1],f2model.w1$fit[order1]-3*f2model.w1$pstd[order1],col=w1.col,lwd=1,lty=3)
		lines(w2$Time[order2],f2model.w2$fit[order2],col=w2.col,lwd=weight)
		lines(w2$Time[order2],f2model.w2$fit[order2]+3*f2model.w2$pstd[order2],col=w2.col,lwd=1,lty=3)
		lines(w2$Time[order2],f2model.w2$fit[order2]-3*f2model.w2$pstd[order2],col=w2.col,lwd=1,lty=3)

		if(csv)
		{
			table <- unique(data.frame(time=w1$Time[order1],spline=f2model.w1$fit[order1],upper=f2model.w1$fit[order1]+3*f2model.w1$pstd[order1],lower=f2model.w1$fit[order1]-3*f2model.w1$pstd[order1]));
			write.csv(table,file=paste(word.1,'f2.csv',sep='_'),row.names=FALSE,quote=FALSE);

			table <- unique(data.frame(time=w2$Time[order2],spline=f2model.w2$fit[order2],upper=f2model.w2$fit[order2]+3*f2model.w2$pstd[order2],lower=f2model.w2$fit[order2]-3*f2model.w2$pstd[order2]));
			write.csv(table,file=paste(word.2,'f2.csv',sep='_'),row.names=FALSE,quote=FALSE);

		}
	}

	if(f3)
	{
		w1<-f3model$data[f3model$data$word=="word1",]
		w2<-f3model$data[f3model$data$word=="word2",]
		n1<-nrow(w1)
		ntot<-nrow(w1)+nrow(w2)
		f3model.pred<-predict(f3model)
		f3model.w1<-data.frame(fit=f3model.pred$fit[1:n1],pstd=f3model.pred$pstd[1:n1])
		f3model.w2<-data.frame(fit=f3model.pred$fit[n1+1:ntot],pstd=f3model.pred$pstd[n1+1:ntot])
		order1<-order(w1$Time)
		order2<-order(w2$Time)
		par(mfrow=c(1,1))
		if(isplotted)
		{
			lines(f3model$data$Time,f3model$data$F3,type="n")
		}
		else
		{
			plot(f3model$data$Time,f3model$data$F3,type="n",xlab="Time (s)",ylab="Frequency (Hz)",ylim=yrange)
			isplotted<-TRUE;
		}
		lines(w1$Time[order1],f3model.w1$fit[order1],col=w1.col,lwd=weight)
		lines(w1$Time[order1],f3model.w1$fit[order1]+3*f3model.w1$pstd[order1],col=w1.col,lwd=1,lty=3)
		lines(w1$Time[order1],f3model.w1$fit[order1]-3*f3model.w1$pstd[order1],col=w1.col,lwd=1,lty=3)
		lines(w2$Time[order2],f3model.w2$fit[order2],col=w2.col,lwd=weight)
		lines(w2$Time[order2],f3model.w2$fit[order2]+3*f3model.w2$pstd[order2],col=w2.col,lwd=1,lty=3)
		lines(w2$Time[order2],f3model.w2$fit[order2]-3*f3model.w2$pstd[order2],col=w2.col,lwd=1,lty=3)

		if(csv)
		{
			table <- unique(data.frame(time=w1$Time[order1],spline=f3model.w1$fit[order1],upper=f3model.w1$fit[order1]+3*f3model.w1$pstd[order1],lower=f3model.w1$fit[order1]-3*f3model.w1$pstd[order1]));
			write.csv(table,file=paste(word.1,'f3.csv',sep='_'),row.names=FALSE,quote=FALSE);

			table <- unique(data.frame(time=w2$Time[order2],spline=f3model.w2$fit[order2],upper=f3model.w2$fit[order2]+3*f3model.w2$pstd[order2],lower=f3model.w2$fit[order2]-3*f3model.w2$pstd[order2]));
			write.csv(table,file=paste(word.2,'f3.csv',sep='_'),row.names=FALSE,quote=FALSE);
		}
	}

	

	# this is just regular plot stuff
	title(paste(word.1,"vs",word.2))
	coords <- par("usr");
	legend(coords[1]+abs(0.05*(coords[2]-coords[1])),coords[4]-abs(0.05*(coords[4]-coords[3])),c(word.1,word.2),lwd=weight,col=c(w1.col,w2.col))
}

comp.plotbci<-function(f1model,f2model,f3model,word.1,word.2,yrange,bw,weight,f1,f2,f3){
	if(bw==TRUE){
		w1.col <- 1;
		w2.col <- 3;
	} else
	{
		w1.col <- 2;
		w2.col <- 3;
	}

	nplots <- 0;
	if(f1){nplots <- nplots+1;}
	if(f2){nplots <- nplots+1;}
	if(f3){nplots <- nplots+1;}
	par(mfrow=c(nplots,2))

	# F1
	if(f1)
	{
		w1<-f1model$data[f1model$data$word=="word1",]
		w2<-f1model$data[f1model$data$word=="word2",]
		n1<-nrow(w1)
		ntot<-nrow(w1)+nrow(w2)
		f1model.pred<-predict(f1model,terms=c(0,1,0,1,0,1))
		f1model.w1<-data.frame(fit=f1model.pred$fit[1:n1],pstd=f1model.pred$pstd[1:n1])
		f1model.w2<-data.frame(fit=f1model.pred$fit[n1+1:ntot],pstd=f1model.pred$pstd[n1+1:ntot])
		order1<-order(w1$Time)
		order2<-order(w2$Time)

		ylimits<-c(min(-f1model.w1$fit[order1]-3*f1model.w1$pstd[order1]),
			   max(-f1model.w1$fit[order1]+3*f1model.w1$pstd[order1]))
		plot(f1model$data$Time,-1*(f1model.pred$fit),type="n",xlab="X",ylab="Y",ylim=ylimits)
		abline(0,0)
		lines(w1$Time[order1],-f1model.w1$fit[order1],col=w1.col,lwd=5)
		lines(w1$Time[order1],-f1model.w1$fit[order1]+3*f1model.w1$pstd[order1],col=w1.col,lwd=1,lty=3)
		lines(w1$Time[order1],-f1model.w1$fit[order1]-3*f1model.w1$pstd[order1],col=w1.col,lwd=1,lty=3)
		title(paste("Interaction effects w/BCI for F1 of",word.1))
		ylimits<-c(min(-f1model.w2$fit[order2]-3*f1model.w2$pstd[order2]),
			   max(-f1model.w2$fit[order2]+3*f1model.w2$pstd[order2]))
		plot(f1model$data$Time,-f1model.pred$fit,type="n",xlab="X",ylab="Y",ylim=ylimits)
		abline(0,0)
		lines(w2$Time[order2],-f1model.w2$fit[order2],col=w2.col,lwd=5)
		lines(w2$Time[order2],-f1model.w2$fit[order2]+3*f1model.w2$pstd[order2],col=w2.col,lwd=1,lty=3)
		lines(w2$Time[order2],-f1model.w2$fit[order2]-3*f1model.w2$pstd[order2],col=w2.col,lwd=1,lty=3)
		title(paste("Interaction effects w/BCI for F1 of",word.2))
	#	list(get.int=f1model.pred)
	}

	if(f2)
	{
		# F2
		w1<-f2model$data[f2model$data$word=="word1",]
		w2<-f2model$data[f2model$data$word=="word2",]
		n1<-nrow(w1)
		ntot<-nrow(w1)+nrow(w2)
		f2model.pred<-predict(f2model,terms=c(0,1,0,1,0,1))
		f2model.w1<-data.frame(fit=f2model.pred$fit[1:n1],pstd=f2model.pred$pstd[1:n1])
		f2model.w2<-data.frame(fit=f2model.pred$fit[n1+1:ntot],pstd=f2model.pred$pstd[n1+1:ntot])
		order1<-order(w1$Time)
		order2<-order(w2$Time)

		ylimits<-c(min(-f2model.w1$fit[order1]-3*f2model.w1$pstd[order1]),
			   max(-f2model.w1$fit[order1]+3*f2model.w1$pstd[order1]))
		plot(f2model$data$Time,-1*(f2model.pred$fit),type="n",xlab="X",ylab="Y",ylim=ylimits)
		abline(0,0)
		lines(w1$Time[order1],-f2model.w1$fit[order1],col=w1.col,lwd=5)
		lines(w1$Time[order1],-f2model.w1$fit[order1]+3*f2model.w1$pstd[order1],col=w1.col,lwd=1,lty=3)
		lines(w1$Time[order1],-f2model.w1$fit[order1]-3*f2model.w1$pstd[order1],col=w1.col,lwd=1,lty=3)
		title(paste("Interaction effects w/BCI for F2 of",word.1))
		ylimits<-c(min(-f2model.w2$fit[order2]-3*f2model.w2$pstd[order2]),
			   max(-f2model.w2$fit[order2]+3*f2model.w2$pstd[order2]))
		plot(f2model$data$Time,-f2model.pred$fit,type="n",xlab="X",ylab="Y",ylim=ylimits)
		abline(0,0)
		lines(w2$Time[order2],-f2model.w2$fit[order2],col=w2.col,lwd=5)
		lines(w2$Time[order2],-f2model.w2$fit[order2]+3*f2model.w2$pstd[order2],col=w2.col,lwd=1,lty=3)
		lines(w2$Time[order2],-f2model.w2$fit[order2]-3*f2model.w2$pstd[order2],col=w2.col,lwd=1,lty=3)
		title(paste("Interaction effects w/BCI for F2 of",word.2))
	#	list(get.int=f2model.pred)
	}

	if(f3)
	{
		# F3
		w1<-f3model$data[f3model$data$word=="word1",]
		w2<-f3model$data[f3model$data$word=="word2",]
		n1<-nrow(w1)
		ntot<-nrow(w1)+nrow(w2)
		f3model.pred<-predict(f3model,terms=c(0,1,0,1,0,1))
		f3model.w1<-data.frame(fit=f3model.pred$fit[1:n1],pstd=f3model.pred$pstd[1:n1])
		f3model.w2<-data.frame(fit=f3model.pred$fit[n1+1:ntot],pstd=f3model.pred$pstd[n1+1:ntot])
		order1<-order(w1$Time)
		order2<-order(w2$Time)

		ylimits<-c(min(-f3model.w1$fit[order1]-3*f3model.w1$pstd[order1]),
			   max(-f3model.w1$fit[order1]+3*f3model.w1$pstd[order1]))
		plot(f3model$data$Time,-1*(f3model.pred$fit),type="n",xlab="X",ylab="Y",ylim=ylimits)
		abline(0,0)
		lines(w1$Time[order1],-f3model.w1$fit[order1],col=w1.col,lwd=5)
		lines(w1$Time[order1],-f3model.w1$fit[order1]+3*f3model.w1$pstd[order1],col=w1.col,lwd=1,lty=3)
		lines(w1$Time[order1],-f3model.w1$fit[order1]-3*f3model.w1$pstd[order1],col=w1.col,lwd=1,lty=3)
		title(paste("Interaction effects w/BCI for F3 of",word.1))
		ylimits<-c(min(-f3model.w2$fit[order2]-3*f3model.w2$pstd[order2]),
			   max(-f3model.w2$fit[order2]+3*f3model.w2$pstd[order2]))
		plot(f3model$data$Time,-f3model.pred$fit,type="n",xlab="X",ylab="Y",ylim=ylimits)
		abline(0,0)
		lines(w2$Time[order2],-f3model.w2$fit[order2],col=w2.col,lwd=5)
		lines(w2$Time[order2],-f3model.w2$fit[order2]+3*f3model.w2$pstd[order2],col=w2.col,lwd=1,lty=3)
		lines(w2$Time[order2],-f3model.w2$fit[order2]-3*f3model.w2$pstd[order2],col=w2.col,lwd=1,lty=3)
		title(paste("Interaction effects w/BCI for F3 of",word.2))
	#	list(get.int=f3model.pred)
	}
}

library(assist)
options(memory=1000000000)