# Estimating the mark mingling and mark differentiation function (Pommerening et al., 2011), Ap - 30.01.2013. Last updated on 18.02.2013.

rm(list = ls())										

cppFile <- "/Users/afs80b/Dropbox/Rcpp/"
# cppFile <- "/Users/arnepommerening/Dropbox/Rcpp/"
# cppFile <- "/Rcpp/"

library(Rcpp)
sourceCpp(paste(cppFile, "MinglingMarkDifferentiation.cpp", sep = ""))

# Randomize
set.seed(-5)

# Load data
dFile <- paste(cppFile, "Clg6.txt", sep = "")
TreeList <- read.table(dFile, header = TRUE)

# How many trees per species?
(xspecies <- tapply(TreeList$Number, TreeList$Species, length))
sum(xspecies)
length(xspecies)

# Create vector with species frequencies
speciesFreq <- 0
for (i in 1 : length(xspecies)) {
  if(i == 1)
    speciesFreq <- xspecies[[i]]
  else
    speciesFreq <- c(speciesFreq, xspecies[[i]])
}

# Sort dbhs ascendingly
dbh <- TreeList$dbh[order(TreeList$dbh, decreasing = FALSE)]

# Settings for mingling and mark differentiation functions
m <- 140 # number of steps
dd <- 0.25 # step width
h <- 2.0 # bandwidth
m * dd
kernel <- 0 # 0: Epanechnikov, 1: box kernel

# Check coordinates
min(TreeList$x)
min(TreeList$y)
max(TreeList$x)
max(TreeList$y)

# Observation window
xmax <- 98 # width
ymax <- 107 # length

# Calculate mingling and mark differentiation functions for observed data
md <- calcMingDiffFunctions(TreeList, xmax, ymax, m, dd, h, speciesFreq, dbh, kernel)

# Simulate random labelling for mark differentiation and L function
rep <- 1000
mdrep <- matrix(0.0, ncol = m, nrow = rep)
ltrep <- matrix(0.0, ncol = m, nrow = rep)
for (k in 1 : rep) {
  cat("Random labelling simulation #", k, "out of", rep, "\n")
  # Reset trees
  if(k > 1)
    rm(SimTreeList)
  SimTreeList <- read.table(dFile, header = TRUE)
  simmd <- NULL
  if(k > 1)
    rm(simmd)
  simmd <- NULL
  while(is.null(simmd) || is.null(simmd$tau)) {
    # Random permutation of tree diameters
    SimTreeList$dbh <- sample(SimTreeList$dbh)
    # Calculate mark differentiation and L function
    simmd <- calcMingDiffFunctions(SimTreeList, xmax, ymax, m, dd, h, speciesFreq, dbh, kernel)
  }
  # Store mark differentiation and L function
  mdrep[k, ] <- simmd$tau
  ltrep[k, ] <- simmd$Lt
}

# Calculate envelopes
MdLow0.05 <- apply(mdrep, MARGIN = 2, FUN = quantile, probs = 0.025) #, na.rm = TRUE)
MdUp0.05 <- apply(mdrep, MARGIN = 2, FUN = quantile, probs = 0.975) #, na.rm = TRUE)
md <- cbind(md, MdLow0.05)
md <- cbind(md, MdUp0.05)

LtLow0.05 <- apply(ltrep, MARGIN = 2, FUN = quantile, probs = 0.025) #, na.rm = TRUE)
LtUp0.05 <- apply(ltrep, MARGIN = 2, FUN = quantile, probs = 0.975) #, na.rm = TRUE)

md <- cbind(md, LtLow0.05)
md <- cbind(md, LtUp0.05)
LTmean <- apply(ltrep, MARGIN = 2, FUN = mean) #, na.rm = TRUE)
md <- cbind(md, LTmean)

# Check figure of Lt functions
plot(md$r, md$Lt - md$L, type = "l", ylim = c(-2, 2)) 
for(i in 1 : (rep)) 
  lines(md$r, ltrep[i,] - md$L, col = grey(0.7))
lines(md$r, LTmean - md$L, lty = 2, col = "red") 

# Release memory
rm(mdrep)
rm(SimTreeList)
rm(MdLow0.05)
rm(MdUp0.05)
rm(LTmean)
rm(LtLow0.05)
rm(LtUp0.05)

# Calculate deviation test (Grabarnik et al., 2011 and Diggle, 2003, p. 14)
# Include observed data as one of the simulation replications
Lts <- rbind(md$Lt, ltrep) # Lts[1,] represent the data, Lt[k,] for k = 2,..., rep + 1 are the simulations

# Calculate for all Lts[k,], k = 1, 2,..., rep + 1 (i.e. for the data and all simulations)
Lts_MeanRep <- matrix(0.0, ncol = m, nrow = rep + 1) 
for(k in 1 : (rep + 1))
  Lts_MeanRep[k, ] <- apply(Lts[-k,], MARGIN = 2, FUN = mean) # (Lts_MeanRep[1, ] is for the observed data)


# Simulate random superposition for mark mingling function
nurep <- matrix(0.0, ncol = m, nrow = rep)
lnrep <- matrix(0.0, ncol = m, nrow = rep)
# Code(s) of species group to be shifted (separated by commas) so that roughly half of the trees are moved
m1 <- c(16, 59, 60)
for (k in 1 : rep) {
  cat("Random superposition Simulation #", k, "out of", rep, "\n")
  # Random uniform point in the observation window for shifting species group
  zx <- runif(1, min = 0, max = xmax)
  zy <- runif(1, min = 0, max = ymax)
  # Reset trees
  if( k > 1)
    rm(SimTreeList)
  SimTreeList <- read.table(dFile, header = TRUE)
  simnu <- NULL
  if(k > 1)
    rm(simnu)
  simnu <- NULL
  while(is.null(simnu) || is.null(simnu$nu)) {
    # Shifting species group
    if( k > 1)
      rm(xy)
    xy <- shiftSpeciesGroup(SimTreeList, m1, zx, zy, xmax, ymax)
    SimTreeList$x <- xy$x
    SimTreeList$y <- xy$y
    # Calculate mark mingling function and L function
    simnu <- calcMingDiffFunctions(SimTreeList, xmax, ymax, m, dd, h, speciesFreq, dbh, kernel)
  }
  # Store mark differentiation and L function
  nurep[k, ] <- simnu$nu
  lnrep[k, ] <- simnu$Ln
}

# Calculate envelopes
NuLow0.05 <- apply(nurep, MARGIN = 2, FUN = quantile, probs = 0.025) #, na.rm = TRUE)
NuUp0.05 <- apply(nurep, MARGIN = 2, FUN = quantile, probs = 0.975) #, na.rm = TRUE)
md <- cbind(md, NuLow0.05)
md <- cbind(md, NuUp0.05)

LnLow0.05 <- apply(lnrep, MARGIN = 2, FUN = quantile, probs = 0.025) #, na.rm = TRUE)
LnUp0.05 <- apply(lnrep, MARGIN = 2, FUN = quantile, probs = 0.975) #, na.rm = TRUE)

md <- cbind(md, LnLow0.05)
md <- cbind(md, LnUp0.05)
LNmean <- apply(lnrep, MARGIN = 2, FUN = mean) #, na.rm = TRUE)
md <- cbind(md, LNmean)

# Check figure of Ln functions
plot(md$r, md$Ln - md$r, type = "l", ylim = c(-2, 2)) 
for(i in 1 : (rep)) 
  lines(md$r, lnrep[i,] - md$r, col = grey(0.7))
lines(md$r, LNmean - md$r, lty = 2, col = "red") 

# Release memory
rm(dbh)
rm(nurep)
rm(SimTreeList)
rm(NuLow0.05)
rm(NuUp0.05)
rm(LNmean)
rm(LnLow0.05)
rm(LnUp0.05)
rm(speciesFreq)
rm(xspecies)

# Calculate deviation test (Grabarnik et al., 2011 and Diggle, 2003, p. 14)
# Include observed data as one of the simulation replications
Lns <- rbind(md$Ln, lnrep) # Lns[1,] represent the data, Ln[k,] for k = 2,..., rep + 1 are the simulations

# Calculate for all Lns[k,], k = 1, 2,..., rep + 1 (i.e. for the data and all simulations)
Lns_MeanRep <- matrix(0.0, ncol = m, nrow = rep + 1) 
for(k in 1 : (rep + 1))
  Lns_MeanRep[k, ] <- apply(Lns[-k,], MARGIN = 2, FUN = mean) # (Lns_MeanRep[1, ] is for the observed data)

# Calculate minimum distance and cut mark functions at that point
(dm <- calcMinDistance(TreeList$x, TreeList$y))
md <- md[md$r > dm, ]
ltrep <- ltrep[md$r > dm, ]
Lts <- Lts[md$r > dm, ]
Lts_MeanRep <- Lts_MeanRep[md$r > dm, ]
lnrep <- lnrep[md$r > dm, ]
Lns <- Lns[md$r > dm, ]
Lns_MeanRep <- Lns_MeanRep[md$r > dm, ]
md$r[1]


# Calculate the deviation measure both for the observed data (k = 1) and for all simulations
dev5t <- 0
dev6t <- 0
for(k in 1 : (rep + 1)) {
  dev5t[k] <- max(abs(Lts[k,] - Lts_MeanRep[k, ]))
  dev6t[k] <- sum((Lts[k,] - Lts_MeanRep[k, ])^2)
}

(p5t.hat <- 1 - 1 / (rep + 1) * length(dev5t[dev5t[1] > dev5t])) # (or length(dev5t[dev5t[1] > dev5t[-1]]), the same since dev5t[1] is never > dev5t[1])
(p6t.hat <- 1 - 1 / (rep + 1) * length(dev6t[dev6t[1] > dev6t]))

dev5n <- 0
dev6n <- 0
for(k in 1 : (rep + 1)) {
  dev5n[k] <- max(abs(Lns[k,] - Lns_MeanRep[k, ]))
  dev6n[k] <- sum((Lns[k,] - Lns_MeanRep[k, ])^2)
}

(p5n.hat <- 1 - 1 / (rep + 1) * length(dev5n[dev5n[1] > dev5n])) # (or length(dev5n[dev5n[1] > dev5n[-1]]), the same since dev5n[1] is never > dev5n[1])
(p6n.hat <- 1 - 1 / (rep + 1) * length(dev6n[dev6n[1] > dev6n]))

# Release memory
rm(ltrep)
rm(Lts_MeanRep)
rm(lnrep)
rm(Lns_MeanRep)

par(mfrow=c(2,2))
# pdf(file = paste(cppFile, "MarkDifferentiation.pdf", sep = ""))
par(mar = c(2.5, 4, 1, 1))
plot(md$r, md$tau, las = 1, main = "",  xlab = "", ylab = "", xlim = c(0, 35), ylim = c(0, 1.5), cex.axis = 1.7, type = "l") #  ylab = expression(gamma(r)), xlab = "r [m]",
polygon(c(md$r, rev(md$r)), c(md$MdLow0.05, rev(md$MdUp0.05)), col = "lightgray", border = "lightgray")
lines(md$r, md$tau, lwd = 2, col = "red")
# abline(h = 1, lty = 2)
curve(1 + x * 0, from = md$r[1], to = 35, lty = 2, add = TRUE)
box(lwd = 2)
# dev.off()

# pdf(file = paste(cppFile, "L_tau.pdf", sep = ""))
par(mar = c(2.5, 3, 1, 1))
plot(md$r, md$Lt, las = 1, main = "",  xlab = "", ylab = "", xlim = c(0, 35), ylim = c(0, 35), lty = 1, type = "l",
cex.axis = 1.7) 
polygon(c(md$r, rev(md$r)), c(md$LtLow0.05, rev(md$LtUp0.05)), col = "lightgray", border = "lightgray")
lines(md$r, md$Lt, lwd = 2, col = "red")
lines(md$r, md$LTmean, lwd = 1, col = "green")
# abline(0, 1, lty = 2)
curve(1 * x, from = md$r[1], to = 35, lty = 2, add = TRUE)
box(lwd = 2)
# dev.off()

# pdf(file = paste(cppFile, "MarkMingling.pdf", sep = ""))
par(mar = c(2.5, 4, 1, 1))
plot(md$r, md$nu, las = 1, main = "",  xlab = "", ylab = "", xlim = c(0, 35), ylim = c(0, 1.5), cex.axis = 1.7, type = "l") #  ylab = expression(gamma(r)), xlab = "r [m]",
polygon(c(md$r, rev(md$r)), c(md$NuLow0.05, rev(md$NuUp0.05)), col = "lightgray", border = "lightgray")
lines(md$r, md$nu, lwd = 2, col = "red")
# abline(h = 1, lty = 2)
curve(1 + x * 0, from = md$r[1], to = 35, lty = 2, add = TRUE)
box(lwd = 2)
# dev.off()

# pdf(file = paste(cppFile, "L_nu.pdf", sep = ""))
par(mar = c(2.5, 3, 1, 1))
plot(md$r, md$Ln, las = 1, main = "",  xlab = "", ylab = "", xlim = c(0, 35), ylim = c(0, 35), lty = 1, type = "l",
cex.axis = 1.7) 
polygon(c(md$r, rev(md$r)), c(md$LnLow0.05, rev(md$LnUp0.05)), col = "lightgray", border = "lightgray")
lines(md$r, md$Ln, lwd = 2, col = "red")
lines(md$r, md$r, lwd = 1, col = "green")
# abline(0, 1, lty = 2)
curve(1 * x, from = md$r[1], to = 35, lty = 2, add = TRUE)
box(lwd = 2)
# dev.off()

rm(TreeList)