expts = list(map(lambda fp: path.basename(fp.rstrip('/')), glob('data/*/')))
print("Expts:", *expts[:10], "...")

def process_expt(expt):
    expt_results = []
    
    def extract_info(filename):
        return re.search(r'kwip/(\d\.?\d*)x-(0\.\d+)-(wip|ip).dist', filename).groups()
    
    # dict of scale: distance matrix, populated as we go
    truths = {}
    
    for distfile in glob("data/{}/kwip/*.dist".format(expt)):
        cov, scale, metric = extract_info(distfile)
        if scale not in truths:
            genome_dist_path = 'data/{ex}/all_genomes-{sc}.dist'.format(ex=expt, sc=scale)
            truths[scale] = load_sample_matrix_to_runs(genome_dist_path)
        exptmat = DistanceMatrix.read(distfile)
        rho = spearmans_rho_distmats(exptmat, truths[scale])
        expt_results.append({
            "coverage": cov,
            "scale": scale,
            "metric": metric,
            "rho": rho,
            "seed": expt,
        })
    return expt_results

#process_expt('3662')

results = []
for res in map(process_expt, expts):
    results.extend(res)
results = pd.DataFrame(results)

"""
Explanation: Calculate performance of kWIP
The next bit of python code calculates the performance of kWIP against the distance between samples calulcated from the alignments of their genomes.
This code caluclates spearman's $\rho$ between the off-diagonal elements of the triagnular distance matrices.
End of explanation
"""


%%R -i results

results$coverage = as.numeric(as.character(results$coverage))
results$scale = as.numeric(as.character(results$scale))

print(summary(results))
str(results)

"""
Explanation: Statistical analysis
Is done is R, as that's easier.
Below we see a summary and structure of the data
End of explanation
"""


%%R

ggplot(results, aes(x=coverage, y=scale)) +
    geom_point() +
    scale_x_log10() +
    scale_y_log10() +
    theme_bw()

"""
Explanation: Experiment design
Below we see the design of the experiment in terms of the two major variables.
We have a series (vertically) that, at 30x coverage, looks at the effect of genetic variation on performance. There is a second series that examines the effect of coverage at an average pairwise genetic distance of 0.001.
There are 100 replicates for each data point, performed as a separate bootstrap across the random creation of the tree and sampling of reads etc.
End of explanation
"""


%%R

dat = results %>%
        filter(scale==0.001, coverage<=30) %>%
        select(rho, metric, coverage)
    
dat$coverage = as.factor(dat$coverage)
ggplot(dat, aes(x=coverage, y=rho, fill=metric)) +
    geom_boxplot(aes(fill=metric))
    

%%R

# AND AGAIN WITHOUT SUBSETTING
dat = results %>%
        filter(scale==0.001) %>%
        select(rho, metric, coverage)
    
dat$coverage = as.factor(dat$coverage)
ggplot(dat, aes(x=coverage, y=rho, fill=metric)) +
    geom_boxplot(aes(fill=metric))
    

%%R

dat = subset(results, scale==0.001 & coverage <=15, select=-scale)
ggplot(dat, aes(x=coverage, y=rho, colour=seed, linetype=metric)) +
    geom_line()

%%R
summ = results %>%
    filter(scale==0.001, coverage <= 50) %>%
    select(-scale) %>%
    group_by(coverage, metric) %>%
    summarise(rho_av=mean(rho), rho_err=sd(rho))
 
ggplot(summ, aes(x=coverage, y=rho_av, ymin=rho_av-rho_err, ymax=rho_av+rho_err, group=metric)) +
    geom_line(aes(linetype=metric)) +
    geom_ribbon(aes(fill=metric), alpha=0.2) +
    xlab('Genome Coverage') +
    ylab(expression(paste("Spearman's ", rho, " +- SD"))) +
    scale_x_log10()+
    ggtitle("Performance of WIP & IP") +
    theme_bw()

%%R
sem <- function(x) sqrt(var(x,na.rm=TRUE)/length(na.omit(x)))
summ = results %>%
    filter(scale==0.001) %>%
    select(-scale) %>%
    group_by(coverage, metric) %>%
    summarise(rho_av=mean(rho), rho_err=sem(rho))
 
ggplot(summ, aes(x=coverage, y=rho_av, ymin=rho_av-rho_err, ymax=rho_av+rho_err, group=metric)) +
    geom_line(aes(linetype=metric)) +
    geom_ribbon(aes(fill=metric), alpha=0.2) +
    xlab('Genome Coverage') +
    ylab(expression(paste("Spearman's ", rho))) +
    scale_x_log10()+
    theme_bw()

%%R
cov_diff = results %>%
    filter(scale==0.001) %>%
    select(rho, metric, coverage, seed) %>%
    spread(metric, rho) %>%
    mutate(diff=wip-ip) %>%
    select(coverage, seed, diff) 
print(summary(cov_diff))

p = ggplot(cov_diff, aes(x=coverage, y=diff, colour=seed)) +
    geom_line() +
    scale_x_log10() +
    ggtitle("Per expt difference in performance (wip - ip)")
print(p)

summ = cov_diff %>%
    group_by(coverage) %>%
    summarise(diff_av=mean(diff), diff_sd=sd(diff))

ggplot(summ, aes(x=coverage, y=diff_av, ymin=diff_av-diff_sd, ymax=diff_av+diff_sd)) +
    geom_line() + 
    geom_ribbon(alpha=0.2) +
    xlab('Genome Coverage') +
    ylab(expression(paste("Improvment in Spearman's ", rho, " (wip - IP)"))) +
    scale_x_log10() +
    theme_bw()

%%R
var = results %>%
        filter(coverage == 30) %>%
        select(-coverage)
var$scale = as.factor(var$scale)
ggplot(var, aes(x=scale, y=rho, fill=metric)) +
    geom_boxplot() +
    xlab('Mean pairwise variation') +
    ylab(expression(paste("Spearman's ", rho))) +
    #scale_x_log10()+
    theme_bw()


%%R
summ = results %>%
        filter(coverage == 30) %>%
        select(-coverage) %>%
        group_by(scale, metric) %>%
        summarise(rho_av=mean(rho), rho_sd=sd(rho))
ggplot(summ, aes(x=scale, y=rho_av, ymin=rho_av-rho_sd, ymax=rho_av+rho_sd, group=metric)) +
    geom_line(aes(linetype=metric)) +
    geom_ribbon(aes(fill=metric), alpha=0.2) +
    xlab('Mean pairwise variation') +
    ylab(expression(paste("Spearman's ", rho))) +
    scale_x_log10()+
    theme_bw()


"""
Explanation: Effect of Coverage
Here we show the spread of data across the 100 reps as boxplots per metric and covreage level.
I note that the weighted product seems slightly more variable, particularly at higher coverage. Though the median is nearly always higher
End of explanation
"""


