library("raster")
library("rgdal")
bio1 <- raster("wc2.0_10m_prec_01_europe_agg_fact_5.tif")
bio1
table(is.na(getValues(bio1)))
plot(bio1)

library(ggplot2)
library(plyr)

file <- "results_10_000_fast.txt"
df <- read.csv(file, sep="\t")

# Multiple plot function
#
# ggplot objects can be passed in ..., or to plotlist (as a list of ggplot objects)
# - cols:   Number of columns in layout
# - layout: A matrix specifying the layout. If present, 'cols' is ignored.
#
# If the layout is something like matrix(c(1,2,3,3), nrow=2, byrow=TRUE),
# then plot 1 will go in the upper left, 2 will go in the upper right, and
# 3 will go all the way across the bottom.
#
multiplot <- function(..., plotlist=NULL, file, cols=1, layout=NULL) {
  library(grid)
  
  # Make a list from the ... arguments and plotlist
  plots <- c(list(...), plotlist)
  
  numPlots = length(plots)
  
  # If layout is NULL, then use 'cols' to determine layout
  if (is.null(layout)) {
    # Make the panel
    # ncol: Number of columns of plots
    # nrow: Number of rows needed, calculated from # of cols
    layout <- matrix(seq(1, cols * ceiling(numPlots/cols)),
                     ncol = cols, nrow = ceiling(numPlots/cols))
  }
  
  if (numPlots==1) {
    print(plots[[1]])
    
  } else {
    # Set up the page
    grid.newpage()
    pushViewport(viewport(layout = grid.layout(nrow(layout), ncol(layout))))
    
    # Make each plot, in the correct location
    for (i in 1:numPlots) {
      # Get the i,j matrix positions of the regions that contain this subplot
      matchidx <- as.data.frame(which(layout == i, arr.ind = TRUE))
      
      print(plots[[i]], vp = viewport(layout.pos.row = matchidx$row,
                                      layout.pos.col = matchidx$col))
    }
  }
}


# Prior predictive sampling distribution
p1 <- ggplot(df, aes(x=r))  + 
  geom_density(kernel="gaussian") +
  stat_function(fun=dunif, args = list(1,20), geom = "area", fill= "green", alpha = 0.35) +
  stat_function(fun=dunif, args = list(1,20), lty = 2)

p2 <- ggplot(df, aes(x=k)) +
  geom_density(kernel="gaussian") +
  stat_function(fun=dunif, args = list(1,500), geom = "area", fill= "green", alpha = 0.35) +
  stat_function(fun=dunif, args = list(1,500), lty = 2)

p3 <- ggplot(df, aes(x=N0)) + 
  geom_density(kernel="gaussian") +
  stat_function(fun=dunif, args = list(3,15), geom = "area", fill= "green", alpha = 0.35) +
  stat_function(fun=dunif, args = list(3,15), lty = 2)


p4 <- ggplot(df, aes(x=a)) +
  geom_density(kernel="gaussian") +
  stat_function(fun=dunif, args = list(100,1000), geom = "area", fill= "green", alpha = 0.35) +
  stat_function(fun=dunif, args = list(100,1000), lty = 2)

multiplot(p1,p2,p3,p4)

loci_names <- names(df)[-c(1,2,3,4)]
ftds <- df[,loci_names]
stacked_ftds <- stack(ftds)
names(stacked_ftds) <- c("FTD", "marker")
ggplot(stacked_ftds, aes(x=FTD)) + geom_density(aes(group=marker, colour=marker))

# Rejection
quantiles <- lapply(ftds, quantile, 0.3)
library (plyr)
qt <- ldply (quantiles, data.frame)
sample <- df[apply(df [loci_names], 1, function(x) all(x < qt[,2])),]
head(sample)

p5 <- ggplot(sample, aes(x=r))  + 
  geom_density(kernel="gaussian") +
  stat_function(fun=dunif, args = list(1,20), geom = "area", fill= "green", alpha = 0.35) +
  stat_function(fun=dunif, args = list(1,20), lty = 2)

p6 <- ggplot(sample, aes(x=k)) +
  geom_density(kernel="gaussian") +
  stat_function(fun=dunif, args = list(1,500), geom = "area", fill= "green", alpha = 0.35) +
  stat_function(fun=dunif, args = list(1,500), lty = 2)

p7<- ggplot(sample, aes(x=N0)) + 
  geom_density(kernel="gaussian") +
  stat_function(fun=dunif, args = list(3,15), geom = "area", fill= "green", alpha = 0.35) +
  stat_function(fun=dunif, args = list(3,15), lty = 2)


p8 <- ggplot(sample, aes(x=a)) +
  geom_density(kernel="gaussian") +
  stat_function(fun=dunif, args = list(100,1000), geom = "area", fill= "green", alpha = 0.35) +
  stat_function(fun=dunif, args = list(100,1000), lty = 2)

multiplot(p5,p6,p7,p8)