# Make sure you have the pacman library installed
require(pacman)
# once pacman  is installed you can run the following lines
# if the packages are not downloaded, they will be downloaded 
# and installed
p_load(adehabitat)
p_load(ks)
p_load(ggplot2)

###

# Calculate 3D volumes and overlap between 3 fish

# set your seed
set.seed(125)

# number of data points to simulate
n = 1000

# simulate 3 different "fish" pathways
a = simm.crw(date=1:n, h = 1, r = 0,
                x0=c(0,0), id="A", burst=1,
                typeII=TRUE)

fish1 = data.frame(x = a[[1]]$x,
                   y = a[[1]]$y,
                   z = rnorm(n,-3,0.5))

b = simm.crw(date=1:n, h = 0.5, r = 0,
             x0=c(0,0), id="B", burst=1,
             typeII=TRUE)

fish2 = data.frame(x = b[[1]]$x,
                   y = b[[1]]$y,
                   z = rnorm(n,-8,0.5))

c = simm.crw(date=1:n, h = 1, r = 0,
             x0=c(0,0), id="C", burst=1,
             typeII=TRUE)

fish3 = data.frame(x = c[[1]]$x,
                   y = c[[1]]$y,
                   z = rnorm(n,-3,0.5))
rm(a,b,c)


ggplot() + geom_path(data = fish1, aes(x,y),col = 1) +
  geom_path(data = fish2, aes(x,y),col = 2) + 
  geom_path(data = fish3, aes(x,y),col = 3) + 
  theme_bw()

# or if you prefer base graphics with points
# undo comments by highlighting area and pressing 'Ctrl + Shift + c'
# plot(range(fish1$x,fish2$x,fish3$x),range(fish1$y,fish2$y,fish3$y),pch="")
# points(fish1$x,fish1$y,pch = 19)
# points(fish2$x,fish2$y,pch = 19,col="blue")
# points(fish3$x,fish3$y,pch = 19,col="red")


#####
#
### This section is adapted from Simpfendorfer et al. 2012
###
# Simpfendorfer, C. A., Olsen, E. M., Heupel, M. R., & Moland, E. (2012).
# Three-dimensional kernel utilization distributions improve estimates
# of space use in aquatic animals. Canadian Journal of Fisheries
# and Aquatic Sciences, 69(3), 565-572.
###
#
#####

Ha = Hpi(fish1, binned=TRUE)
Hb = Hpi(fish2, binned=TRUE)
Hc = Hpi(fish3, binned=TRUE)
# 

# Set the minimum and maximum x,y,and z coordinates so that all
# "animals" are on the same scale
minX = min(fish1$x, fish2$x, fish3$x)
minY = min(fish1$y, fish2$y, fish3$y)
minZ = min(fish1$z, fish2$z, fish3$z)

maxX = max(fish1$x, fish2$x, fish3$x)
maxY = max(fish1$y, fish2$y, fish3$y)
maxZ = max(fish1$z, fish2$z, fish3$z)

fhata = kde(fish1, H = Ha, xmin = c(minX,minY,minZ),
            xmax = c(maxX,maxY,maxZ))
fhatb = kde(fish2, H = Hb, xmin = c(minX,minY,minZ),
            xmax = c(maxX,maxY,maxZ))
fhatc = kde(fish3, H = Hc, xmin = c(minX,minY,minZ),
            xmax = c(maxX,maxY,maxZ))


# Pull out the 'heights' using the contourLevels function
# see ?contourLevels for more info
cta = contourLevels(fhata, cont = 50, approx = TRUE)
ctb = contourLevels(fhatb, cont = 50, approx = TRUE)
ctc = contourLevels(fhatc, cont = 50, approx = TRUE)

# In a 2D world areas are calculated in pixels (width * length)
# In a 3D world volumes are calculated in voxels (height * length * height)

# volume of fish1
# calculate the volume of the voxel
vol.voxela = prod(sapply(fhata$eval.points, diff)[1,])
# calculate how many of the voxels exceed the value returned 
# from the countourLevels function above
no.voxela = sum(fhata$estimate > cta)
# Calvulate the discrete volume of space use
vol50a = no.voxela * vol.voxela
vol50a

# volume of fish2
vol.voxelb = prod(sapply(fhatb$eval.points, diff)[1,])
no.voxelb = sum(fhatb$estimate > ctb)
vol50b = no.voxelb * vol.voxelb
vol50b

# volume of fish2
vol.voxelc = prod(sapply(fhatc$eval.points, diff)[1,])
no.voxelc = sum(fhatc$estimate > ctc)
vol50c = no.voxelc * vol.voxelc
vol50c


# Since we set our xmin and xmax values to be the same
# among all fish, we can directly compare the overlap 
# between individuals

# Overlap of fish1 with fish2
vol.voxela = prod(sapply(fhata$eval.points, diff)[1,])
no.voxela = sum(fhata$estimate > cta & fhatb$estimate > ctb)
overlap12 = no.voxela * vol.voxela
# no overlap
overlap12


# Overlap of fish2 with fish3
vol.voxelb = prod(sapply(fhatb$eval.points, diff)[1,])
no.voxelb = sum(fhatb$estimate > ctb & fhatc$estimate > ctc)
overlap23 = no.voxela * vol.voxela
#no overlap
overlap23


# Overlap of fish1 with fish3
vol.voxela = prod(sapply(fhata$eval.points, diff)[1,])
no.voxela = sum(fhata$estimate > cta & fhatc$estimate > ctc)
overlap13 = no.voxela * vol.voxela
# some overlap
overlap13


# and plot the points
p_load(rgl)
plot3d(fish1$x, fish1$y, fish1$z, alpha = 0.4, xlab = "", ylab = "",
       zlab = "",aspect = c(1,1,0.5))
plot3d(fish2$x, fish2$y, fish2$z, alpha = 0.4, add = T, col="blue")
plot3d(fish3$x, fish3$y, fish3$z, alpha = 0.4, add = T, col="red")



# plot the volumes using rgl and stretch the depth a bit
# to show the overlap a bit more clearly
# You can change the cont value to 95 to see what 
# the 95% contours would look like in 3D space
cont = 50
plot(fhata, cont = cont, colors = "black", drawpoints = F,
     xlab = "", ylab = "", zlab = "", 
     axes = F, box = F, alphavec = 0.3,
     xlim = range(fish1$x, fish2$x, fish3$x),
     ylim = range(fish1$y, fish2$y, fish3$y),
     zlim = range(fish1$z, fish2$z, fish3$z),
     aspect = c(1,1,1.5))
plot(fhatb, cont = cont, colors = "red", drawpoints = F,
     xlab = "", ylab = "", zlab = "", size = 2, 
     add = T, axes = F, box = F, alphavec = 0.3)
plot(fhatc, cont = cont, colors = "blue", drawpoints = F,
     xlab = "", ylab = "", zlab = "", size = 2, 
     add = T, axes = F, box = F, alphavec = 0.3)




#####
#
### Randomization Test
#
#####

randomization.test  = function(var1,  var2,  reps = 10000){
  s = sqrt((((length(var1)-1)*sd(var1)^2) + ((length(var2)-1)*sd(var2)^2)) / (length(var1)+length(var2) - 2))
  # calculate the T statistic
  T = (mean(var1) - mean(var2)) / (s * sqrt(1/length(var1) + 1/length(var2)))
  # the null hypothesis is that there is no difference between RT overlap with LT
  # The alternative hypothesis is that RT overlap with RT is larger than RT with LT
  p.value = pt(-abs(T), df = length(var1)+length(var2)-2)
  D = mean(var1) - mean(var2)
  allVars = c(var1, var2)
  for(i in 2:reps){
    newall = sample(allVars)
    newVar1 = newall[1:length(var1)]
    newVar2 = newall[(length(var1) + 1):length(newall)]
    s = sqrt((((length(newVar1)-1)*sd(newVar1) ^ 2)+((length(newVar2)-1)*sd(newVar2)^2))/(length(newVar1)+length(newVar2)-2))
    T = (mean(newVar1) - mean(newVar2))/(s * sqrt(1/length(newVar1) + 1/length(newVar2)))
    p.value[i] = pt(-abs(T), df = length(allVars) - 2)
    D[i] = mean(newVar1) - mean(newVar2)
  }
  par(mfrow = c(1,1))
  hist(D, breaks = 50)
  abline(v = D[1], col = "red")
  return(as.data.frame(list(Difference = length(D[D >= D[1]])/length(D),
                            p.value = p.value[1],
                            numberVar1 = length(var1),
                            numberVar2 = length(var2),
                            meanVar1 = mean(var1),
                            meanVar2 = mean(var2))))
}


set.seed(3)
var1 = rnorm(15, 100, 5)
var2 = rnorm(17, 97, 3)

#check out means of var1 and var2
mean(var1)
# 99.04699
mean(var2)
# 96.5791

ggplot() + geom_density(aes(var1,fill = "black", lty = NA),alpha = 0.5) +
  geom_density(aes(var2,fill = "blue", lty = NA), alpha = 0.5) + xlab("Values")

x = randomization.test(var1, var2)

x$Difference
# Difference between means are statistically significant p-value < 0.05


set.seed(3)
var1 = rnorm(15, 100, 5)
var2 = rnorm(17, 99, 3)
ggplot() + geom_density(aes(var1,fill = "black",lty = NA),alpha = 0.5) +
  geom_density(aes(var2,fill = "blue",lty = NA),alpha = 0.5) + xlab("Values")

x = randomization.test(var1, var2)
x$Difference
# Difference between means are not statistically significant p-value > 0.05



#####
#
### Simulate 2D buffer around cage
#
#####

cageCenter = data.frame(x = 0, y = 0)
circle = seq(0,2*pi,length = 100)
radius = 25
coords = as.data.frame(t(rbind(cageCenter$x + sin(circle) * radius,
                               cageCenter$y + cos(circle) * radius)))
names(coords) = c("x","y")
plot(coords)
# and plot the center of the cage
points(cageCenter$x, cageCenter$y, pch = 19, cex = 3)


#####
#
### Simulate 3D near cage habitat
#
#####

createPointsX = runif(10000, min = cageCenter$x - radius, max = cageCenter$x + radius)
createPointsY = runif(10000, min = cageCenter$y - radius, max = cageCenter$y + radius)
createPointsZ = runif(10000, min = -15, max = 0)
cageSim = data.frame(x = createPointsX, y = createPointsY, depth = createPointsZ)


# use the inpoly function in the GEOmap library to retain points within limits of created circle
p_load(GEOmap)
inp = inpoly(cageSim$x, cageSim$y, coords)

# plot the x & y locs of the simulated data
plot(cageSim$x, cageSim$y, pch = 19)

# change the colour of the points inside the polygon
plot(cageSim$x, cageSim$y, pch = 19, col = inp + 1)

# see the points in 3d
plot3d(cageSim$x, cageSim$y, cageSim$depth, col = inp + 1)

# remove the points that lie outside the polygon
cageSim = cageSim[inp == 1,]

# and re-plot the points
plot(cageSim$x, cageSim$y, pch = 19)

plot3d(cageSim$x, cageSim$y, cageSim$depth)

# You can use the new data to create a new object and measure overlap