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

## First, let's get ask() and rNeymanScott() working:

ask = function(message = "Type in datum"){
    eval(parse(prompt = paste(message, ": ", sep="")))
}


x = ask(1)  ## now it works!

rNeymanScott = function(lambda, rmax, rcluster, 
win = owin(c(0, 1), c(0, 1)), ..., lmax = NULL) {
    win <- as.owin(win)
    frame <- bounding.box(win)
    dilated <- owin(frame$xrange + c(-rmax, rmax), frame$yrange + 
        c(-rmax, rmax))
    if (is.im(lambda) && !is.subset.owin(as.owin(lambda), dilated)) 
        stop("The window in which the image 'lambda' is defined\n\nis not large enough to contain the dilation of the window 'win'")
    parents <- rpoispp(lambda, lmax = lmax, win = dilated)
    result <- ppp(numeric(0), numeric(0), window = win)
    if (parents$n == 0) 
        return(result)
    for (i in seq(parents$n)) {
        cluster <- rcluster(parents$x[i], parents$y[i], ...)
        cluster <- ppp(cluster$x, cluster$y, window = frame)
        cluster <- cluster[, win]
        result <- superimpose(result, cluster)
    }
    return(result)
}

runifdisc = function(n, r=1, x = 0, y = 0){
theta = runif(n, min = 0, max = 2*pi)
s = sqrt(runif(n,min=0,max=r^2))
return(list(x=x+s*cos(theta),y=y+s*sin(theta)))
}


#### Simulate a Neyman-Scott model
nclust =  function(x0, y0, radius, n) {
                                  return(runifdisc(n, radius, x0, y0))
                                }
z = rNeymanScott(10, 0.2, nclust, radius=0.05, n=20)
par(mfrow=c(1,1))
plot(z,pch=".")
z1 = cbind(z$x,z$y)



####  Now, estimating variograms for geostatistical data sampled at 
####  locations (x,y)

http://www.weather.com/maps/news/forecastsummary/uscurrenttemperatures_large.html?from=wxcenter_maps
z = c(52, 52, 43, 45, 
60, 57, 33, 34, 38, 38, 
56, 60, 53, 46, 25, 26, 36, 27,
52, 35, 20, 24, 30, 30, 27, 5, 23)

x = c(.05, .05, .02, .08, 
.2, .17, .18, .19, .16, .2, 
.6, .56, .4, .38, .45, .42, .43, .42,
.8, .82, .75, .76, .63, .62, .61, .69, .81)

y = c(.30, .35, .40, .60, 
.10, .15, .38, .70, .76, .87, 
.4, .9, .17, .23, .43, .50, .62, .82,
.01, .21, .31, .42, .52, .63, .85, .99, .76)

n = length(x)

plot(c(0,1),c(0,1),type="n",xlab="x-coordinate",ylab="y-coordinate")
for(i in 1:n) text(x[i],y[i],as.character(z[i]))

zmin1 = -60
zmax1 = 140

par(mfrow=c(2,2))
data2 = data.frame(x=x,y=y,z=z)
my.ls = surf.ls(3,data2)
## my.ls is the least squares fit of a degree 3 polynomial surface to the data
my.matrix.ls = trmat(my.ls,0,1,0,1,17)  ## you can use whatever you like here, for 17
### it's going to output the results on a 17 x 17 grid
image(my.matrix.ls,zlim=c(zmin1,zmax1),col=grey(c(64:20)/64),main="LS fit, degree 3")
## the next 9 lines are just to make a legend:
zrng = zmax1 - zmin1
zmid = zmin1 + zrng/2
plot(c(0,10),c(zmid-2*zrng/3,zmid+2*zrng/3),type="n",axes=F,xlab="",ylab="")
zgrid = seq(zmin1,zmax1,length=100)
image(c(-1:1),zgrid,matrix(rep(zgrid,2),ncol=100,byrow=T),add=T,col=gray((64:20)/64))
text(2.5,zmin1,as.character(signif(zmin1,2)),cex=1)
text(2.5,zmax1,as.character(signif(zmax1,2)),cex=1)
text(2.5,zmid,as.character(signif(zmid,2)),cex=1)
text(4.5,zmid,"degrees F",srt=-90)
## GLS fit, degree 3:
my.gls = surf.gls(3,expcov,data2,d=0.1)
## my.gls is the fit of a degree 3 surface to the data by GLS,
## estimating the variogram exponentially with range d=0.1. (change the 0.1 if you'd like)
my.matrix.gls = trmat(my.gls,0,1,0,1,17)  ## again, feel free to change 17 
image(my.matrix.gls, zlim=c(zmin1,zmax1), col=grey(c(64:20)/64),
    main="GLS fit, degree 3\n with expcov(0.1)")
my.krig = prmat(my.gls,0,1,0,1,17)  ## Kriging estimates
image(my.krig,zlim=c(zmin1,zmax1), col=grey(c(64:20)/64),
    main="Kriging of deg. 3 fit\n with expcov(.1)")


### DIFFERENCES:
## to get an idea of the range of sizes of the differences, so you can 
## scale zlim in the following plots appropriately:
median(my.krig$z-my.matrix.gls$z)
sd(my.krig$z-my.matrix.gls$z)
max(abs(my.krig$z-my.matrix.gls$z))
zmin2 = -20
zmax2 = 20
par(mfrow=c(2,2))
image(my.krig$x,my.krig$y,my.krig$z-my.matrix.gls$z,zlim=c(zmin2,zmax2), col=grey(c(64:20)/64),
    main="Kriging minus GLS fit")
image(my.krig$x,my.krig$y,my.krig$z-my.matrix.ls$z,zlim=c(zmin2,zmax2), col=grey(c(64:20)/64),
    main="Kriging minus LS fit")
image(my.krig$x,my.krig$y,my.matrix.gls$z-my.matrix.ls$z,zlim=c(zmin2,zmax2), col=grey(c(64:20)/64),
    main="GLS minus LS fit")
## the next 9 lines are just to make a legend:
zrng = zmax2 - zmin2
zmid = zmin2 + zrng/2
plot(c(0,10),c(zmid-2*zrng/3,zmid+2*zrng/3),type="n",axes=F,xlab="",ylab="")
zgrid = seq(zmin2,zmax2,length=100)
image(c(-1:1),zgrid,matrix(rep(zgrid,2),ncol=100,byrow=T),add=T,col=gray((64:20)/64))
text(2.5,zmin2,as.character(signif(zmin2,2)),cex=1)
text(2.5,zmax2,as.character(signif(zmax2,2)),cex=1)
text(2.5,zmid,as.character(signif(zmid,2)),cex=1)
text(4.5,zmid,"degrees F",srt=-90)







####### VARIOGRAMS AND COVARIOGRAMS:
sill1 = var(c(z))
my.ls0 = surf.ls(0,data2)
c1 = correlogram(my.ls0,100,plot=F)  ## 100 lags -- you can change this
plot(c1$x, c1$y * sill1, xlab="distance", ylab="covariance",type="l") ## covariogram
m1 = max((1:length(c1$x))[c1$x<.5]) ## to only go up to distance 0.5
plot(c1$x[1:m1], c1$y[1:m1] * sill1, xlab="distance", ylab="covariance",type="l") 
plot(c1$x[1:m1], c1$y[1:m1], xlab="distance", ylab="correlation",type="l") ## correlogram
c2 = variogram(my.ls0,100,plot=F)  ## again, 100 lags -- you may change this.
plot(c2$x[1:m1], c2$y[1:m1], xlab="distance, h", ylab="gamma(h)",type="l",
main="semi-variogram of raw data") 
c3 = variogram(my.ls,100,plot=F)  ## fits semi-variogram to LS residuals
plot(c3$x[1:m1], c3$y[1:m1], xlab="distance, h", ylab="gamma(h)",type="l", 
    main="semi-variogram of LS residuals")
c4 = variogram(my.gls,100,plot=F) ## fits semi-variogram to GLS residuals
plot(c4$x[1:m1], c4$y[1:m1], xlab="distance, h", ylab="gamma(h)",type="l",
    main="semi-variogram of GLS residuals")



### VARIOGRAM MODELS:
par(mfrow=c(2,2))
a1 = seq(0.01,0.5,length=50)
plot(c1$x[1:m1], 2*sill1 - 2*c1$y[1:m1] * sill1, 
    xlab="distance, h", ylab="2gamma(h)",type="l",
    lty=2,main="exponential, d=5") 
lines(a1, 2*sill1 - 2*expcov(a1,5)*sill1)
plot(c1$x[1:m1], 2*sill1 - 2*c1$y[1:m1] * sill1, xlab="distance, h", ylab="2gamma(h)",type="l",
    lty=2,main="exponential, d=0.5") 
lines(a1, 2*sill1 - 2*expcov(a1,0.5)*sill1)
plot(c1$x[1:m1], 2*sill1 - 2*c1$y[1:m1] * sill1, xlab="distance, h", ylab="2gamma(h)",type="l",
    lty=2,main="exponential, d=0.1") 
lines(a1, 2*sill1 - 2*expcov(a1,0.1)*sill1)
plot(c1$x[1:m1], 2*sill1 - 2*c1$y[1:m1] * sill1, xlab="distance, h", ylab="2gamma(h)",type="l",
    lty=2,main="exponential, d=0.01") 
lines(a1, 2*sill1 - 2*expcov(a1,0.01)*sill1)

par(mfrow=c(2,2))
plot(c1$x[1:m1], 2*sill1 - 2*c1$y[1:m1] * sill1, xlab="distance, h", ylab="2gamma(h)",type="l",
    lty=2,main="Gaussian, d=5") 
lines(a1, 2*sill1 - 2*gaucov(a1,5)*sill1)
plot(c1$x[1:m1], 2*sill1 - 2*c1$y[1:m1] * sill1, xlab="distance, h", ylab="2gamma(h)",type="l",
    lty=2,main="Gaussian, d=0.5") 
lines(a1, 2*sill1 - 2*gaucov(a1,0.5)*sill1)
plot(c1$x[1:m1], 2*sill1 - 2*c1$y[1:m1] * sill1, xlab="distance, h", ylab="2gamma(h)",type="l",
    lty=2,main="Gaussian, d=0.1") 
lines(a1, 2*sill1 - 2*gaucov(a1,0.1)*sill1)
plot(c1$x[1:m1], 2*sill1 - 2*c1$y[1:m1] * sill1, xlab="distance, h", ylab="2gamma(h)",type="l",
    lty=2,main="Gaussian, d=0.01") 
lines(a1, 2*sill1 - 2*gaucov(a1,0.01)*sill1)

par(mfrow=c(2,2))
plot(c1$x[1:m1], 2*sill1 - 2*c1$y[1:m1] * sill1, xlab="distance, h", ylab="2gamma(h)",type="l",
    lty=2,main="spherical, d=5") 
lines(a1, 2*sill1 - 2*sphercov(a1,5)*sill1)
plot(c1$x[1:m1], 2*sill1 - 2*c1$y[1:m1] * sill1, xlab="distance, h", ylab="2gamma(h)",type="l",
    lty=2,main="spherical, d=0.5") 
lines(a1, 2*sill1 - 2*sphercov(a1,0.5)*sill1)
plot(c1$x[1:m1], 2*sill1 - 2*c1$y[1:m1] * sill1, xlab="distance, h", ylab="2gamma(h)",type="l",
    lty=2,main="spherical, d=0.1") 
lines(a1, 2*sill1 - 2*sphercov(a1,0.1)*sill1)
plot(c1$x[1:m1], 2*sill1 - 2*c1$y[1:m1] * sill1, xlab="distance, h", ylab="2gamma(h)",type="l",
    lty=2,main="spherical, d=0.01") 
lines(a1, 2*sill1 - 2*sphercov(a1,0.01)*sill1)

## Compare all 3 with d = 0.1:
par(mfrow=c(1,1))
plot(c1$x[1:m1], 2*sill1 - 2*c1$y[1:m1] * sill1, xlab="distance, h", ylab="2gamma(h)",type="l",
    lty=2,main="covariograms, d=0.1") 
lines(a1, 2*sill1 - 2*expcov(a1,0.1)*sill1,col=3)
lines(a1, 2*sill1 - 2*gaucov(a1,0.1)*sill1,col=4)
lines(a1, 2*sill1 - 2*sphercov(a1,0.1)*sill1,col=5)
legend(0.3,12,c("nonparametric","exponential", "Gaussian", "spherical"),
    col=c(1,3,4,5),lty=c(2,1,1,1))








http://www.montereybaywhalewatch.com/sightings/smap0303.htm

x = 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)

y = 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)

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(x,y,pch=1)


### Pseudo-Log-Likelihood.  The model is lambda_p ( z | z_1, ..., z_k) = 
## mu + alpha x + beta y + gamma SUM_{i = 1 to k} exp{-a1 D(z_i,z)} 
## where z = (x,y), and where D means distance.
## So, if a1 is positive, then there is clustering.
d1 = as.matrix(dist(cbind(x,y))) ## matrix of distances between pts

f = function(p){  
## 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)
    p0 = exp(- (p[1] + p[2]/2 + p[3]/2))
    lam = p[1] + p[2] * x + p[3] * y
    for(i in 1:n){
	for(j in c(1:n)[-i]){
	    lam[i] = lam[i] + p[4] * exp(-p[5] * d1[i,j])
	}
    }
###
    lam2 = p[1] + p[2] * simx + p[3] * simy
    for(i in 1:100){
	for(j in c(1:n)){
	    lam2[i] = lam2[i] + p[4] * exp(-p[5] * 
		sqrt((simx[i]-x[j])^2+(simy[i]-y[j])^2))
	}
    }
    cat(mean(lam2),"  ",p,"\n") ## the 1st column should be roughly n when it's done
    return(mean(lam2)-sum(log(lam)))
}
    
    
simx = runif(100)
simy = runif(100)
pstart = c(10, 10, 10, 10, 10)
z1 = optim(pstart,f,control=list(maxit=50))
pend = z1$par

### TO CHECK, COMPARE THE FOLLOWING:
f(pstart)  ## -135
f(pend)    ## -145
### the latter one should be less.
pend


### Plot 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)
z2 = 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){
z2[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(z2,z3))
zmax = max(c(z2,z3))
image(x2,y2,z2,col=gray((64:20)/64),zlim=c(zmin,zmax),add=T)
points(x,y)
######### 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(x,y)


### PLOT LAMBDA
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] * 
		sqrt((x2[i]-x[k])^2+(y2[j]-y[k])^2))
	    zz3[i,j] = zz3[i,j] + pstart[4] * exp(-pstart[5] * 
		sqrt((x2[i]-x[k])^2+(y2[j]-y[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(x,y)
######### 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(x,y)