# Estimating the partial pair and the mark connection function (Illian et al., 2008), Ap - 25.10.2013. Last updated on 25.10.2013.
# Species codes should be numeric.

rm(list = ls())										

cppFile <- "/Users/pommeren/Dropbox/Rcpp/"

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

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

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

(xspecies <- xspecies[order(xspecies, decreasing = TRUE)])

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

# Create vector with species codes of the 3 most abundant species
(speciesCodes <- as.numeric(row.names(xspecies[1:3])))

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

# Observation window
xmax <- 100 # width
ymax <- 100 # length

# Settings for partial pair and mark connection functions
m <- 200 # number of steps
dd <- 0.25 # step width
h <- 3.0 # bandwidth
m * dd
kernel <- 0 # 0: Epanechnikov, 1: box kernel

# Calculate partial pair and the mark connection functions for observed data
ppmc <- calcPartialPairMarkConnection(TreeList, xmax, ymax, speciesFreq, speciesCodes, m, dd, h, kernel)

par(mfrow = c(1, 2))
# pdf(file = paste(cppFile, "gmm.pdf", sep = ""))
par(mar = c(2.5, 2.5, 1, 1))
plot(ppmc$r, ppmc$g11, las = 1, main = "",  xlab = "", ylab = "", xlim = c(0, m * dd), ylim = c(0, 4), cex.axis = 1.7, type = "l", col = "black") 
lines(ppmc$r, ppmc$g12, type = "l", col = "red")
lines(ppmc$r, ppmc$g13, type = "l", col = "blue")
lines(ppmc$r, ppmc$g22, type = "l", col = "green")
lines(ppmc$r, ppmc$g23, type = "l", col = "purple")
lines(ppmc$r, ppmc$g33, type = "l", col = "orange")
curve(1 + x * 0, from = ppmc$r[1], to = m * dd, lty = 2, add = TRUE)
box(lwd = 2)
# dev.off()

par(mar = c(2.5, 4, 1, 1))
plot(ppmc$r, ppmc$p11, las = 1, main = "",  xlab = "", ylab = "", xlim = c(0, 5), ylim = c(0, 0.08), cex.axis = 1.7, type = "l", col = "black") 
lines(ppmc$r, ppmc$p12, type = "l", col = "red")
lines(ppmc$r, ppmc$p13, type = "l", col = "blue")
lines(ppmc$r, ppmc$p22, type = "l", col = "green")
lines(ppmc$r, ppmc$p23, type = "l", col = "purple")
lines(ppmc$r, ppmc$p33, type = "l", col = "orange")
# curve(1 + x * 0, from = ppmc$r[1], to = m * dd, lty = 2, add = TRUE)
box(lwd = 2)


rm(TreeList)