DS.Rmd

---
title: "DS[1,2,3,3s,4,5]"
date: '`r format(Sys.time(), "%d %B, %Y")`'
output: pdf_document
toc: true
---

```{r setup, include=FALSE}
library(reshape2)
library(ggplot2)
library(knitr)
library(tidyr)
library(dplyr)
library(kableExtra)
library(ape)
knitr::opts_chunk$set(echo = TRUE, fig.path='Figs/')
```
\newpage

# The dubious methods

- SRF: posterior distribution cacluated with golden runs
- NMC: naive Monte Carlo with 10,000 samples
- ELBO: evidence lower bound (mean field Gaussian)
- HM: harmonic mean
- SHM: stabilized harmonic mean
- SS: stepping stone with 50 power posteriors
- GSS: generalized stepping stone with 50 power posterior. Reference distribution is independent gamma distributions.
- PS: path sampling with beta(0.3,1) and 50 power posteriors
- PS2: modified path sampling with beta(0.3,1) and 50 power posteriors
- ML: maximum likleihood
- MAP: maximum a posteriori using exp(10) as prior on each branch
- BS: bridge sampling
- NS: nested sampling
- CPO: conditional predictive ordinate
- LPPD: log pointwise predictive density
- GL: gamma Laplace
- LL: lognormal Laplace
- BL: beta' Laplace
- GLIS: gamma Laplace with importance sampling (10,000 samples)
- VBIS: Importance sampling using a variational distribution as the importance distribution 10,000 samples

```{r}
log.sum.exp<- function(x) {
  # Computes log(sum(exp(x))
  # Uses offset trick to avoid numeric overflow: http://jblevins.org/notes/log-sum-exp
  # if ( max(abs(x)) > max(x) )
  #   offset <- min(x)
  # else
    offset <- max(x)
  log(sum(exp(x - offset))) + offset
}

kl.log.pq <- function(P,Q, normalize=FALSE){
  if(normalize){
    Q = Q - log.sum.exp(Q)
    P = P - log.sum.exp(P)
  }
  sum(exp(P)*(P-Q))
}

normalize <- function(P){
  return(P - log.sum.exp(P))
}

h.log.pq <-function(P, Q, normalize=FALSE){
  if(normalize){
    Q = Q - log.sum.exp(Q)
    P = P - log.sum.exp(P)
  }
  sqrt(sum( (sqrt(exp(P)) - sqrt(exp(Q)))^2 ))/sqrt(2)
}

rmsd <- function(a,b){
  sqrt(mean((a-b)^2))
}

rmsd2 <- function(tree, P, Q, ptrees, pp){
  Q = Q - log.sum.exp(Q)
  P = P - log.sum.exp(P)
  
  Q = exp(Q[order(tree)])
  P = exp(P[order(tree)])
  
  msd = 0
  for(i in 1:length(pp)){
    msdi = 0
    for(j in 1:length(ptrees)){
      if(list(pp[[i]]) %in% ptrees[[j]]){
        msdi = msdi + Q[j] - P[j]
      }
    }
    msd = msd + msdi*msdi
  }
  sqrt(msd/length(pp))
}

splits <- function(data, ptrees){
  Q = data$marginal - log.sum.exp(data$marginal)
  P = data$SRF - log.sum.exp(data$SRF)
  
  Q = exp(Q[order(data$tree)])
  P = exp(P[order(data$tree)])
  splits = list()
  f = c()
  fhat = c()
  count = c()
  # iterate over every treee
  for(i in 1:length(ptrees)){
    # iterate over every splits
    for(k in 1:length(ptrees[[i]])){
      found = F
      bipart = ptrees[[i]][[k]]
      if(length(splits) > 0){
        for(j in 1:length(splits)){
            if(list(bipart) %in% list(splits[[j]])){
              found = T
              f[j] = f[j]+P[i]
              fhat[j] = fhat[j]+Q[i]
              count[j] = count[j] + 1
              break
            }
        }
      }
      else
        j = 0
      
      if(!found){
        j = j + 1
        f[j] = P[i]
        fhat[j] = Q[i]
        splits[[j]] = bipart
        count[j] = 1
      }
    }
  }
  list(data=data.frame(f=f, fhat=fhat,count=count), splits=splits)
}

read.ds <- function(dataset, force=F){
  #trees.file  = file.path('data', paste0(dataset, ".trees", sep=""))
  trees.file = paste0('/Users/mathieu/Desktop/marginal-experiments/data/JC_no_gamma_credible_set_ds', substr(dataset, nchar(dataset), nchar(dataset)))
  if(file.exists(paste0(dataset, '.Rdata')) & force!=T){
    load(paste0(dataset, '.Rdata'))
    return(res)
  }
  
  df = read.csv(file.path(dataset, paste(dataset, ".csv", sep="")), sep='\t')
  
  # total time
  total = sum(df$time)

  srf = read.csv(file.path(dataset, "data.csv"), sep='\t')
  srf = srf %>% mutate(SRF=log(SRF))
  df = mutate(df, rep=rep+1)

  dfff = df %>% filter(algorithm %in% c('mcmc', 'mmcmc', 'mmcmc-gss'))%>%
    spread(algorithm,time) %>% select(-marginal)
  df = filter(df, !(algorithm %in% c('mcmc', 'mmcmc', 'mmcmc-gss')))
  dfff = dfff %>% right_join(df, by=c('rep','tree'))
  #MMCMC
  for(m in c('SS', 'PS', 'PS2')){
    dfff[dfff$algorithm==m,]$time  = dfff[dfff$algorithm==m,]$time + 
      dfff[dfff$algorithm==m,]$mmcmc
  }
  #MCMC
  for(m in c('BS', 'CPO', 'LPPD', 'HM', 'SHM')){
    dfff[dfff$algorithm==m,]$time  = dfff[dfff$algorithm==m,]$time +
      dfff[dfff$algorithm==m,]$mcmc
  }
  dfff[dfff$algorithm=='GSS',]$time  = dfff[dfff$algorithm=='GSS',]$time +
    dfff[dfff$algorithm=='GSS',]$`mmcmc-gss`
  
  df.time = dfff %>% group_by(algorithm) %>%
    summarise(median=median(time),mean=mean(time), sd=sd(time), cv=sd(time)/mean(time)) %>%
    arrange(algorithm) %>% mutate(dataset=dataset)
  
  # KL
  df.kl = dfff %>% full_join(srf, by='tree') %>% group_by(rep, algorithm) %>%
    summarize(KL=kl.log.pq(SRF, marginal, TRUE), H=h.log.pq(SRF, marginal, TRUE)) %>%
    arrange(algorithm) %>% mutate(dataset=dataset)
  
  #RMSD
  trees = read.tree(trees.file, keep.multi =TRUE)
  trees = .compressTipLabel(trees)
  ntree <- length(trees)
  class(trees) <- NULL
  for (i in 1:ntree) storage.mode(trees[[i]]$Nnode) <- "integer"
  class(trees) <- "multiPhylo"
  trees <- reorder(trees, "postorder")
  ptrees = lapply(trees, prop.part)
  
  
  dfff = dfff %>% group_by(rep, algorithm) %>% mutate(Posterior=normalize(marginal)) %>% 
    ungroup() %>% mutate(dataset=dataset)
  
  temp = dfff %>% full_join(srf, by='tree')
  split.freqs = NULL
  replicates = unique(as.character(temp$rep))
  for(replicate in replicates){
    for(a in unique(as.character(temp$algorithm))){
      data = temp %>% filter(algorithm==a&rep==replicate) %>% select(tree, SRF, marginal)
      res = splits(data, ptrees)
      split.freq = res[['data']]
      bips = res[['splits']]
      split.freq = cbind(split.freq, algorithm=rep(a, nrow(split.freq)), rep=rep(replicate, nrow(split.freq)))
      if(is.null(split.freqs))
        split.freqs = split.freq
      else
        split.freqs = rbind(split.freq,split.freqs)
    }
  }
  df.rmsd = split.freqs %>% group_by(algorithm, rep) %>% summarize(RMSD=rmsd(f,fhat)) %>%
    arrange(algorithm) %>% mutate(dataset=dataset)
  
  # plot
  myorder = df.kl %>% filter(rep==1) %>% arrange(KL) %>% ungroup() %>%
    select(algorithm) %>% unlist(use.names = FALSE) %>% as.vector
  
  df.plot = dfff %>% full_join(srf, by='tree') %>%
    mutate(algorithm=factor(algorithm, levels=myorder), dataset=dataset)

  res = list(time=df.time, kl=df.kl, plot=df.plot, marginals=dfff,rmsd=df.rmsd,total=total,splits=split.freqs)
  save(res, file=paste0(dataset, '.Rdata'))
  res
}

f = function(.o,.f,.kl,.rmsd) paste0('atop(atop(',
                                     'textstyle(',deparse(.f),'),',
                                     'textstyle(RMSD == "',formatC(.rmsd, format='e', digits=1),'")),',
                                     'textstyle(KL == "',formatC(.kl, format='e', digits=1),'"))')

datasets = paste0('DS', 1:5)

origin=c('GSS', 'GLIS', 'VBIS', 'BS', 'SS', 'PS', 'PS2', 'LL', 'ML', 'MAP', 'GL',
                             'ELBO', 'LPPD', 'BL', 'CPO', 'SHM', 'HM', 'NS', 'NMC')
final=c("GSS", "GLIS", "VBIS", "BS", "SS", "PS", "MPS", "LL", "ML", "MAP", "GL",
                             "ELBO", "PPD", "BL", "CPO", "SHM", "HM", "NS", "NMC")

# final=c("GSS", expression("LG"["IS"]), expression("VB"["IS"]), "BS", "SS", "PS", "PSm", "LL", "ML", "MAP", "LG",
#                              "ELBO", "PPD", "LB", "CPO", "SHM", "HM", "NS", "NMC")
# final2=c("GSS", bquote("LG"["IS"]), bquote("VB"["IS"]), "BS", "SS", "PS", "PSm", "LL", "ML", "MAP", "LG",
#                              "ELBO", "PPD", "LB", "CPO", "SHM", "HM", "NS", "NMC")

color_datasets = c(rgb(141,160,203,max=255), rgb(252,141,98,max=255), rgb(102,194,165,max=255),
                   rgb(225,198,47,max=255),#rgb(175,175,175,max=255),
                   rgb(204,121,167,max=255))
shape_datasets <- c(3,17,18,19,4)

theme_set(theme_bw(16) + theme(strip.background = element_blank()))
```



\newpage

# DS1 to DS5 together

```{r load-all}
df = read.ds(datasets[1])
df = lapply(df, as.data.frame)
all.kl = as.data.frame(df[['kl']])
all.rmsd = df[['rmsd']]
all.time = df[['time']]
all.plot = df[['plot']]
all.marginal = df[['marginals']]
total = df[['total']]
for(i in 2:length(datasets)){
  df = read.ds(datasets[i])
  total = total + df[['total']]
  df = lapply(df, as.data.frame)
  all.kl = rbind(all.kl, df[['kl']])
  all.rmsd = rbind(all.rmsd, df[['rmsd']])
  all.time = rbind(all.time, df[['time']])
  all.plot = rbind(all.plot, df[['plot']])
  all.marginal = rbind(all.marginal, df[['marginals']])
}
myrange<-function(x){max(x)-min(x)}
df.plot = all.plot %>% filter(algorithm=="GSS") %>% group_by(rep,dataset) %>%  summarize(range=myrange(marginal))  %>% spread(dataset, range)
kable(df.plot)
```


\newpage

## RMSD vs time
```{r RMSD_vs_time}
all.rmsd.time = all.rmsd %>% filter(!(algorithm %in% c('MAP', 'PS', 'PS2', 'BL'))) %>%
    group_by(algorithm, dataset) %>% summarise(medianRMSD=median(RMSD),meanRMSD=mean(RMSD)) %>% ungroup() %>%
    left_join(all.time, by=c('dataset','algorithm'))
all.rmsd.time = all.rmsd.time %>% mutate(d=substr(dataset, 3,3)) %>% mutate(label=paste0(algorithm, d))
label_formated= Map(.a=as.character(all.rmsd.time$algorithm), .b=all.rmsd.time$d, f = function(.a,.b) bquote(list(.(.a)^.(.b))))
all.rmsd.time$label2 = sapply(label_formated, deparse)

ggplot(all.rmsd.time, aes(x=mean, y=meanRMSD, col=algorithm,label=label2)) +
  scale_x_log10()+
  scale_y_log10() +
  geom_hline(yintercept = 0.01) +
  geom_hline(yintercept = 0.05, linetype="dashed") +
  geom_text(size=3, show.legend = FALSE, parse=TRUE) +
  xlab("Running time (seconds) per tree") +
  ylab("RMSD")
```


\newpage

## KL vs time
```{r kl_vs_time}
all.kl.time = all.kl %>% filter(!(algorithm %in% c('MAP', 'PS', 'PS2', 'BL'))) %>%
    group_by(algorithm, dataset) %>% summarise(medianKL=median(KL),meanKL=mean(KL)) %>% ungroup() %>%
    left_join(all.time, by=c('dataset','algorithm'))
all.kl.time = all.kl.time %>% mutate(d=substr(dataset, 3,3)) %>% mutate(label=paste0(algorithm, d))
label_formated= Map(.a=as.character(all.kl.time$algorithm), .b=all.kl.time$d, f = function(.a,.b) bquote(list(.(.a)^.(.b))))
all.kl.time$label2 = sapply(label_formated, deparse)
ggplot(all.kl.time, aes(x=mean, y=meanKL, col=algorithm,label=label2)) +
  scale_x_log10() +
  scale_y_log10() +
  geom_text(size=3, show.legend = FALSE, parse=TRUE) +
  xlab("Running time (seconds) per tree") +
  ylab("KL divergence")
```

\newpage

## RMSD of 10 replicates
```{r RMSD_by_method_and_dataset_with_replicates}
all.rmsd$rep = factor(all.rmsd$rep, levels=unique(all.rmsd$rep))
all.rmsd.determ = filter(all.rmsd, algorithm %in% c('ML', 'MAP', 'GL', 'LL', 'BL') & rep==1)
all.rmsd = filter(all.rmsd, !(algorithm %in% c('ML', 'MAP', 'GL', 'LL', 'BL')))
all.rmsd = rbind(all.rmsd, all.rmsd.determ)
myorder = all.rmsd %>% group_by(algorithm) %>% summarize(medianRMSD=median(RMSD),meanRMSD=mean(RMSD)) %>%
  arrange(meanRMSD) %>% select(algorithm) %>% unlist(use.names = FALSE) %>% as.vector

all.rmsd$Method = factor(all.rmsd$algorithm, levels=myorder)
legend = final[match(myorder, origin)]
ggplot(all.rmsd, aes(x=Method, y=RMSD,col=dataset)) +#, shape=dataset)) +
  theme(legend.position = c(0.9, 0.2),
        legend.title=element_blank(), legend.background=element_blank(),
        legend.key = element_rect(fill = NA, colour = NA, size = 0.25),
        axis.title.x=element_blank(),
        axis.text.x = element_text(angle = 30, hjust = 1)) +
  geom_hline(yintercept = 0.01) +
  geom_hline(yintercept = 0.05, linetype="dashed") +
  geom_jitter(width = 0.3, size=1) +
  scale_x_discrete(labels=legend) +
  scale_colour_manual(name="dataset", values=color_datasets) +
  #scale_shape_manual(name="dataset", values=shape_datasets) +
  scale_y_log10()
```

\newpage

## KL of 10 replicates
```{r kl_by_method_and_dataset_with_replicates}
all.kl$rep = factor(all.kl$rep, levels=unique(all.kl$rep))
all.kl.determ = filter(all.kl, algorithm %in% c('ML', 'MAP', 'GL', 'LL', 'BL') & rep==1)
all.kl = filter(all.kl, !(algorithm %in% c('ML', 'MAP', 'GL', 'LL', 'BL')))
all.kl = rbind(all.kl, all.kl.determ)
myorder = all.kl %>% group_by(algorithm) %>% summarize(medianKL=median(KL),meanKL=mean(KL)) %>%
  arrange(meanKL) %>% select(algorithm) %>% unlist(use.names = FALSE) %>% as.vector

all.kl$Method = factor(all.kl$algorithm, levels=myorder)
legend = final[match(myorder, origin)]
ggplot(all.kl, aes(x=Method, y=KL,col=dataset)) +#, shape=dataset)) +
  theme(legend.position = c(0.9, 0.2),
        legend.title=element_blank(), legend.background=element_blank(),
        legend.key = element_rect(fill = NA, colour = NA, size = 0.25),
        axis.title.x=element_blank(),
        axis.text.x = element_text(angle = 30, hjust = 1)) +
  geom_jitter(width = 0.3, size=1) +
  scale_x_discrete(labels=legend) +
  scale_colour_manual(name="dataset", values=color_datasets) +
  #scale_shape_manual(name="dataset", values=shape_datasets) +
  scale_y_log10() +
  ylab("KL divergence")
```

## Standard error
```{r marginal-se, fig.height=10, fig.width=10}
margs = all.marginal %>% group_by(dataset,algorithm, tree) %>%
  summarize(SE=sd(marginal), mean=mean(marginal),CV=sd(marginal)/mean(marginal)) %>%
  filter(SE!=0) %>% as.data.frame

 myorder = margs %>% group_by(algorithm) %>% summarize(meanSE=mean(SE)) %>% as.data.frame
 margs$Method = factor(margs$algorithm, levels=myorder$algorithm[order(myorder$meanSE,decreasing=F)])

mylabs = final[match(levels(margs$Method), origin)]

 ggplot(margs, aes( x=Method, y=SE)) +
   theme(axis.title.x=element_blank(),
        axis.text.x = element_text(angle = 30, hjust = 1)) +
   facet_wrap(~dataset, ncol=1, scales = "free_y") + 
   geom_boxplot() +
   scale_y_log10() +
   scale_x_discrete(labels=mylabs) +
   ylab("Standard error (log scale)")
```

\newpage

```{r}

plot.posterior.srf <- function(ds, replicate, breaks=NULL){
  df.plot = all.plot %>% filter(dataset==ds&rep==replicate) %>% select(-rep)
  df.rmsd = all.rmsd %>% filter(dataset==ds&rep==replicate) %>% select(-rep)
  df.kl = all.kl %>% filter(dataset==ds&rep==replicate) %>% mutate(algorithm=droplevels(algorithm)) %>% select(-rep)
  temp = df.kl %>% full_join(df.rmsd, by=c('algorithm', 'dataset','Method')) %>% arrange(RMSD)
  df.plot = df.plot %>% mutate(algorithm=factor(algorithm, levels=temp$algorithm))

  origin_ordered = origin[match(temp$algorithm, origin)]
  final_ordered = final[match(temp$algorithm, origin)]
  label_formated= Map(.o=origin_ordered, .f=final_ordered, .kl=temp$KL, .rmsd=temp$RMSD, f)
  levels(df.plot$algorithm) = unlist(label_formated)
  p = ggplot(df.plot, aes(x=SRF, y=Posterior)) +
    facet_wrap( ~algorithm, scales = "free_y", ncol=4,labeller=label_parsed) +
    theme(aspect.ratio = 1,strip.text.x = element_text(size = 9)) +
    geom_point(color='#66666670')+
    xlab("Sample relative frequency (log scale)") +
    ylab("Approximate posterior probability (log scale)")
  if(!is.null(breaks))
    p = p + scale_x_continuous(breaks = breaks)
  print(p)
}

plot.splits <- function(ds, replicate){
  df = read.ds(ds)
  df.splits = df[['splits']] %>% filter(rep==replicate)
  df.rmsd = df[['rmsd']] %>% filter(rep==replicate)
  df.kl = df[['kl']] %>% filter(rep==replicate)
  df.splits = df.splits %>%
    full_join(df.rmsd, by=c('algorithm')) %>%
    full_join(df.kl, by=c('algorithm'))
  colors = c(rgb(225,198,47,max=255),rgb(102,194,165,max=255),rgb(252,141,98,max=255))
  df.splits$color = colors[1]
  df.splits[df.splits$RMSD<0.01,]$color = colors[2]
  df.splits[df.splits$RMSD>0.05,]$color = colors[3]
  df.splits$color = factor(df.splits$color, levels=colors)

  temp = df.kl %>% full_join(df.rmsd, by=c('algorithm'))
  temp = temp[match(origin,temp$algorithm),]

  df.splits = df.splits %>%
    mutate(algorithm=factor(algorithm, levels=temp$algorithm[order(temp$RMSD,decreasing=F)]))
  label_formated= Map(.o=origin, .f=final, .kl=temp$KL, .rmsd=temp$RMSD, f)
  label_formated = label_formated[match(levels(df.splits$algorithm),names(label_formated))]
  levels(df.splits$algorithm) = unlist(label_formated)
  ggplot(df.splits, aes(x=f, y=fhat, color=color)) +
    facet_wrap( ~algorithm, scales = "free_y", ncol=4,labeller=label_parsed) +
    theme(aspect.ratio = 1,strip.text.x = element_text(size = 9), legend.position="none") +
    geom_point()+
    scale_color_manual(values=colors) +
    scale_x_continuous(breaks = c(0, 0.5, 1)) +
    labs(y="Approximate split posterior probability", x="SRF split posterior probabitlity")
}
```

\newpage

# DS1
## Approximate posterior vs. SRF posteriors of replicate 1
```{r DS1_scatterplot, fig.height=10, fig.width=20}
plot.posterior.srf('DS1', 1, c(-6, -4, -2))
```

## Split posterior probabilities approximate vs. SRF (replicate 1)
```{r DS1_split_probs, fig.height=10, fig.width=20}
plot.splits('DS1', 1)
```

\newpage

# DS2
## Approximate posterior vs. SRF posteriors of replicate 1
```{r DS2_scatterplot, fig.height=10, fig.width=20}
plot.posterior.srf('DS2', 1, c(-3, -2, -1))
```

## Split posterior probabilities approximate vs. SRF (replicate 1)
```{r DS2_split_probs, fig.height=10, fig.width=20}
plot.splits('DS2', 1)
```

\newpage

# DS3
## Approximate posterior vs. SRF posteriors of replicate 1
```{r DS3_scatterplot, fig.height=10, fig.width=20}
plot.posterior.srf('DS3', 1, c(-5, -3, -1))
```

## Split posterior probabilities approximate vs. SRF (replicate 1)
```{r DS3_split_probs, fig.height=10, fig.width=20}
plot.splits('DS3', 1)
```

\newpage

# DS4
## Approximate posterior vs. SRF posteriors of replicate 1
```{r DS4_scatterplot, fig.height=10, fig.width=20}
plot.posterior.srf('DS4', 1)#, c(-8, -6, -4, -2))
```

## Split posterior probabilities approximate vs. SRF (replicate 1)
```{r DS4_split_probs, fig.height=10, fig.width=20}
plot.splits('DS4', 1)
```

\newpage

# DS5
## Approximate posterior vs. SRF posteriors of replicate 1
```{r DS5_scatterplot, fig.height=10, fig.width=20}
plot.posterior.srf('DS5', 1, c(-9, -8, -7))
```

## Split posterior probabilities approximate vs. SRF (replicate 1)
```{r DS5_split_probs, fig.height=10, fig.width=20}
plot.splits('DS5', 1)
```