for.multi.R

for.multi <- function(comm_data,site1,site2,trait_data,res,bandwith.fac = 1, kde.bandwidth.type = "default"){

  first <- comm_data[site1, comm_data[site1,] > 0]
  second <- comm_data[site2, comm_data[site2,] > 0]
  trait1 <- trait_data[names(first),]
  trait2 <- trait_data[names(second),]

  if(kde.bandwidth.type=="default"){
    kde.bandwidth1 <-  estimate_bandwidth(trait1)
    kde.bandwidth2 <-  estimate_bandwidth(trait2)
  }else if(kde.bandwidth.type=="uniform"){
    trait12 <- rbind(trait1,trait2)
    kde.bandwidth1 <- kde.bandwidth2 <- estimate_bandwidth(trait12)/2
  }
  
  ##compute hypervolumes with package hypervolume
  hypervolume1 = hypervolume_gaussian(trait1,name='trait1', verbose = FALSE, kde.bandwidth = kde.bandwidth1)
  hypervolume2 = hypervolume_gaussian(trait2,name='trait2', verbose = FALSE, kde.bandwidth = kde.bandwidth2)
  
  lims <- c(range(c(hypervolume1@RandomPoints[,1],hypervolume2@RandomPoints[,1])))
  for(i in 2:ncol(hypervolume1@RandomPoints))
    lims <- c(lims,range(c(hypervolume1@RandomPoints[,i],hypervolume2@RandomPoints[,i])))

  ##2dimensions - compute the KDEs from the random points generated by hypervolume_gaussian()
  hvtest1.kkde.multi <- kde(x=hypervolume1@RandomPoints,H=diag(estimate_bandwidth(hypervolume1@Data)/bandwith.fac),xmin=lims[c(1,3)],xmax=lims[c(2,4)],gridsize=c((lims[2]-lims[1])/res,(lims[4]-lims[3])/res))
  hvtest2.kkde.multi <- kde(x=hypervolume2@RandomPoints,H=diag(estimate_bandwidth(hypervolume2@Data)/bandwith.fac),xmin=lims[c(1,3)],xmax=lims[c(2,4)],gridsize=c((lims[2]-lims[1])/res,(lims[4]-lims[3])/res))
  
  ##rescale the KDEs
  hvtest1.kkde.multi$estimate=hvtest1.kkde.multi$estimate/max(hvtest1.kkde.multi$estimate)
  hvtest2.kkde.multi$estimate=hvtest2.kkde.multi$estimate/max(hvtest2.kkde.multi$estimate)

  ##Compute the indices using the kernel integrals
  hv.bothT.multi <- array(c(hvtest1.kkde.multi$estimate,hvtest2.kkde.multi$estimate),c(dim(hvtest1.kkde.multi$estimate),2))
  hv.bothT.min.multi <- apply(hv.bothT.multi,1:length(dim(hvtest1.kkde.multi$estimate)),min)
  hv.bothT.max.multi <- apply(hv.bothT.multi,1:length(dim(hvtest1.kkde.multi$estimate)),max)
  Jac <- 1-(sum(hv.bothT.min.multi)/sum(hv.bothT.max.multi))
  Turn <- 2*(min(sum(hvtest1.kkde.multi$estimate),sum(hvtest2.kkde.multi$estimate))-sum(hv.bothT.min.multi))/sum(hv.bothT.max.multi)
  Nest <- Jac-Turn
  
  output <- list()
  output$indices <- c(Jac,Turn,Nest)
  output$kernel1 <- hvtest1.kkde.multi
  output$kernel2 <- hvtest2.kkde.multi
  output$points1 <- hypervolume1@RandomPoints
  output$points2 <- hypervolume2@RandomPoints
  
  return(output)

}