# R script to plot output data from Spatial Invasion Simulator
# WS Harpole 9/1/2007

dat1<-read.csv("./log.txt", header=F, skip=4 ) 
dat2<-read.table("./time.txt", header=T, sep=" " ) #time,species1, species2
mat1<-read.csv("./initial.txt", header=T )  
mat2<-read.csv("./final.txt", header=T )

# Percent Abundances
dat2$p1<- 100*(dat2$spp1/(dat2$spp1+dat2$spp2))
dat2$p2<- 100*(dat2$spp2/(dat2$spp1+dat2$spp2))

pdf(file = "sisout.pdf", width = 8, height = 8, onefile = T) # open output pdf file

def.par <- par(no.readonly = TRUE) # save default, for resetting...

palette(c("dark orange","blue")) # define colors: species 1, species 2

layout(matrix(c(1,2,3,4),2,2,byrow=TRUE), TRUE) # divide device into 4 panels

# Panel 1, Write out the Parameters
sp1dis<-dat1[1,2]
sp2dis<-dat1[2,2]
sp1fb<-dat1[3,2]
sp2fb<-dat1[4,2]
sp1m<-dat1[5,2]
sp2m<-dat1[6,2]
scenario<-dat1[7,2]
initabund<-dat1[8,2]
seed<-dat1[9,2]
edge<-dat1[10,2]
dim<-dat1[11,2]
tmax<-dat1[12,2]
pid<-dat1[13,2]

par( mar=c(0,1,3,0), mgp=c(1,0,0), adj=0 )
plot(1:10,1:10, type="n", xaxt="n", yaxt="n", xlab="", ylab="", bty="n", main="Spatial Invasion Simulator",cex.main=2 )
text(1,10,paste("process ID =",pid,sep=" "), cex=.75 )
text(2,9,"Parameters", cex=1.2 )
text(5,9,"exotic", cex=1.2, col=1)
text(8,9,"native", cex=1.2, col=2)
text(2,8, "dispersal:")
	if(sp1dis==1) text(5,8,"global") else text(5,8,"local")
	if(sp2dis==1) text(8,8,"global") else text(8,8,"local")
text(2,7, "feedback:")
	text(5.3,7,sp1fb)
	text(8.3,7,sp2fb)
text(2,6, "mortality:")
	text(5.3,6,sp1m)
	text(8.3,6,sp2m)
text(2,4.5, paste("exotic/native seed ratio: ",seed, ":1", sep="" ))
if(scenario==1){ text(2,4,"initial feedback state: exotic") 
	} else if(scenario==2){ text(2,4,"initial feedback state: native")
	} else { text(2,4,"initial feedback state: neutral") 
	}
text(2,3.5, paste("initial native abundance: ", initabund, "%", sep="" ))
if(edge==0){ text(2,3,"NO external exotic seed source") 
	} else { text(2,3,"YES external exotic seed source")
	}
text(2,2.5, paste("lattice:",dim,"x",dim,",    time steps:", tmax, sep=" " ))
lines(c(1,10),c(8.5,8.5), col="gray50" )
rect<-matrix(c(2,4),1,2)
symbols(5.5,7.5,rectangles=rect, fg=1, add=T, inches=F, lwd=4 )
symbols(8.5,7.5,rectangles=rect, fg=2, add=T, inches=F, lwd=4 )


# Panel 2, Time Series Plot
par( mar=c(4,4,3,3), mgp=c(2,1,0), adj=0.5 )
plot(dat2$time, dat2$p2, type="l" , xlim=c(0,tmax), ylim=c(0,100), ylab="% Native", xlab="Years", col=2, lwd=4, cex.axis=1, cex.lab=1 , bty="o" )
abline(h=100, col="gray30", lty=3)
abline(h=0, col="gray30", lty=3)

# Panel 3, Initial matrix
par( mar=c(0,0,1,0), mgp=c(-1,0,0) )
symbols( mat1$i, mat1$j, squares=rep( 1, (max(mat1$i)+1)*(max(mat1$j)+1) ), inches=F, fg=mat1$x, bg=mat1$x, xlab="", ylab="", xaxt="n", yaxt="n", add=F, bty="n" )
text(dim*1.1,dim/2,"->",cex=2)
text(dim/2,dim*1.1,"Initial",cex=1.8)

# Panel 4, Final matrix (after tmax time steps)
par( mar=c(0,0,1,0), mgp=c(-1,0,0) )
symbols( mat2$i, mat2$j, squares=rep( 1, (max(mat2$i)+1)*(max(mat2$j)+1) ), inches=F, fg=mat2$x, bg=mat2$x, xlab="", ylab="", xaxt="n", yaxt="n", add=F, bty="n" )
text(dim/2,dim*1.1,"Final",cex=1.8)

palette("default")
par(def.par)#- reset to default
#dev.off()