library(Hmsc)
library(ggplot2)
library(abind)
library(dplyr)


i=1

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


########################################################################
#
#	Simple prediction per year AB data - Finnish grid Sites
#			Habitat scenario
#
#######################################################################

# Get model information Finnish grid Sites - HMSC run

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

# 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)

dim(mapDataALL)
mapcoord1 = mapDataALL[mapDataALL$Year == 1998, 14:15]
mapDataALL$SiteID = paste(mapDataALL$X, '', mapDataALL$Y)
SiteFIn_un = unique(mapDataALL$SiteID)
MapdataALL.list = split(mapDataALL, f = mapDataALL$Year) 
Fin_coord = MapdataALL.list[[1]][, 14:15]


# Prepare data
#........................................................................................
## 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]


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

## Hab change scenario : Clim Stay the same
ClimY1.list= list(mapDataALL[Sini[1]:Send[1], c(2:6)], mapDataALL[Sini[2]:Send[2], c(2:6)],
                  mapDataALL[Sini[3]:Send[3], c(2:6)], mapDataALL[Sini[4]:Send[4], c(2:6)],
                  mapDataALL[Sini[5]:Send[5], c(2:6)])

ClimY1 = Reduce('+', ClimY1.list)/length(ClimY1.list)
summary(ClimY1)

year.un = unique(mapDataALL$Year)

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

# cli  does not change
Xdata_FIN_small[, 1:5] = ClimY1

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()

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


predYnewAB1000 = predict(modelyear[[2]], 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 # 10.5 min

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.pred1000 = Arr.predPA1000*exp(Arr.predAB1000)


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

#predYnew2current = predYnew2
predYnew2AB = abind(Arr.pred1000, along=3)  # take too much space  # dim site 14484  sp 80  mcmc 100

rm(Arr.predPA1000, Arr.predAB1000, Arr.pred1000)
#save(predYnew2AB, file = 'predABcurentY1_1508.RDATA')



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


library(vegan)
library(betapart)

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)

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

