# Martha's graphs

mm <- read.delim("~/data/marthasmodel.txt")
postscript("~/data/platy/svg/marthamodel.eps",
	height = 2, width = 3, horizontal = FALSE,
	onefile = FALSE, paper = "special")
xyplot(platyproportion ~ distance,
	data = mm,
	groups = release,
	pch = c(49,50,51,16),
	col = "black",
	cex = 0.5)
dev.off()

rs <- read.delim("~/data/releasesummary.txt")
postscript("~/data/svg/releasesummary.eps",
	height = 4, width = 6, horizontal = FALSE,
	onefile = FALSE, paper = "special")
xyplot(parasitismrate ~ distance | scale,
	data = rs,
	groups = release,
	lty = 1,
	pch = c(49,50,51),
	col = "black",
	layout = c(1,3),
	type = "b")
dev.off()

foraging <- read.delim("~/data/foraging.txt")
foraging.plot <- function() {
	attach(foraging)
	x <- log(eggs + 1)
	y <- timeonplant/eggs
	plot(y ~ x)
	abline(coef(lm(y ~ x)))
	n <- data.frame(x = seq(min(x)-2, max(x)+2, .1))
	yhat <- predict(lm(y ~ x), n, interval = "confidence")
	lines(n$x, yhat[,2], col = "gray")
	lines(n$x, yhat[,3], col = "gray")
	detach(foraging)
	print(paste("R2:", summary(lm(y ~ x))$r.squared))
	print(paste("y = ", as.numeric(summary(lm(y ~ x))$coef[1]), 
		"±", as.numeric(summary(lm(y ~ x))$coef[3]) * qt(.975, 
			as.numeric(summary(lm(y ~ x, data = foraging))$df[2])),
		"+", as.numeric(summary(lm(y ~ x))$coef[2]), 
		"±", as.numeric(summary(lm(y ~ x))$coef[4]) * qt(.975, 
			as.numeric(summary(lm(y ~ x, data = foraging))$df[2])),
		"x"))
	}
postscript("~/data/foraging.eps",
	height = 3, width = 3, horizontal = FALSE,
	onefile = FALSE, paper = "special")
foraging.plot()
dev.off()


fp <- data.frame(eggs = c(16, 40, 8, 12, 14, 52, 69, 131, 132, 62, 92),
	parasitism = c(0.3125, 0.125, 0.625, 0.583333333, 0.785714286,
		0, 0.072463768, 0.106870229, 0.212121212, 0.548387097, 0.358695652))
fp.plot <- function() {
	attach(fp)
	x <- eggs
	y <- parasitism
	plot(y ~ x)
	abline(coef(lm(y ~ x)))
	n <- data.frame(x = seq(min(x)-10, max(x)+10, .1))
	yhat <- predict(lm(y ~ x), n, interval = "confidence")
	lines(n$x, yhat[,2], col = "gray")
	lines(n$x, yhat[,3], col = "gray")
	detach(fp)
	print(paste("R2:", summary(lm(y ~ x))$r.squared))
	print(paste("y = ", as.numeric(summary(lm(y ~ x))$coef[1]), 
		"±", as.numeric(summary(lm(y ~ x))$coef[3]) * qt(.975, 
			as.numeric(summary(lm(y ~ x, data = fp))$df[2])),
		"+", as.numeric(summary(lm(y ~ x))$coef[2]), 
		"±", as.numeric(summary(lm(y ~ x))$coef[4]) * qt(.975, 
			as.numeric(summary(lm(y ~ x, data = fp))$df[2])),
		"x"))
	}
postscript("~/data/foragingparasitism.eps",
	height = 3, width = 3, horizontal = FALSE,
	onefile = FALSE, paper = "special")
fp.plot()
dev.off()

# Parasitization rate from release experiment

load(url("http://anthony.darrouzet-nardi.net/data/platyrelease.Rdata"))
attach(platy.df)
midge <- deadmidge + emergedmidge
platy <- deadplaty + emergedplaty
para <- platy / (platy + midge)
para[is.nan(para)] <- 0
meanpara <- as.numeric(tapply(para, interaction(distance, release), mean))
releasemean.df <- data.frame(meanpara,
	distance = rep(c(.75, 1.5, 3, 6, 12), 3),
	release = rep(c(1,2,3), each = 5))

xyplot(meanpara ~ distance, groups = release,
	data = releasemean.df,
	t = "b", pch = c(49, 50, 51), col = "black", cex = 2, lty = 1)


load(url("http://anthony.darrouzet-nardi.net/data/platyrelease.Rdata"))
attach(platy.df)
midge <- deadmidge 
platy <- deadplaty
para <- platy / (platy + midge)
para[is.nan(para)] <- 0
meanpara <- as.numeric(tapply(para, interaction(distance, release), mean))
releasemean.df <- data.frame(meanpara,
	distance = rep(c(.75, 1.5, 3, 6, 12), 3),
	release = rep(c(1,2,3), each = 5))

xyplot(meanpara ~ distance, groups = release,
	data = releasemean.df,
	t = "b", pch = c(49, 50, 51), col = "black", cex = 2, lty = 1)


# Diffusion Model Behavior
load(url("http://anthony.darrouzet-nardi.net/data/gc.Rdata"))
gc("100", up = 1, down = 0, left = 0, right = 12)
for(d in c(1, 10*1:10))
	gc("exp(-x^2/d)", add = T)
text(1, .2, "d = 1")
text(3.1, .3, "d = 10")
text(6, .75, "d = 100")

# our data
source("~/r/gc")
gc("1-exp(-1*(x^2/9.96))", col = "red", lwd = 3,
	up = 1, down = 0, left = 0, right = 12)
	
# AND, the whole model from scratch
load(url("http://anthony.darrouzet-nardi.net/data/platyrelease.Rdata"))
platy <- platy.df$deadplaty + platy.df$emergedplaty
distance <- unique(platy.df$distance)
diffusion <- 1 - exp(-1*(distance^2/9.96))[1]  # 9.96 (d) was fit by Martha
for (x in 2:5) exp(-1*(distance^2/9.96))[x - 1] -
	exp(-1*(distance^2/9.96))[x] -> diffusion[x]
platyperplant <- tapply(platy, interaction(platy.df$distance,
	platy.df$release), sum) / rep(2^(1:5), 3)
platydata <- platyperplant / 
	rep(tapply(platyperplant, rep(1:3, each = 5), sum), each = 5)
xyplot(c(diffusion, platydata) ~ rep(distance, 4), groups = rep(1:4, each = 5),
	pch = c(16,49:51), col = "black", cex = 3)