
########################################################################
#
#	Simple prediction per year PA data - Sampling Sites
#			Climate scenario
#
#######################################################################

library(Hmsc)
library(tidyr)
library(dplyr)
library(abind)
library(scales)
library(vegan)
library(betapart)

## Model results
#..................................................................................
setwd("E:/Helsinki/ResultsAug23/For submission/1.Data/3.ModelRun")
load('test10908cent_yearPA.RDATA')
modelyear = list()
modelyear[[1]] = models[[1]]
load('test10908cent_yearAB.RDATA')
modelyear[[2]] = models[[2]]


## Trait info
#..................................................................................
setwd("E:/Helsinki/ResultsAug23/For submission/1.Data/2.ObsInfo")
Traits = read.csv('Traits_Moths.csv', header=T)
head(Traits)

Traits.sub= Traits %>%
  dplyr::select(Genus.and.species, female.size..average.wingspan..mm.,
                HostPlantForm)

Traits.sub$HostPlantForm = as.factor(Traits.sub$HostPlantForm)
levels(Traits.sub$HostPlantForm) <- c('herbaceous', 'graminoids', 'lichens-algae', 'Other', 'tree-stemmed', 
                                      'decideoustrees-shrub', 'Conifer-shrubs')

head(Traits.sub)
colnames(Traits.sub)[1] = 'species'
colnames(Traits.sub)[2] = 'WingspanF'

Traits.sub$WingGP = NA
Traits.sub$WingGP[Traits.sub$WingspanF <28] = 'W<28'
Traits.sub$WingGP[Traits.sub$WingspanF > 28 & Traits.sub$WingspanF < 34] = '28<W<34'
Traits.sub$WingGP[Traits.sub$WingspanF > 34] = 'W>34'


# cluster info
#..................................................................................
setwd("E:/Helsinki/ResultsAug23/For submission/1.Data/1.EnvironmentalInfo")
load('MothclusterFIN.RDATA')
dim(Moth_knn)  # 109    13
head(Moth_knn)
Moth_knn2 = Moth_knn[, 7:8]
colnames(Moth_knn2)[2] = 'sample'


STDES = modelyear[[1]]$studyDesign
STDES_clus = merge(STDES, Moth_knn2, by='sample')
head(STDES_clus)

STDES_clus = STDES_clus[order(STDES_clus$sample),]
dim(STDES_clus)

modelyear[[1]]$XData
Ypa = modelyear[[1]]$Y
dim(Ypa)


## get enviro info all sites all years
#..................................................................................
setwd("E:/Helsinki/ResultsAug23/For submission/1.Data/1.EnvironmentalInfo")
load('Site_env2508.RDATA')
head(Siteenv)

Siteenv = Siteenv[order(Siteenv$Year),]

year = unique(Siteenv$Year)

i=1


Site.unique = unique(modelyear[[1]]$studyDesign[, c(1,3)])
StudyDesignA = data.frame(sample = as.factor(rep(Site.unique[,1], 23)),
                          year = as.factor(rep(year, each = dim(Site.unique)[1])),
                          region = as.factor(rep(Site.unique[,2], 23)))
dim(StudyDesignA)

#StudyDesignALL[order(StudyDesignALL$year, StudyDesignALL$sample),]


StudyDesignN = StudyDesignA[StudyDesignA$year == year[i],]
dim(StudyDesignN)

sample_all = c(modelyear[[1]]$studyDesign$sample, as.factor(StudyDesignN$sample))
rLsample = HmscRandomLevel(units = as.factor(sample_all))
region_all = c(modelyear[[1]]$studyDesign$region, as.factor(StudyDesignN$region))    
rLregion = HmscRandomLevel(units = region_all)
year_all = c(modelyear[[1]]$studyDesign$year, as.factor(StudyDesignN$year))
rLyear= HmscRandomLevel(units = as.factor(year_all))

XData = Siteenv[, -c(1,2,4, 10)]
XData$Effort = mean(modelyear[[1]]$XData$Effort)
dim(XData)
colnames(XData) = colnames(modelyear[[1]]$XData)[-c(2,8)]

seqin = seq(1, by=109, length.out = 23)
seqend = seq(109, by=109, length.out = 23)

XDatasub = (XData[seqin[i]:seqend[i] ,])
head(XDatasub)
dim(XDatasub)

# Hab does not change
HabY1 = XData[seqin[1]:seqend[1], 5:10]
XDatasub[, 5:10] = HabY1


## PA predictions
#..................................................................................
start_time <- Sys.time()
predYnew = predict(modelyear[[1]], studyDesign=StudyDesignN, XData=XDatasub, 
                   ranLevels=list('sample'=rLsample, 'year'=rLyear, 'region'= rLregion),
                   expected=TRUE, predictEtaMean=TRUE)
end_time <- Sys.time()
end_time - start_time

length(predYnew)
predYnew2 = abind(predYnew, along=3)
Raw1_Y1 = apply(predYnew2, c(1,2), mean) # mean over MCMC
dim(Raw1_Y1)



S_Y1 = rowSums(Raw1_Y1)

# group species
Raw1H1_Y1 = Raw1_Y1[, which(colnames(Raw1_Y1) %in% Traits.sub[Traits.sub$HostPlantForm == 'herbaceous',]$species)]
Raw1H2_Y1 = Raw1_Y1[, which(colnames(Raw1_Y1) %in% Traits.sub[Traits.sub$HostPlantForm == 'graminoids',]$species)]
Raw1H4_Y1 = Raw1_Y1[, which(colnames(Raw1_Y1) %in% Traits.sub[Traits.sub$HostPlantForm %in% c('lichens-algae', 'Other'),]$species)]
Raw1H31_Y1 = Raw1_Y1[, which(colnames(Raw1_Y1) %in% Traits.sub[Traits.sub$HostPlantForm == 'tree-stemmed',]$species)]
Raw1H32_Y1 = Raw1_Y1[, which(colnames(Raw1_Y1) %in% Traits.sub[Traits.sub$HostPlantForm == 'decideoustrees-shrub',]$species)]
Raw1H33_Y1 = Raw1_Y1[, which(colnames(Raw1_Y1) %in% Traits.sub[Traits.sub$HostPlantForm == 'Conifer-shrubs',]$species)]
dim(Raw1H1_Y1)
dim(Raw1H2_Y1)

Raw1H1_Y1 = rowSums(Raw1H1_Y1)#/dim(Raw1H1_Y1)[2]
Raw1H2_Y1 = rowSums(Raw1H2_Y1)#/dim(Raw1H2_Y1)[2]    
Raw1H4_Y1 = rowSums(Raw1H4_Y1)#/dim(Raw1H4_Y1)[2]    
Raw1H31_Y1 = rowSums(Raw1H31_Y1)#/dim(Raw1H31_Y1)[2]    
Raw1H32_Y1 = rowSums(Raw1H32_Y1)#/dim(Raw1H32_Y1)[2]    
Raw1H33_Y1 = rowSums(Raw1H33_Y1)#/dim(Raw1H33_Y1)[2]    
length(Raw1H1_Y1)



Raw1HW1_Y1 = Raw1_Y1[, which(colnames(Raw1_Y1) %in% Traits.sub[Traits.sub$WingGP == 'W<28',]$species)]
Raw1HW2_Y1 = Raw1_Y1[, which(colnames(Raw1_Y1) %in% Traits.sub[Traits.sub$WingGP == '28<W<34',]$species)]
Raw1HW3_Y1 = Raw1_Y1[, which(colnames(Raw1_Y1) %in% Traits.sub[Traits.sub$WingGP %in% c('W>34'),]$species)]


Raw1HW1_Y1 = rowSums(Raw1HW1_Y1)#/dim(Raw1HH1_Y1)[2]
Raw1HW2_Y1 = rowSums(Raw1HW2_Y1)#/dim(Raw1HH2_Y1)[2]    
Raw1HW3_Y1 = rowSums(Raw1HW3_Y1)#/dim(Raw1HH4_Y1)[2]    



setwd("/home/local/guilbaul/tsclient/E/Helsinki/ResultsAug23/RDATA/4.Predictions")

save(S_Y1, Raw1H1_Y1, Raw1H2_Y1, Raw1H4_Y1, Raw1H31_Y1, Raw1H32_Y1,
     Raw1H33_Y1, Raw1HW1_Y1, Raw1HW2_Y1, Raw1HW3_Y1, file='predPArichSITclim_Y1.RDATA')


## AB predictions
#..................................................................................

start_time <- Sys.time()

predYnew = predict(modelyear[[1]], studyDesign=StudyDesignN, XData=XDatasub, 
                   ranLevels=list('sample'=rLsample, 'year'=rLyear, 'region'= rLregion),
                   expected=FALSE, predictEtaMean=TRUE)
predYnew1000 = predYnew
rm(predYnew)

predYnewAB = predict(modelyear[[2]], studyDesign=StudyDesignN, XData=XDatasub, 
                     ranLevels=list('sample'=rLsample, 'year'=rLyear, 'region'= rLregion),
                     expected=TRUE, predictEtaMean=TRUE)

predYnewAB1000 = predYnewAB

end_time <- Sys.time()
end_time - start_time # 25 sec

rm(predYnewAB)

Arr.predPA1000 = array((unlist(predYnew1000)), dim = c(dim(predYnew1000[[1]]), length(predYnew1000)))
Arr.predAB1000 = array((unlist(predYnewAB1000)), dim = c(dim(predYnewAB1000[[1]]), length(predYnewAB1000)))
Arr.pred = Arr.predPA1000*exp(Arr.predAB1000)

predYnew2AB = abind(Arr.pred, along=3)  # 
#save(predYnew2AB, file = 'SitepredABY1.RDATA')


setwd("E:/Helsinki/ResultsAug23/For submission/1.Data/4.Predictions")

load('colnamesp.RDATA')

ScurrentW1 = eve.currentW1 = list()
ScurrentW2 = eve.currentW2 = list()
ScurrentW3 = eve.currentW3 = list()

ScurrentH1 = eve.currentH1 = list()
ScurrentH2 = eve.currentH2 = list()
ScurrentH4 = eve.currentH4 = list()
ScurrentH31 = eve.currentH31 = list()
ScurrentH32 = eve.currentH32 = list()
ScurrentH33 = eve.currentH33 = list()
ScurrentAll = eve.currentAll = list()
for(j in 1:1000){
  
  predY.j = predYnew2AB[,,j]
  colnames(predY.j) = col.namesp
  
  # separate per trait
  predYRW1_Y1 = predY.j[, which(colnames(predY.j) %in% Traits.sub[Traits.sub$WingGP == 'W<28',]$species)]
  predYRW2_Y1 = predY.j[, which(colnames(predY.j) %in% Traits.sub[Traits.sub$WingGP == '28<W<34',]$species)]
  predYRW3_Y1 = predY.j[, which(colnames(predY.j) %in% Traits.sub[Traits.sub$WingGP %in% c('W>34'),]$species)]
  
  ScurrentW1[[j]] <- apply(predYRW1_Y1>0, 1, sum)
  eve.currentW1[[j]] = (diversity(predYRW1_Y1, index="invsimpson"))/(ScurrentW1[[j]])
  
  ScurrentW2[[j]] <- apply(predYRW2_Y1>0, 1, sum)
  eve.currentW2[[j]] = (diversity(predYRW2_Y1, index="invsimpson"))/(ScurrentW2[[j]])
  
  ScurrentW3[[j]] <- apply(predYRW3_Y1>0, 1, sum)
  eve.currentW3[[j]] = (diversity(predYRW3_Y1, index="invsimpson"))/(ScurrentW3[[j]])
  
  # separate per trait
  predYRH1_Y1 = predY.j[, which(colnames(predY.j) %in% Traits.sub[Traits.sub$HostPlantForm == 'herbaceous',]$species)]
  predYRH2_Y1 = predY.j[, which(colnames(predY.j) %in% Traits.sub[Traits.sub$HostPlantForm == 'graminoids',]$species)]
  predYRH4_Y1 = predY.j[, which(colnames(predY.j) %in% Traits.sub[Traits.sub$HostPlantForm %in% c('lichens-algae', 'Other'),]$species)]
  predYRH31_Y1 = predY.j[, which(colnames(predY.j) %in% Traits.sub[Traits.sub$HostPlantForm == 'tree-stemmed',]$species)]
  predYRH32_Y1 = predY.j[, which(colnames(predY.j) %in% Traits.sub[Traits.sub$HostPlantForm == 'decideoustrees-shrub',]$species)]
  predYRH33_Y1 = predY.j[, which(colnames(predY.j) %in% Traits.sub[Traits.sub$HostPlantForm == 'Conifer-shrubs',]$species)]
  dim(predYRH1_Y1)
  
  ScurrentH1[[j]] <- apply(predYRH1_Y1>0, 1, sum)
  eve.currentH1[[j]] = (diversity(predYRH1_Y1, index="invsimpson"))/(ScurrentH1[[j]])
  
  ScurrentH2[[j]] <- apply(predYRH2_Y1>0, 1, sum)
  eve.currentH2[[j]] = (diversity(predYRH2_Y1, index="invsimpson"))/(ScurrentH2[[j]])
  
  ScurrentH4[[j]] <- apply(predYRH4_Y1>0, 1, sum)
  eve.currentH4[[j]] = (diversity(predYRH4_Y1, index="invsimpson"))/(ScurrentH4[[j]])
  
  ScurrentH31[[j]] <- apply(predYRH31_Y1>0, 1, sum)
  eve.currentH31[[j]] = (diversity(predYRH31_Y1, index="invsimpson"))/(ScurrentH31[[j]])
  
  ScurrentH32[[j]] <- apply(predYRH32_Y1>0, 1, sum)
  eve.currentH32[[j]] = (diversity(predYRH32_Y1, index="invsimpson"))/(ScurrentH32[[j]])
  
  ScurrentH33[[j]] <- apply(predYRH33_Y1>0, 1, sum)
  eve.currentH33[[j]] = (diversity(predYRH33_Y1, index="invsimpson"))/(ScurrentH33[[j]])
  
  predYRAll_Y1 = predY.j
  ScurrentAll[[j]] <- apply(predYRAll_Y1>0, 1, sum)
  eve.currentAll[[j]] = (diversity(predYRAll_Y1, index="invsimpson"))/(ScurrentAll[[j]])
  
}

eveH1_Y1 = rowMeans(do.call(cbind, eve.currentH1)*is.finite(do.call(cbind, eve.currentH1)), 
                    na.rm = TRUE)#Reduce("+", eve.currentH1)/length(eve.currentH1)
eveH2_Y1 = rowMeans(do.call(cbind, eve.currentH2)*is.finite(do.call(cbind, eve.currentH2)), 
                    na.rm = TRUE)#Reduce("+", eve.currentH2)/length(eve.currentH2)
eveH4_Y1 = rowMeans(do.call(cbind, eve.currentH4)*is.finite(do.call(cbind, eve.currentH4)), 
                    na.rm = TRUE)#Reduce("+", eve.currentH4)/length(eve.currentH4)
eveH31_Y1 = rowMeans(do.call(cbind, eve.currentH31)*is.finite(do.call(cbind, eve.currentH31)), 
                     na.rm = TRUE)#Reduce("+", eve.currentH31)/length(eve.currentH31)
eveH32_Y1 = rowMeans(do.call(cbind, eve.currentH32)*is.finite(do.call(cbind, eve.currentH32)), 
                     na.rm = TRUE)#Reduce("+", eve.currentH32)/length(eve.currentH32)
eveH33_Y1 = rowMeans(do.call(cbind, eve.currentH33)*is.finite(do.call(cbind, eve.currentH33)), 
                     na.rm = TRUE)#Reduce("+", eve.currentH33)/length(eve.currentH33)
eveAll_Y1 = rowMeans(do.call(cbind, eve.currentAll)*is.finite(do.call(cbind, eve.currentAll)), na.rm = TRUE)


eveW1_Y1 = rowMeans(do.call(cbind, eve.currentW1)*is.finite(do.call(cbind, eve.currentW1)), 
                    na.rm = TRUE)#Reduce("+", eve.currentH1)/length(eve.currentH1)
eveW2_Y1 = rowMeans(do.call(cbind, eve.currentW2)*is.finite(do.call(cbind, eve.currentW2)), 
                    na.rm = TRUE)#Reduce("+", eve.currentH2)/length(eve.currentH2)
eveW3_Y1 = rowMeans(do.call(cbind, eve.currentW3)*is.finite(do.call(cbind, eve.currentW3)), 
                    na.rm = TRUE)#Reduce("+", eve.currentH4)/length(eve.currentH4)


setwd("E:/Helsinki/ResultsAug23/For submission/1.Data/4.Predictions")

save(eveW3_Y1, eveW2_Y1, eveW1_Y1, eveAll_Y1, eveH33_Y1, eveH32_Y1, eveH31_Y1,
     eveH4_Y1, eveH2_Y1, eveH1_Y1, file = 'predABeveSITclim_Y1.RDATA')
