#### Make character maps of locomotory regimes for forelimb traits ####

#### libraries ####
library(tidyverse)
library(plyr)
library(plotrix)
library(geiger)
library(phytools)
library(phangorn)

#### read data ####
data  <- read.csv("data_fore_loc3b.csv") #read species mean data
rownames(data) <- data[,1]

#Read multitree
multitree <- read.nexus(file = "Upham_tree_batforelimbs_V2.nex")
multitree <- lapply(multitree, force.ultrametric)
multitree <- lapply(multitree, ladderize)

multitree_prune <- list(length=length(multitree))
for(i in 1:length(multitree)){
  multitree_prune[[i]] <- treedata(phy = multitree[[i]], data = data, warnings=TRUE)$phy
}
class(multitree_prune) <- "multiPhylo"

#### obtain mcc tree from 1000 trees #### 
tree_mcc <- read.tree("Upham_mcctree_batforelimbs_V2.tre") #read mcc tree
tree_prune <- treedata(tree_mcc, data, sort=TRUE)$phy #check that all species are in tree (nothing should appear)

data <- data[tree_mcc$tip.label,] #Sorts the rows of data so that they match the order of the tip labels. SUPER IMPORTANT

#### variables ####
forelimb_lsr <- as.matrix(select(data, ends_with("_lsr_fore"))) #select log shape ratio traits

loc3 <- data$Loc_mode_3groups
names(loc3) <- data[,1]

loc3b <- data$Loc_mode_3bgroups
names(loc3b) <- data[,1]

loc3c <- data$Loc_mode_3cgroups
names(loc3c) <- data[,1]

loc4 <- data$Loc_mode_4groups
names(loc4) <- data[,1]

#### simmaps ####
multitree_prune_500 <- sample(multitree_prune[1:1000],500) #subsample 500 of 1k  trees

### loc3 (generalists, arbs, gliders+flyers)  ###
simmap_loc3_5k <- make.simmap(multitree_prune_500,loc3, nsim=10) #10 simmap simulations across 500 trees
save(simmap_loc3_5k, file="simmap_loc3_5k_V2.Rdata")
simmap_loc3_5k_s250 <- sample(simmap_loc3_5k[1:5000],250) #subsample 250 of 5k simmap trees
save(simmap_loc3_5k_s250, file="simmap_loc3_5k_s250.Rdata")
simmap_loc3_mcc_s50_fore <- make.simmap(tree_prune,loc3, nsim=50) #10 simmap simulations across 1k trees
save(simmap_loc3_mcc_s50_fore, file="simmap_loc3_mcc_s50_fore.Rdata")

### loc3b (generalists, arb+gliders, flyers) ###
simmap_loc3b_5k <- make.simmap(multitree_prune_500,loc3b, nsim=10) #10 simmap simulations across 500 trees
save(simmap_loc3b_5k, file="simmap_loc3b_5k_V2.Rdata")
simmap_loc3b_5k_s250 <- sample(simmap_loc3b_5k[1:5000],250) #subsample 250 of 5k simmap trees
save(simmap_loc3b_5k_s250, file="simmap_loc3b_5k_s250.Rdata")
simmap_loc3b_mcc_s50_fore <- make.simmap(tree_prune,loc3b, nsim=50) #10 simmap simulations across 1k trees
save(simmap_loc3b_mcc_s50_fore, file="simmap_loc3b_mcc_s50_fore.Rdata")

### loc3c (generalists+arb, gliders, flyers) ###
simmap_loc3c_5k <- make.simmap(multitree_prune_500,loc3c, nsim=10) #10 simmap simulations across 500 trees
save(simmap_loc3c_5k, file="simmap_loc3c_5k_V2.Rdata")
simmap_loc3c_5k_s250 <- sample(simmap_loc3c_5k[1:5000],250) #subsample 250 of 5k simmap trees
save(simmap_loc3c_5k_s250, file="simmap_loc3c_5k_s250.Rdata")
simmap_loc3c_mcc_s50_fore <- make.simmap(tree_prune,loc3c, nsim=50) #10 simmap simulations across 1k trees
save(simmap_loc3c_mcc_s50_fore, file="simmap_loc3c_mcc_s50_fore.Rdata")


### loc 4 (already run) ###
simmap_loc4_5k <- make.simmap(multitree_prune_500,loc4, nsim=10) #10 simmap simulations across 500 trees
save(simmap_loc4_5k, file="simmap_loc4_5k_V2.Rdata")
simmap_loc4_5k_s250 <- sample(simmap_loc4_5k[1:5000],250) #subsample 250 of 5k simmap trees
save(simmap_loc4_5k_s250, file="simmap_loc4_5k_s250.Rdata")
simmap_loc4_mcc_s50_fore <- make.simmap(tree_prune,loc4, nsim=50) #10 simmap simulations across 1k trees
save(simmap_loc4_mcc_s50_fore, file="simmap_loc4_mcc_s50_fore.Rdata")