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

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

ncells<-50

path.in<-path.out<-"D:/Data/LPJmL-FIT/validation/"

sites<-c("seitseminen","kalkalpen","bialowieza","hainich","laegeren","dundo")
lons<-c(23.25,14.25,23.75,10.25,8.25,14.75)
lats<-c(61.75,47.75,52.75,51.25,47.25,44.75)

startyear_simulation<-1901

startyear<-2004
endyear<-2013
nyears<-113

istart<-startyear-startyear_simulation+1
iend<-endyear-startyear_simulation+1



path.in.clm<-"D:/Data/LPJmL-Fit/Input/"
load(file=paste(path.in.clm,"climate_EU.Rdata",sep=""))
### load grid ###
headersize<-51
file  <- file(paste(path.in.clm,"grid_europe.clm",sep=""),"rb") # read binary mode
seek(file, where = headersize, origin="start")
x<-readBin(file,integer(),n = 2*5095, size=2)/100
lon.clm <- x[c(1:5095)*2-1]
lat.clm <- x[c(1:5095)*2]
close(file)


calc.divInd<-function(data,weight=F){
  if(nrow(data)>10){
    #extract data
    dat<-cbind(data$SLA,data$Wooddens,data$LL)
    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
    
    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])
    
    ### 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),3))
    
    hull[,1]<-dat[vertices,1]
    hull[,2]<-dat[vertices,2]
    hull[,3]<-dat[vertices,3]    
    
    cofg<-apply(hull,2,mean)#center of gravity (mean for each trait)
    dGi<-sqrt(apply(((dat-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<-NA
  }
  
  return(c(richness,divergence,evenness))
}

calc.divInd.4d<-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
    H<-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]-H[1])/(H[2]-H[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<-NA
  }
  
  return(c(richness,divergence,evenness))
}

filenames<-list.files(path.in) #get all .nc filenames in path.output.FIT

data_table<-array(NA,c(length(sites),10))
data_divInd_pft_even<-data_divInd_pft_weight<-array(NA,c(length(sites),3,4))
for(isite in 1:length(sites)){
  
  print(sites[isite])
  #### MCWD & MAT ####
  
  cell.clm<-intersect(which(lon.clm==lons[isite]),which(lat.clm==lats[isite]))
  
  MCWD_site<-mean(MCWD[cell.clm,istart:iend])
  MAT_site<-mean(temp.annual.mean[cell.clm,istart:iend])
  
  
  #### GPP ####
  
  #filter filename that contain outputname
  filename<-grep('mgpp', filenames, value = TRUE) 
  filename<-grep(sites[isite], filename, value = TRUE) 
  
  file  <- file(paste(path.in,filename,sep=""),"rb") # read binary mode
  data<-readBin(file,double(),n = ncells*12*113, size=4)
  close(file)
  str(data)
  
  dim(data)<-c(ncells,12,113)
  
  data<-apply(data[,,istart:iend],c(2),mean)
  GPP_sum_site<-sum(data)
  
  #### VEGC ####
  
  #filter filename that contain outputname
  filename<-grep('vegc', filenames, value = TRUE) 
  filename<-grep(sites[isite], filename, value = TRUE) 
  
  file  <- file(paste(path.in,filename,sep=""),"rb") # read binary mode
  data<-readBin(file,double(),n = ncells*1*113, size=4)
  close(file)
  str(data)
  
  dim(data)<-c(ncells,113)
  
  VEGC_site<-mean(data[,istart:iend])/1000
  
  
  #### HEIGHT ####
  
  #filter filename that contain outputname
  filename<-grep('height', filenames, value = TRUE) 
  filename<-grep(sites[isite], filename, value = TRUE) 
  
  file  <- file(paste(path.in,filename,sep=""),"rb") # read binary mode
  data<-readBin(file,double(),n = 50*100*nyears*4, size=4)
  close(file)
  
  dim(data)<-c(50,100,4,nyears)
  bins<-seq(0, 50,length.out = 100)
  
  str(data)
  data<-apply(data[,,,istart:iend],c(2,3),mean)
  data<-apply(data,1,sum)
  
  HEIGHT_site<-sum(data*bins/sum(data))
  
  
  #### Diversity Indices ####
  
  filename<-grep('trait_ind', filenames, value = TRUE) 
  filename<-grep(sites[isite], filename, value = TRUE) 
  data<-fread(paste(path.in,"/",filename,sep=""),blank.lines.skip=T)
  names(data)<-c("Year","Cell","SLA","Wooddens","LL","PFT","Biomass","Height")
  
  data$SLA<-data$SLA*455
  data$Wooddens<-data$Wooddens*(10^-6)/0.455
  
  ### remove unphysical values
  data[ which(data[,3:5]==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),]
  
  years<-2004:2013
  cells<-(1:50)-1
  
  weight_bool<-T
  
  data.array<-array(NA,c(10,50,3))
  data.array_weight<-data.array_even<-array(NA,c(10,50,3,4))
  
  for(iyear in 1:10){
    dat_year<-subset(data,Year==years[iyear])
    for(icell in 1:50){ 
      dat<-subset(dat_year,Cell==cells[icell])
      Ind<-calc.divInd.4d(data=dat,weight=weight_bool)
      #print(paste(icell,Ind))
      data.array[iyear,icell,]<-Ind
    for(ipft in 1:4){
        datpft<-subset(dat,PFT==ipft-1)
        Ind<-calc.divInd.4d(data=datpft,weight=weight_bool)
        #print(paste(icell,Ind))
        data.array_weight[iyear,icell,,ipft]<-Ind
      }
    }
  }
  
  Rich_weight_site<-mean(data.array[,,1])
  Div_weight_site<-mean(data.array[,,2])
  Even_weight_site<-mean(data.array[,,3])
  
  for(ipft in 1:4){
    data_divInd_pft_weight[isite,,ipft]<-apply(data.array_weight,c(3,4),mean,na.rm=T)[,ipft]
  }
  
  
  weight_bool<-F
  
  data.array<-array(NA,c(10,50,3))
  
  for(iyear in 1:10){
    dat_year<-subset(data,Year==years[iyear])
    for(icell in 1:50){ 
      dat<-subset(dat_year,Cell==cells[icell])
      Ind<-calc.divInd.4d(data=dat,weight=weight_bool)
      #print(paste(icell,Ind))
      data.array[iyear,icell,]<-Ind
      
      for(ipft in 1:4){
        datpft<-subset(dat,PFT==ipft-1)
        Ind<-calc.divInd.4d(data=datpft,weight=weight_bool)
        #print(paste(icell,Ind))
        data.array_even[iyear,icell,,ipft]<-Ind
      }
    }
  }
  
  for(ipft in 1:4){
    data_divInd_pft_even[isite,,ipft]<-apply(data.array_even,c(3,4),mean,na.rm=T)[,ipft]
  }
  
  Rich_even_site<-mean(data.array[,,1])
  Div_even_site<-mean(data.array[,,2])
  Even_even_site<-mean(data.array[,,3])
  
  data_table[isite,1]<-MCWD_site
  data_table[isite,2]<-MAT_site
  data_table[isite,3]<-GPP_sum_site
  data_table[isite,4]<-VEGC_site
  data_table[isite,5]<-HEIGHT_site
  data_table[isite,6]<-Rich_even_site
  data_table[isite,7]<-Div_even_site
  data_table[isite,8]<-Even_even_site
  data_table[isite,9]<-Div_weight_site
  data_table[isite,10]<-Even_weight_site
    
  print(apply(data.array_weight,c(3,4),mean,na.rm=T))
  print(apply(data.array_even,c(3,4),mean,na.rm=T))
}

data_table

data_divInd_pft_even
data_divInd_pft_weight