
library(ggplot2)
library(ggridges)
library(grid)
library(data.table)
library(reshape2)

library(Hmsc)
library(reshape2)
library(abind)
library(tidyr)
library(dplyr)
#library(berryFunctions)
#library(lemon)
#library(cowplot)
library(patchwork)
library(ggh4x)
library(ggpubr)

setwd("D:/Helsinki/ResultsAug23/RDATA/3.ModelRun")
load('VPcondAFIN100cdd_PA_1408.RDATA')
VPa_PA = VP.a
VP.cond_PA = VP.cond

load('VPcondAFIN100cdd_AB_1408.RDATA')
VPa_AB = VP.a
VP.cond_AB = VP.cond


# yearly model
#####  species
# VP all ----------------------------------------------------------------------------------------------------------------

VP.a = list(VPa_PA, VPa_AB)

MeanA.vp = list()
conf.up = conf.low = list()
for (i in 1:2) {
 
  # mean over species and MCMC samples
  dim(VP.a[[i]]$vals)  # covar / species /  / MCMC sample
  MeanA.vp[[i]] = apply(VP.a[[i]]$vals, c(1, 2), mean)
  dim(MeanA.vp[[i]])  # covar / btw-with
  
  conf.up[[i]] = apply(VP.a[[i]]$vals, c(1, 2), function(x){quantile(x, 0.975)})
  conf.low[[i]] = apply(VP.a[[i]]$vals, c(1, 2), function(x){quantile(x, 0.025)})
}


lapply(conf.up, rowMeans)
lapply(conf.low, rowMeans)
Mean.Fin = lapply(MeanA.vp, rowMeans)
Mean.Fin

modelyear = list()
load('test10908cent_yearPA.RDATA')
modelyear[[1]] = models[[1]]


# Try to change order of species per trait
modelyear[[1]]$TrData[order(modelyear[[1]]$TrData$HostPlantForm),]

spT_order = modelyear[[1]]$TrData[order(modelyear[[1]]$TrData$HostPlantForm, 
                            modelyear[[1]]$TrData$WingspanF) ,]

lab1 = list('A - Presence Absence', 'B - Abundance')

# VP all ----------------------------------------------------------------------------------------------------------------
group1 = c(1,1,1, 1,1, 1,1, 1,1, 2,2, 2,2,2,2,2,2, 2,2, 3)
groupnames1 = c('Climate','Habitat', 'Effort')

multistack.bar <- function(x=list(x), betweenspace = 2, withinspace=0, group.names = NULL,...){
  # note that number of rows in each component matrix of the list
  # should be equal
  # if you have missing value use "NA", to make it complete 
  # number of column can differ
  mylist <- x
  space = list()
  space[[1]] <- c(rep(withinspace, ncol(mylist[[1]])), betweenspace) 
  
  for ( i in 2:length(mylist)){
    if(i == length(mylist)){
      space[[i]] <- c(rep(withinspace, ncol(mylist[[i]])-1)) 
    } else{
      if(is.null(ncol(mylist[[i]]))){
        space[[i]] <- c(betweenspace) 
      }else{
        space[[i]] <- c(rep(withinspace, ncol(mylist[[i]])-1), betweenspace) 
      }
    }
    
  }
  un.space <- c(unlist(space))
  newdata <- do.call("cbind", mylist)
  b.pos= barplot(newdata, space= un.space, ...)
  
  if(is.null(group.names)){}else{
    tg = list()
    
    for (i in 1:length(unique(group.names))) {
      t.i = group.names[group.names == unique(group.names)[i]]
      if((length(t.i) %% 2) == 0){
        tg[[i]] = c(rep('', ceiling(length(t.i)-1)/2), rep(as.character(unique(group.names)[i]), 1), rep('', ceiling(length(t.i))/2))
      }else{
        tg[[i]] = c(rep('', floor(length(t.i)-1)/2), rep(as.character(unique(group.names)[i]), 1), rep('', ceiling(length(t.i)-1)/2))
        
      }
    }
    
    lab = unlist(tg)
    
    par(mgp=c(2,1,0))
    #axis(4, at=b.pos, labels = lab, tick = F, cex.axis=0.6, las=2)
    mapply(axis, side = 4, at=b.pos, labels = lab, tick = F, cex.axis=1.1, las=2)
    par(mgp=c(3,1,0))
  }
  
}


coul = c("#66C2A5", 'purple',
                     "#E5C494", 'grey78','blue3')
                     

head(spT_order)

levels(spT_order$HostPlantForm)[levels(spT_order$HostPlantForm) == '5'] = '4'


levels(spT_order$HostPlantForm)[levels(spT_order$HostPlantForm) == '1'] = 'Herbaceous'
levels(spT_order$HostPlantForm)[levels(spT_order$HostPlantForm) == '2'] = 'Graminoids'
levels(spT_order$HostPlantForm)[levels(spT_order$HostPlantForm) == '4'] = 'Others'
levels(spT_order$HostPlantForm)[levels(spT_order$HostPlantForm) == '31'] = 'Dwarf shrubs'
levels(spT_order$HostPlantForm)[levels(spT_order$HostPlantForm) == '32'] = 'Deciduous'
levels(spT_order$HostPlantForm)[levels(spT_order$HostPlantForm) == '33'] = 'Coniferous'




#####
####   ----------------------------------------------------------------------------
## with traits info
setwd("D:/Helsinki/Moth_project/ResultsAug23/RDATA/2.ObsInfo")
load('prep_centCLC_A23.RDATA')

head(spT)
colnames(spT)[1] = 'species'
head(spT_order)
levels(spT_order$HostPlantForm)[5] = 'Deciduous'
spT_order$species = rownames(spT_order)

spt2 = merge(spT_order, spT[,-3], by='species', all.x=T)
head(spt2)
spt2 = spt2[order(spt2$HostPlantForm, spt2$WingspanF.x),]
head(spt2)


par(mfrow=c(1,3), mar=c(1,3,2,1), oma=c(1,0,1,0))
layout(matrix(c(1, 2, 3), ncol=3), widths=c(2,1,2))

par(mar = c(1,1,6,2))
##..................................................................................PA
# Finland
VP2_order = MeanA.vp[[1]][, c(spT_order$species)]

ng = dim(VP2_order)[1]
leg = rownames(VP2_order)

mylist2 <- list(VP2_order[,which(spT_order$HostPlantForm == 'Herbaceous')], VP2_order[,which(spT_order$HostPlantForm == 'Graminoids')], 
                VP2_order[,which(spT_order$HostPlantForm == 'Others')], #VP2_order[,which(spT_order$HostPlantForm == '5')],
                VP2_order[,which(spT_order$HostPlantForm == 'Dwarf shrubs')],
                VP2_order[,which(spT_order$HostPlantForm == 'Deciduous')], VP2_order[,which(spT_order$HostPlantForm == 'Coniferous')])
mylist2neg = lapply(mylist2, function(x){(-1)*x})
#mylist2neg = lapply(mylist2, function(y) apply(y, 2, function(x) (-1)*x / sum(x, na.rm = TRUE)))


mn = Mean.Fin[[1]]
means2 = round(100 *mn, 1)

leg2 = c()
for (j in 1:ng) {
  leg2[j] = paste(leg[j], " (mean = ", toString(means2[j]),
                  ")", sep = "")
}
multistack.bar(mylist2neg, betweenspace=1.2, withinspace=0.1, las=1, group.names = '',
               col = coul, xlab="", ylab="", xaxt='n', legend = leg2,  horiz=T, cex.lab =1.7, cex.axis=1.7, yaxt='n',
               main = '', cex.main = 1.7, cex.names = 2,
               args.legend=list(x ='topleft', bty='n', cex=1.4, inset=c(0,-0.08), cex=1.6))
axis(1, at = c(-1, -0.8, -0.6, -0.4, -0.2, 0), labels = c(1, 0.8, 0.6, 0.4, 0.2, 0), pos=0, cex.axis=1.7)
title('a) VP Presence Absence', adj = 0.95, line = 0, cex.main = 1.7)

#### middle plot
par(mar = c(1,1,6,7))

spt2$row = 1:78

# barplot(spt2$WingspanF~spt2$row, at= c(1:22, 24:30, 32:38, 40:55, 57:78, 80:83),
#         horiz=TRUE)

wing.list = list(spt2[which(spt2$HostPlantForm == 'Herbaceous'),]$WingspanF.x, 
                 spt2[which(spt2$HostPlantForm == 'Graminoids'),]$WingspanF.x, 
                 spt2[which(spt2$HostPlantForm == 'Others'),]$WingspanF.x, 
                 spt2[which(spt2$HostPlantForm == 'Dwarf shrubs'),]$WingspanF.x,
                 spt2[which(spt2$HostPlantForm == 'Deciduous'),]$WingspanF.x, 
                 spt2[which(spt2$HostPlantForm == 'Coniferous'),]$WingspanF.x)

load('space_barplot.RDATA')
space.saved = space.save

multistack.bar2 <- function(x2=list(x), betweenspace = 2, withinspace=0, group.names = NULL, cex.axis=1.3,...){
  # note that number of rows in each component matrix of the list
  # should be equal
  # if you have missing value use "NA", to make it complete 
  # number of column can differ
  mylist <- x2
  space = space.saved
  un.space <- c(unlist(space))
  newdata <-unlist(mylist)
  b.pos= barplot(newdata, space= un.space, ...)
  
  if(is.null(group.names)){}else{
    tg = list()
    
    for (i in 1:length(unique(group.names))) {
      t.i = group.names[group.names == unique(group.names)[i]]
      if((length(t.i) %% 2) == 0){
        tg[[i]] = c(rep('', ceiling(length(t.i)-1)/2), rep(as.character(unique(group.names)[i]), 1), rep('', ceiling(length(t.i))/2))
      }else{
        tg[[i]] = c(rep('', floor(length(t.i)-1)/2), rep(as.character(unique(group.names)[i]), 1), rep('', ceiling(length(t.i)-1)/2))
        
      }
    }
    
    lab = unlist(tg)
    
    par(mgp=c(2,1,0))
    #axis(4, at=b.pos, labels = lab, tick = F, cex.axis=0.6, las=2)
    mapply(axis, side = 4, at=b.pos, labels = lab, tick = F, cex.axis=cex.axis, las=2, adj=0.1)
    par(mgp=c(3,1,0))
  }
  
}
multistack.bar2(x2 = wing.list, betweenspace=1.2, withinspace=0.1, las=1, 
                group.names = as.character(spT_order$HostPlantForm),
               col = 'burlywood3', border= 'burlywood4', xlab="", ylab="", xaxt='n', legend = '',  
               horiz=T, cex.lab =1.7, cex.axis=1.4, yaxt='n',
               main = '', cex.main = 1.7, cex.names = 3, 
               args.legend=list(x ='topleft', bty='n', cex=1.15, inset=c(0,-0.54)))
axis(1, at = c(0, 20, 40, 60), labels = c(0, 20, 40, 60), pos=0, cex.axis=1.7)
title('b) Wingsize (mm)', adj = 0.3, line = 0, cex.main = 1.7)


##..................................................................................AB
## Finland
par(mar = c(1,2,6,1))

VPAB_order = MeanA.vp[[2]][, c(spT_order$species)]

ng = dim(VPAB_order)[1]
leg2 = rownames(VPAB_order)

mylist2AB <- list(VPAB_order[,which(spT_order$HostPlantForm == 'Herbaceous')], VPAB_order[,which(spT_order$HostPlantForm == 'Graminoids')], 
                  VPAB_order[,which(spT_order$HostPlantForm == 'Others')], #VPAB_order[,which(spT_order$HostPlantForm == '5')],
                  VPAB_order[,which(spT_order$HostPlantForm == 'Dwarf shrubs')],
                  VPAB_order[,which(spT_order$HostPlantForm == 'Deciduous')], VPAB_order[,which(spT_order$HostPlantForm == 'Coniferous')])


#mylist2ABneg = lapply(mylist2AB, function(x){(-1)*x})
#mylist2ABneg = lapply(mylist2AB, function(y) apply(y, 2, function(x) (-1)*x / sum(x, na.rm = TRUE)))


mn = Mean.Fin[[2]]
means2 = round(100 *mn, 1)
for (j in 1:ng) {
  leg2[j] = paste(leg[j], " (mean = ", toString(means2[j]),
                  ")", sep = "")
}
multistack.bar(mylist2AB, betweenspace=1.2, withinspace=0.1, las=1, group.names = rep('', 77),
               col = coul, xlab="", ylab="", xaxt='n', legend = leg2,  horiz=T, cex.lab =1.7, cex.axis=1.7, yaxt='n',
               main = '', cex.main = 1.7, cex.names = 2,
               args.legend=list(x ='topright', bty='n', cex=1.4, inset=c(0,-0.08), cex=1.6))
axis(1, at = c(0, 0.2, 0.4, 0.6, 0.8, 1), labels = c(0, 0.2, 0.4, 0.6, 0.8, 1), pos=0, cex.axis=1.7)
title('c) VP Abundance', adj = 0.05, line = 0, cex.main = 1.7)






### ----------------------------------------   SITES
setwd("D:/Helsinki/Moth_project/ResultsAug23/RDATA/3.ModelRun")
load('VPall_Sites1808.RDATA')


par(mfrow=c(1,3), mar=c(1,3,2,1), oma=c(1,0,1,0))
layout(matrix(c(1, 2, 3), ncol=3), widths=c(2,1,2))

par(mar = c(1,1,6,2))
##..................................................................................PA
# 
VP2_order_sites = VP2_order_sites[, c(spT_order$species)]

ng = dim(VP2_order_sites)[1]
leg = rownames(VP2_order_sites)

mylist2 <- list(VP2_order_sites[,which(spT_order$HostPlantForm == 'Herbaceous')], 
                VP2_order_sites[,which(spT_order$HostPlantForm == 'Graminoids')], 
                VP2_order_sites[,which(spT_order$HostPlantForm == 'Others')], #VP2_order[,which(spT_order$HostPlantForm == '5')],
                VP2_order_sites[,which(spT_order$HostPlantForm == 'Dwarf shrubs')],
                VP2_order_sites[,which(spT_order$HostPlantForm == 'Deciduous')], 
                VP2_order_sites[,which(spT_order$HostPlantForm == 'Coniferous')])
mylist2neg = lapply(mylist2, function(x){(-1)*x})
#mylist2neg = lapply(mylist2, function(y) apply(y, 2, function(x) (-1)*x / sum(x, na.rm = TRUE)))


mn = Mean.Sites[[1]]
means2 = round(100 * mn, 1)

leg2 = c()
for (j in 1:ng) {
  leg2[j] = paste(leg[j], " (mean = ", toString(means2[j]),
                  ")", sep = "")
}
multistack.bar(mylist2neg, betweenspace=1.2, withinspace=0.1, las=1, group.names = '',
               col = coul, xlab="", ylab="", xaxt='n', legend = leg2,  horiz=T, cex.lab =1.7, cex.axis=1.7, yaxt='n',
               main = '', cex.main = 1.7, cex.names = 2,
               args.legend=list(x ='topleft', bty='n', cex=1.4, inset=c(0,-0.08), cex=1.6))
axis(1, at = c(-1, -0.8, -0.6, -0.4, -0.2, 0), labels = c(1, 0.8, 0.6, 0.4, 0.2, 0), pos=0, cex.axis=1.7)
title('a) VP Presence Absence', adj = 0.95, line = 0, cex.main = 1.7)

#### middle plot
par(mar = c(1,1,6,7))

spt2$row = 1:78

# barplot(spt2$WingspanF~spt2$row, at= c(1:22, 24:30, 32:38, 40:55, 57:78, 80:83),
#         horiz=TRUE)

wing.list = list(spt2[which(spt2$HostPlantForm == 'Herbaceous'),]$WingspanF.x, 
                 spt2[which(spt2$HostPlantForm == 'Graminoids'),]$WingspanF.x, 
                 spt2[which(spt2$HostPlantForm == 'Others'),]$WingspanF.x, 
                 spt2[which(spt2$HostPlantForm == 'Dwarf shrubs'),]$WingspanF.x,
                 spt2[which(spt2$HostPlantForm == 'Deciduous'),]$WingspanF.x, 
                 spt2[which(spt2$HostPlantForm == 'Coniferous'),]$WingspanF.x)

multistack.bar2(x2 = wing.list, betweenspace=1.2, withinspace=0.1, las=1, group.names = as.character(spT_order$HostPlantForm),
                col = 'burlywood3', border= 'burlywood4', xlab="", ylab="", xaxt='n', legend = '',  
                horiz=T, cex.lab =1.7, cex.axis=1.4, yaxt='n',
                main = '', cex.main = 1.7, cex.names = 3,
                args.legend=list(x ='topleft', bty='n', cex=1.15, inset=c(0,-0.14)))
axis(1, at = c(0, 20, 40, 60), labels = c(0, 20, 40, 60), pos=0, cex.axis=1.7)
title('b) Wingsize (mm)', adj = 0.3, line = 0, cex.main = 1.7)

##..................................................................................AB
## 
par(mar = c(1,2,6,1))

VPAB_order_sites = VPAB_order_sites[, c(spT_order$species)]

ng = dim(VPAB_order_sites)[1]
leg2 = rownames(VPAB_order_sites)

mylist2AB <- list(VPAB_order_sites[,which(spT_order$HostPlantForm == 'Herbaceous')], 
                  VPAB_order_sites[,which(spT_order$HostPlantForm == 'Graminoids')], 
                  VPAB_order_sites[,which(spT_order$HostPlantForm == 'Others')], #VPAB_order[,which(spT_order$HostPlantForm == '5')],
                  VPAB_order_sites[,which(spT_order$HostPlantForm == 'Dwarf shrubs')],
                  VPAB_order_sites[,which(spT_order$HostPlantForm == 'Deciduous')], 
                  VPAB_order_sites[,which(spT_order$HostPlantForm == 'Coniferous')])


#mylist2ABneg = lapply(mylist2AB, function(x){(-1)*x})
#mylist2ABneg = lapply(mylist2AB, function(y) apply(y, 2, function(x) (-1)*x / sum(x, na.rm = TRUE)))


mn = Mean.Sites[[2]]
means2 = round(100 * mn, 1)
for (j in 1:ng) {
  leg2[j] = paste(leg[j], " (mean = ", toString(means2[j]),
                  ")", sep = "")
}
multistack.bar(mylist2AB, betweenspace=1.2, withinspace=0.1, las=1, group.names = rep('', 77),
               col = coul, xlab="", ylab="", xaxt='n', legend = leg2,  horiz=T, cex.lab =1.7, cex.axis=1.7, yaxt='n',
               main = '', cex.main = 1.7, cex.names = 2,
               args.legend=list(x ='topright', bty='n', cex=1.4, inset=c(0,-0.08), cex=1.6))
axis(1, at = c(0, 0.2, 0.4, 0.6, 0.8, 1), labels = c(0, 0.2, 0.4, 0.6, 0.8, 1), pos=0, cex.axis=1.7)
title('c) VP Abundance', adj = 0.05, line = 0, cex.main = 1.7)

