library(sm)
library(spatstat)
library(splancs)
library(spatial)
library(geoR)


## 1) Enter Data
## http://www.montereybaywhalewatch.com/sightings/smap0303.htm

x1 = c(0.1,0.77,.97,.99,
.96,.9,.92,.81,.78,.77,
.76,.755,.75,
.757,.754,.745,.74,
.732,.734,.736,
.71,.715,.715,
.68,
.72,.725,.7,.707,.712,.71,
.75,.76,.765,.755,.757,.756,
.77,.78,.73,.73,
.83)

y1 = c(.92,.9,.97,.99,
.7,.66,.59,.65,.7,.64,
.69,.68,.675,
.63,.64,.638,.625,
.65,.645,.648,
.64,.641,.646,
.65,
.61,.6,.6,.595,.59,.583,
.575,.57,.565,.56,.555,.55,
.555,.569,.54,.541,
.35)

z1 = c(7,3,5,8,7,
2,7,5,4,5,
3,3,9,4,5,
3,2,4,1,8,
3,7,6,8,4,
4,4,8,9,8,
4,8,5,7,2,
4,3,5,8,9,
1)

z2 = (z1 - min(z1))/(max(z1)-min(z1)) ## rescaled to [0,1].

par(mfrow=c(1,1))
plot(c(0,1),c(0,1),type="n",xlab="x-coordinate",ylab="y-coordinate",
main="whales and dolphins")
points(x1,y1,pch=1,cex=1+3*z2)

n = length(x1)


### 2) Marked J-function:
b2 = as.ppp(cbind(x1,y1), W = c(0,1,0,1))   
b2$marks = z1
b2$n = n 
## the above convert the points into a "marked ppp" object, 
## using as a window [0,1] x [0,1]
par(mfrow=c(1,1))
jm4 = Jmulti(b2, b2$marks < 4.5, b2$marks > 5.5)
plot(jm4)
## or to control the plot yourself, try:
plot(jm4$r[jm4$r<.5],jm4$rs[jm4$r<.5],xlab="h",ylab="J(h)",type="l",lty=1) 
lines(jm4$r[jm4$r<.5],jm4$theo[jm4$r<.5],lty=2) 
legend(.2,.2,lty=c(1,2),legend=c("data","Poisson"))

## try playing around with this, for different classes of marks instead 
## of < 4.5, and > 5.5


######### 3) Marked K-function & L-function:
km4 = Kmulti(b2, b2$marks < 4.5, b2$marks > 5.5)
plot(km4$r[km4$r<.3],km4$border[km4$r<.3],xlab="h",ylab="K(h)",type="l",lty=1) 
lines(km4$r[km4$r<.3],km4$theo[km4$r<.3],lty=2) 
legend(.2,.2,lty=c(1,2),legend=c("data","Poisson"))
Lm4 = sqrt(km4$border[km4$r<.3]/pi)-km4$r[km4$r<.3]
plot(c(0,.3),range(Lm4),type="n",xlab="lag, h",ylab="L4(h) - h")
points(km4$r[km4$r<.3],Lm4,pch="*")
lines(km4$r[km4$r<.3],Lm4)
abline(h=0,lty=2)


### THEORETICAL BOUNDS for L-function
## bounds = 1.96 * sqrt(2*pi*A) * h / E(N), where 
## A = area of space, and 
## E(N) = expected # of pts of type j in the space  (approximated here using
## the observed # of pts of type j
s = km4$r[km4$r<.3]
L4upper = 1.96 * sqrt(2*pi*1*1) * s / sum(b2$marks > 5.5)
L4lower = -1.0 * L4upper
lines(s,L4upper,lty=3)
lines(s,L4lower,lty=3)



###  4) Kernel smoothing, using bandwidth h
par(mfrow=c(1,2))
h = .05 
## h should usually be a fraction (roughly 1/4 or so) 
## of the range of your x-coordinates or y-coordinates
## you can use bw.nrd0(x1) or bw.nrd0(y1) as a guide
n1 = 20
mygrid1 = seq(0,1,length=n1)
mygrid2 = seq(0,1,length=n1)
a1 = matrix(0,ncol=n1,nrow=n1)
for(i in 1:n1){
for(j in 1:n1){
    a1[i,j] = sum(z1 * dnorm(
	    sqrt((mygrid1[i]-x1)^2 + (mygrid2[j]-y1)^2),sd=h)) 
    # / sum(dnorm(sqrt((mygrid1[i]-x1)^2 + (mygrid2[j]-y1)^2),sd=h))
}
}

image(mygrid1,mygrid2,a1,xlab="x",ylab="y",zlim=range(a1),
    col=grey(c(64:20)/64))
points(x1,y1,pch=1,cex=1+3*z2)

## legend
x = a1
zmin = min(x)
zmax = max(x)
zrng = zmax - zmin
zmid = zmin + zrng/2
plot(c(0,10),c(zmid-2*zrng/3,zmid+2*zrng/3),type="n",axes=F,xlab="",ylab="")
zgrid = seq(zmin,zmax,length=100)
## zgrid = vector of 100 equally-spaced numbers spanning range of the values.
image(c(-1:1),zgrid,matrix(rep(zgrid,2),ncol=100,byrow=T),add=T,
    zlim=range(x),col=grey((64:20)/64))
text(2.5,zmin,as.character(signif(zmin,2)),cex=1)
text(2.5,zmax,as.character(signif(zmax,2)),cex=1)
text(2.5,zmid,as.character(signif(zmid,2)),cex=1)
text(4.5,zmid,"Values",srt=-90)


##### 5) Quadrat totals

x = matrix(0,ncol=10,nrow=10)
for(i in 1:10){
for(j in 1:10){
for(k in 1:n){
if((x1[k]<i/10) && (x1[k] >= (i-1)/10) &&
(y1[k]<j/10) && (y1[k] >= (j-1)/10)) x[i,j] = x[i,j] + z1[k]
}
}}
## can check that sum(x) should = sum(z1)

##### Plot the quadrat counts
par(mfrow=c(1,2))  ## makes a 1x2 grid of plots on the graphic screen
plot(c(0,1),c(0,1),type="n",xlab="x-coordinate",ylab="y-coordinate",
main="Quadrat counts of \n Whales & Dolphins")
image(x=c(0:9)/10+.05,y=c(0:9)/10+.05,z=x,col=grey(c(64:20)/64),add=T)
points(x1,y1,pch=1,cex=1+3*z2)
######### LEGEND:
zmin = min(x)
zmax = max(x)
zrng = zmax - zmin
zmid = zmin + zrng/2
plot(c(0,10),c(zmid-2*zrng/3,zmid+2*zrng/3),type="n",axes=F,xlab="",ylab="")
zgrid = seq(zmin,zmax,length=100)
## zgrid = vector of 100 equally-spaced numbers spanning range of the values.
image(c(-1:1),zgrid,matrix(rep(zgrid,2),ncol=100,byrow=T),add=T,col=gray((64:20)/64))
text(2.5,zmin,as.character(signif(zmin,2)),cex=1)
text(2.5,zmax,as.character(signif(zmax,2)),cex=1)
text(2.5,zmid,as.character(signif(zmid,2)),cex=1)
text(4.5,zmid,"Values",srt=-90)

### 6) Fitting a Pseudo-Likelihood model  
## I'm using the model lambda_p ( s | s_1, ..., s_k) = 
## mu + alpha x + beta y + gamma SUM_{i = 1 to k} exp{-a1 * z_i * D(s_i,s)} 
## where s = (x,y), and where D means distance.
## So, if gamma is positive, then there is clustering; otherwise inhibition
d1 = as.matrix(dist(cbind(x1,y1))) ## matrix of distances between pts

f = function(p){  
## returns the negative pseudo log-likelihood
## p = (mu,alpha,beta,gamma,a1)
    if(p[1] < 0) return(99999)
    if(p[1] + p[2] < 0) return(99999)
    if(p[1] + p[3] < 0) return(99999)
    if(p[1] + p[2] + p[3] < 0) return(99999)
    if(p[4] < 0) return(99999)
    lam = p[1] + p[2] * x1 + p[3] * y1
    for(i in 1:n){
	for(j in c(1:n)[-i]){
	    lam[i] = lam[i] + p[4] * exp(-p[5] * z1[i] * d1[i,j])
	}
    }
    if (min(lam) < 0) return (99999)
    simx = runif(155)
    simy = runif(155)
    lam2 = p[1] + p[2] * simx + p[3] * simy
    for(i in 1:155){
	for(j in c(1:n)){
	    lam2[i] = lam2[i] + p[4] * exp(-p[5] * z1[j] * 
		sqrt((simx[i]-x1[j])^2+(simy[i]-y1[j])^2))
	}
    }
    cat(mean(lam2),mean(lam2)-sum(log(lam)), "  ",p,"\n") ## the 1st column should be roughly n when it's done
    return(mean(lam2)-sum(log(lam)))
}
    
pstart = c(10, .1, .1, .1, 1)
fit1 = optim(pstart,f,control=list(maxit=200))
pend = fit1$par

### TO CHECK, COMPARE THE FOLLOWING:
f(pstart)  
f(pend)    
### the latter one should be less.
pend 
## Run maximum likelihood again, now starting where you left off.
fit2 = optim(pend, f, control=list(maxit=200), hessian=T)
pend = fit2$par ## interpret these parameters!!!
b3 = sqrt(diag(solve(fit2$hess))) ## interpret these standard errors!!!
pend/b3

### 7) Plot the Model's Background Rate
par(mfrow=c(1,3))
plot(c(0,1),c(0,1),type="n",xlab="x-coordinate",ylab="y-coordinate",
main="background rate")
x2 = seq(0.05,0.95,length=10)
y2 = seq(0.05,0.95,length=10)
zz2 = matrix(rep(0,(10*10)),ncol=10)
z3 = matrix(rep(0,(10*10)),ncol=10)
for(i in 1:10){
for(j in 1:10){
zz2[i,j] = pend[1] + pend[2]*x2[i] + pend[3]*y2[j]
z3[i,j] = pstart[1] + pstart[2]*x2[i] + pstart[3]*y2[j]
}}
zmin = min(c(zz2,z3))
zmax = max(c(zz2,z3))
image(x2,y2,zz2,col=gray((64:20)/64),zlim=c(zmin,zmax),add=T)
points(x1,y1,pch=1,cex=1+3*z2)
######### LEGEND:
zrng = zmax - zmin
zmid = zmin + zrng/2
plot(c(0,10),c(zmid-2*zrng/3,zmid+2*zrng/3),type="n",axes=F,xlab="",ylab="")
zgrid = seq(zmin,zmax,length=100)
## zgrid = vector of 100 equally-spaced numbers spanning range of the values.
image(c(-1:1),zgrid,matrix(rep(zgrid,2),ncol=100,byrow=T),add=T,col=gray((64:20)/64))
text(2.5,zmin,as.character(signif(zmin,2)),cex=1)
text(2.5,zmax,as.character(signif(zmax,2)),cex=1)
text(2.5,zmid,as.character(signif(zmid,2)),cex=1)
text(4.5,zmid,"Values",srt=-90)
plot(c(0,1),c(0,1),type="n",xlab="x-coordinate",ylab="y-coordinate",
main="original guess")
image(x2,y2,z3,col=gray((64:20)/64),zlim=c(zmin,zmax),add=T)
points(x1,y1,pch=1,cex=1+3*z2)


### 8) PLOT LAMBDA_p
par(mfrow=c(1,3))
plot(c(0,1),c(0,1),type="n",xlab="x-coordinate",ylab="y-coordinate",
main="lambda_p")
x2 = seq(0.05,0.95,length=10)
y2 = seq(0.05,0.95,length=10)
zz2 = matrix(rep(0,(10*10)),ncol=10)
zz3 = matrix(rep(0,(10*10)),ncol=10) 
for(i in 1:10){
    for(j in 1:10){
	zz2[i,j] = pend[1] + pend[2] * x2[i] + pend[3] * y2[j]
	zz3[i,j] = pstart[1] + pstart[2] * x2[i] + pstart[3] * y2[j]
    for(k in c(1:n)){
	    zz2[i,j] = zz2[i,j] + pend[4] * exp(-pend[5] * z1[k] *
		sqrt((x2[i]-x1[k])^2+(y2[j]-y1[k])^2))
	    zz3[i,j] = zz3[i,j] + pstart[4] * exp(-pstart[5] * z1[k] *
		sqrt((x2[i]-x1[k])^2+(y2[j]-y1[k])^2))
	}
    }
}
zmin = min(c(zz2,zz3))
zmax = max(c(zz2,zz3))
image(x2,y2,zz2,col=gray((64:20)/64),zlim=c(zmin,zmax),add=T)
points(x1,y1,pch=1,cex=1+3*z2)
######### LEGEND:
zrng = zmax - zmin
zmid = zmin + zrng/2
plot(c(0,10),c(zmid-2*zrng/3,zmid+2*zrng/3),type="n",axes=F,xlab="",ylab="")
zgrid = seq(zmin,zmax,length=100)
## zgrid = vector of 100 equally-spaced numbers spanning range of the values.
image(c(-1:1),zgrid,matrix(rep(zgrid,2),ncol=100,byrow=T),add=T,col=gray((64:20)/64))
text(2.5,zmin,as.character(signif(zmin,2)),cex=1)
text(2.5,zmax,as.character(signif(zmax,2)),cex=1)
text(2.5,zmid,as.character(signif(zmid,2)),cex=1)
text(4.5,zmid,"Values",srt=-90)
plot(c(0,1),c(0,1),type="n",xlab="x-coordinate",ylab="y-coordinate",
main="original guess")
image(x2,y2,zz3,col=gray((64:20)/64),zlim=c(zmin,zmax),add=T)
points(x1,y1,pch=1,cex=1+3*z2)