rm(list=ls()) # clears workspace
library(stats)
library(MASS)
library(time)

			#####################################################################
			# If running on a computer, change the directory names as appropriate
			#####################################################################
####  Is this code being submitted to a server or run on a personal computer?
Answer<-'SR'	# Desktop - 'PC'
####  If a PC desktop, specify the following directories:
# Specify working directory
Rdir<-'/Users/marknovak/Documents/Research/R-Codes/LoopAnalysis/SppConSimulations/Output/'
#Specify directory in which Mathematica and the 'LoopAnalysis.m' file are located
Mathdir<-'/Applications/Mathematica.app/Contents/MacOS/'
			#####################################################################
# Set the species richness levels
Spp<-seq(6,26,2)
# Set the connectance levels
EW_cons<- c(0.1,0.15,0.25,0.35)
# Set tolerance level for the connectance of niche webs (+/-%)
tol<-0.025				
#### Note that it may not be possible to creat all Spp X Connectance levels, dependent upon the tolerance level
# On a scale of 0 to 1, what should the cut-off value for a "strong" weighted prediction be?
Cut<-0.5
# Define the number of replicate webs per Richness/Connectance level needed
Reps<-5
# Define the maximum number of times to try making a given web before giving up
MaxTries<-2e5
# Self-limitation for all species?  T - yes, F - no.  Self-limitation almost ensures a positive determinant which is needed to calculated the adjoint matrix
SelfLim<-'T'

						##################################################
						################### Run Code #####################
						##################################################
if (Answer=='PC'){setwd(Rdir)}
Raw.WPreds<-dim(0)
Attempts<-dim(0)
				
for (z in 1:length(Spp)){
for (y in 1:length(EW_cons)){
	S<-Spp[z]
	EW_con<-EW_cons[y]
	
# initialize arrays and looping variable
	chat<-dim(0)
	Niche_web<-array(NA,dim=c(S,S,Reps))
	Adj<-array(NA,dim=c(S,S,Reps))
	AbsFeed<-array(NA,dim=c(S,S,Reps))
	WPred<-array(NA,dim=c(S,S,Reps))
	Body_mass<-array(NA,dim=c(S,Reps))
	Diet_breadth<-array(NA,dim=c(S,2,Reps))
	nac<-array(NA,dim=c(S,Reps))
	mean.WPred<-array(NA,dim=c(1,Reps))
	prop.cut.WPred<-array(NA,dim=c(1,Reps))
	quant.WPred<-array(NA,dim=c(1,Reps))
	num.cut.WPred<-array(NA,dim=c(1,Reps))
	Times<-array(NA,dim=c(1,Reps))
	Links<-dim(0)
	NW_conIs<-dim(0)
	NW_conDs<-dim(0)
	NW_conIs.tries<-dim(0)
	j<-0
	tries<-0
	NAC<-0
	NAC.tries<-0
	
	while(j<Reps&tries<=MaxTries){
		
	# Define species body mass M (position on niche axis)
		M<-runif(S,0,1)
	# # Define the diet range (diet breadth) for each species, 'r'.
		# 1st shape parameter of the Beta function.
		a<-1                          
		# 2nd shape parameters of the Beta function defined by the connectance of the web
		# add noise to EW_con to counter effect of measuring interactive connectance as opposed to directed connectance
		temp.EW_con<-EW_con*runif(1,1,1.5)
		b<-(1-2*temp.EW_con)/(2*temp.EW_con)  	
					#############################
					###  BAGDASSARIAN METHOD  ###
					#############################
# 		x<-rbeta(S,a,b)/2		# Divide by two since E(ri)=2*Connectance
# 		r<-M*x
# 	# Define the center of the diet range, 'chat'.
#  		for (i in 1:S){	chat[i]<-runif(1,r[i]/2,M[i])}
		
					####################################
					###  WILLIAMS & MARTINEZ METHOD  ###
					###  according to 2004 paper	 ###
					####################################
	 	r<-rbeta(S,a,b)/2		# Divide by two since E(ri)=2*Connectance
	 # Define the center of the diet range, 'chat'.
	 	temp.c<-apply(cbind(M,1-(r/2)),1,min)
	 	for (i in 1:S){	chat[i]<-runif(1,r[i]/2,temp.c[i]) 	}
										####################################
	# Reorder all vectors M, C, r by M
		o<-order(M)
		M<-M[o]
		r<-r[o]
		chat<-chat[o]
		
	# Assign species with lowest 'M' an 'r' of zero (from Williams & Martines 2004)
		r[1]<-0
	
	# Define upper and lower limits to r for each species.
		 DB<-cbind(chat-r/2,chat+r/2)               
	
	# Now find which species feed on each other
		A <- matrix(0,S,S)
		for (i in 1:S){
		    # Define binary matrix for food web
		    ind <- which(M>DB[i,1] & M<DB[i,2])
		    A[i,ind] <- 1
		    A[ind,i] <- -1
		}
	
		# Rearrange the matrix into upper and lower triangular form
	    A1 <- A
	    A1[which(A1>0)]<-0
	    sift <- apply(abs(A1),2,sum)	# column sums			
	    o <- order(sift)
	    A <- A[o,o]
	    M <- M[o]
	    chat <- chat[o]
	    DB <- DB[o,]
	    
		# Remove cannabalistic links  
		# This is not done in Williamms & Martinez's orginal forumation, but is required since we will, for the sake of loop analysis, be setting all diagonal elements to -1
	 	A <- A-diag(diag(A))	
# 		diag(A)<-0		# same as above	
		# Niche Web interactive connectance (only correct if diagonal elements are removed)
		NW_conI<-sum(abs(A-diag(diag(A))))/(S*(S-1))
		# count non-diag elements for each species
		nac.tries<-apply(abs(A),2,sum)
		Z<-ginv(abs(0.9*A))
		comps<-length(which(Z==0))
		# Count the number of attempted webs in which not all species were connected
		if(min(nac.tries)==0 & comps>0){NAC.tries<-NAC.tries+1}		

# 		# Determine if any species have identical pred & prey
# 		ident<-dim(0)
# 		for (p in 1:S){ for (o in 1:p-1){
# 			if (identical(A[p,p],A[o,o]) & identical(A[p,-c(p,o)],A[o,-c(p,o)]) & identical(A[-c(p,o),p],A[-c(p,o),o])){
# 			ident<-rbind(ident,c(p,o))
# 		}}}
	###############################
	# Accept this niche web if... 
	###############################
		# its connectance is +/- x% of the target
		if (NW_conI>(EW_con-tol) & NW_conI<(EW_con+tol)){ 	
		# all species are connected
		if(min(nac.tries)!=0 & comps==0){									
		# no two species are identical in their predators and prey
# 		if(length(ident)==0){								
 		# make all diagonal elements -1. This almost ensures stability and a determinant > 0
 		if (SelfLim=='T'){diag(A)<- -1}
 		# Niche Web directed connectance
		NW_conD<-(sum(abs(A)))/(S^2)

 		# if the determinant of the matrix is zero it is 'neutrally stable' and cannot by analysed qualitatively.  The matrix cannot be inverted.
 		if(det(A)>0){ 									
			# count successfuly created webs +1
			j<-j+1										
			Niche_web[,,j] <- A
			Body_mass[,j] <- M
			Diet_breadth[,,j] <- DB
			NW_conIs[j]<-NW_conI
			NW_conDs[j]<-NW_conD
			Links[j]<-sum(abs(A-diag(diag(A))))/2
			# count non-diagonal links per predator
			nac[,j]<-apply(abs(A),2,sum)				
	
			fout1<- file(paste('Output/Math/NicheWeb',S,'spp','ExpConn-', round(EW_con,2),'Web-',j,'NWconn-',round(NW_conI,2),'.csv'),"w")
				write.table(A,file=fout1,sep=',',row.names=F, col.names=F,quote=F)
			close(fout1)

	#############		
	# Redirect to Mathematica in order to calculate adjoint, total feedback, and weighted prediction matrix
	if (Answer=='PC'){setwd(Mathdir)}
		fout1<- file(paste('A.csv'),"w");
		write.table(A,file=fout1,sep=',',row.names=F,col.names=F,quote=F); close(fout1)
# 		time.start<-getTime()
# 			shell('math < LoopAnalysis.m',mustWork=T) # for Windows computers
# 			system('Mathkernel < LoopAnalysis.m',wait=T) # for mac computers
			sys.time<-system.time(system('math < LoopAnalysis.m',wait=T)) # for server
# 		time.end<-getTime()
		Adj[,,j]<-as.matrix(read.csv("Adj.csv",header=FALSE))
		AbsFeed[,,j]<-as.matrix(read.csv("AbsFeed.csv",header=FALSE))
		WPred[,,j]<-as.matrix(read.csv("WPred.csv",header=FALSE))
	if (Answer=='PC'){setwd(Rdir)}
		fout1<- file(paste('Output/Math/NicheWeb',S,'spp','ExpConn-', round(EW_con,2),'Web-',j,'NWconn-', round(NW_conI,2),'Adj.csv'),"w")
		write.table(Adj[,,j],file=fout1,sep=',', row.names=F,col.names=F,quote=F); close(fout1)
		fout1<- file(paste('Output/Math/NicheWeb',S,'spp','ExpConn-', round(EW_con,2),'Web-',j,'NWconn-',round(NW_conI,2), 'AbsFeed.csv'),"w")
		write.table(AbsFeed[,,j],file=fout1,sep=',', row.names=F,col.names=F,quote=F); close(fout1)
		fout1<- file(paste('Output/Math/NicheWeb',S,'spp','ExpConn-', round(EW_con,2),'Web-',j,'NWconn-', round(NW_conI,2),'WPred.csv'),"w")
		write.table(WPred[,,j],file=fout1,sep=',', row.names=F,col.names=F,quote=F); close(fout1)
	##############
	# calculate the mean weighted prediction value for the web
	mean.WPred[j]<-mean(WPred[,,j])
	# calculate the proportion of weighted predictions falling above a given cut-off value
	prop.cut.WPred[j]<-length(which(WPred[,,j]>Cut))/S^2		
	quant.WPred[j]<-quantile(WPred[,,j],0.95)	
	num.cut.WPred[j]<-length(which(WPred[,,j]>Cut))		
# 	time.elapsed<- (time.end-time.start)[1]*60+(time.end-time.start)[2]+(time.end-time.start)[3]/60 #convert to minutes
# 	Times[j]<-time.elapsed
	Times[j]<-sys.time[3]/60 # Elapsed time in seconds converted to minutes
		}
# 		}
		}
		}
		tries<-tries+1		#  count the number of webs that were attempted
		NW_conIs.tries[tries]<-NW_conI
	}
	Tries<-tries
	Succ<-j
	Attempts<-rbind(Attempts,c(S,EW_con,Tries,Succ))

	if (Succ<Reps&Tries==MaxTries){stop("Number of attempts at producing desired replicate webs exceeds set cut-off")}
	if(Succ==Reps){
		PSucc<-round(Succ/Tries,3)
		cols<-nac
		cols[which(nac>0)]<-'black'
		cols[which(nac==0)]<-'red'
		unconBS<-Body_mass[which(nac==0)]
		conBS<-Body_mass[which(nac>1)]
		pdf(paste('Output/Math/NicheModel-Niches',S,'spp','ExpConn-',EW_con,'.pdf'), width=8, height=10)
				plot(c(-0.1,1.1),c(0,Reps+1),type='n',xlab='Niche Axis',ylab='Replicate Web',frame.plot=F)
				segments(0,seq(1:Reps),1,seq(1:Reps),lwd=2)
				for (j in 1:Reps){
					segments(Body_mass[,j],j,Body_mass[,j],j+0.1,lwd=4,col=cols[,j])
				}
				legend('topright',col='red',lwd=4,legend='Unconnected',bty='n')
			par(mfrow=c(3,2))
				if(length(conBS)>0){hist(conBS, main=paste('Body Sizes of All Species\n',S,'Species',EW_con,'connectance'), xlab='Body Size (0 to 1)') }
		dev.off()
		
		colnames(Body_mass)<-paste('Web',seq(1:Reps))
		fout1<-file(paste('Output/Math/NicheModel-NicheValues',S,'spp','ExpConn-', EW_con,'.csv'),"w")
			write.table(Body_mass,file=fout1,sep=',',row.names=F,col.names=T,quote=F)
		close(fout1)
		############
		raw.WPreds<-data.frame(rep(S,Reps),rep(EW_con,Reps), NW_conIs,NW_conDs, Links, t(mean.WPred),t(prop.cut.WPred),t(quant.WPred),t(num.cut.WPred),t(Times))
		Raw.WPreds<-rbind(Raw.WPreds,raw.WPreds)
	}
}}

Attempts<-data.frame(Attempts); colnames(Attempts)<-c('S','EW_con','Tries','Succ')
fout1<-file(paste('Output/Math/a-NicheModel-WPreds-SppVSConnectance-Attempts.csv'),"w")
	write.table(Attempts,file=fout1,sep=',',row.names=F,col.names=T,quote=F)
close(fout1)

if(nrow(Raw.WPreds)>0){
Raw.WPreds<-data.frame(Raw.WPreds)
colnames(Raw.WPreds)<- c('S','EW_con','NW_conI','NW_conD','Links','WPred.Mean','WPred.PropCut', 'WPred.95Quantile','WPred.NumCut', 'RunTime')
fout1<-file(paste('Output/Math/a-NicheModel-WPreds-SppVSConnectance.csv'),"w")
	write.table(Raw.WPreds,file=fout1,sep=',',row.names=F,col.names=T,quote=F)
close(fout1)
}

#############################################################################################
#############################################################################################
#############################################################################################