# to subset a vector of columns
subset_cols <- function(df, col.list){
  idx <- match(col.list, names(df))
  if(is.data.table(df)==T){
    df <- df[,na.omit(idx), with = F]
  }else{
    df <- df[,na.omit(idx)]
  }
  return(df)
}

# to plot an ellipse (dependent function)
veganCovEllipse<-function (cov, center = c(0, 0), scale = 1, npoints = 100) 
{
  theta <- (0:npoints) * 2 * pi/npoints
  Circle <- cbind(cos(theta), sin(theta))
  t(center + scale * t(Circle %*% chol(cov)))
}

# to plot data and color based on group, rename group, and plot ellipse
nmds_ellipse_plot<-function(group1, group2 = NA, MDS, group.names = NA, 
                            xlim = NA, ylim = NA){
  
  if(is.na(group2[1])){
    
    NMDS = data.frame(NMDS1 = MDS$points[,1], 
                      NMDS2 = MDS$points[,2],
                      group = as.factor(group1))
  }else{
    NMDS = data.frame(NMDS1 = MDS$points[,1], 
                      NMDS2 = MDS$points[,2],
                      group1 = as.factor(group1),
                      group2 = as.factor(group2),
                      group = interaction(group1, group2))
    NMDS$group<-droplevels(NMDS$group)
  }
  
  
  ord<-plot.ellipse.invisible(MDS,NMDS$group)
  
  df_ell <- data.frame()
  for(g in levels(NMDS$group)){
    df_ell <- rbind(df_ell, 
                    cbind(as.data.frame(
                      with(NMDS[NMDS$group==g,],
                           veganCovEllipse(ord[[g]]$cov,
                                           ord[[g]]$center,
                                           ord[[g]]$scale))), 
                      group=g))
  }
  
  
  p<-ggplot() + 
    geom_point(data = NMDS, aes(NMDS1, NMDS2), color = "grey", size = 0.75)
  p<-p +
    geom_point(data = subset(NMDS, !is.na(group)), 
               aes(NMDS1, NMDS2, color = group), size = 0.85) +
    geom_path(data=df_ell, 
              aes(x=NMDS1, y=NMDS2,colour=group), 
              size=0.85, linetype=1)+
    theme_classic()+ theme(legend.title=element_blank())
  
  if(!is.na(group.names[1])){
    grps<-levels(NMDS$group)
    p<- p + scale_colour_discrete(breaks=c(grps),
                                  labels=c(group.names)) 
  }
  
  if(!(is.na(ylim[1]))){
    p<- p+ xlim(xlim)+ylim(ylim)
  }
  
  print(p)
  
}

# dependent function for ellipse plotting
plot.ellipse.invisible <- function(MDS, MDSgroup){
  ff <- tempfile()
  png(filename=ff)
  plot.new()
  ord<-ordiellipse(MDS, MDSgroup, display = "sites", draw = "none", 
                   kind = "se", conf = 0.90, label = F)
  dev.off()
  unlink(ff)
  return(ord)
}
