######################################################################
# Exercise 1
######################################################################

#Load package for the analysis of spatial point patterns
library("spatstat")

#Manually create a point pattern with n=30 points by clicking
#into the window.
X <- clickppp(30)

#Test point pattern for complete spatial randomness
plot(quadratcount(X, nx=5,ny=5),add=TRUE,col=2, cex=1.5, lty=2)
quadrat.test(X,nx=5,ny=5)

#Compute p-value by Monte-Carlo simulation and compare
#the p-value with that of the above chisq-approximation
observed <- quadrat.test(X, nx=5,ny=5)$statistic
MCstat <- c(observed,sapply(1:999, function(i) {
  X <- runifpoint(100, win=window(X))
  quadrat.test(X,nx=5,ny=5)$statistic
}))
mean(MCstat >= observed)
hist(MCstat,nclass=50,freq=FALSE)
lines(rep(observed,2),c(-0.001,1),lwd=3,col=3,lty=2)
legend(x="topleft",c("Observed X^2"),col=3,lwd=3,lty=2)


######################################################################
# Exercise 2
#
# Compute the K-function for the Norwegian spruce data
######################################################################

data("spruces")
plot(Kest(spruces,correction="Ripley"),xlab="h (metres)",ylab="K(h)")
legend(x="topleft",c("Kest""CSR"),col=1:2,lty=1:2)

#Manually confirm K(h) of homogeneous Poisson process
u <- seq(0,10,length=1000)
lines(u, <<<<<REPLACE ME WITH CORRECT EXPRESSION>>>>>, lwd=3,col=2)

#Simulate under null model. 
env <- envelope(spruces,Kest)
plot(env)


######################################################################
# Exercise 3
# Simulate simple instance of Thomas process
######################################################################

library("MASS")
sigma <- 0.05
mu <- 10
nParents <- rpois(1, lambda=10 )
nOffspring <- rpois(nParents,lambda=mu)

#Simulate parents
parents <- runifpoint(nParents, win=unit.square())

#Simulate offspring. Described whats happening here?
offspring <- sapply(1:nParents, function(i) {
  p <- parents[i,]
  t(mvrnorm(nOffspring[i],mu=c(p$x,p$y),Sigma=sigma^2*diag(2)))
})
offspring <- data.frame(matrix(unlist(offspring),ncol=2,byrow=TRUE,dimnames=list(NULL,c("x","y"))))

nspall <- superimpose(parents=parents, offspring=offspring)
plot(nspall,cols=c(1,2))

#<<ACTUALLY THE ABOVE IS NOT A CLUSTER PROCESS. WHAT IS WRONG?>>


#R has a function to the above
#X <- rThomas(kappa, sigma, mu, win = owin(c(0,1),c(0,1)))