rm(list=ls())
setwd('/Users/marknovak/Documents/Research/R-Codes/LoopAnalysis-SppConSimulations/Output/Math-FullRun-Final/')
dat<-read.csv("a-NicheModel-WPreds-SppVSConnectance.csv",header=T)

se<-function(x){se<-sd(x)/sqrt(length(x));se}

Means<-aggregate(dat,list(S=dat$S,EW_con=dat$EW_con),mean)[,-c(3,4)]
SE<-aggregate(dat,list(S=dat$S,EW_con=dat$EW_con),se)[,-c(3,4)]

##########################################
EW_cons<-unique(Means$EW_con)

pt.symb<-c(25,22,23,21)
pdf(paste('NicheModels-WPreds-xSppxCon-Stacked-Means_Cutoff-Prop.pdf'),width=3, height=6.5)
par(mfrow=c(3,1),cex=0.5,cex.axis=1,cex.lab=1.3,mar=c(5,5,2,1),tck=-0.03,mgp=c(2.5,0.4,0))
	plot(Means$WPred.Mean~Means$S,type='n', xaxt='n', ylim=c(0,1),xlab='Food Web Species Richness',ylab='Grand Mean Prediction Weight',axes=F)
		axis(1,seq(6,24,2),at=seq(6,24,2),mgp=c(3,0.8,0))
		axis(2,mgp=c(3,0.8,0))
		box(lwd=1.5)
		for (k in 1:length(EW_cons)){
			arrows(Means$S[Means$EW_con==EW_cons[k]], Means$WPred.Mean[Means$EW_con==EW_cons[k]]+ SE$WPred.Mean[SE$EW_con==EW_cons[k]], Means$S[Means$EW_con==EW_cons[k]], Means$WPred.Mean[Means$EW_con==EW_cons[k]]- SE$WPred.Mean[SE$EW_con==EW_cons[k]], lwd=1,length=0.02,angle=90,code=3, col='black')
			points(Means$S[Means$EW_con==EW_cons[k]], Means$WPred.Mean[Means$EW_con==EW_cons[k]], type='o',pch=pt.symb[k], bg='grey', cex=1.6,lwd=1)
		}
# 	legend('bottomleft',legend=EW_cons,cex=1,pch=pt.symb,pt.bg='grey',pt.cex=1.5, title='Connectance',inset=0.01,bty='n')	
	mtext('a',side=3,line=0.5,at=2.1)

# plot(1,1,xlab='',ylab='',axes=F,type='n')
# legend('center',legend=EW_cons,cex=1.5,pch=pt.symb,pt.bg='grey',pt.cex=1.6, title='Connectance')	

plot(Means$WPred.PropCut~Means$S, type='n', xaxt='n', ylim=c(0,1), xlab='Food Web Species Richness', ylab=paste('Mean Proportion of W > 0.5'),axes=F)
		axis(1,seq(6,24,2),at=seq(6,24,2),mgp=c(3,0.8,0))
		axis(2,mgp=c(3,0.8,0))
		box(lwd=1.5)
		for (k in 1:length(EW_cons)){
			arrows(Means$S[Means$EW_con==EW_cons[k]], Means$WPred.PropCut[Means$EW_con==EW_cons[k]]+ SE$WPred.PropCut[SE$EW_con==EW_cons[k]], Means$S[Means$EW_con==EW_cons[k]], Means$WPred.PropCut[Means$EW_con==EW_cons[k]]- SE$WPred.PropCut[SE$EW_con==EW_cons[k]], lwd=1,length=0.02,angle=90,code=3, col='black')
			
			points(Means$S[Means$EW_con==EW_cons[k]], Means$WPred.PropCut[Means$EW_con==EW_cons[k]], type='o',pch=pt.symb[k],bg='grey',cex=1.6,lwd=1)
		}
# 	legend('bottomleft',legend=EW_cons,cex=1,pch=pt.symb,pt.bg='grey',pt.cex=1.5, title='Connectance',inset=0.01,bty='n')	
	mtext('b',side=3,line=0.5,at=2.1)

plot(dat$WPred.PropCut~dat$L, type='n',  ylim=c(0,1), xlim=c(0,105),xlab='Predator-Prey Links', ylab=paste('Mean Proportion of W > 0.5'),axes=F)
		axis(1,seq(0,100,20),at=seq(0,100,20),mgp=c(3,0.8,0))
		axis(2,mgp=c(3,0.8,0))
		box(lwd=1.5)
		for (k in 1:length(EW_cons)){
			points(dat$L[dat$EW_con==EW_cons[k]], dat$WPred.PropCut[dat$EW_con==EW_cons[k]], pch=pt.symb[k],bg='grey',cex=1.6,lwd=1)
		}
	legend('topright',legend=EW_cons,cex=1,pch=pt.symb,pt.bg='grey',pt.cex=1.5, title='Connectance',inset=0.05,bty='n')	
	mtext('c',side=3,line=0.5,at=-23)
dev.off()
	

#########################
pdf(paste('NicheModels-WPreds-xSppxCon-Stacked-Means_Cutoff-Num.pdf'),width=3, height=4.33)
par(mfrow=c(2,1),cex=0.5,cex.axis=1,cex.lab=1.3,mar=c(5,5,2,1),tck=-0.03,mgp=c(2.5,0.4,0))
plot(Means$WPred.NumCut~Means$S, type='n', xaxt='n', ylim=c(0,300), xlab='Food Web Species Richness', ylab=paste('Mean Number of W > 0.5'),axes=F)
		axis(1,seq(6,24,2),at=seq(6,24,2),mgp=c(3,0.8,0))
		axis(2,mgp=c(3,0.8,0))
		box(lwd=1.5)
		for (k in 1:length(EW_cons)){
			arrows(Means$S[Means$EW_con==EW_cons[k]], Means$WPred.NumCut[Means$EW_con==EW_cons[k]]+ SE$WPred.NumCut[SE$EW_con==EW_cons[k]], Means$S[Means$EW_con==EW_cons[k]], Means$WPred.NumCut[Means$EW_con==EW_cons[k]]- SE$WPred.NumCut[SE$EW_con==EW_cons[k]], lwd=1,length=0.02,angle=90,code=3, col='black')
			
			points(Means$S[Means$EW_con==EW_cons[k]], Means$WPred.NumCut[Means$EW_con==EW_cons[k]], type='o',pch=pt.symb[k],bg='grey',cex=1.6,lwd=1)
		}
# 	legend('topleft',legend=EW_cons,cex=1,pch=pt.symb,pt.bg='grey',pt.cex=1.5, title='Connectance',inset=0.05,bty='n')	
	mtext('a',side=3,line=0.5,at=2.1)

plot(dat$WPred.NumCut~dat$L, type='n', xlim=c(0,105),ylim=c(0,300),xlab='Predator-Prey Links', ylab=paste('Mean Number of W > 0.5'),axes=F)
		axis(1,seq(0,100,20),at=seq(0,100,20),mgp=c(3,0.8,0))
		axis(2,mgp=c(3,0.8,0))
		box(lwd=1.5)
		for (k in 1:length(EW_cons)){
			points(dat$L[dat$EW_con==EW_cons[k]], dat$WPred.NumCut[dat$EW_con==EW_cons[k]], pch=pt.symb[k],bg='grey',cex=1.6,lwd=1)
		}
# 	legend('topright',legend=EW_cons,cex=1,pch=pt.symb,pt.bg='grey',pt.cex=1.5, title='Connectance',inset=0.05,bty='n')	
	mtext('e',side=3,line=0.5,at=-23)

dev.off()

#################################################
# Runtime regression for webs
Conn<-0.15
regn.data<-Means[Means$S>8 & Means$EW_con==Conn,]
regn<-lm(log10(regn.data$RunTime)~regn.data$S)
time.min.100sp<-10^(regn$coef[1]+regn$coef[2]*100)
time.year.100sp<-time.min.100sp/(60*24*365)

# Plot runtimes for one connectance level
time.means<-Means[Means$EW_con==Conn,]
time.sd<-SE[SE$EW_con==Conn,]
ylim=c(4*10^-3,1.2*max(time.means$RunTime+time.sd$RunTime))
pdf(paste('NicheModels-WPreds-xSppxCon-RunTimes-1con-Means.pdf'),width=3, height=2.25)
par(cex=0.5,cex.axis=1,cex.lab=1.3,mar=c(5,5,2,1))
	plot(time.means$RunTime~time.means$S,type='n',xlab='Food Web Species Richness', log='y',ylab='Run Time (minutes)',ylim=ylim,axes=F)
		box(lwd=1.5)
		axis(1,seq(6,24,2),at=seq(6,24,2),mgp=c(3,0.8,0),tck=-0.03)
		tck1<-c(1e-3,1e-2,1e-1,1,1e1,1e2,1e3)
		axis(2,labels=tck1,at=tck1,las=2,mgp=c(3,0.8,0),tck=-0.03)
		arrows(time.means$S, time.means$RunTime+time.sd$RunTime, time.means$S, time.means$RunTime-time.sd$RunTime,				   lwd=1,length=0.02,angle=90,code=3, col='black')
		points(time.means$S,time.means$RunTime, type='p',pch=21, bg='grey',cex=1.6,lwd=1)
dev.off()

ylim=c(8*10^-3,1.5*max(time.means$RunTime+time.sd$RunTime))
pdf(paste('NicheModels-WPreds-xSppxCon-RunTimes-1con-Raw.pdf'),width=3, height=2.25)
par(cex=0.5,cex.axis=1,cex.lab=1.3,mar=c(5,5,2,1))
time.dat<-dat[dat$EW_con==Conn,]
plot(time.dat$RunTime~time.dat$S,type='n',xlab='Food Web Species Richness', log='y',ylab='Run Time (min)',axes=F,ylim=ylim)
		box(lwd=1.5)
		axis(1,seq(6,24,2),at=seq(6,24,2),mgp=c(3,0.8,0),tck=-0.03)
		tck1<-c(1e-3,1e-2,1e-1,1,1e1,1e2,1e3)
		axis(2,labels=tck1,at=tck1,las=2,mgp=c(3,0.8,0),tck=-0.03)
		points(time.dat$S,time.dat$RunTime, type='p',pch=21, bg='grey',cex=1.6,lwd=1)
dev.off()


# # Plot runtimes for two connectance levels
# con<-c(0.1,0.3)
# time.means<-Means[Means$EW_con==con[1] | Means$EW_con==con[2],]
# time.sd<-SE[SE$EW_con==con[1] | SE$EW_con==con[2],]
# ylim=c(5*10^-4,1.1*max(time.means$RunTime+time.sd$RunTime))
# 
# pdf(paste('NicheModels-WPreds-xSppxCon-RunTimes-2con-Means.pdf'),width=3, height=2.25)
# par(cex=0.5,cex.axis=1,cex.lab=1.3,mar=c(5,5,2,1))
# 	plot(time.means$RunTime~time.means$S,type='n',xlab='Food Web Species Richness', log='y',ylab='Run Time',ylim=ylim,axes=F)
# 		box(lwd=1.5)
# 		axis(1,seq(6,22,2),at=seq(6,22,2))
# 		tck1<-c(1e-3,1e-2,1e-1,1,1e1,1e2,1e3)
# 		axis(2,labels=tck1,at=tck1,las=2)
# 	for (i in con){
# 		temp.means<-time.means[time.means$EW_con==i,]
# 		temp.sd<-time.sd[time.sd$EW_con==i,]
# 		arrows(temp.means$S, temp.means$RunTime+temp.sd$RunTime, temp.means$S, temp.means$RunTime-temp.sd$RunTime,				   lwd=1,length=0.02,angle=90,code=3, col='black')
# 		points(temp.means$S,temp.means$RunTime, type='p',pch=21, bg='grey',cex=1.6,lwd=1)
# 		}
# dev.off()