library(Hmsc)
library(ggplot2)
library(abind)
library(dplyr)
library(ggpubr)
library(patchwork)

########################################################################
#
#	Simple prediction per year PA data - Finnish grid Sites
#			Climate scenario
#
#######################################################################

# Get model information Finnish grid Sites - HMSC run

setwd("E:/Helsinki/ResultsAug23/For submission/1.Data/4.Predictions")
load('Predsmall2408.RDATA') # 

Sini = seq(1, by=12340, length.out=23)
Send = seq(12340, by=12340, length.out=23)

i = 23

# Grid sites info global
setwd("E:/Helsinki/ResultsAug23/For submission/1.Data/1.EnvironmentalInfo")
load('CurrentSmallFIN_1408.RDATA')
tail(mapData2)
tail(mapData2hab)

mapDataALL = cbind(mapData2[, c(1:2, 6, 15, 10:11)], mapData2hab[, c(2:10)])
head(mapDataALL)


# Prepare data
#........................................................................................
## Clim change scenario : hab Stay the same
HabY1 = mapDataALL[Sini[1]:Send[1], c(7:13)]
summary(HabY1)

year.un = unique(mapDataALL$Year)

Xdata_FIN_small =  mapDataALL[Sini[i]:Send[i], -c(1, 14, 15)]
summary(Xdata_FIN_small)

# Hab does not change
Xdata_FIN_small[, 6:12] = HabY1

Xdata_FIN_small$Effort = mean(modelyear[[1]]$XData$Effort)

colnames(Xdata_FIN_small) = colnames(modelyear[[1]]$XData)

#sample.id
StudyDesignNew = data.frame(sample=as.factor(sample.id[Sini[i]:Send[i]]),
                            year=as.factor(year[Sini[i]:Send[i]]),
                            region=as.factor(region[Sini[1]:Send[1]]))  # because it is always the same
StudyDesignAll = rbind(modelyear[[1]]$studyDesign, StudyDesignNew)

XData_pred = rbind(modelyear[[1]]$XData, Xdata_FIN_small)
sample_all = c(modelyear[[1]]$studyDesign$sample, as.factor(sample.id[Sini[i]:Send[i]]))
rLsample=HmscRandomLevel(units = as.factor(sample_all))
region_all = c(modelyear[[1]]$studyDesign$region, as.factor(region[Sini[1]:Send[1]]))    # because it is always the same
rLregion=HmscRandomLevel(units = region_all)
year_all = c(modelyear[[1]]$studyDesign$year, as.factor(year[Sini[i]:Send[i]]))
rLyear= HmscRandomLevel(units = as.factor(year_all))



#object.size(Xdata_FIN_small)


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

rm(Xdata_FIN_small, mapData2)
rm(mapDataALL, mapDataALL.sd, Xdata_FIN_small, modelyear, StudyDesignAll, StudyDesignNew)

#head(predYnew) 
length(predYnew) ## should be mcmc posterior samples 6 000
dim(predYnew[[1]]) ## site species 12340 77

predYnew2 = abind(predYnew, along=3)  # take too much space


#save(predYnew2, file = 'predPAcurentY230_1408.RDATA')


# Prepare predictions
#........................................................................................
## 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('clusterdata2.RDATA')
dim(mapData3)

mapData3 = mapData3[order(mapData3$X.x) ,] # same order as model
cluster = mapData3[, 20]


predYR_Y23 = apply(predYnew2, c(1,2), mean) # mean over MCMC     file='ABtrend_SIT_1306.RDATA')

dim(predYR_Y23)

#predYR_Y23sd = apply(predYnew2, c(1,2), sd)

# Species richness
#........................................................................................
dim(predYR_Y23)  # ## site species 12340 80
S_Y23 = rowSums(predYR_Y23)


predYRH1_Y23 = predYR_Y23[, which(colnames(predYR_Y23) %in% Traits.sub[Traits.sub$HostPlantForm == 'herbaceous',]$species)]
predYRH2_Y23 = predYR_Y23[, which(colnames(predYR_Y23) %in% Traits.sub[Traits.sub$HostPlantForm == 'graminoids',]$species)]
predYRH4_Y23 = predYR_Y23[, which(colnames(predYR_Y23) %in% Traits.sub[Traits.sub$HostPlantForm %in% c('lichens-algae', 'Other'),]$species)]
predYRH31_Y23 = predYR_Y23[, which(colnames(predYR_Y23) %in% Traits.sub[Traits.sub$HostPlantForm == 'tree-stemmed',]$species)]
predYRH32_Y23 = predYR_Y23[, which(colnames(predYR_Y23) %in% Traits.sub[Traits.sub$HostPlantForm == 'decideoustrees-shrub',]$species)]
predYRH33_Y23 = predYR_Y23[, which(colnames(predYR_Y23) %in% Traits.sub[Traits.sub$HostPlantForm == 'Conifer-shrubs',]$species)]
dim(predYRH1_Y23)
dim(predYRH2_Y23)

predYRH1_Y23 = rowSums(predYRH1_Y23)#/dim(predYRH1_Y23)[2]
predYRH2_Y23 = rowSums(predYRH2_Y23)#/dim(predYRH2_Y23)[2]    
predYRH4_Y23 = rowSums(predYRH4_Y23)#/dim(predYRH4_Y23)[2]    
predYRH31_Y23 = rowSums(predYRH31_Y23)#/dim(predYRH31_Y23)[2]    
predYRH32_Y23 = rowSums(predYRH32_Y23)#/dim(predYRH32_Y23)[2]    
predYRH33_Y23 = rowSums(predYRH33_Y23)#/dim(predYRH33_Y23)[2]    




predYRW1_Y23 = predYR_Y23[, which(colnames(predYR_Y23) %in% Traits.sub[Traits.sub$WingGP == 'W<28',]$species)]
predYRW2_Y23 = predYR_Y23[, which(colnames(predYR_Y23) %in% Traits.sub[Traits.sub$WingGP == '28<W<34',]$species)]
predYRW3_Y23 = predYR_Y23[, which(colnames(predYR_Y23) %in% Traits.sub[Traits.sub$WingGP %in% c('W>34'),]$species)]


predYRW1_Y23 = rowSums(predYRW1_Y23)#/dim(predYRH1_Y23)[2]
predYRW2_Y23 = rowSums(predYRW2_Y23)#/dim(predYRH2_Y23)[2]    
predYRW3_Y23 = rowSums(predYRW3_Y23)#/dim(predYRH4_Y23)[2]    



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

save(S_Y23, predYRH1_Y23, predYRH2_Y23, predYRH4_Y23, predYRH31_Y23, predYRH32_Y23,
     predYRH33_Y23, predYRW1_Y23, predYRW2_Y23, predYRW3_Y23, file='predPArich_Y23clim.RDATA')
