### Geneland ####
library(Geneland)
library(PBSmapping)
library(RgoogleMaps)
coords <- read.table(file = "Coords.txt",sep = "\t")
colnames(coords) <- c("Y","X")
genetable <- read.csv(file = "Last_N_pops.stru",header = TRUE,sep="\t")
attr(coords, "projection") <- "LL"
attr(coords, "zone") <- NA
colnames(coords) <- c("Y","X")
coord.utm <- PBSmapping::convUL(coords) # converting Long/Lat to UTM

write.table(cbind(coords, pop = Dzik@pop), file = "pops.csv",
            sep = ",",quote = FALSE,col.names = NA)

gene <- genetable[,9:32]
folder <- "Geneland"
MCMC(coordinates = coord.utm, geno.dip.codom = gene, varnpop = TRUE,
               npopmax = 2, npopmin = 2, spatial = TRUE, freq.model = "Correlated",
               nit = 10000000, thinning = 1000,path.mcmc = folder)
PostProcessChain(coordinates = coord.utm, path.mcmc = paste0(folder,"/"),
                           nxdom = 100, nydom = 100, burnin = 1000)
Plotnpop(path.mcmc = paste0(folder,"/"), burnin = 1000)
GenlandRes <- cbind(coords,Clust = read.csv(file = paste0(folder,"/modal.pop.indiv.txt"),sep = " ",header = FALSE)[,3])
write.table(GenlandRes,file = paste0("Genland_opt.csv"), sep = ",",quote = FALSE,col.names = NA)
Cluster = read.csv(file = paste0(folder,"/modal.pop.indiv.txt"),sep = " ",header = FALSE)[,3]
Geneland_opt = cbind(coords,Cluster)
table(Cluster)

pdf("Geneland_K_opt.pdf", width = 10, height = 6)
ggmap::ggmap(N_map)+
  scale_fill_manual(values = as.character(glasbey(10))) +
  geom_point(data = data.frame(Geneland_opt),shape = 21, size = 4, aes(X, Y, fill = factor(Cluster))) +
  labs(fill = "Cluster")+
  coord_equal()
dev.off()


####### ####K2
folder <- "Geneland_K2"
MCMC(coordinates = coord.utm, geno.dip.codom = gene, varnpop = TRUE,
     npopmax = 2, npopmin = 2, spatial = TRUE, freq.model = "Correlated",
     nit = 10000000, thinning = 1000,path.mcmc = folder)
PostProcessChain(coordinates = coord.utm, path.mcmc = paste0(folder,"/"),
                 nxdom = 100, nydom = 100, burnin = 1000)
Plotnpop(path.mcmc = paste0(folder,"/"), burnin = 1000)
GenlandRes <- cbind(coords,Clust = read.csv(file = paste0(folder,"/modal.pop.indiv.txt"),sep = " ",header = FALSE)[,3])
write.table(GenlandRes,file = paste0("Genland_K2.csv"), sep = ",",quote = FALSE,col.names = NA)
Cluster = read.csv(file = paste0(folder,"/modal.pop.indiv.txt"),sep = " ",header = FALSE)[,3]
Geneland_opt = cbind(coords,Cluster)

pdf("Geneland_K2.pdf", width = 10, height = 6)
ggmap::ggmap(N_map)+
  scale_fill_manual(values = as.character(glasbey(2))) +
  geom_point(data = data.frame(Geneland_opt),shape = 21, size = 4, aes(X, Y, fill = factor(Cluster))) +
  labs(fill = "Cluster")+
  coord_equal()
dev.off()

###### ####K3
folder <- "Geneland_K3"
MCMC(coordinates = coord.utm, geno.dip.codom = gene, varnpop = TRUE, 
     npopmax = 3, npopmin = 3, spatial = TRUE, freq.model = "Correlated", 
     nit = 10000000, thinning = 1000,path.mcmc = folder)
PostProcessChain(coordinates = coord.utm, path.mcmc = paste0(folder,"/"), 
                 nxdom = 100, nydom = 100, burnin = 1000)
Plotnpop(path.mcmc = paste0(folder,"/"), burnin = 1000)
GenlandRes <- cbind(coords,Clust = read.csv(file = paste0(folder,"/modal.pop.indiv.txt"),sep = " ",header = FALSE)[,3])
write.table(GenlandRes,file = paste0("Genland_K3.csv"), sep = ",",quote = FALSE,col.names = NA)
Cluster = read.csv(file = paste0(folder,"/modal.pop.indiv.txt"),sep = " ",header = FALSE)[,3]
Geneland_opt = cbind(coords,Cluster)

pdf("Geneland_K3.pdf", width = 10, height = 6)
ggmap::ggmap(N_map)+
  scale_fill_manual(values = as.character(glasbey(3))) +
  geom_point(data = data.frame(Geneland_opt),shape = 21, size = 4, aes(X, Y, fill = factor(Cluster))) +
  labs(fill = "Cluster")+
  coord_equal()
dev.off()