#!/usr/bin/env python3

# make R test for whether R (i) is found, (ii) runs, and (iii) can generate svgs
# switch to pathlib
# --outdir-exists=move,remove,fail

import shutil
import argparse
import re
import os
from os import path
from pathlib import Path
import subprocess
import json
import itertools
import glob
import sys
import logging
import shlex

def add_suffix(p,s):
    p = Path(p)
    if not p.name:
        return p / s
    else:
        return p.with_suffix(p.suffix+s)

def print_model(model):
    main = model['main']

    assign = '='
    if main[0:1] == '~':
        main = main[1:]
        assign = '~'

    section = main + '\n'
    extracted = model['extracted']

    if extracted:
        section += '<table class="model">\n'
        for arg in extracted:
            name = arg[0]
            submodel = arg[1]
            (subassign, submodel) = print_model(submodel)
            section += f'<tr><td>{name}</td><td>{subassign}</td><td>{submodel}</td></tr>\n'
        section += '</table>'
    return (assign, section)


def print_models(name,models):
    section = '<table class="models">\n'
    i = 1
    for model in models:
        model_name =  f"{name}{i}"
        (assign, model_html) = print_model(model)
        section += f'<tr><td><a name="{model_name}" class="anchor"></a><span class="modelname">{model_name}</span></td><td>{assign}</td><td>\n'
        section += model_html
        section += '</td></tr>\n'
        i += 1
    section += '</table>\n'
    return section

def flatten(l):
    return [item for sublist in l for item in sublist]

def file_is_empty(filename):
    return path.getsize(filename) == 0

def hsv_to_rgb(H,S,V):
    # decompose color range [0,6) into a discrete color (i) and fraction (f)
    h = (H * 6);
    i = int(h);
    f = h - i;
    i = i % 6;

    p = V*(1.0 - S);
    q = V*(1.0 - (S*f));
    t = V*(1.0 - (S*(1.0-f)));

    if i==0:
        return (V,t,p)
    elif i==1:
        return (q,V,p)
    elif i==2:
        return (p,V,t)
    elif i==3:
        return (p,q,V)
    elif i==4:
        return (t,p,V)
    elif i==5:
        return (V,p,q)

def rgb_to_color(rgb):
    (r,g,b) = rgb
    return f"{r} {g} {b}"

def interpolate(i,total,start,end):
    return start + (end-start)*(i/(total-1))

def get_x3d_of_mds(point_string):
    points = []
    x = []
    y = []
    z = []
    num = dict()

    for line in point_string.strip().split('\n'):
        line.strip()
        point = line.split('\t')
        x.append(float(point[0]))
        y.append(float(point[1]))
        z.append(float(point[2]))
        group = 1
        if len(point) > 3:
            group = int(point[3])
        if group not in num:
            num[group]=0
        num[group] += 1
        points.append((float(point[0]),float(point[1]),float(point[2]),group))

    xmin = min(x)
    xmax = max(x)
    xw = abs(xmax - xmin)

    ymin = min(y)
    ymax = max(y)
    yw = abs(ymax - ymin)

    zmin = min(z)
    zmax = max(z)
    zw = abs(zmax - zmin)

    x3d = ""
    seen = dict()
    for group in num.keys():
        seen[group] = 0

    lines = ""
    spheres = ""
    lastpoint=None
    lastgroup=None
    for (x,y,z,group) in points:
        xx = (x - xmin)/xw*5.0 - 2.5
        yy = (y - ymin)/yw*5.0 - 2.5
        zz = (z - zmin)/zw*5.0 - 2.5

        size = 0.04
        scale = f"{size} {size} {size}"

        if group == 1:
            color = rgb_to_color(hsv_to_rgb(interpolate(seen[group], num[group], 11/12, 1),
                                            interpolate(seen[group], num[group], 0.4,      1.0),
                                            1)
                                )
        elif group == 2:
            color = rgb_to_color(hsv_to_rgb(interpolate(seen[group], num[group], 7/12     , 8/12),
                                            interpolate(seen[group], num[group], 0.4,      1.0),
                                            1)
                                )
        elif group == 3:
            color = rgb_to_color(hsv_to_rgb(interpolate(seen[group], num[group], 3/12, 4/12),
                                            interpolate(seen[group], num[group], 0.4,      1.0),
                                            1)
                                )
        elif group == 4:
            color = rgb_to_color(hsv_to_rgb(interpolate(seen[group], num[group], 9/12, 10/12),
                                            interpolate(seen[group], num[group], 0.4,      1.0),
                                            1)
                                )
        else:
            color = "1 0 0"

        sphere = f"<transform translation='{xx} {yy} {zz}' scale='{scale}'><shape><appearance><material specularColor='{color}' shininess='0.8' ambientIntensity='0.8' diffuseColor='{color}'/></material></appearance><sphere></sphere></shape></transform>\n"
        if lastgroup is not None and lastgroup == group:
            tmp=lastpoint+[xx,yy,zz]
            xyzs=' '.join([str(t) for t in tmp])
            line = f"<transform><shape><appearance><material ambientIntensity='0.9' diffuseColor='{color}' specularColor='{color}' transparency='0.95'/></appearance><lineset vertexCount='2'><coordinate point='{xyzs}'/></lineset></shape></transform>\n"
            lines += line
        spheres += sphere
        seen[group] += 1
        lastpoint=[xx,yy,zz]
        lastgroup=group

    x3d = f"{lines}\n{spheres}"
    return x3d

def more_recent_than(f1,f2):
    if not path.exists(f1):
        return False
    return path.getmtime(f1) > path.getmtime(f2)

def more_recent_than_all_of(f1,f2s):
    for f2 in f2s:
        if not more_recent_than(f1,f2):
            return False;
    return True

def find_exe(name,message=None):
    exe = shutil.which(name)
    if exe is None:
        if message is not None:
            print(f"Program '{name}' not found: {message}")
        else:
            print(f"Program '{name}' not found.")
    return exe

def get_n_lines(filename):
    count = 0
    with open(filename) as ifile:
        for line in ifile:
            count += 1
    return count

def get_tree_counts(filename):
    counts = []
    with open(filename, encoding='utf-8') as countfile:
        for line in countfile:
            m = re.match(r'Kept (\d+) trees', line)
            if m:
                counts.append(m.group(1))
    return counts

def get_json_from_file(filename):
    if filename is None:
        return None;
    with open(filename, encoding='utf-8') as json_file:
        return json.load(json_file)

def make_extracted(x):
    if not "main" in x:
        return {"main": x, "extracted": []}
    else:
        return x

def all_same(xs):
    if len(xs) < 2:
        return True
    for x in xs:
        if x != xs[0]:
            return False
    return True

def first_all_same(xs):
    if not all_same(xs):
        raise Exception("Not all the same!")
    return xs[0]

def get_unused_dir_name(dirbase):
    for i in itertools.count(1):
        dirname = Path(f"{dirbase}.{i}")
        if not dirname.exists():
            return dirname

def get_value_from_file(filename,attribute):
    with open(filename,'r',encoding='utf-8') as file:
        for line in file:
            m = re.match(f'{attribute} ([^ ]*)($| )', line)
            if m:
                return m.group(1)
    return None

def print_or(filenames):
    qfilenames = [f'"{filename}"' for filename in filenames]
    result = qfilenames[0]
    for i in range(1,len(filenames)-1):
        if i == len(filenames)-1:
            result += " or " + qfilenames[i]
        else:
            result += ", " + qfilenames[i]
    return result

def check_one_file_exists(filenames, required=True):
    for filename in filenames:
        if path.exists(filename):
            return Path(filename)
    if required:
        print(f"Error: I cannot find file {print_or(filenames)}")
        exit(1)
    else:
        return None

def check_file_exists(filename,required=True):
    return check_one_file_exists([filename], required)

def arg_max(xs):
    arg = None
    val = None
    for i in range(len(xs)):
        x = xs[i]
        if val is None or x > val:
            arg = i
            val = x
    return arg

def arg_min(xs):
    arg = None
    val = None
    for i in range(len(xs)):
        x = xs[i]
        if val is None or x < val:
            arg = i
            val = x
    return arg

class MCMCRun(object):
    def __init__(self, mcmc_output,**kwargs):
        self.mcmc_output = Path(mcmc_output)
        self.out_file = None
        self.trees_file = None
        self.log_file = None
        self.alignments_files = None
        self.input_files = None
        self.command = None
        self.version = None
        self.code = None

        if self.mcmc_output.is_dir():
            self.dir = self.mcmc_output
            self.prefix = None
        else:
            self.dir = self.mcmc_output.parent
            self.prefix = self.mcmc_output.name

    def personality(self):
        return None

    def has_trees(self):
        return self.trees_file

    def has_parameters(self):
        return self.log_file

    def has_alignments(self):
        return self.alignments_files

    def has_model(self):
        return False

    def has_code(self):
        if self.code:
            return True;
        else:
            return False;

    def get_parent_dir(self):
        return self.dir.resolve().parent

    def get_version(self):
        return self.version

    def get_dir(self):
        return self.dir

    def get_prefix(self):
        return self.prefix
    
    def get_trees_file(self):
        return self.trees_file

    def get_input_files(self):
        return self.input_files

    def get_alignments_files(self):
        return self.alignments_files

    def get_log_file(self):
        return self.log_file

    def get_command(self):
        return self.command

    def input_files(self):
        return self.input_files

    def n_partitions(self):
        return 1

    def get_n_sequences(self,p):
        return None

    def get_alphabets(self):
        return None

    def n_iterations(self):
        return None

#    def compute_initial_alignments(self):
#        return []

    def get_smodels(self):
        return None

    def get_smodel_indices(self):
        return None

    def get_smodel_for_partition(self,p):
        smodels = self.get_smodels()
        if smodels is None:
            return None

        indices = self.get_smodel_indices()
        if indices is None:
            return None

        index = indices[p]
        if index is None:
            return None

        return smodels[index]

    def get_imodels(self):
        return None

    def get_imodel_for_partition(self,p):
        imodels = self.get_imodels()
        if imodels is None:
            return None

        indices = self.get_imodel_indices()
        if indices is None:
            return None

        index = indices[p]
        if index is None:
            return None

        return imodels[index]

    def get_imodel_indices(self):
        return None

    def get_scale_models(self):
        return None

    def get_scale_model_indices(self):
        return None

    def get_scale_model_for_partition(self,p):
        scale_models = self.get_scale_models()
        if scale_models is None:
            return None

        indices = self.get_scale_model_indices()
        if indices is None:
            return None

        index = indices[p]
        if index is None:
            return None

        return scale_models[index]

    def get_topology_prior(self):
        return None

    def get_branch_length_prior(self):
        return None

    def get_tree_prior(self):
        return None

    def get_code(self):
        return []

class LogFileRun(MCMCRun):
    def __init__(self,mcmc_output,**kwargs):
        super().__init__(mcmc_output, **kwargs)
        self.log_file = mcmc_output
        logging.debug(f"LogFileRun: {mcmc_output}")

    def n_iterations(self):
        n = get_n_lines(self.get_log_file())-1
        if n <= 0:
            print(f"Error: Log file '{self.get_log_file()}' has no samples!")
            exit(1)
        return n


class TreeFileRun(MCMCRun):
    def __init__(self, mcmc_output, **kwargs):
        super().__init__(mcmc_output, **kwargs)
        self.trees_file = mcmc_output
        logging.debug(f"TreeFileRun: {mcmc_output}")

    def n_iterations(self):
        n = get_n_lines(self.get_trees_file())-1
        if n <= 0:
            print(f"Error: Tree file '{self.get_trees_file()}' has no samples!")
            exit(1)
        return n

def translate_none(indices):
    return [None if index == '--' or index == '-1' else int(index) for index in indices]

def add_range(filename,r):
    name = filename
    if r:
        name = filename + ":" + r
    return name

class BAliPhyRun(TreeFileRun,LogFileRun):

    def __init__(self,mcmc_output,**kwargs):
        super().__init__(mcmc_output,**kwargs)
        self.prefix = 'C1'
        # What I should actually do is to create the run_info from C1.out if C1.run.json does not exist.
        self.run_file = check_file_exists( self.get_dir() / 'C1.run.json', required=False)
        if self.run_file:
            self.run_info = get_json_from_file(self.run_file)
            self.out_file = None
        else:
            self.out_file = check_file_exists( self.get_dir() / 'C1.out', required=False)
            if not self.out_file:
                print("Error: cannot find either 'C1.out' or 'C1.run.json")
                exit(1)
            self.run_info = translate_to_dict()

        self.log_file = check_one_file_exists( [self.get_dir() / 'C1.log.json',
                                                self.get_dir() / 'C1.log'] )

        self.input_files = [ add_range(p['filename'],p['range']) for p in self.run_info['partitions']]
        self.trees_file = check_file_exists( self.get_dir() / 'C1.trees')
        self.alignments_files = self.get_alignment_files()

        self.smodel_indices = [p["smodel"] for p in self.run_info["partitions"]]
        self.imodel_indices = [p["imodel"] for p in self.run_info["partitions"]]
        self.scale_model_indices = [p["scale"] for p in self.run_info["partitions"]]
        self.alphabets = [p["alphabet"] for p in self.run_info["partitions"]]

        self.smodels = self.find_models("subst models", "smodels")
        self.imodels = self.find_models("indel models", "imodels")
        self.scale_models = self.find_models("scale model", "scales")

        self.command = ' '.join([shlex.quote(word) for word in self.run_info["command"]])
        self.version = self.run_info["program"]["version"]
        if (isinstance(self.run_info["tree"],str)):
            self.topology_prior = None
            self.branch_prior = None
            self.tree_prior = self.run_info["tree"]
        else:
            self.topology_prior = self.run_info["tree"]["topology"]
            self.branch_prior = self.run_info["tree"]["lengths"]
            self.tree_prior = None

        self.code = []
        for file_path in os.listdir( self.get_dir() ):
            if Path(file_path).suffix == '.hs':
                self.code.append( self.get_dir() / file_path )

        logging.debug(f"BAliPhyRun: {mcmc_output}")

    def translate_to_dict():
        run_info = dict()
        run_info["command"] = self.find_header_attribute("command")
        run_info["program"]["version"] = self.find_header_attribute("VERSION").split(r'\s+')[0]
        run_info["tree"]["topology"] = self.find_in_header('T:topology ')
        run_info["tree"]["lengths"] = self.find_in_header('T:lengths ')

        smodel_indices = translate_none(self.find_partition_values('smodel-index'))
        imodel_indices = translate_none(self.find_partition_values('imodel-index'))
        scale_indices = translate_none(self.find_partition_values('scale-index'))
        alphabets = self.find_partition_values('alphabet')

        input_filenames = []
        with open(outfile,'r',encoding='utf-8') as outf:
            for line in outf:
                m = re.match(r'data(.+) = (.+)', line)
                if m:
                    input_filenames.append(m.group(2))

        for i in range(0,n_partitions()-1):
            p["smodel"] = smodel_indices[i]
            p["imodel"] = imodel_indices[i]
            p["scale"] = scale_indices[i]
            p["alphabet"] = alphabets[i]
            p["filename"] = input_filenames[i]
            p["range"] = ""
            if ':' in input_filenames[i]:
                p["filename"] = input_filenames.split(':',1)[0]
                p["range"] = input_filenames.split(':',1)[1]
            run_info["partitions"] = p

        return run_info

    def find_models(self, name1, name2):
        return [make_extracted(model) for model in self.run_info[name2]]

    def get_alphabets(self):
        return self.alphabets

    def has_model(self):
        return True

    def get_smodels(self):
        return self.smodels

    def get_smodel_indices(self):
        return self.smodel_indices

    def get_imodels(self):
        return self.imodels

    def get_imodel_indices(self):
        return self.imodel_indices

    def get_scale_models(self):
        return self.scale_models

    def get_scale_model_indices(self):
        return self.scale_model_indices

    def get_topology_prior(self):
        return self.topology_prior

    def get_branch_length_prior(self):
        return self.branch_prior

    def get_tree_prior(self):
        return self.tree_prior

    def get_n_sequences(self,p):
        features = self.get_alignment_info( self.dir / "C1.P{p+1}.initial.fasta")
        return features["n_sequences"]

    def get_alignment_files(self):
        filenames = []
        for p in range(1,self.n_partitions()+1):
            filename = self.get_dir() / f'C1.P{p}.fastas'
            if not path.exists(filename):
                filename = None
            filenames.append(filename)
        return filenames

    def find_header_attribute(self, attribute):
        return self.find_in_header(attribute + ': ')

    def find_in_header(self, key):
        with open(self.out_file,'r',encoding='utf-8') as file:
            reg = key + '(.*)$'
            for line in file:
                m = re.match(reg,line)
                if m:
                    return m.group(1)
                if line.startswith("iterations"):
                    break
        return None

    def n_partitions(self):
        return len(self.get_input_files())

    def find_partition_values(self,prefix):
        indices = []
        with open(self.out_file,encoding='utf-8') as file:
            regexp = prefix+"(.+) = (.+)"
            for line in file:
                m = re.match(regexp, line)
                if m:
                    indices.append(m.group(2))
                if line.startswith("iterations"):
                    break
        return indices

    def get_code(self):
        return self.code


class TreeLogFileRun(TreeFileRun,LogFileRun):
    def __init__(self,mcmc_output, **kwargs):
        super().__init__(mcmc_output, **kwargs)
        logging.debug(f"TreeLogFile Run: {mcmc_output}")

# Since PhyloBayes doesn't make new run files for each dir, maybe we need to file the different prefixes, and pass them to PhyloBayes?
# A prefix could be like dir/prefix (e.g. dir/prefix.trace) or prefix (e.g. prefix.trace)
class PhyloBayesRun(MCMCRun):
    def __init__(self,mcmc_output, **kwargs):
        super().__init__(mcmc_output, **kwargs)
        self.log_file = check_file_exists(f"{prefix}.trace")
        self.out_file = check_file_exists(f"{prefix}.log")
        logging.debug(f"PhyloBayesRun: {mcmc_output}")

    def personality(self):
        return "phylobayes"

    def get_input_file_names_for_outfile(self,outfile):
        return get_value_from_file(outfile,"data file :")

    def get_n_sequences(self, p):
        return get_value_from_file(self.out_file, "number of taxa:")

class BEASTRun(MCMCRun):
    def __init__(self,mcmc_output, **kwargs):
        super().__init__(mcmc_output, **kwargs)
        for tree_file in glob.glob( self.get_dir() / '*.trees'):
            check_file_exists(tree_file)
            self.prefix = get_file_prefix(tree_file)
            self.log_file = check_file_exists(f"{prefix}.log")

def get_bali_phy_runs(mcmc_args):
    runs = []
    good = []
    bad = []
    for mcmc_dir in mcmc_args:
        if mcmc_dir.exists() and \
           (Path( mcmc_dir , 'C1.out').exists() or Path( mcmc_dir , 'C1.run.json').exists()) and \
           (Path( mcmc_dir , 'C1.log' ).exists() or Path (mcmc_dir, 'C1.log.json')):
            runs.append(BAliPhyRun(mcmc_dir))
            good.append(mcmc_dir)
        else:
            bad.append(mcmc_dir)
    if len(bad) > 0 and len(good) > 0:
        raise Exception(f"{good} appear to be BAli-Phy 3/4 runs, but {bad} do not!")
    return runs

def get_phylobayes_runs(mcmc_args):
    runs = []
    good = []
    bad = []
    for run_name in mcmc_args:
        if path.exists(add_suffix(run_name,'.treelist')):
            runs.append(PhyloBayesRun(run_name))
            good.append(run_name)
        else:
            bad.append(run_name)
    if len(bad) > 0 and len(good) > 0:
        raise Exception(f"{good} appear to be phylobayes runs, but {bad} do not!")
    return runs

def get_tree_log_runs(mcmc_args):
    runs = []
    good = []
    bad = []
    for run_name in mcmc_args:
        if path.exists(str(run_name)+".trees") and path.exists(str(run_name)+".log"):
            runs.append(TreeLogFileRun(run_name))
            good.append(run_name)
        else:
            bad.append(run_name)
    if len(bad) > 0 and len(good) > 0:
        raise Exception(f"{good} appear to be generic tree/parameters runs, but {bad} do not!")
    return runs

def get_tree_runs(mcmc_args):
    runs = []
    good = []
    bad = []
    for run_name in mcmc_args:
        if run_name.exists() and str(run_name).endswith(".trees"):
            runs.append(TreeFileRun(run_name))
            good.append(run_name)
        elif add_suffix(run_name,".trees").exists() and not add_suffix(run_name,'.log').exists():
            runs.append(TreeFileRun(add_suffix(run_name,".trees")))
            good.append(run_name)
        else:
            bad.append(run_name)
    if len(bad) > 0 and len(good) > 0:
        raise Exception(f"{good} appear to be generic tree-only runs, but {bad} do not!")
    return runs

def get_log_runs(mcmc_args):
    runs = []
    good = []
    bad = []
    for run_name in mcmc_args:
        if path.exists(run_name) and str(run_name).endswith(".log"):
            runs.append(LogFileRun(run_name))
            good.append(run_name)
        elif not add_suffix(run_name,'.trees').exists() and add_suffix(run_name,'.log').exists():
            runs.append(LogFileRun(add_suffix(run_name,".log")))
            good.append(run_name)
        else:
            bad.append(run_name)
    if len(bad) > 0 and len(good) > 0:
        raise Exception(f"{good} appear to be generic parameter-only runs, but {bad} do not!")
    return runs

def construct_runs(mcmc_outputs):
    for mcmc_output in mcmc_outputs:
        if not mcmc_output.exists() and not glob.glob(str(mcmc_output)+('*')):
            print(f"Error: No file or directory '{mcmc_output}' or '{str(mcmc_output)+'*'}'")
            exit(1)

    run_groups = []

    run_groups.append( get_bali_phy_runs(mcmc_outputs) )

    run_groups.append( get_phylobayes_runs(mcmc_outputs) )

    run_groups.append( get_tree_log_runs(mcmc_outputs) )

    run_groups.append( get_tree_runs(mcmc_outputs) )

    run_groups.append( get_log_runs(mcmc_outputs) )

    for run_group in run_groups:
        if run_group:
            return run_group

    print(f"Error: MCMC runs '{mcmc_outputs}' are not BAli-Phy, BEAST, or PhyloBayes run.", file=sys.stderr)
    exit(3)
        


class Analysis(object):

    def __init__(self,args,mcmc_outputs,outdir):
        self.mcmc_runs = construct_runs(mcmc_outputs)

        self.trees_consensus_exe = find_exe('trees-consensus', message="See the main for adding the bali-phy programs to your PATH.")
        if self.trees_consensus_exe is None:
            exit(1)

        self.draw_tree_exe = find_exe('draw-tree', message="Tree pictures will not be generated.\n")
        # FIXME - maybe switch from gnuplot to R?
        self.gnuplot_exe = find_exe('gnuplot', message='Some graphs will not be generated.\n')
        self.R_exe = find_exe('R', message='Some mixing graphs will not be generated.\n')
        self.fourmolu_exe = shutil.which("fourmolu")
        self.ormolu_exe = shutil.which("ormolu")
        self.hindent_exe = shutil.which("hindent")
        if not self.fourmolu_exe and not self.ormolu_exe and not self.hindent_exe:
            print("Can't find fourmolu, ormolu or hindent -- not formatting Haskell code.")

        self.subpartitions = args.subpartitions
        self.subsample = args.subsample
        self.prune = args.prune
        self.burnin = args.skip
        self.until = args.until
        self.verbose = args.verbose
        self.speed = 1
        self.outdir = Path(outdir)

        # map from alignment name to alphabet name
        self.alignments = list()

        prefix=path.dirname(path.dirname(self.trees_consensus_exe))
        self.libexecdir = Path(prefix,"lib","bali-phy","libexec")
        if not path.exists(self.libexecdir):
            print("Can't find bali-phy libexec path '{self.libexecdir}'")

        for i in range(len(self.mcmc_runs)):
            if self.mcmc_runs[i].get_command() != self.mcmc_runs[0].get_command():
                print(f"WARNING: Commands differ!\n   {self.mcmc_runs[0].get_command()}\n   {self.mcmc_runs[i].get_command()}\n")

        self.get_input_files()

        self.tree_consensus_levels = [0.5,0.66,0.8,0.9,0.95,0.99,1.0]
        self.alignment_consensus_values = [0.1,0.25,0.5,0.75]

        self.determine_burnin()
        self.initialize_results_directory()

        self.log_shell_cmds = (self.outdir / "commands.log").open('a+',encoding='utf-8')
        print("-----------------------------------------------------------",file=self.log_shell_cmds)

        if self.has_parameters():
            self.summarize_numerical_parameters()
        if self.has_trees():
            self.summarize_topology_distribution()
            self.draw_trees()
            self.compute_tree_mixing_diagnostics()
            self.compute_srq_plots()
            self.compute_tree_MDS()
        if self.has_alignments():
            self.compute_initial_alignments()
            self.compute_wpd_alignments()
            self.compute_ancestral_states()
            self.draw_alignments()
            self.compute_and_draw_AU_plots()
        if self.has_code():
            self.copy_code()
        self.print_index_html(self.outdir / "index.html")

    def n_chains(self):
        return(len(self.mcmc_runs))

    def has_trees(self):
        return self.mcmc_runs[0].has_trees()

    def has_parameters(self):
        return self.mcmc_runs[0].has_parameters()

    def has_alignments(self):
        return self.mcmc_runs[0].has_alignments()

    def has_model(self):
        return self.mcmc_runs[0].has_model()

    def has_code(self):
        return self.mcmc_runs[0].has_code()

    def get_code(self):
        # We should check that all the runs have the SAME code!
        return self.mcmc_runs[0].get_code()

    def run(self,p):
        return self.mcmc_runs[p]

    def get_libexec_script(self,file):
        filename = self.libexecdir / file
        if not path.exists(filename):
            print(f"Can't find script '{file}'!")
        return filename

    def Rexec(self,script,args=[]):
        if not self.R_exe:
            return
        return self.exec_show([self.R_exe,'--slave','--vanilla','--args']+args,infile=script)

    def run_gnuplot(self,script):
        if not self.gnuplot_exe:
            return
        subprocess.Popen(self.gnuplot_exe,stdin=subprocess.PIPE).communicate(script.encode('utf-8'))

    def get_input_files(self):
        return first_all_same([run.get_input_files() for run in self.mcmc_runs])

    def n_partitions(self):
        return first_all_same([run.n_partitions() for run in self.mcmc_runs])

    def get_n_sequences(self, p):
        return first_all_same([run.get_n_sequences(p) for run in self.mcmc_runs])

    def get_imodels(self):
        return first_all_same([run.get_imodels() for run in self.mcmc_runs])

    def get_imodel_indices(self):
        return first_all_same([run.get_imodel_indices() for run in self.mcmc_runs])

    def get_imodel_for_partition(self,p):
        return first_all_same([run.get_imodel_for_partition(p) for run in self.mcmc_runs])

    def get_smodels(self):
        return first_all_same([run.get_smodels() for run in self.mcmc_runs])

    def get_smodel_indices(self):
        return first_all_same([run.get_smodel_indices() for run in self.mcmc_runs])

    def get_smodel_for_partition(self,p):
        return first_all_same([run.get_smodel_for_partition(p) for run in self.mcmc_runs])

    def get_scale_models(self):
        return first_all_same([run.get_scale_models() for run in self.mcmc_runs])

    def get_scale_model_indices(self):
        return first_all_same([run.get_scale_model_indices() for run in self.mcmc_runs])

    def get_scale_model_for_partition(self,p):
        return first_all_same([run.get_scale_model_for_partition(p) for run in self.mcmc_runs])

    def get_tree_prior(self):
        return first_all_same([run.get_tree_prior() for run in self.mcmc_runs])

    def get_topology_prior(self):
        return first_all_same([run.get_topology_prior() for run in self.mcmc_runs])

    def get_branch_length_prior(self):
        return first_all_same([run.get_branch_length_prior() for run in self.mcmc_runs])

    def find_models(self, name1, name2):
        return first_all_same([run.find_models(name1,name2) for run in self.mcmc_runs])

    def get_alphabets(self):
        return first_all_same([run.get_alphabets() for run in self.mcmc_runs])

    def get_alignments_files(self):
        return [run.get_alignments_files() for run in self.mcmc_runs]

    def get_log_files(self):
        return [run.get_log_file() for run in self.mcmc_runs]

    def get_trees_files(self):
        return [run.get_trees_file() for run in self.mcmc_runs]

    def get_alignments_for_partition(self,p):
        return [run.get_alignments_files()[p] for run in self.mcmc_runs]

    def exec_show_result(self,cmd,**kwargs):
        showcmd = ' '.join([f"'{word}'" for word in cmd])

        subargs = dict()
        if "infile" in kwargs:
            infile = kwargs["infile"]
            subargs["stdin"] = open(infile,encoding='utf-8')
            showcmd += f" < '{infile}'"
        elif "stdin" in kwargs:
            subargs["stdin"] = kwargs["stdin"]

        if "outfile" in kwargs:
            outfile = kwargs["outfile"]
            subargs["stdout"] = open(outfile,'w+',encoding='utf-8')
            showcmd += f" > '{outfile}'"
        elif "stdout" in kwargs:
            subargs["stdout"] = kwargs["stdout"]
        else:
            subargs["stdout"] = subprocess.PIPE

        if "errfile" in kwargs:
            errfile = kwargs["errfile"]
            subargs["stderr"] = open(errfile,'w+',encoding='utf-8')
            showcmd += f" 2> '{errfile}'"
        elif "stderr" in kwargs:
            subargs["stderr"] = kwargs["stderr"]
        else:
            subargs["stderr"] = subprocess.PIPE

        if "cwd" in kwargs:
            subargs["cwd"] = kwargs["cwd"]

        print(showcmd,file=self.log_shell_cmds)
        result = subprocess.run(cmd,**subargs)
        if result.returncode != 0:
            print(f"command: {showcmd}",file=sys.stderr)
            if "outfile" in kwargs and path.exists(kwargs["outfile"]):
                os.remove(kwargs["outfile"])
            if "errfile" in kwargs and path.exists(kwargs["errfile"]):
                os.remove(kwargs["errfile"])
            if "handler" in kwargs:
                handler = kwargs["handler"]
                handler(result.returncode)
        elif self.verbose:
            print(f"\n\t{showcmd}\n")
        return result

    def exec_show(self,cmd,**kwargs):
        result = self.exec_show_result(cmd,**kwargs)

        out_message = None
        if "stdout" not in kwargs and "outfile" not in kwargs:
            out_message = result.stdout.decode('utf-8')

        err_message = None
        if "stderr" not in kwargs and "errfile" not in kwargs:
            err_message = result.stderr.decode('utf-8')

        # Always record error messages in the log file.
        if err_message:
            print(f"  err: {err_message}", file=self.log_shell_cmds)
        if self.verbose and out_message:
            print(f"  out: {out_message}", file=self.log_shell_cmds)

        code = result.returncode
        if code != 0:
            print(f" exit: {code}", file=sys.stderr)

            print(f" exit: {code}", file=self.log_shell_cmds)
            if out_message:
                print(f"  out: {out_message}", file=self.log_shell_cmds)
                print(f"  out: {out_message}", file=sys.stdout)
            if err_message:
                print(f"  err: {err_message}", file=self.log_shell_cmds)

            exit(code)
        return out_message

    def load_analysis_properties(self):
        filename = self.outdir / "properties.json"
        if not filename.exists():
            return dict()
        return get_json_from_file(filename)

    def save_analysis_properties(self, properties):
        filename = self.outdir / "properties.json"
        with filename.open('w') as outfile:
            json.dump(properties, outfile, indent=2)

    def read_analysis_property(self,key):
        return self.load_analysis_properties()[key]

    def check_analysis_property(self,key,value):
        properties = self.load_analysis_properties()
        return key in properties and value == properties[key]

    def write_analysis_property(self,key,value):
        properties = self.load_analysis_properties()
        properties[key] = value
        self.save_analysis_properties(properties)
        return properties
        
    def write_input_file_names(input_file_names):
        filename = self.outdir / "input_files"
        with filename.open("w",encoding='utf-8') as out:
            print(input_file_names.join('\n'),file=out)

    def initialize_results_directory(self):
        reuse = self.check_analysis_property("burnin", self.burnin)
        reuse = reuse and self.check_analysis_property("subsample", self.subsample)
        reuse = reuse and self.check_analysis_property("until", self.until)
        reuse = reuse and self.check_analysis_property("prune", self.prune)
        reuse = reuse and self.check_analysis_property("alignment_file_names",[str(file) for file in self.get_alignments_files()])
        reuse = reuse and self.check_analysis_property("input_files", self.get_input_files())

        if not reuse and self.outdir.exists():
            new_dir_name = get_unused_dir_name(self.outdir)
            if self.outdir.exists():
                print(f"Renaming '{self.outdir}' to '{new_dir_name}'\n")
                os.rename(self.outdir,new_dir_name)

        if not self.outdir.exists():
            os.mkdir(self.outdir)
            os.mkdir(self.outdir / "Work")
            print(f"Creating new directory '{self.outdir}' for summary files.")

        self.write_analysis_property("burnin", self.burnin)
        self.write_analysis_property("subsample", self.subsample)
        self.write_analysis_property("until", self.until)
        self.write_analysis_property("prune", self.prune)
        self.write_analysis_property("input_files", self.get_input_files())
        self.write_analysis_property("alignment_file_names", [str(file) for file in self.get_alignments_files()])

    def determine_burnin(self):
        if self.burnin is not None:
            for run in self.mcmc_runs:
                if self.burnin > run.n_iterations():
                    print(f"MCMC run {run.mcmc_output} has only {run.n_iterations()} iterations.")
                    print(f"Error!  The burnin (specified as {self.burnin}) cannot be higher than this.")
                    exit(1)
            return

        iterations = [run.n_iterations() for run in self.mcmc_runs]

        max_iterations = max(iterations)
        max_i = arg_max(iterations)

        min_iterations = min(iterations)
        min_i = arg_min(iterations)

        if max_iterations > 3*min_iterations:
            print( "The variation in run length between MCMC chains is too great.\n")
            print(f"  Chain #{max_i} has {max_iterations} iterations.")
            print(f"  Chain #{min_i} has {min_iterations} iterations.")
            print( "Not going to guess a burnin: please specify one.")
            exit(1)

        self.burnin = int(0.1*min_iterations)

    def summarize_numerical_parameters(self):
        log_files = self.get_log_files()
        if log_files is None or None in log_files:
            return

        print("Summarizing distribution of numerical parameters: ",end='')
        if not more_recent_than_all_of(self.outdir / "Report",log_files):
            cmd = ["statreport"] + log_files
            if self.subsample is not None and self.subsample != 1:
                cmd.append(f"--subsample={self.subsample}")
            if self.burnin is not None:
                cmd.append(f"--skip={self.burnin}")
            if self.until is not None:
                cmd.append(f"--until={self.until}")
            self.exec_show(cmd,outfile=self.outdir / "Report")
        print("done.")

        print("Analyzing scalar variables: ",end='',flush=True)
        self.median = dict()
        self.CI_low = dict()
        self.CI_high = dict()

        self.ACT = dict()
        self.ESS = dict()
        self.Burnin = dict()
        self.PSRF_CI80 = dict()
        self.PSRF_RCF = dict()

        with open(self.outdir / "Report",encoding='utf-8') as report:
            lines = report.readlines()
            i = 0
            while i < len(lines):
                line = lines[i].strip()

                if not line:
                    i += 1
                    continue
                m = re.search(r'([^\s]+) ~ ([^\s]+)\s+\((.+),(.+)\)',line)
                if m:
                    var = m.group(1)
                    self.median[var] = m.group(2)
                    self.CI_low[var] = m.group(3)
                    self.CI_high[var] = m.group(4)

                    i += 1
                    m = re.search(r't @ (.+)\s+Ne = ([^ ]+)\s+burnin = (Not Converged!|[^ ]+)', lines[i])
                    if m:
                        self.ACT[var] = float(m.group(1))
                        self.ESS[var] = float(m.group(2))
                        self.Burnin[var] = m.group(3)
                        if self.Burnin[var] != "Not Converged!":
                            self.Burnin[var] = float(self.Burnin[var])

                    i += 1
                    m = re.search(r'PSRF-80%CI = ([^ ]+)\s+PSRF-RCF = ([^ ]+)',lines[i])
                    if m:
                        self.PSRF_CI80[var] = float(m.group(1))
                        self.PSRF_RCF[var] = float(m.group(2))
                    i += 1
                    continue

                m = re.search(r'\s+(.+) = (.+)', line)
                if m:
                    var = m.group(1)
                    self.median[var] = m.group(2)
                i += 1

        self.min_ESS = min(self.ESS.values())
        print("done.")


    def summarize_topology_distribution(self):
        assert(self.has_trees())
        print("\nSummarizing topology distribution: ",end='')
        logging.debug(f"tree files = {self.get_trees_files()}")
        cmd = ['trees-consensus',f'--map-tree={self.outdir / "MAP.PP.tree"}',f'--greedy-consensus={self.outdir / "greedy.PP.tree"}','--report',self.outdir / 'consensus']+self.get_trees_files()
        cmd.append(f"--support-levels={self.outdir / 'c-levels.plot'}")
        if self.subpartitions:
            cmd.append("--sub-partitions")
            cmd.append(f"--extended-support-levels={self.outdir / 'extended-c-levels.plot'}")
        if self.prune is not None:
            cmd.append(f"--ignore={self.prune}")
        if self.subsample is not None and self.subsample != 1:
            cmd.append(f"--subsample={self.subsample}")
        if self.burnin is not None:
            cmd.append(f"--skip={self.burnin}")
        if self.until is not None:
            cmd.append(f"--until={self.until}")

        self.trees = [("greedy","greedy"),("MAP","MAP")]
        tree_names=[self.outdir / 'greedy.PP.tree',self.outdir / 'MAP.PP.tree']
        consensus_trees=[]
        for level in self.tree_consensus_levels:
            value = int(level*100)
            name = f"c{value}"
            filename = self.outdir / f"{name}.PP.tree"
            tree_names.append(filename)
            self.trees.append((name,f'{value}% consensus'))
            consensus_trees.append(f"{level}:{filename}")
        cmd.append("--consensus={}".format(','.join(consensus_trees)))

        extended_consensus_trees=[]
        for level in self.tree_consensus_levels:
            filename = self.outdir / f"c{int(level*100)}.mtree"
            extended_consensus_trees.append(f"{level}:{filename}")
        if self.subpartitions:
            cmd.append("--extended-consensus={}".format(','.join(extended_consensus_trees)))

        extended_consensus_L=[]
        for level in self.tree_consensus_levels:
            filename = self.outdir / f"c{int(level*100)}.mlengths"
            extended_consensus_L.append(f"{level}:{filename}")
        if self.subpartitions:
            cmd.append("--extended-consensus-L={}".format(','.join(extended_consensus_L)))

        if not more_recent_than_all_of(self.outdir / "consensus", self.get_trees_files()):
            self.exec_show(cmd)
        for tree in tree_names:
            if not path.exists(tree) or file_is_empty(tree):
                raise Exception(f"Tree '{tree}' not found!")
            assert(str(tree).endswith('.PP.tree'))

            # Replace suffix ".PP.tree" with just ".tree"
            tree2 = tree
            while tree2.suffix:
                tree2 = tree2.with_suffix('')
            tree2 = tree2.with_suffix('.tree')

            if not more_recent_than(tree2,tree):
                self.exec_show(['tree-tool',tree,'--strip-internal-names','--name-all-nodes'],outfile=tree2)
        print(" done.")

    def get_alignment_info(self,filename):
        output = self.exec_show(['alignment-info',filename])
        features = dict()
        for line in output.split('\n'):
            m = re.match(r'Alphabet: (.*)$', line)
            if m:
                features["alphabet"] = m.group(1)
                continue

            m = re.match(r'Alignment: (.+) columns of (.+) sequences', line)
            if m:
                features["length"] = m.group(1)
                features["n_sequences"] = m.group(2)
                continue

            m = re.search(r'sequence lengths: ([^ ]+)-([^ ]+) ', line)
            if m:
                features["min_length"] = m.group(1)
                features["max_length"] = m.group(2)

            m = re.search(r' const.: ([^ ]+) \(([^ ]+)\%\)', line)
            if m:
                features["n_const"] = m.group(1)
                features["p_const"] = m.group(2)

            m = re.search(r'inform.: ([^ ]+) \(([^ ]+)\%\)', line)
            if m:
                features["n_inform"] = m.group(1)
                features["p_inform"] = m.group(1)

            m = re.search(r' ([^ ]+)% minimum sequence identity', line)
            if m:
                features["min_p_identity"] = m.group(1)

        return features

    def draw_trees(self):
        if not self.draw_tree_exe:
            return
        print("Drawing trees: ",end='')
        for level in self.tree_consensus_levels:
            value = int(level*100)
            tree = f"c{value}"

            filename1 = self.outdir / f"{tree}.tree"
            filename2 = self.outdir / f"{tree}.mtree"

            if self.speed < 2 and path.exists(filename2):
                cmd = ['draw-tree', self.outdir / f'{tree}.mlengths', '--out', self.outdir / f'{tree}-mctree', '--draw-clouds=only']
                if not more_recent_than(self.outdir / f"{tree}-mctree.svg",filename2):
                    self.exec_show(cmd + ['--output=svg'])
                if not more_recent_than(self.outdir / f"{tree}-mctree.pdf",filename2):
                    self.exec_show(cmd + ['--output=pdf'])

            if not more_recent_than(self.outdir / f"{tree}-tree.pdf", filename1):
                cmd = ['draw-tree',tree+".tree",'--layout=equal-daylight']
                self.exec_show(cmd,cwd=self.outdir)
            if not more_recent_than(self.outdir / f"{tree}-tree.svg", filename1):
                cmd = ['draw-tree',tree+".tree",'--layout=equal-daylight','--output=svg']
                self.exec_show(cmd,cwd=self.outdir)
            print(tree+' ',end='',flush=True)

        for tree in ['greedy','MAP']:
            filename = self.outdir / f"{tree}.tree"
            if path.exists(filename):
                outname = self.outdir / f"{tree}-tree"
                if not more_recent_than(outname.with_suffix(".pdf"), filename):
                    self.exec_show(['draw-tree', tree+'.tree','--layout=equal-daylight'],cwd=self.outdir)
                if not more_recent_than(outname.with_suffix(".svg"), filename):
                    self.exec_show(['draw-tree', tree+'.tree','--layout=equal-daylight','--output=svg'],cwd=self.outdir)
            print(f'{tree} ', end='')

        print(". done.")

    def compute_tree_mixing_diagnostics(self):
        print("\nGenerate mixing diagnostics for topologies ...",end='')

        if not more_recent_than(self.outdir / "partitions",self.outdir / "consensus"):
            self.exec_show(['pickout','--no-header','--large','pi'],
                           infile=self.outdir / "consensus",
                           outfile=self.outdir / "partitions")

        # This just adds blank lines between the partitions.
        if not more_recent_than(self.outdir / "partitions.pred",self.outdir / "partitions"):
            with open(self.outdir / "partitions",encoding='utf-8') as infile:
                with open(self.outdir / "partitions.pred","w+",encoding='utf-8') as outfile:
                    for line in infile:
                        print(line,file=outfile)

        if not more_recent_than_all_of(self.outdir / "partitions.bs",self.get_trees_files()):
            cmd = ['trees-bootstrap',
                   '--pred',self.outdir / 'partitions.pred',
                   '--LOD-table',self.outdir / 'LOD-table',
                   '--pseudocount=1']
            cmd += self.get_trees_files()
            if self.prune is not None:
                cmd.append(f"--ignore={self.prune}")
            if self.subsample is not None and self.subsample != 1:
                cmd.append(f"--subsample={self.subsample}")
            if self.burnin is not None:
                cmd.append(f"--skip={self.burnin}")
            if self.until is not None:
                cmd.append(f"--until={self.until}")
            self.exec_show(cmd,outfile=self.outdir / "partitions.bs")

        if self.n_chains() < 2:
            print(" done.")
            return

        if not more_recent_than_all_of(self.outdir / "convergence-PP.pdf",self.get_trees_files()):
            script = self.get_libexec_script("compare-runs.R")
            self.Rexec(script,[self.outdir / "LOD-table",self.outdir / "convergence-PP.pdf"])

        if (not more_recent_than_all_of(self.outdir / "convergence1-PP.svg",self.get_trees_files()) or
            not more_recent_than_all_of(self.outdir / "convergence2-PP.svg",self.get_trees_files())):
            script = self.get_libexec_script("compare-runs2.R")
            self.Rexec(script,[self.outdir / "LOD-table",self.outdir / "convergence1-PP.svg",self.outdir / "convergence2-PP.svg"])
        print(" done.")


    def compute_srq_plots(self):
        self.srq = []
        print("Generate SRQ plot for partitions: ",end='')

        if not more_recent_than_all_of(self.outdir / "partitions.SRQ",self.get_trees_files()):
            cmd = ['trees-to-SRQ',self.outdir / 'partitions.pred','--max-points=1000']
            if self.subsample is not None and self.subsample != 1:
                cmd.append(f"--subsample={self.subsample}")
            if self.burnin is not None:
                cmd.append(f"--skip={self.burnin}")
            if self.until is not None:
                cmd.append(f"--until={self.until}")
            cmd += self.get_trees_files()
            self.exec_show(cmd,outfile=self.outdir / "partitions.SRQ")
        print("done.");
        self.srq.append("partitions")

        print("Generate SRQ plot for c50 tree: ", end='')
        if not more_recent_than_all_of(self.outdir / "c50.SRQ", self.get_trees_files()):
            cmd = ['trees-to-SRQ',self.outdir / 'c50.tree','--max-points=1000']
            if self.subsample is not None and self.subsample != 1:
                cmd.append(f"--subsample={self.subsample}")
            if self.burnin is not None:
                cmd.append(f"--skip={self.burnin}")
            if self.until is not None:
                cmd.append(f"--until={self.until}")
            cmd += self.get_trees_files()
            self.exec_show(cmd,outfile=self.outdir / "c50.SRQ")
        print("done.");
        self.srq.append("c50")

        for srq in self.srq:
            cmd = ['gnuplot']
            self.run_gnuplot(f"""\
set terminal png size 300,300
set output "{self.outdir / f'{srq}.SRQ.png'}"
set key right bottom
set xlabel "Regenerations (fraction)"
set ylabel "Time (fraction)"
set title "Scaled Regeneration Quantile (SRQ) plot: $srq"
plot "{self.outdir / f'{srq}.SRQ'}" title "{srq}" with linespoints lw 1 lt 1, x title "Goal" lw 1 lt 3
""")
        if self.subpartitions:
            self.run_gnuplot(f"""\
set terminal svg
set output "{self.outdir / 'c-levels.svg'}"
set xlabel "Log10 posterior Odds (LOD)"
set ylabel "Supported Splits"
set style data lines
plot [0:][0:] '{self.outdir / "c-levels.plot"}' title 'Full Splits','{self.outdir / "extended-c-levels.plot"}' title 'Partial Splits'""")
        else:
            self.run_gnuplot(f"""\
set terminal svg
set output "{self.outdir / 'c-levels.svg'}"
set xlabel "Log10 posterior Odds (LOD)"
set ylabel "Supported Splits"
plot [0:][0:] '{self.outdir / "c-levels.plot"}' with lines notitle
""")
    def compute_tree_MDS(self):
        if not self.R_exe:
            return
        print ("\nGenerate MDS plots of topology burnin: ", end='',flush=True)

        C = min(self.n_chains(),4)
        N = int(800/self.n_chains())
        dist_cmd = ['trees-distances','matrix',f'--max={N}', '--jitter=0.3', '-V']

        if self.subsample is not None and self.subsample != 1:
            dist_cmd.append(f"--subsample={self.subsample}")

        if self.burnin is not None:
            dist_cmd.append(f"--skip={self.burnin}")

        if self.until is not None:
            dist_cmd.append(f"--until={self.until}")

        tree_files = self.get_trees_files()[0:C]
        matfile = self.outdir / "tree-MDS.M"
        countfile = self.outdir / "tree-MDS.count"

        script = self.get_libexec_script(f"tree-plot.R")
        outfile = self.outdir / "tree-MDS.svg"

        script3d = self.get_libexec_script(f"tree-plot-3D.R")
        outfile3d = self.outdir / "tree-MDS.points.html"

        if not more_recent_than_all_of(matfile, tree_files):
            self.exec_show(dist_cmd+tree_files, outfile=matfile,errfile=countfile)

        tree_counts = get_tree_counts(countfile)

        if not more_recent_than(outfile, matfile):
            self.Rexec(script,[matfile,outfile]+tree_counts)

        if not more_recent_than(outfile3d,matfile):
            point_string = self.Rexec(script3d, [matfile]+tree_counts)
            self.write_x3d_file(self.outdir, "tree-MDS.points", point_string)

        print(" done.")

    def write_x3d_file(self,dir,filename,point_string):
        assert(point_string is not None)

        if dir:
            filename = dir / (filename+".html")
        with open(filename,"w+",encoding='utf-8') as x3d_file:
            print("""\
<html>
 <head>
   <title>MDS 3D Plot</title>
   <script type='text/javascript' src='http://www.x3dom.org/download/x3dom.js'> </script>
   <link rel='stylesheet' type='text/css' href='http://www.x3dom.org/download/x3dom.css'></link>
    <style>
      x3d { border:2px solid darkorange; }
    </style>
  </head>
  <body>
    <x3d width='1000px' height='1000px'>
      <scene>""",file=x3d_file)
            print(get_x3d_of_mds(point_string),file=x3d_file)
            print("""\
      </scene>
    </x3d>
  </body>
</html>""",file=x3d_file)

    def compute_initial_alignments(self):
        print("Computing initial alignments: ",end='',flush=True)
        alignment_names = []
        for i in range(self.run(0).n_partitions()):
            name = f"P{i+1}.initial"
            source = Path(self.run(0).dir , f"C1.{name}.fasta")
            dest= Path(self.outdir, "Work", name+"-unordered.fasta")
            if not more_recent_than(dest,source):
                # strip empty columns
                cmd = ['alignment-cat', '-p', '-e', source]
                self.exec_show(cmd, outfile=dest)
            alignment_names.append(name)
        print(" done.")

        for i in range(len(alignment_names)):
            name = alignment_names[i]
            self.alignments.append((name, self.get_alphabets()[i], "Initial"))
            self.make_ordered_alignment(name)


    def compute_wpd_alignments(self):
        print("\nComputing WPD alignments: ", end='',flush=True)
        for i in range(self.n_partitions()):
            if self.get_imodel_for_partition(i) is None:
                continue
            afiles = self.get_alignments_for_partition(i)
            name = f"P{i+1}.max"
            self.alignments.append((name, self.get_alphabets()[i], "Best (WPD)"))
            if not more_recent_than_all_of(self.outdir / f"Work/{name}-unordered.fasta",afiles):
                cut_cmd=['cut-range']+afiles
                if self.burnin is not None:
                    cut_cmd.append(f"--skip={self.burnin}")
                if self.until is not None:
                    cut_cmd.append(f"--until={self.until}")
                p1 = subprocess.Popen(cut_cmd,  stdout=subprocess.PIPE)

                chop_cmd=['alignment-chop-internal','--tree',self.outdir / 'MAP.tree']
                p2 = subprocess.Popen(chop_cmd, stdout=subprocess.PIPE, stdin=p1.stdout)

                max_cmd=['alignment-max']
                self.exec_show(max_cmd,stdin=p2.stdout, outfile=self.outdir / f"Work/{name}-unordered.fasta")

                p1.wait()
                p2.wait()
                self.make_ordered_alignment(name)
        print(" done.")

    def reorder_alignment_by_tree(self,alignment,tree,outfile):
        if not path.exists(tree):
            print(f"Can't reorder alignment by tree '{tree}': tree file does not exist!")
        if file_is_empty(tree):
            print(f"Can't reorder alignment by tree '{tree}': tree file is empty!")
        if not path.exists(alignment):
            print(f"Can't reorder alignment '{tree}' by tree: alignment file does not exist!")
        if file_is_empty(alignment):
            print(f"Can't reorder alignment '{tree}' by tree: alignment file is empty!")
        if not more_recent_than_all_of(outfile,[tree,alignment]):
            cmd = ['alignment-cat', alignment, f'--reorder-by-tree={tree}']
            self.exec_show(cmd, outfile=outfile)
        assert(path.exists(outfile))

    def make_ordered_alignment(self,alignment):
        ufilename = self.outdir / f"Work/{alignment}-unordered.fasta"
        filename = self.outdir / f"{alignment}.fasta"
        if not more_recent_than_all_of(filename, [ufilename,self.outdir / "c50.tree"] ):
            self.reorder_alignment_by_tree(ufilename,self.outdir / "c50.tree",filename)
        assert(path.exists(filename))
        return filename

    def color_scheme_for_alphabet(self,alphabet):
        if alphabet == "Amino-Acids":
            return "AA+contrast"
        elif re.match(r'Numeric\[\d+\]',alphabet):
            return "plain+contrast"
        else:
            return "DNA+contrast"

    def draw_alignment(self,filename,**kwargs):
        cmd = ['alignment-draw',filename]

        if "color_scheme" in kwargs:
            color_scheme = kwargs["color_scheme"]
            if color_scheme is not None:
                cmd += ['--color-scheme',color_scheme]

        if "ruler" in kwargs:
            if kwargs["ruler"]:
                cmd += ['--show-ruler']

        if "AU" in kwargs:
            AU = kwargs["AU"]
            if AU is not None:
                cmd += ['--AU',AU]

        if "outfile" in kwargs:
            outfile = kwargs["outfile"]
        else:
            outfile = None
        if outfile is None:
            outfile = os.path.splitext(filename)[0]+'.html'

        self.exec_show(cmd,outfile=outfile)
        return outfile

    def draw_alignments(self):
        if not self.alignments:
            return

        print("Drawing alignments: ", end='', flush=True)
        for (alignment,alphabet,name) in self.alignments:
            filename = self.outdir / f"{alignment}.fasta"
            color_scheme = self.color_scheme_for_alphabet(alphabet)
            self.draw_alignment(filename, color_scheme=color_scheme, ruler=True)
            print('*',end='',flush=True)

        for i in range(self.n_partitions()):
            if self.get_imodel_for_partition(i) is None:
                continue

            outfile = self.outdir / f"P{i+1}.initial-diff.AU"
            initial = self.outdir / f"P{i+1}.initial.fasta"
            wpd = self.outdir / f"P{i+1}.max.fasta"
            self.exec_show(['alignments-diff',initial,wpd],outfile=outfile)

            outhtml = self.outdir / f"P{i+1}.initial-diff.html"
            self.exec_show(['alignment-draw',initial,'--scale=identity','--AU',outfile,'--show-ruler','--color-scheme=diff[1]+contrast'],outfile=outhtml)
            print("*",end='',flush=True);

        print(" done.")


    def compute_ancestral_states(self):
        print("Computing ancestral state alignment: ",end='',flush=True)
        for i in range(self.n_partitions()):
            afiles = self.get_alignments_for_partition(i)
            name = f'P{i+1}.ancestors'
            if self.get_imodel_for_partition(i) is None:
                bad = 0
                for afile in afiles:
                    if afile is None or not path.exists(afile):
                        bad += 1
                if bad > 0:
                    continue
            assert(len(afiles) == self.n_chains())

            self.alignments.append((name,self.get_alphabets()[i], "Ancestral"))

            template = self.outdir / f"P{i+1}.max.fasta"
            if self.get_imodel_for_partition(i) is None:
                template = self.outdir / f"P{i+1}.initial.fasta"
            tree = self.outdir / "greedy.tree"
            cmd = ['summarize-ancestors',template,'-n',tree,'--nodes-min=0.1']
#           QUESTION: Should we make a flag to log ancestors at the end of branches too?
            if False:
                cmd += ['-g',tree,'--groups-min=0.1']
            for dir in [run.dir for run in self.mcmc_runs]:
                cmd += ['-A',dir / f'C1.P{i+1}.fastas',
                        '-T',dir / 'C1.trees']
            output = self.outdir / f"P{i+1}.ancestors.fasta"
            if not more_recent_than_all_of(output,self.get_trees_files()+list(filter(lambda x:x is not None,flatten(self.get_alignments_files())))):
                self.exec_show(cmd,outfile=output)
        print(" done.")

    def compute_and_draw_AU_plots(self):
        for (alignment,alphabet,name) in self.alignments:
            # Don't try and compute AU plots for ancestral sequences
            m = re.match(r'^P([0-9]+).ancest.*',alignment)
            if m:
                continue

            m = re.match(r'^P([0-9]+).*',alignment)
            if m:
                i = int(m.group(1))-1
                afile = self.outdir / f"{alignment}.fasta"
                afiles = self.get_alignments_for_partition(i)
                if None in afiles:
                    continue
                print(f"Generating AU values for '{alignment}'...", end='',flush=True)
                AUfile = self.outdir / f"{alignment}-AU.prob"
                if not more_recent_than_all_of(AUfile, afiles):
                    cut_cmd=['cut-range']+afiles
                    if self.burnin is not None:
                        cut_cmd.append(f"--skip={self.burnin}")
                    if self.until is not None:
                        cut_cmd.append(f"--until={self.until}")
                    p1 = subprocess.Popen(cut_cmd,  stdout=subprocess.PIPE)

                    chop_cmd=['alignment-chop-internal','--tree',self.outdir / 'MAP.tree']
                    p2 = subprocess.Popen(chop_cmd, stdout=subprocess.PIPE, stdin=p1.stdout)

                    gild_cmd = ['alignment-gild',afile,self.outdir / 'MAP.tree','--max-alignments=500']
                    self.exec_show(gild_cmd, stdin=p2.stdout, outfile=AUfile)

                    p1.wait()
                    p2.wait()

                html_file = self.outdir / f"{alignment}-AU.html"
                if not more_recent_than_all_of(html_file,[afile,AUfile]):
                    color_scheme=self.color_scheme_for_alphabet(alphabet)
                    color_scheme += "+fade+fade+fade+fade"
                    self.draw_alignment(afile,ruler=True,AU=AUfile,color_scheme=color_scheme,outfile=html_file)
                print(' done.')



    def print_index_html(self,filename):
        with filename.open("w+",encoding='utf-8') as index:
            title = "MCMC Post-hoc Analysis"
            output  = self.html_header(title)

            #FIXME: We should only write sections that are not empty for this analysis!
            output += self.topbar()
            output += '<div class="content">\n'
            output += f"<h1>{title}</h1>\n"
            if self.has_model():
                output += self.section_data_and_model()
            if self.has_parameters():
                output += self.section_scalar_variables()
            if self.has_trees():
                output += self.section_phylogeny_distribution()
            if self.has_alignments():
                output += self.section_alignment_distribution()
            output += self.section_ancestral_sequences() # FIXME!
            output += self.section_mixing()
            output += self.section_analysis()
            if self.has_model():
                output += self.section_model_and_priors()
            if self.has_code():
                output += self.section_code()
            output += self.section_end()
            print(output, file=index)

        print(f"\nReport written to '{filename}")

    def html_header(self,title):
        section = """\
<!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Transitional//EN" "http://www.w3.org/TR/xhtml1/DTD/xhtml1-transitional.dtd">
<html xmlns="http://www.w3.org/1999/xhtml">
  <head>
"""
        section += f"<title>BAli-Phy: {title}</title>\n"
        section += """\
      <style type="text/css">
      body {margin:0; padding:0}
      ol li {padding-bottom:0.5em}

      th {padding-left: 0.5em; padding-right:0.5em}
      td {padding: 0.1em;}
      .backlit2 td {padding-left: 0.5em;}
      .backlit2 td {padding-right: 1.0em;}

#topbar {
	background-color: rgb(201,217,233);
	margin: 0;
	padding: 0.5em;
	display: table;
        width: 100%;
        position: fixed;
        top: 0;
}

#topbar #menu {
//	font-size: 90%;
	display: table-cell;
	text-align: right;
        // Leave space for the menubar
        padding-right: 0.75em;
}
#topbar #path {
	font-weight: bold;
	display: table-cell;
	text-align: left;
	white-space: nowrap;
}

      .content {margin:1em; margin-top: 3em}
      .backlit td {background: rgb(220,220,220);}

      :target:before {
        content:"";
        display:block;
        height:2em; /* fixed header height*/
        margin:-2em 0 0; /* negative fixed header height */
      }

      h1 {font-size: 150%;}
      h2 {font-size: 130%; margin-top:2.5em; margin-bottom: 1em}
      h3 {font-size: 110%; margin-bottom: 0.2em}// margin-top: 0.3em}

      ul {margin-top: 0.4em;}

      .center td {text-align: center};

      *[title] {cursor:help;}
      a[title] {vertical-align:top;color:#00f;font-size:70%; border-bottom:1px dotted;}
      .floating_picture {float:left;padding-left:0.5em;padding-right:0.5em;margin:0.5em;}
      .r_floating_picture {float:right;padding-left:0.5em;padding-right:0.5em;margin:0.5em;}
      table.backlit2 {border-collapse: collapse; }
      table.backlit2 tr:nth-child(odd) td { background-color: #F0F0F0; margin:0}
      .clear {clear:both;}

      .modelname {font-weight: bold; color: blue}
      table.model {margin-left: 2em}
      table.model td {vertical-align: top}
      table.models td {vertical-align: top}

    </style>
</head>
<body>
"""
        return section

    def topbar(self):
        return """\
  <p id="topbar">
    <span id="path">Sections:</span>
    <span id="menu">
      [<a href="#data">Data+Model</a>]
      [<a href="#parameters">Parameters</a>]
      [<a href="#topology">Phylogeny</a>]
      [<a href="#alignment">Alignment</a>]
      [<a href="#mixing">Mixing</a>]
      [<a href="#analysis">Analysis</a>]
      [<a href="#models">Models+Priors</a>]
      [<a href="#code">Code</a>]&nbsp;
    </span>
"""

    def section_data_and_model(self):
        section = """\
<h2 style="clear:both"><a class="anchor" name="data"></a>Data &amp; Model</h2>
<table class="backlit2">
  <tr><th>Partition</th><th>Sequences</th><th>Lengths</th><th>Alphabet</th><th>Substitution&nbsp;Model</th><th>Indel&nbsp;Model</th><th>Scale&nbsp;Model</th></tr>
"""
        for i in range(self.n_partitions()):
            section +=  '<tr>\n'
            section += f'  <td>{i+1}</td>\n'
            section += f'  <td>{self.get_input_files()[i]}</td>\n'

            features = self.get_alignment_info(self.outdir / f"P{i+1}.initial.fasta")
            minlength = features['min_length']
            maxlength = features['max_length']
            section += f'<td>{minlength} - {maxlength}</td>\n'

            alphabet = self.get_alphabets()[i]
            section += f'<td>{alphabet}</td>\n'

            smodel = self.get_smodel_for_partition(i)
            if smodel:
                smodel = smodel['main']
                index = self.get_smodel_indices()[i]+1
                target = f"S{index}"
                link = f'<a href="#{target}">{target}</a>'
                section += f'<td>{link} = {smodel}</td>\n'
            elif self.get_smodels() is not None:
                section += '<td>none</td>\n'
            else:
                section += '<td></td>\n'

            imodel = self.get_imodel_for_partition(i)
            if imodel:
                imodel = imodel['main']
                index = self.get_imodel_indices()[i]+1
                target = f"I{index}"
                link = f'<a href="#{target}">{target}</a>'
                section += f'<td>{link} = {imodel}</td>\n'
            elif self.get_imodels() is not None:
                section += '<td>none</td>\n'
            else:
                section += '<td></td>\n'

            scale_model = self.get_scale_model_for_partition(i)
            if scale_model:
                scale_model = scale_model['main']
                index = self.get_scale_model_indices()[i]+1
                target = f"scale{index}"
                link = f'<a href="#{target}">{target}</a>'
                assign = '='
                if scale_model[0] == '~':
                    scale_model = scale_model[1:]
                    assign = '~'
                section += f'<td>{link} {assign} {scale_model}</td>\n'
            elif self.get_scale_models() is not None:
                section += '<td>none</td>\n'
            else:
                section += '<td></td>\n'

            section += '</tr>\n'
        section += '</table>\n'
        return section


    def format_hs(self,infile,outfile):
        shutil.copyfile(infile,outfile)
        if self.fourmolu_exe:
            self.exec_show([self.fourmolu_exe, '-i', outfile,'--column-limit','120'])
        elif self.ormolu_exe:
            self.exec_show([self.ormolu_exe, '-i', outfile])
        elif self.hindent_exe:
            self.exec_show([self.hindent_exe, outfile])


    def create_viewable_source(self,inpath,outpath):
        infile = open(inpath)
        code = infile.read()
        infile.close()
        code = re.sub('<','&lt;',code)

        with open(outpath, 'w') as outfile:
            print(f"""
<html>
<head>
<title>{inpath.name}</title>
<link rel="stylesheet" href="https://cdnjs.cloudflare.com/ajax/libs/highlight.js/11.8.0/styles/atom-one-light.min.css">
<script src="https://cdnjs.cloudflare.com/ajax/libs/highlight.js/11.8.0/highlight.min.js"></script>
<script src="https://cdnjs.cloudflare.com/ajax/libs/highlight.js/11.8.0/languages/haskell.min.js"></script>
<script>hljs.highlightAll();</script>
</head>
<body>
<pre><code>{code}
</code>
</pre>
</body>
</html>
""", file=outfile)

    def copy_code(self):
        for file_path in self.get_code():
            out_file_path = self.outdir / file_path.name
            self.format_hs(file_path, out_file_path)
            self.create_viewable_source(out_file_path, out_file_path.with_suffix('.html'))

    def section_code(self):
        if not self.has_code():
            return None

        section = """\
<h2 class="clear"><a class="anchor" name="code"></a>Model Code</h2>
    <ul>
"""
        for file_path in self.get_code():
            url = file_path.with_suffix(".html").name
            section += f"        <li><a href=\"{url}\">{file_path.name}</a></li>\n"
        section += "    </ul>\n"
        section += '<br/><hr/>\n'

        return section


    def section_scalar_variables(self):
        log_files = self.get_log_files()
        if log_files is None or None in log_files:
            return

        section = """\
<h2 class="clear"><a class="anchor" name="parameters"></a>Scalar variables</h2>
    <table class="backlit2">
        <tr><th>Statistic</th><th>Median</th><th title="95% Bayesian Credible Interval">95% BCI</th><th title="Auto-Correlation Time">ACT</th><th title="Effective Sample Size">ESS</th><th>burnin</th><th title="Potential Scale Reduction Factor based on width of 80% credible interval">PSRF-CI80%</th><th>PSRF-RCF</th></tr>
"""
        for var in self.median:
            if var == "iter":
                continue
            if var == "time" and self.personality() == "phylobayes":
                continue
            if var == "#treegen" and self.personality() == "phylobayes":
                continue

            row_style = "" if var in self.CI_low else 'style="font-style:italic"'
            section += f'<tr {row_style}>\n'
            section += f'  <td>{var}</td>\n'
            section += f'  <td>{self.median[var]}</td>\n'
            if var in self.CI_low:
                section += f'  <td>({self.CI_low[var]}, {self.CI_high[var]})</td>\n'

                ACT_style = ' style="color:red"' if self.ESS[var] <= self.min_ESS else ""
                section += f'  <td {ACT_style}>{self.ACT[var]}</td>\n'

                if self.ESS[var] < 100:
                    ESS_style = ' style="color:red"'
                elif self.ESS[var] < 300:
                    ESS_style = ' style="color:orange"'
                else:
                    ESS_style = ''
                section += f'  <td {ESS_style}>{self.ESS[var]}</td>\n'

                burnin_style = ' style="color:red"' if self.Burnin[var] == "Not Converged!" else ""
                section += f'  <td {burnin_style}>{self.Burnin[var]}</td>\n'

                if var in self.PSRF_CI80:
                    if self.PSRF_CI80[var] >= 1.2:
                        CI80_style = ' style="color:red"'
                    elif self.PSRF_CI80[var] >= 1.05:
                        CI80_style = ' style="color:orange"'
                    else:
                        CI80_style = ''
                    section += f'  <td {CI80_style}>{self.PSRF_CI80[var]}</td>\n'
                else:
                    section += '  <td>NA</td>\n'

                if var in self.PSRF_RCF:
                    if self.PSRF_RCF[var] >= 1.2:
                        RCF_style = ' style="color:red"'
                    elif self.PSRF_RCF[var] >= 1.05:
                        RCF_style = ' style="color:orange"'
                    else:
                        RCF_style = ''
                    section += f'  <td {RCF_style}>{self.PSRF_RCF[var]}</td>\n'
                else:
                    section += '  <td>NA</td>\n'

            else:
                section += '  <td></td>\n'
                section += '  <td></td>\n'
                section += '  <td></td>\n'
                section += '  <td></td>\n'
                section += '  <td></td>\n'
                section += '  <td></td>\n'

            section += '</tr>\n'
        section += '</table>\n'

        return section

    def section_phylogeny_distribution(self):
        section = '<h2><a class="anchor" name="topology"></a>Phylogeny Distribution</h2>\n'
        section += self.html_svg('c-levels.svg','35%','',['r_floating_picture'])
        section += self.html_svg('c50-tree.svg','25%','',['floating_picture'])

        section += """\
<table>
  <tr>
    <td>Partition support: <a href="consensus">Summary</a></td>
    <td><a href="partitions.bs">Across chains</a></td>
  </tr>
</table>
"""
        section += '<table>\n'
        for (tree,name) in self.trees:
            section += f"""\
  <tr>
    <td>{name}</td>
    <td><a href="{tree}.tree">Newick</a></td>
    <td>(<a href="{tree}.PP.tree">+PP</a>)</td>
    <td><a href="{tree}-tree.pdf">PDF</a></td>
    <td><a href="{tree}-tree.svg">SVG</a></td>
"""
            if self.subpartitions and (path.exists(self.outdir / f"{tree}.mtree") or
                                       path.exists(self.outdir / f"{tree}-mctree.svg") or
                                       path.exists(self.outdir / f"{tree}-mctree.pdf")):
                section += '<td>MC Tree:</td>\n'
            else:
                section +='<td></td>'

            if self.subpartitions and path.exists(self.outdir / f"{tree}.mtree"):
                section += f'<td><a href="{tree}.mtree">-L</a></td>\n'
            else:
                section +='<td></td>'

            if self.subpartitions and path.exists(self.outdir / f"{tree}-mctree.pdf"):
                section += f'<td><a href="{tree}-mctree.pdf">PDF</a></td>\n'
            else:
                section +='<td></td>'

            if self.subpartitions and path.exists(self.outdir / f"{tree}-mctree.svg"):
                section += f'<td><a href="{tree}-mctree.svg">SVG</a></td>\n'
            else:
                section +='<td></td>'

            section += '  </tr>\n'
        section += '</table>'
        return section

    def section_alignment_distribution(self):
        if not self.n_partitions():
            return

        section = '<h2 class="clear"><a class="anchor" name="alignment"></a>Alignment Distribution</h2>\n'
        for i in range(self.n_partitions()):
            section += f'<h3>Partition {i+1}</h3>\n'
            section += """\
  <table>
    <tr>
       <th style="padding-right:3em"></th>
       <th></th>
       <th></th>
       <th title ="Comparison of this alignment (top) to the WPD alignment (bottom)">Diff</th>
       <th></th>
       <th style="padding-right:0.5em;padding-left:0.5em" title="Percent identity of the most dissimilar sequences">Min. %identity</th>
       <th style="padding-right:0.5em;padding-left:0.5em" title="Number of columns in the alignment"># Sites</th>
       <th style="padding-right:0.5em;padding-left:0.5em" title="Number of invariant columns">Constant</th>
       <th title="Number of parsimony-informative columns.">Informative</th>
    </tr>
"""
            for (alignment,alphabet,name) in self.alignments:
                if not alignment.startswith(f'P{i+1}.'):
                    continue
                section +=  '    <tr>\n'
                section += f'      <td style="padding-right:3em">{name}</td>\n'
                section += f'      <td><a href="{alignment}.fasta">FASTA</a></td>\n'

                if path.exists(self.outdir / f"{alignment}.html" ):
                    section += f'      <td><a href="{alignment}.html">HTML</a></td>\n'
                else:
                    section +=  '      <td></td>\n'

                if path.exists(self.outdir / f"{alignment}-diff.html" ):
                    section += f'      <td><a href="{alignment}-diff.html">Diff</a></td>\n'
                else:
                    section +=  '      <td></td>\n'

                if path.exists(self.outdir / f"{alignment}-AU.html" ):
                    section += f'      <td><a href="{alignment}-AU.html">AU</a></td>\n'
                else:
                    section +=  '      <td></td>\n'

                features = self.get_alignment_info(self.outdir / f"{alignment}.fasta")
                section += f'<td style="text-align: center">{features["min_p_identity"]}%</td>\n'
                section += f'<td style="text-align: center">{features["length"]}</td>\n'
                section += f'<td style="text-align: center">{features["n_const"]} ({features["p_const"]}%)</td>\n'
                section += f'<td style="text-align: center">{features["n_inform"]} ({features["p_inform"]}%)</td>\n'

                section += '    </tr>\n'
            section += '  </table>\n'
        return section

    def get_value_from_file(self,filename,attribute):
        with open(filename,encoding='utf-8') as file:
            for line in file:
                m = re.search(attribute + ' ([^ ]*)($| )', line)
                if m:
                    return m.group(1)
        return None

    def html_svg(self,url,width=None,height=None,classes=[],extra=None):
        classes.append('svg-image')
        class_=' '.join(classes)
        svg = f'<object data="{url}" type="image/svg+xml" class="{class_}"'
        if width:
            svg += f" width={width}"
        if height:
            svg += f" height={height}"
        if extra:
            svg += " "+extra
        svg += ' ></object>'
        return svg

    def section_ancestral_sequences(self):   #REMOVE!
        section = ""
        return section

    def section_tree_mixing2(self):
        section = ""
        assert(self.has_trees())
        MDS_figure = "tree-MDS.svg"
        MDS_figure_3d = "tree-MDS.points"
        C = min(self.n_chains(),4)
        if self.n_chains() > 4:
            MDS_title = "the first 4 chains"
        else:
            MDS_title = f"{C} chains"

        section += f"""\
<table style="width:100%;clear:both">
  <tr>
    <td style="width:40%;vertical-align:top">
      <h4 style=\"text-align:center\">Projection of RF distances for {MDS_title} (<a href=\"https://doi.org/10.1080/10635150590946961\">Hillis et al 2005</a>)</h4>
"""
        if path.exists( self.outdir / MDS_figure ):
            section += self.html_svg(MDS_figure,"90%","",[])
            section += f'<a href="{MDS_figure_3d}.html">3D</a>'
        else:
            if not self.R_exe:
                section += "<p>Not generated: can't find R</p>"
        section += '</td>'
        section += '<td style="width:40%;vertical-align:top">'
        section += '<h4 style="text-align:center">Variation of split PPs across chains (<a href="https://doi.org/10.1080/10635150600812544">Beiko et al 2006</a>)</h4>'
        if path.exists(self.outdir / "convergence1-PP.svg"):
            section += self.html_svg("convergence1-PP.svg","90%","",[])
        else:
            if not self.R_exe:
                section += "<p>Not generated: can't find R.</p>"
            elif self.n_chains() == 1:
                section += "<p>Not generated: multiple chains needed.</p>"
            else:
                section += "<p>Not generated.</p>"

        if path.exists(self.outdir / "convergence2-PP.svg"):
            section += "<br/><br/><br/><br/>"
            section += self.html_svg("convergence2-PP.svg","90%","",[])
        section += "</td>"
        section += "</tr></table>"
        return section

    def section_mixing(self):
        section = """\
<h2><a class="anchor" name="mixing"></a>Mixing</h2>
  <table style="width:100%;clear:both">
    <tr>
      <td style="vertical-align:top">
        <p><b>Statistics:</b></p>
        <table class="backlit2">
"""
        burnin_before = "NA"
        min_NE = "NA"
        if self.has_parameters():
            if self.get_log_files():
                burnin_before = self.get_value_from_file(self.outdir / "Report", 'min burnin <=')
                if burnin_before == "Not":
                    burnin_before = "Not Converged!"
                min_NE = self.get_value_from_file(self.outdir / "Report", 'Ne  >=')
            section += f"<tr><td><b>scalar burnin</b></td><td>{burnin_before}</td></tr>\n"
            section += f"<tr><td><b>scalar ESS</b></td><td>{min_NE}</td></tr>\n"


        if self.has_trees():
            min_NE_partition = self.get_value_from_file(self.outdir / 'partitions.bs','min Ne =')
            section += f'<tr><td><b title=\"Effective Sample Size for bit-vectors of partition support (smallest)\">topological ESS</b></td><td>{min_NE_partition}</td></tr>\n'

            asdsf = self.get_value_from_file(self.outdir / 'partitions.bs',r'ASDSF\[min=0.100\] =')
            if asdsf is None:
                asdsf = "NA"
            else:
                asdf = asdsf.rstrip()
            section += f'<tr><td><b title=\"Average Standard Deviation of Split Frequencies\">ASDSF</b></td><td>{asdsf}</td></tr>'

            msdsf = self.get_value_from_file(self.outdir / 'partitions.bs','MSDSF =')
            if msdsf is None:
                msdsf = "NA"
            else:
                msdsf = msdsf.rstrip()
            section += f'<tr><td><b title=\"Maximum Standard Deviation of Split Frequencies\">MSDSF</b></td><td>{msdsf}</td></tr>'

        psrf_80 = 'NA'
        psrf_rcf = 'NA'
        if self.has_parameters():
            if self.get_log_files() and len(self.get_log_files()) >= 2:
                psrf_80 = self.get_value_from_file(self.outdir / 'Report','PSRF-80%CI <=')
                psrf_rcf = self.get_value_from_file(self.outdir / 'Report','PSRF-RCF <=')
            section += f'<tr><td><b>PSRF CI80%</b></td><td>{psrf_80}</td><tr>\n'
            section += f'<tr><td><b>PSRF RCF</b></td><td>{psrf_rcf}</td><tr>\n'

        section += '</table>\n'
        section += '</td>\n'

        section += '    <td>\n'
        if self.has_trees():
            section += '      <img src="c50.SRQ.png" alt="SRQ plot for supprt of 50% consensus tree."/>'
            section += '      <img src="partitions.SRQ.png" alt="SRQ plot for supprt of each partition."/>'
        section += '    </td>\n'
        section += '  </tr></table>\n'

        if self.has_trees():
            section += self.section_tree_mixing2()

        print("")
        if burnin_before:
            print(f"NOTE: burnin (scalar) <= {burnin_before}")
        if min_NE:
            print(f"NOTE: min_ESS (scalar)    = {min_NE}")
        if self.has_trees():
            if min_NE_partition:
                print(f"NOTE: min_ESS (partition)    = {min_NE_partition}")
            if asdsf:
                print(f"NOTE: ASDSF = {asdsf}")
            if msdsf:
                print(f"NOTE: MSDSF = {msdsf}")
        if psrf_80:
            print(f"NOTE: PSRF-80%CI = {psrf_80}")
        if psrf_rcf:
            print(f"NOTE: PSRF-RCF = {psrf_rcf}")
        return section

    def section_analysis(self):
        section = ""

        commands = [run.get_command() for run in self.mcmc_runs]
        versions = [run.get_version() for run in self.mcmc_runs]
        parent_dirs = [run.get_parent_dir() for run in self.mcmc_runs]

        something_common = all_same(commands) or all_same(versions) or all_same(parent_dirs)

        section += '<br/><hr/><br/>\n'
        section += '<h2><a class="anchor" name="analysis"></a>Analysis</h2>\n'
        if something_common:
            section += "<p>"
        if all_same(commands):
            section += f"<b>command line</b>: {commands[0]}</br>\n"
        if all_same(parent_dirs):
            section += f"<b>directory</b>: {parent_dirs[0]}</br>\n"
        if all_same(versions):
            section += f"<b>version</b>: {versions[0]}\n"
        if something_common:
            section += "</p>"
        #section += '<table style="width:100%">'."\n"

        #section += '<table class="backlit2 center" style="width:100%">'."\n"
        section += '<table class="backlit2 center">\n'
        section += "<tr><th>chain #</th>"
        if not all_same(versions):
            section += "<th>version</th>" 
        section += "<th>burnin</th><th>subsample</th><th>samples</th>"
        if not all_same(commands):
            section += "<th>command line</th>"
        section += "<th>subdirectory</th>"
        if not all_same(parent_dirs):
            section += "<th>directory</th>"
        section += "</tr>\n"

        print()
        for i in range(self.n_chains()):
            section += "<tr>\n"

            section += f"  <td>{i+1}</td>\n"
            if not all_same(versions):
                section += f"  <td>{versions[i]}</td>\n"

            section += f"  <td>{self.burnin}</td>\n"
            section += f"  <td>{self.subsample}</td>\n"

            remaining = (self.run(i).n_iterations() - self.burnin)/self.subsample
            print(f"RUN {i+1}: directory = {self.run(i).get_dir()}     iterations = {self.run(i).n_iterations()}     burnin = {self.burnin}")
            section += f"  <td>{remaining}</td>\n"

            if not all_same(commands):
                section += f"  <td>{commands[i]}</td>\n"
            section += f"  <td>{self.run(i).get_dir()}</td>\n"
            if not all_same(parent_dirs):
                section += f"  <td>{parent_dirs[i]}/td>\n"

            section += "</tr>\n"

        section += "</table>\n"
        return section

    def section_model_and_priors(self):
        section = '<h2 style="clear:both"><a class="anchor" name="models"></a>Model and priors</h2>\n'

        section += '<h3 style="clear:both"><a class="anchor" name="tree"></a>Tree (+priors)</h3>\n'
#    $section .= "<table class=\"backlit2\">\n";
        section += "<table>\n"
        if self.get_topology_prior():
            section += f'<tr><td class="modelname">topology</td><td>{self.get_topology_prior()}</td></tr>'
        if self.get_branch_length_prior():
            section += f'<tr><td class="modelname">branch lengths</td><td>{self.get_branch_length_prior()}</td></tr>'
        if self.get_tree_prior():
            section += f'<tr><td class="modelname">tree</td><td>{self.get_tree_prior()}</td></tr>'
        section += '</table>\n'

        section += '<h3 style="clear:both"><a class="anchor" name="smodel"></a>Substitution model (+priors)</h3>\n'
        section += print_models("S",self.get_smodels())

        section += '<h3 style="clear:both"><a class="anchor" name="imodel"></a>Indel model (+priors)</h3>\n'
        section += print_models("I",self.get_imodels())

        if (self.get_scale_models()):
            section += '<h3 style="clear:both"><a class="anchor" name="scales"></a>Scales (+priors)</h3>\n'
            section += print_models("scale",self.get_scale_models());

        return section

    def section_end(self):
        return """\
     </div>
   </body>
</html>
"""


if __name__ == '__main__':
    parser = argparse.ArgumentParser(description="Generate an HTML report summarizing MCMC runs for BAli-Phy and other software.",
                                     epilog= "Examples:\n   bp-analyze analysis-dir-1/ analysis-dir-2/\n   bp-analyze trees1.trees trees2.trees\n   bp-analyze log1.log log2.log",
                                     formatter_class=argparse.RawTextHelpFormatter)

    parser.add_argument("mcmc_outputs", default=['.'], help="Subdirectories with MCMC runs",nargs='*')
    parser.add_argument("--clean", default=False,help="Delete generated files",action='store_true')
    parser.add_argument("--verbose",default=0, help="Be verbose",action='store_true')
    parser.add_argument("--skip", metavar='NUM',type=int,default=None,help="Skip NUM iterations as burnin")
    parser.add_argument("--subsample",metavar='NUM',type=int,default=1,help="Keep only every NUM iterations")
    parser.add_argument("--until",type=int,default=None)
    parser.add_argument("--prune", default=None,help="Taxa to remove")
    parser.add_argument("--outdir", default="Results",help="Output directory")
#    parser.add_argument("--muscle")
#    parser.add_argument("--probcons")
#    parser.add_argument("--mafft")
    parser.add_argument("--subpartitions",default=False,action='store_true')
    args = parser.parse_args()

    if args.verbose:
        logging.basicConfig(level=logging.DEBUG)

    mcmc_outputs = []
    for p in args.mcmc_outputs:
        mcmc_outputs.append(Path(p))

    analysis = Analysis(args, mcmc_outputs, args.outdir)