###################################################
###################################################
###################################################
##### FUNCTION TO FORMAT DATA FOR THE WORKSHOP ANALYSIS

#########################
## 
##' .. Compute the basal area of competitor in a plot..
##'
##' .. content for \details{} ..
##' @title 
##' @param diam diameters of each individuals in cm
##' @param weights weight to compute the basal area in cm^2/m^2
##' @return basal area in cm^2/m^2 
##' @author Kunstler
BA.fun <- function(diam,weights){
((diam/2)^2)*pi*weights
}



#### 
##' .. function to compute the basal area of neighborood tree in plots ..
##'
##' .. content for \details{} ..
##' @title 
##' @param id.tree plot identifier
##' @param diam diam of tree in cm
##' @param sp species name or code
##' @param id.plot identifier of the plot
##' @param weights weights to compute basal area in cm^2/m^2 SO THE 1/AREA of teh plot (or subplot) with area in m^2
##' @param weights.full.plot weights for the whole plot to compute basal area in cm^2/m^2 or if NA use weights of the individuals (for simple plots)
##' @return data frame with tree.id and one column per species with basal area of the species (without the target tree)
##' @author Kunstler
BA.SP.FUN <-  function(id.tree,diam,sp,id.plot,weights,weight.full.plot){
# compute BA tot per species per plot
if(!(length(id.tree)==length(diam) & length(id.tree)==length(sp) & length(id.tree)==length(id.plot) & length(id.tree)==length(weights)))
    stop("length of id.tree diam,sp id.plot & weights need to be the same")
sp <- as.character(sp)
BASP <- tapply(BA.fun(diam,weights),INDEX=list(id.plot,sp),FUN=sum,na.rm=T)
print(dim(BASP))
DATA.BASP <-  data.frame(id.plot= rownames(BASP),BASP)
names( DATA.BASP) <- c("id.plot",colnames(BASP)) 

#### MERGE with indivudal tree
data.indiv <- data.frame(id.tree,sp,id.plot,diam,weights)
data.merge <-  merge(data.indiv,DATA.BASP,by="id.plot", sort = FALSE)

## substracte the BA of target species in the column of the good species (with transpose of the matrix )
t.arbre.BASP <- t(as.matrix(data.merge[,-(1:5)]))
if(!is.na(weight.full.plot)){
  t.arbre.BASP[t(outer(data.merge$sp,names(data.merge[,-(1:5)]),FUN="=="))] <-
  t.arbre.BASP[t(outer(data.merge$sp,names(data.merge[,-(1:5)]),FUN="=="))]  -BA.fun(data.merge$diam,rep(weight.full.plot,length(diam)))
  }else{
  t.arbre.BASP[t(outer(data.merge$sp,names(data.merge[,-(1:5)]),FUN="=="))] <-
  t.arbre.BASP[t(outer(data.merge$sp,names(data.merge[,-(1:5)]),FUN="=="))]  -BA.fun(data.merge$diam,data.merge$weights)
  }
data.res <- data.frame(data.merge$id.tree,t(t.arbre.BASP))
names(data.res) <- c("id.tree",colnames(BASP))
return( data.res)
}


#### function with X Y coordinates based on a neighborhood of radius R

BA.SP.FUN.XY <-  function(id.tree,xy.table,diam,sp,Rlim){
dist.mat <- as.matrix(dist(xy.table,upper=TRUE,diag=TRUE))
dist.mat[dist.mat <Rlim] <- 1
dist.mat[dist.mat >Rlim] <- 0
diag(dist.mat) <- 0
BA <- BA.fun(diam,weights=10000/(pi*Rlim^2))
BA.mat <- matrix(rep(BA,length(BA)),nrow=length(BA),byrow=TRUE)
fun.sum.sp <- function(x,sp) tapply(x,INDEX=sp,FUN=sum,na.rm=TRUE)
return(t(apply(dist.mat*BA.mat,MARGIN=1,FUN=fun.sum.sp ,sp)))
}
 

############################
## FUNCTION remove trailing white space
trim.trailing <- function (x) sub("\\s+$", "", x)

## clean species.tab
fun.clean.species.tab <- function(species.tab){
species.tab2 <- species.tab[!is.na(species.tab$Latin_name),]
 
### species IFN reformat names
## clean species names and synonyme names
species.tab2$Latin_name <- (gsub("_", " ", species.tab2$Latin_name))
species.tab2$Latin_name_syn<- (gsub("_", " ", species.tab2$Latin_name_syn))
## remove trailing white space
species.tab2$Latin_name_syn<- trim.trailing(species.tab2$Latin_name_syn)

species.clean <- species.tab2[!duplicated(species.tab2$Latin_name),
                              c("code","Latin_name","Exotic_Native_cultivated")]
return(species.clean)}


### compute quantile 99% and sd with a bootstrap

library(boot)

f.quantile <- function (x,ind,probs){quantile(x[ind],probs=probs,na.rm=TRUE)}

f.quantile.boot <-  function(i,x,fac,R,probs=0.99){
    require(boot)
    if(length(na.exclude(x[fac==i]))>0){
quant.boot <- boot(x[fac==i],f.quantile,R=R,probs=probs)
return(as.matrix(c(mean=mean(quant.boot$t),sd=sd(quant.boot$t),nobs=length(na.exclude(x[fac==i]))),ncol=3,nrow=1))
}else{
return(as.matrix(c(mean=NA,sd=NA,nobs=NA),ncol=3,nrow=1))
}
}

#######################
### function to compute number of dead per plot 
function.perc.dead <-  function(dead){
 sum(dead)/length(dead)}