#rm(list=ls(all=TRUE))
#dev.off()

library(data.table)
library(ggplot2)
library(ape)
library(vegan)
library(geometry)

#setwd("D:/Data/LPJmL-FIT/Europe/Historical")

filename<-"trait_ind_europe_1951_20"
weight_bool<-F
quant<-0.9
quant_bool<-T
#cells<-1:4
cells<-1:5095

path.in<-path.out<-getwd()

data<-fread(paste(path.in,"/",filename,".txt",sep=""),blank.lines.skip=T,select="Year")

numlines<-table(data$Year)
print(numlines)
years<-as.numeric(names(numlines))

rm(data)

calc.divInd<-function(data,weight=F){
  if(nrow(data)>10){
    #extract data
    dat<-cbind(data$SLA,data$Wooddens,data$LL,data$Height)
    if(weight){
      wi<-data$Biomass
      wi<-wi/sum(wi)
     }
    #normalization
    SLA<-c(0.005,0.07)*455
    LL<-c(0.2,12)
    WD<-c(0.7e5,6.5e5)*(10^-6)/0.455
    Height<-c(5,50)

    dat[,1]<-(dat[,1]-SLA[1])/(SLA[2]-SLA[1])
    dat[,2]<-(dat[,2]-WD[1])/(WD[2]-WD[1])
    dat[,3]<-(dat[,3]-LL[1])/(LL[2]-LL[1])
    dat[,4]<-(dat[,4]-Height[1])/(Height[2]-Height[1])

 ### Functional Richness ###

    richness<-convhulln(dat,options = "FA")$vol #calculate volume

   ### Functional Divergence ###
    hull<-convhulln(dat)
    vertices<-unique(as.vector(hull))
    hull<-array(NA,dim=c(length(vertices),4))

    hull[,1]<-dat[vertices,1]
    hull[,2]<-dat[vertices,2]
    hull[,3]<-dat[vertices,3]
    hull[,4]<-dat[vertices,4]

    cofg<-apply(hull,2,mean)#center of gravity (mean for each trait)
    dGi<-sqrt(apply(((sweep(dat, 2, cofg, `-`))^2), 1,sum)) #euklidian distance of each point to cofg
    dG<-mean(dGi)

    if(weight){
       delta_d<-sum(wi*(dGi-dG))
       delta_d_abs<-sum(wi*abs(dGi-dG))
    }else{ #no weighting
       delta_d<-sum(dGi-dG)/length(data$SLA)
       delta_d_abs<-sum(abs(dGi-dG))/length(data$SLA)
    }
    divergence<-(delta_d+dG)/(delta_d_abs+dG)
    #divergence

    ### Functional Evenness ###
    S<-d<-min.span.tree<-EW1<-evenness<-0
    S<-length(dat[,1]) # Number of individuals
    d <- dist(dat) #distance matrix
    min.span.tree<-spantree(d, toolong = 0)#calculate minimim spanning tree
 if(weight){
      kid<-min.span.tree$kid
      dist<-min.span.tree$dist
      wi_i<-wi[2:length(wi)]
      wi_j<-wi[kid]
      EWl<-min.span.tree$dist/(wi_i+wi_j)
    }else{ #no weighting
      EWl<-min.span.tree$dist
    }

     PEWl<-EWl/sum(EWl) #see: NEW MULTIDIMENSIONAL FUNCTIONAL DIVERSITY INDICES FOR A MULTIFACETED FRAMEWORK IN FUNCTIONAL ECOLOGY Sébastien Villéger
     PEWl.min<-which(PEWl>(1/(S-1)))
     PEWl[PEWl.min]<-1/(S-1)# if PEWl is bigger than 1/(S-1) replace it with 1/(S-1)

     evenness<-(sum(PEWl)-1/(S-1))/(1-1/(S-1)) #formula for evenness
   }else{
     richness<-divergence<-evenness<--1
   }

  return(c(richness,divergence,evenness))
}

data.array<-array(NA, c(length(cells),3,length(years),5)) #cell,divInd,year,PFT
quant.array<-array(NA, c(length(cells),length(years)))

options(warn=2)
for(iyear in (length(years)-30):length(years)){
#for(iyear in 1:2){
  print(years[iyear])
  nskip <- if(iyear>1){sum(numlines[1:(iyear-1)])+iyear}else{1}
  #print(numlines[iyear])
  #print(paste(path.in,"/",filename,".txt",sep=""))
  data<-fread(paste(path.in,"/",filename,".txt",sep=""),blank.lines.skip=T, skip=nskip,nrows =numlines[iyear],select=2:8)
  names(data)<-c("Cell","SLA","Wooddens","LL","PFT","Biomass","Height")
  #print(head(data))
  data$SLA<-data$SLA*455
  data$Wooddens<-data$Wooddens*(10^-6)/0.455
  #head(data)
  ### remove unphysical values
  data[ which(data[,2:4]==0,arr.ind = T)[,1],]<-NA #trait is zero
  data[ which(data$Biomass<10,arr.ind = T),]<-NA #Biomass smaller than 10g
  #print(head(data))
  data[ which(data$LL>20),]<-NA #longevity bigger than 20a
  #print(head(data))
  data<-data[complete.cases(data),]
  #print(head(data))
  #cells<-unique(data$Cell)
  #cells<-sort(cells)
  #print(cells)
    for(icell in 1:length(cells)){
     dat<-subset(data,Cell==icell-1)

     if(quant_bool){
       quant_height<-quantile(dat$Height,probs = quant)
       dat<-subset(dat,Height>quant_height)
       quant.array[icell,iyear]<-quant_height
       }

     Ind<-calc.divInd(data=dat,weight=weight_bool)
     #print(paste(icell,Ind))
     data.array[icell,,iyear,1]<-Ind

     for(ipft in 1:4){
     datpft<-subset(dat,PFT==ipft-1)
     Ind<-calc.divInd(data=datpft,weight=weight_bool)
     #print(paste(icell,Ind))
     data.array[icell,,iyear,ipft+1]<-Ind
     }
    }
}

if(quant_bool){

  if(weight_bool){
    save(data.array,file=paste(path.out,"/","DivInd_",filename,'_weight_4d_quant.Rdata', sep=''))
  }else{
    save(data.array,file=paste(path.out,"/","DivInd_",filename,'_4d_quant.Rdata', sep=''))
  }
  save(quant.array,file=paste(path.out,"/","quantile_height_",filename,'.Rdata', sep=''))

  }else{

  if(weight_bool){
    save(data.array,file=paste(path.out,"/","DivInd_",filename,'_weight_4d.Rdata', sep=''))
  }else{
  save(data.array,file=paste(path.out,"/","DivInd_",filename,'_4d.Rdata', sep=''))
  }
  }