#    Copyright Katja Schiffers,  Frank M. Schurr, Katja Tielboerger, 
#    Carsten Urbach, Kirk Moloney, and Florian Jeltsch 2008, 2010

#    This program is free software: you can redistribute it and/or modify
#    it under the terms of the GNU General Public License as published by
#    the Free Software Foundation, either version 3 of the License, or
#    (at your option) any later version.

#    This program is distributed in the hope that it will be useful,
#    but WITHOUT ANY WARRANTY; without even the implied warranty of
#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#    GNU General Public License for more details.

#   For details of the GNU General Public License
#   see <http://www.gnu.org/licenses/>.


# developed under R-version 2.5.1
# updated for R-version 2.10.1
# spatstat and spatial - package needed



library(spatstat)
library(spatial)


# function for numerically calculating derivatives 
nderiv <- function(x, y)
  {
    r <- rep(NA, length(x))
       
    for(i in 2:(length(x)-3))
      {
        r[i] <- (y[i+1]-y[i-1])/(x[i+1]-x[i-1])
      }
    
    return(invisible(data.frame(x=x, dy=r)))
  }



# function for estimating the K2-index of the pattern
K2 <- function(X, nsim=399, CI=95, stoyan=0.15, stepwidth=0)
{
  lengthx <- sqrt(area.owin(X$window))/4
  if(stepwidth==0) stepwidth <- lengthx/200
  r <- seq(from=0, to=lengthx, by=stepwidth)

  res <- pcf.ppp(X, r=r, stoyan=stoyan, correction="translate")
  D <- nderiv(res$r, res$trans)
  D$dy[1:2] <- 0


  lambda <- X$n/area.owin(X$window)

  lo <- rep(0, length(res$r))
  hi <- rep(0, length(res$r))

  nrank <- ceiling(((nsim)*(1-CI/100))/2)

  simvals <- matrix(0, nrow=nsim, ncol=length(res$r))

  cat("Generating", nsim, "simulations", "\n")
  for(i in 1:nsim)
    {
      XSim <- rpoispp(lambda, win=X$window)
      SimRes <- pcf.ppp(XSim, r=r, stoyan=stoyan, correction="translate")
      DSim <- nderiv(SimRes$r, SimRes$trans)
      simvals[i,] <- DSim$dy
      cat(i, " ")
      if(i%%15 == 0) cat("\n")
    }
  
  cat("\n", "Done \n")
  
  for(R in 1:length(res$r))
    {
      sortSim <- sort(simvals[,R])
      lo[R] <- sortSim[nrank]
      hi[R] <- sortSim[nsim - (nrank-1)]
    }

  
  K2 <- data.frame(r=res$r, theo=rep(0, length(res$r)))
  
  desc <- c("distance argument r", "theoretical Poisson K2(r)")
  K2 <- fv(K2, "r", substitute(K2(r), NULL), "theo", , c(0, max(res$r)), c("r", "K2"), desc, fname = "K2")

  K2 <- bind.fv(K2, data.frame(hi=hi), "hi", "higher envelope", "hi")
  K2 <- bind.fv(K2, data.frame(lo=lo), "lo", "lower envelope", "lo")
  K2 <- bind.fv(K2, data.frame(obs=D$dy), "obs", "observations", "obs")
  attr(K2, "fmla") <- . ~ r
  unitname(K2) <- unitname(X)
  return(K2)
}



# function for creating an example pattern with overdispersion at smaller scales and attraction at larger scales
twoscalepp <- function(n, attsc=attsc, attint=attint, repsc=repsc, repint=repint)
  {

    loc <- matrix(NA, nrow=n, ncol=2)
    
    set <- FALSE
    
    loc[1,1]<- runif(1, min=0, max=20)
    loc[1,2]<- runif(1, min=0, max=20)
    
    N <- 1
    att <- runif(1, min=0, max=1)
    rep <- runif(1, min=0, max=1)
    
    while(N < n)  {
      x <- runif(1, min=0, max=20)
      y <- runif(1, min=0, max=20)
      
      
      for(i in 1:N){
        if(set==FALSE){
          dist <- sqrt(((loc[i,1]-x)^2)+((loc[i,2]-y)^2))
          
          if((dist < attsc) && (att < attint)) set <- TRUE
          if((dist > attsc) && (att > attint)) set <- TRUE
        }
      }
      for(j in 1:N){
        dist <- sqrt(((loc[j,1]-x)^2)+((loc[j,2]-y)^2))
        
        if((dist < repsc) && (rep < repint)) set <- FALSE
      }
      
      if(set==TRUE){
        att <- runif(1, min=0, max=1)
        rep <- runif(1, min=0, max=1)
        
        
        N <- N+1
        loc[N,1] <- x
        loc[N,2] <- y
        cat(N, " ")
        if(N%%10==0) cat("\n")
        set <- FALSE
      }
      
    }
    pp <- ppp(loc[,1], loc[,2], window=owin(c(0,20), c(0,20)))
    return(pp) 
  }



# generating a pattern with overdispersion at smaller scales and attraction at larger scales
# arguments:
# n = number of points
# attsc = scale of attraction
# attint = intensity of attraction
# resps = scale of repulsion
# repint = intensity of repulsion

two <- twoscalepp(150, attsc=0.9, attint=0.95, repsc=0.6, repint=0.75)
plot(two, pch=20, cex=0.8, xlab="", ylab="")


# analysing pattern with Ripley's K-function and plotting the L-function

twoK <- envelope(two, fun=Kest)
plot(twoK, sqrt(./pi)-r ~ r, xlab="scale(r)",ylab="L(r)", main="Ripley's K")

# analysing pattern with pair correlation function with kernel width parameter 'stoyan' = 0.15

twoPCF <- envelope(two, fun=pcf.ppp, stoyan=0.15)
plot(twoPCF, xlab="scale(r)",  ylab="g(r)", main="Pair correlation function")

# analysing pattern with K2-index with kernel width parameter 'stoyan' = 0.15

twoK2 <- K2(two, nsim=99, stoyan=0.15)
plot(twoK2, main="K2-index", xlab="scale(r)", ylab="K2(r)", legendpos="bottomright")