# load libraries
library(tidyverse)
library(RColorBrewer)
library(ozmaps)
library(doParallel)
library(sf)
library(usdm)
library(hypervolume)
library(rgdal)
library(rgeos)
library(raster)
library(sf)
library(alphahull)
library(sp)
library(foreach)
library(parallel)

# set function for converting alpha hulls to polygon
require(alphahull)
require(sp)
ashape_to_SPLDF <- function(x, proj4string=NA)
{
  if(class(x) != 'ashape')
    stop('this function only works with `ashape` class objects')
  
  # convert ashape edges to DF
  x.as.df <- as.data.frame(x$edges)
  
  # convert each edge to a line segment
  l.list <- list()
  for(i in 1:nrow(x.as.df))
  {
    # extract line start and end points as 1x2 matrices
    p1 <- cbind(x.as.df$x1[i], x.as.df$y1[i])
    p2 <- cbind(x.as.df$x2[i], x.as.df$y2[i])
    # row-bind into 2x3 matrix
    l.list[[i]] <- Line(rbind(p1, p2))
  }
  
  # promote to Lines class, then to SpatialLines class
  l <- Lines(l.list, ID=1)
  
  # copy over CRS data from original point data
  l.spl <- SpatialLines(list(l), proj4string=CRS(as.character(NA)))
  
  # promote to SpatialLinesDataFrame, required for export to GRASS / OGR
  l.spldf <- SpatialLinesDataFrame(l.spl, data=data.frame(id=1), match.ID=FALSE)
  
  return(l.spldf)
}

# set working directory
#setwd('/Users/jarrodsopniewski/Library/CloudStorage/OneDrive-AustralianNationalUniversity/OneDrive Folder - ANU/RA/NEE-submission/Updated_run/') 

setwd('')

# read in occs
occs <- read.csv('Data/final_frog_data.csv')
species_list <- c('Litoria_aurea',
                  'Litoria_booroolongensis',
                  'Litoria_dayi',
                  'Litoria_littlejohni',
                  'Litoria_nannotis',
                  'Litoria_raniformis',
                  'Litoria_rheocola',
                  'Litoria_spenceri',
                  'Litoria_subglandulosa',
                  'Mixophyes_balbus',
                  'Mixophyes_fleayi',
                  'Mixophyes_iteratus',
                  'Philoria_frosti',
                  'Taudactylus_eungellensis',
                  'Adelotus_brevis',
                  'Geocrinia_victoriana',
                  'Litoria_lesueurii',
                  'Litoria_nudidigita',
                  'Litoria_pearsoniana',
                  'Litoria_serrata',
                  'Litoria_verreauxii',
                  'Litoria_wilcoxii',
                  'Pseudophryne_bibronii',
                  'Pseudophryne_dendyi',
                  'Taudactylus_liemi',
                  'Litoria_brevipalmata',
                  'Litoria_cooloolensis',
                  'Litoria_freycineti',
                  'Litoria_olongburensis',
                  'Philoria_kundagungan',
                  'Philoria_loveridgei',
                  'Philoria_sphagnicola',
                  'Pseudophryne_australis',
                  'Litoria_burrowsae',
                  'Litoria_chloris',
                  'Litoria_citropa',
                  'Litoria_dentata',
                  'Litoria_ewingii',
                  'Litoria_fallax',
                  'Litoria_gracilenta',
                  'Litoria_infrafrenata',
                  'Litoria_jervisiensis',
                  'Litoria_latopalmata',
                  'Litoria_paraewingi',
                  'Litoria_peronii',
                  'Litoria_phyllochroa',
                  'Litoria_revelata',
                  'Litoria_tyleri',
                  'Litoria_xanthomera',
                  'Mixophyes_fasciolatus',
                  'Mixophyes_schevilli',
                  'Pseudophryne_coriacea',
                  'Pseudophryne_major',
                  'Pseudophryne_semimarmorata',
                  'Litoria_nigrofrenata')

occs <- read.csv('final_frog_data.csv')

# world borders available at https://koordinates.com/layer/7354-tm-world-borders-03/
world <- readOGR(dsn = '',
                 layer = 'TM_WORLD_BORDERS-0.3')
aus <- world[world@data$name=='Australia',]

# read in 30 second resolution raster layers, available from http://www.worldclim.com/version2
hv_stack <- stack(raster('wc2.1_30s_bio/wc2.1_30s_bio_1.tif'),
                  raster('wc2.1_30s_bio/wc2.1_30s_bio_2.tif'),
                  raster('wc2.1_30s_bio/wc2.1_30s_bio_12.tif'),
                  raster('wc2.1_30s_bio/wc2.1_30s_elev.tif'))
hv_stack <- crop(hv_stack, aus)
hv_stack <- mask(hv_stack, aus)
box <- extent(hv_stack)
box@xmin <- 141
hv_stack <- crop(hv_stack, box)

# get xy coordinates for Australia
xy_aus <- xyFromCell(hv_stack[[1]], cell = 1:ncell(hv_stack[[1]]),
                     spatial = TRUE)

# set number of repetitions
reps_null <- 1000 # (how many null distributions)
reps_base <- 200 # (how many replications of the base hypervolumes)
# set number of repetitions
reps_null <- 1000 # (how many null distributions)
reps_base <- 200 # (how many replications of the base hypervolumes)

# set up parallel
n.cores = detectCores()
n.cores = round(n.cores*0.7)
my.cluster <- parallel::makeCluster(n.cores, 
                                    type = "PSOCK")
registerDoParallel(cl = my.cluster)

# do loop
for(i in 1:length(species_list)) {
  tryCatch({
    soi <- species_list[i]
    
    print(paste('Commencing Species ',
                i, 
                ' (',
                soi,
                ') at ',
                Sys.time(),
                sep = ''))
    # isolate data and create hulls (data already spatially thinned)
    species_data_historical <- subset(occs, Species == soi)
    coordinates(species_data_historical) <- ~Lon+Lat
    cells <- cellFromXY(hv_stack, species_data_historical)
    dups <- duplicated(cells)
    species_data_historical <- as.data.frame(species_data_historical[!dups,])
    
    ashape_data_historical <- ashape(x = species_data_historical$Lon, 
                                     y = species_data_historical$Lat, 
                                     alpha = 2)
    hull_data_historical <- ashape_to_SPLDF(ashape_data_historical, 
                                            proj4string = "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")
    hull_data_historical <- st_as_sf(hull_data_historical)
    hull_data_historical <- st_polygonize(hull_data_historical)
    hull_data_historical <- st_cast(hull_data_historical, warn = FALSE)
    hull_data_historical <- as(hull_data_historical, "Spatial")
    proj4string(hull_data_historical) <-"+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
    hull_data_historical <- crop(hull_data_historical, aus)
    area_hull_historical <- sum(raster::area(hull_data_historical))
    
    species_data_contemporary <- subset(species_data_historical, Period == 'Contemporary')
    ashape_data_contemporary <- ashape(x = species_data_contemporary$Lon, 
                                       y = species_data_contemporary$Lat, 
                                       alpha = 2)
    hull_data_contemporary <- ashape_to_SPLDF(ashape_data_contemporary, 
                                              proj4string = "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")
    hull_data_contemporary <- st_as_sf(hull_data_contemporary)
    hull_data_contemporary <- st_polygonize(hull_data_contemporary)
    hull_data_contemporary <- st_cast(hull_data_contemporary, warn = FALSE)
    hull_data_contemporary <- as(hull_data_contemporary, "Spatial")
    proj4string(hull_data_contemporary) <-"+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
    hull_data_contemporary <- crop(hull_data_contemporary, aus)
    area_hull_contemporary <- sum(raster::area(hull_data_contemporary))
    
    
    # extract data for points
    values <- cbind(as.data.frame(raster::extract(hv_stack,
                                                  species_data_historical[,c('Lon', 'Lat')])),
                    species_data_historical$Period)
    values <- values[complete.cases(values),]
    
    # scale variables
    values_scaled <- cbind(as.data.frame(scale(values[,1:4], scale = TRUE, center = TRUE)), 
                           values[,5])
    colnames(values_scaled) <- c('bio1', 'bio2', 'bio12', 'elevation', 'period')
    # construct hypervolumes
    # construct contemporary hypervolumes
    print(paste('Commencing contemporary HV construction at',
                Sys.time()))
    contemporary_hypervolume <- foreach(j = 1:reps_base,
                                        .combine = 'c',
                                        .packages = 'hypervolume',
                                        .inorder = FALSE,
                                        .verbose = T) %dopar% {
                                          hypervolume_svm(data = subset(values_scaled, period == 'Contemporary')[1:4],
                                                          svm.nu = 0.01,
                                                          svm.gamma = 2.5,
                                                          verbose = F)
                                        }
    
    # construct null distribution of hypervolumes
    print(paste('Commencing null distribution HV construction at',
                Sys.time()))
    null_hypervolume <- foreach(j = 1:reps_null,
                                .combine = 'c',
                                .packages = 'hypervolume',
                                .inorder = FALSE,
                                .verbose = T) %dopar% {
                                  hypervolume_svm(data = values_scaled[sample(1:nrow(values_scaled),
                                                                              size = nrow(subset(values_scaled, period == 'Contemporary'))),1:4],
                                                  svm.nu = 0.01,
                                                  svm.gamma = 2.5,
                                                  verbose = F)
                                }                                
    
    # save HV objects
    save(contemporary_hypervolume,
         null_hypervolume,
         file = paste0('Results/HV_objects/',
                       soi,
                       '_hv_objects_3march.RData'))
    
    # get results
    # contemporary results
    contemporary_results <- list()
    for(j in 1:length(contemporary_hypervolume)) {
      centroids <- get_centroid(contemporary_hypervolume[[j]])
      contemporary_results[[j]] <- data.frame(Species = soi,
                                              Area_EOO = area_hull_contemporary,
                                              Niche_Volume = get_volume(contemporary_hypervolume[[j]]),
                                              bio1 = centroids[1],
                                              bio2 = centroids[2],
                                              bio12 = centroids[3],
                                              elevation = centroids[4]) 
    }
    contemporary_results <- bind_rows(contemporary_results)
    contemporary_results$Period = 'Contemporary'
    
    # contemporary results
    null_results <- list()
    for(j in 1:length(null_hypervolume)) {
      centroids <- get_centroid(null_hypervolume[[j]])
      null_results[[j]] <- data.frame(Species = soi,
                                      Area_EOO = area_hull_historical,
                                      Niche_Volume = get_volume(null_hypervolume[[j]]),
                                      bio1 = centroids[1],
                                      bio2 = centroids[2],
                                      bio12 = centroids[3],
                                      elevation = centroids[4]) 
    }
    null_results <- bind_rows(null_results)
    null_results$Period = 'Null'
    results <- rbind(contemporary_results,
                     null_results)
    write.csv(results,
              paste0('Results/Numeric_results/',
                     soi,
                     '_amalgamated_results_3march.csv'),
              row.names = FALSE)
  },error=function(e){cat("mistake:", conditionMessage(e),"\n")})
}