#This code accompanies the paper Quicke, Jones & Epstein (2007) Correcting the problem of false incongruence due to noise imbalance in the incongruence length difference (ILD) test 
#Systematic Biology, Volume 56, Issue 3 May 2007 , pages 496 - 503
#The code was written by Owen R. Jones (owen.jones@imperial.ac.uk)
#
#
#This code takes the outputs from PAUP and parses them to generate Figure 1


#parsing the linearisation outputs
#import thedata
setwd("/Users/orj/Documents/ILD Test MS/Programming/Simulation_files/Type_6_noise/linearisation/PAUP outputs/")
thedata<-list.files(pattern="OUT4")

ntax<-c(8,16,32,64,128,256,512)
all.lineardata<-NULL

for (i in 1:length(thedata)){
	
	path<-thedata[i]
	x1<-read.delim(path)
	
	#standardise the length
	xmax<-subset(x1,pseudorep==999)
	x<-subset(x1,pseudorep!=999)
	x$noiselev<-x$noiselev+1
	
	nnoise<-length(unique(x$noiselev))
	noiselevs<-unique(x$noiselev)
	pseudoreps<-unique(x$pseudorep)
	numpseudo<-length(pseudoreps)
	
	#calculate the standardised lengths and 1-RI
	#standardise partition 1
	xnoisy1<-subset(xmax,part==1)
	meanmaxnoisyL1<-mean(xnoisy1$lengths)
	meanmaxnoisyOMRI1<-mean(xnoisy1$OneMinusRI)
		
	x$standardisedlength[x$part==1]<-(x$lengths[x$part==1]-min(x$lengths[x$part==1]))/(meanmaxnoisyL1-min(x$lengths[x$part==1]))
	x$standardOneMinusRI[x$part==1]<-(x$OneMinusRI[x$part==1]-0)/(meanmaxnoisyOMRI1-0)
	x$standardisedlength[x$part==1][is.nan(x$standardisedlength[x$part==1])==TRUE]<-0
	x$standardOneMinusRI[x$part==1][is.nan(x$standardOneMinusRI[x$part==1])==TRUE]<-0
	
	#standardise partition 2
	xnoisy2<-subset(xmax,part==2)
	meanmaxnoisyL2<-mean(xnoisy2$lengths)
	meanmaxnoisyOMRI2<-mean(xnoisy2$OneMinusRI)
		
	x$standardisedlength[x$part==2]<-(x$lengths[x$part==2]-min(x$lengths[x$part==2]))/(meanmaxnoisyL2-min(x$lengths[x$part==2]))
	x$standardOneMinusRI[x$part==2]<-(x$OneMinusRI[x$part==2]-0)/(meanmaxnoisyOMRI2-0)
	x$standardisedlength[x$part==2][is.nan(x$standardisedlength[x$part==2])==TRUE]<-0
	x$standardOneMinusRI[x$part==2][is.nan(x$standardOneMinusRI[x$part==2])==TRUE]<-0
	
	x$asinOneMinusRI[x$standardOneMinusRI<=1]<-asin(x$standardOneMinusRI[x$standardOneMinusRI<=1])
	x$asinOneMinusRI[x$standardOneMinusRI>1]<-asin(1)
	
	x$asinstandardisedlength[x$standardisedlength<=1]<-asin(x$standardisedlength[x$standardisedlength<=1])
	x$asinstandardisedlength[x$standardisedlength>1]<-asin(1)
	
	x<-subset(x,part==2)
	x$taxa<-ntax[i]
	
	all.lineardata<-rbind(all.lineardata,x)
	}



#generate the plots
quartz(h=14,w=10.5)
par(mfrow=c(6,2))

for (i in 1:6){
	df1<-subset(all.lineardata,taxa==ntax[i])
	par(mar=(c(5, 4, 4, 6) + 0.1))
	par(cex.axis=1.3,cex.lab=1.3)
	
	tax<-df1$tax[i]
	rs<-2.5
	plot(c(0,1),c(0,rs),xlab="Proportion noisy characters",ylab="Standardised length",type="n",axes=F)
	mtext("Arcsine standardised length",4,3,cex=0.8)
	title(paste("Length:",tax,"taxa"))
	axis(1,at=seq(0,1,0.1))
	axis(2,at=seq(0,1,0.2))
	abline(0,rs,lty=3)
	tickpoints<-c(seq(0,asin(1),0.2),asin(1))
	
	axis(4,at=rescale(tickpoints,c(0,rs)),labels=c(tickpoints[1:length(tickpoints)-1],"1.57"))
	points(jitter(df1$propnoise),df1$standardisedlength,pch=1,cex=2)
	asL<-rescale(df1$asinstandardisedlength,c(0,rs))
	points(jitter(df1$propnoise),asL,pch=16,cex=2)
	
	
	plot(c(0,1),c(0,rs),xlab="Proportion noisy characters",ylab="Standardised 1-RI",type="n",axes=F)
	mtext("Arcsine standardised 1-RI",4,3,cex=0.8)
	title(paste("1-RI:",tax,"taxa"))
	axis(1,at=seq(0,1,0.1))
	axis(2,at=seq(0,1,0.2))
	abline(0,rs,lty=3)
	tickpoints<-c(seq(0,asin(1),0.2),asin(1))
	
	axis(4,at=rescale(tickpoints,c(0,rs)),labels=c(tickpoints[1:length(tickpoints)-1],"1.57"))
	points(jitter(df1$propnoise),df1$standardOneMinusRI,pch=1,cex=2)
	
	asOMRI<-rescale(df1$asinOneMinusRI,c(0,rs))
	points(jitter(df1$propnoise),asOMRI,pch=16,cex=2)
	}

dev.copy2eps(file="Figure_1_ILD_linearisation.eps")