#!/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 csv
import re
import os
from os import path
from pathlib import Path
import subprocess
import json
import itertools
import glob
import html
import sys
import logging
import shlex
import tempfile
from contextlib import ExitStack

scores = {'multiply': 'very short/log',
          'wsum2':    'short/EA',
          'wsum':     'medium/AMA',
          'sum':      'long/TC'}

# Read ordinary and conditioned statreport output without assuming diagnostics follow each summary.
# Each call owns its dictionaries, so conditional estimates cannot overwrite the main scalar report.
def parse_scalar_report(filename):
    results = {name: {} for name in (
        "mean", "stddev", "median", "CI_low", "CI_high", "ACT", "ESS", "Burnin",
        "PSRF_CI80", "PSRF_RCF", "model_support_statistics", "matching_samples")}
    results["constants"] = set()
    var = None
    for raw in Path(filename).read_text(encoding="utf-8").splitlines():
        line = raw.strip()
        if not line or line.startswith("Condition:"):
            continue
        m = re.fullmatch(r'Matching samples(?: \[([0-9]+)\])? = ([0-9]+)', line)
        if m:
            results["matching_samples"][m[1] or "pooled"] = int(m[2])
            continue
        m = re.fullmatch(r'(.+):\s+posterior-probability = ([^\s]+)\s+\[log-odds = ([^\]]+)\]', line)
        if m:
            results["model_support_statistics"][m[1]] = {
                "posterior_probability": m[2], "log_odds": m[3]}
            continue
        m = re.fullmatch(r'E ([^\s]+) = ([^\s]+)\s+\[\+- ([^\s]+)\]', line)
        if m:
            results["mean"][m[1]], results["stddev"][m[1]] = m[2], m[3]
            continue
        m = re.search(r'([^\s]+) ~ ([^\s]+)\s+\((.+),(.+)\)', line)
        if m:
            var = m[1]
            results["median"][var] = m[2]
            results["CI_low"][var], results["CI_high"][var] = m[3], m[4]
            continue
        m = re.search(r't @ (.+)\s+Ne = ([^ ]+)\s+burnin = (Not Converged!|[^ ]+)', line)
        if m and var is not None:
            results["ACT"][var], results["ESS"][var] = float(m[1]), float(m[2])
            results["Burnin"][var] = m[3] if m[3] == "Not Converged!" else float(m[3])
            continue
        m = re.search(r'PSRF-80%CI = ([^ ]+)\s+PSRF-RCF = ([^ ]+)', line)
        if m and var is not None:
            results["PSRF_CI80"][var], results["PSRF_RCF"][var] = float(m[1]), float(m[2])
            continue
        m = re.fullmatch(r'([^=]+?) = ([^=]+)', line)
        if m:
            var, value = m[1], m[2].strip()
            results["median"][var] = value
            try:
                float(value)
            except ValueError:
                pass
            else:
                results["constants"].add(var)
    return results

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_string(model):
    m = model
    assign = '='
    if m[0] == '~':
        assign = '~'
        m = m[1:]
    return f'<code class="model-expression">{assign} {html.escape(m)}</code>'

def print_model(model):
    main = model['main']
    section = print_model_string(main) + '\n'
    extracted = model['extracted']

    if extracted:
        section += '<table class="model" role="presentation">\n<tbody>\n'
        for arg in extracted:
            name = arg[0]
            submodel = arg[1]
            submodel = print_model(submodel)
            section += (
                f'<tr><td class="model-variable"><code class="model-identifier">{html.escape(name)}</code></td>'
                f'<td>{submodel}</td></tr>\n'
            )
        section += '</tbody>\n</table>'
    return section


def print_models(name,models):
    section = '<table class="models" role="presentation">\n<tbody>\n'
    i = 1
    for model in models:
        model_name =  f"{name}{i}"
        model_html = print_model(model)
        section += (
            f'<tr><td><span id="{model_name}" class="modelname">{model_name}</span></td>'
            f'<td>\n'
        )
        section += model_html
        section += '</td></tr>\n'
        i += 1
    section += '</tbody>\n</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):
    # A singleton has no position along a gradient, so use the midpoint rather than an endpoint.
    if total <= 1:
        return (start + end) / 2
    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:
        # A flat MDS dimension carries no displacement; center it while scaling varying axes.
        xx = 0 if xw == 0 else (x - xmin)/xw*5.0 - 2.5
        yy = 0 if yw == 0 else (y - ymin)/yw*5.0 - 2.5
        zz = 0 if zw == 0 else (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]
    if len(qfilenames) == 1:
        return qfilenames[0]
    if len(qfilenames) == 2:
        return " or ".join(qfilenames)
    return ", ".join(qfilenames[:-1]) + ", or " + qfilenames[-1]

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.character_property_files = None
        self.column_name_map = 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_character_property_files(self):
        return self.character_property_files

    def get_column_name_map(self):
        return self.column_name_map

    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_subst_rates(self):
        return None

    def get_indel_rates(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())
        if n <= 0:
            print(f"Error: Tree file '{self.get_trees_file()}' has no samples!")
            exit(1)
        return n

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'
        self.run_file = check_file_exists(self.get_dir() / 'C1.run.json')
        self.run_info = get_json_from_file(self.run_file)

        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.character_property_files = self.find_character_property_files()
        column_name_map_file = self.get_dir() / "C1.log.column-map.json"
        if column_name_map_file.exists():
            self.column_name_map = get_json_from_file(column_name_map_file)

        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.subst_rates = self.run_info["subst rates"] if "subst rates" in self.run_info else None
        self.indel_rates = self.run_info["indel rates"] if "indel rates" in self.run_info else 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 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_subst_rates(self):
        return self.subst_rates

    def get_indel_rates(self):
        return self.indel_rates

    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

    # Return the optional category-state/property log associated with each partition.
    def find_character_property_files(self):
        filenames = []
        for partition in range(1, self.n_partitions() + 1):
            filename = self.get_dir() / f"C1.P{partition}.site-property-samples.jsonl"
            filenames.append(filename if filename.exists() else None)
        return filenames

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

    def get_code(self):
        return self.code


class TreeLogFileRun(TreeFileRun,LogFileRun):
    def __init__(self,mcmc_output, **kwargs):
        MCMCRun.__init__(self, mcmc_output, **kwargs)
        self.trees_file = add_suffix(mcmc_output, ".trees")
        self.log_file = add_suffix(mcmc_output, ".log")
        logging.debug(f"TreeLogFile Run: {mcmc_output}")

    def n_iterations(self):
        # A paired tree file may be Nexus; the tabular log has exactly one header and one row per sample.
        return LogFileRun.n_iterations(self)

# 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.run.json').exists() and \
           (Path(mcmc_dir, 'C1.log').exists() or Path(mcmc_dir, 'C1.log.json').exists()):
            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")
        self.R_exe = find_exe('R', message='Some report plots 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.sub_partitions
        self.node_bubbles = args.node_bubbles
        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(f"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()
            self.summarize_branch_models()
        if self.has_trees():
            self.sampled_topology_count = self.summarize_topology_distribution()
            if self.node_bubbles:
                self.compute_mean_branch_lengths()
            self.draw_trees()
            self.compute_tree_mixing_diagnostics()
            self.compute_srq_plots()
            self.compute_tree_MDS()
        if self.has_alignments():
            self.summarize_character_properties()
            self.compute_initial_alignments()
            # Selection-table coordinates require the posterior display alignments first.
            self.compute_wpd_alignments()
            self.summarize_positive_selection()
            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_topology_variation(self):
        return self.sampled_topology_count > 1

    def has_mixing_results(self):
        return self.has_parameters() or (self.has_trees() and self.has_topology_variation())

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

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

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

    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 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_subst_rates(self):
        return first_all_same([run.get_subst_rates() for run in self.mcmc_runs])

    def get_indel_rates(self):
        return first_all_same([run.get_indel_rates() 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_character_property_files(self):
        return [run.get_character_property_files() for run in self.mcmc_runs]

    # Use the shared complete TSV name mapping only when every chain provides the same mapping.
    def get_column_name_map(self):
        mappings = [run.get_column_name_map() for run in self.mcmc_runs]
        if None in mappings:
            return None
        return first_all_same(mappings)

    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]

    # Collect one partition's optional raw property log from every chain.
    def get_character_properties_for_partition(self, p):
        return [
            run.get_character_property_files()[p]
            if run.get_character_property_files() is not None else None
            for run in self.mcmc_runs
        ]

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

        with ExitStack() as opened_files:
            subargs = dict()
            if "infile" in kwargs:
                infile = kwargs["infile"]
                subargs["stdin"] = opened_files.enter_context(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"] = opened_files.enter_context(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"] = opened_files.enter_context(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

    # Run a streaming command pipeline, checking every stage and removing incomplete output.
    def exec_pipeline(self, commands, outfile):
        showcmd = " | ".join(shlex.join([str(word) for word in command]) for command in commands)
        showcmd += f" > {shlex.quote(str(outfile))}"
        print(showcmd, file=self.log_shell_cmds)

        processes = []
        errors = []
        with ExitStack() as resources:
            output = resources.enter_context(open(outfile, "wb"))
            previous_stdout = None
            for index, command in enumerate(commands):
                error_output = resources.enter_context(tempfile.TemporaryFile(mode="w+", encoding="utf-8"))
                stdout = output if index == len(commands) - 1 else subprocess.PIPE
                process = subprocess.Popen(command, stdin=previous_stdout, stdout=stdout, stderr=error_output)
                if previous_stdout is not None:
                    previous_stdout.close()
                previous_stdout = process.stdout
                processes.append(process)
                errors.append(error_output)

            for process in reversed(processes):
                process.wait()

            error_messages = []
            for error_output in errors:
                error_output.seek(0)
                error_messages.append(error_output.read())

        for index, error_message in enumerate(error_messages):
            if error_message:
                print(f"  stage {index + 1} err: {error_message}", file=self.log_shell_cmds)

        failed_stage = next((index for index, process in enumerate(processes) if process.returncode != 0), None)
        if failed_stage is not None:
            if path.exists(outfile):
                os.remove(outfile)
            code = processes[failed_stage].returncode
            print(f"command pipeline stage {failed_stage + 1} failed: {showcmd}", file=sys.stderr)
            print(f" exit: {code}", file=self.log_shell_cmds)
            raise SystemExit(code)

        if self.verbose:
            print(f"\n\t{showcmd}\n")

    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):
        resolved_outdir = self.outdir.resolve()
        for run in self.mcmc_runs:
            input_directory = run.get_dir().resolve()
            try:
                input_directory.relative_to(resolved_outdir)
            except ValueError:
                pass
            else:
                print(
                    f"Error: output directory '{self.outdir}' is an input directory or contains one.",
                    file=sys.stderr,
                )
                exit(1)

        character_property_file_names = [
            [str(filename) if filename is not None else None for filename in files]
            if files is not None else None
            for files in self.get_character_property_files()
        ]
        if self.outdir.exists() and not self.outdir.is_dir():
            print(f"Error: output path '{self.outdir}' is not a directory.", file=sys.stderr)
            exit(1)
        existing_entries = list(self.outdir.iterdir()) if self.outdir.exists() else []
        if existing_entries:
            properties = self.load_analysis_properties()
            identifying_properties = {"burnin", "subsample", "input_files"}
            recognized_report = (self.outdir / "Work").is_dir() and identifying_properties <= properties.keys()
            if not recognized_report:
                print(
                    f"Error: existing directory '{self.outdir}' is not a recognized bpy-summarize report; "
                    "leaving it unchanged.",
                    file=sys.stderr,
                )
                exit(1)

        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("alphabets", self.get_alphabets())
        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("character_property_file_names", character_property_file_names)
        reuse = reuse and self.check_analysis_property("input_files", self.get_input_files())

        if not reuse and existing_entries:
            new_dir_name = get_unused_dir_name(self.outdir)
            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)
            print(f"Creating new directory '{self.outdir}' for summary files.")
        (self.outdir / "Work").mkdir(exist_ok=True)

        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("alphabets", self.get_alphabets())
        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()])
        self.write_analysis_property("character_property_file_names", character_property_file_names)

    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 + 1} has {max_iterations} iterations.")
            print(f"  Chain #{min_i + 1} 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='')
        report_file = self.outdir / "Report"
        report_is_current = more_recent_than_all_of(report_file, log_files)
        if report_is_current:
            report_text = report_file.read_text(encoding='utf-8')
            # Timestamp-current reports made by older versions still need mean and SD records.
            if re.search(r'^\s*[^\s]+ ~ ', report_text, re.MULTILINE) and not re.search(
                r'^\s*E [^\s]+ = ', report_text, re.MULTILINE
            ):
                report_is_current = False

        if not report_is_current:
            cmd = ["statreport", "--mean", "--median"] + 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=report_file)
        print("done.")

        print("Analyzing scalar variables: ",end='',flush=True)
        results = parse_scalar_report(report_file)
        for name, values in results.items():
            setattr(self, name, values)

        print("done.")

    # Match the outer displayed constructor, not arbitrary occurrences inside a composite model.
    def branch_model_kinds(self):
        models = self.get_smodels() or []
        kinds = {}
        for index, model in enumerate(models, 1):
            match = re.search(r'(?:^|\+>\s*)(BranchModel(?:_test)?)\s*$', model["main"])
            if match:
                kinds[index] = match[1]
        return kinds

    # Only matching per-chain metadata permits automatic omission; absent metadata is not an error.
    def branch_category_usage(self, smodel, hypotheses):
        usages = []
        for run in self.mcmc_runs:
            filename = run.get_dir() / "C1.branch-category-usage.json"
            if not filename.exists():
                usages.append(None)
                continue
            data = get_json_from_file(filename)
            if not isinstance(data, dict):
                raise ValueError(f"Invalid branch-category metadata in {filename}")
            models = data.get("substitution_models")
            usage = data.get("used_categories")
            if (type(data.get("version")) is not int or data["version"] != 1 or not isinstance(models, list)
                    or any(type(n) is not int or n < 1 or n > len(self.get_smodels()) for n in models)
                    or len(set(models)) != len(models)
                    or not isinstance(usage, list) or not usage):
                raise ValueError(f"Invalid branch-category metadata in {filename}")
            if smodel not in models:
                usages.append(None)
                continue
            if len(usage) != len(hypotheses):
                raise ValueError(f"Hypothesis count differs for S{smodel} in {filename}")
            for categories in usage:
                if (not isinstance(categories, list)
                        or any(type(c) is not int or c < 0 for c in categories)
                        or categories != sorted(set(categories))):
                    raise ValueError(f"Invalid categories for S{smodel} in {filename}")
            usages.append(usage)
        if not usages or any(u is None for u in usages) or not all_same(usages):
            return None
        return usages[0]

    # Reuse unconditional summaries and regenerate only the hypothesis-filtered reports.
    # Full logged paths keep shared parameters and model instances distinct.
    def summarize_branch_models(self):
        self.branch_model_reports = []
        mapping = self.get_column_name_map() or {}
        for smodel, kind in self.branch_model_kinds().items():
            prefix = f"S{smodel}/{kind}:"
            omegas = {}
            support = {}
            for full, short in mapping.items():
                if not full.startswith(prefix):
                    continue
                field = full[len(prefix):]
                match = re.fullmatch(r'omegas\[([1-9][0-9]*)\]', field)
                if match:
                    omegas[int(match[1]) - 1] = short
                match = re.fullmatch(r'LogOddsBranchHypothesis\[([0-9]+)\]', field)
                if match:
                    support[int(match[1])] = self.model_support_statistics.get(short, {})
            selector = mapping.get(prefix + "hypothesis")
            hypotheses = sorted(support)
            if hypotheses and hypotheses != list(range(len(hypotheses))):
                raise ValueError(f"Nonconsecutive hypotheses for S{smodel}")
            usage = (self.branch_category_usage(smodel, hypotheses)
                     if kind == "BranchModel_test" and hypotheses else None)
            report = {"smodel": smodel, "kind": kind, "omegas": omegas,
                      "support": support, "usage": usage, "selector": selector, "conditional": {}}
            self.branch_model_reports.append(report)
            if kind != "BranchModel_test" or not omegas or selector is None:
                continue
            for hypothesis in hypotheses:
                used = usage[hypothesis] if usage is not None else sorted(omegas)
                selected = [c for c in used if c in omegas]
                filename = self.outdir / f"S{smodel}.hypothesis-{hypothesis}.Report"
                command = ["statreport", "--mean", "--median", f"--condition={selector}={hypothesis}"]
                command += [f"--select={omegas[c]}" for c in selected]
                # Even without logged used categories, request the selector for matching counts
                # instead of letting an empty --select list request every log column.
                if not selected:
                    command.append(f"--select={selector}")
                if self.burnin is not None:
                    command.append(f"--skip={self.burnin}")
                if self.subsample is not None and self.subsample != 1:
                    command.append(f"--subsample={self.subsample}")
                if self.until is not None:
                    command.append(f"--until={self.until}")
                command += self.get_log_files()
                self.exec_show(command, outfile=filename)
                report["conditional"][hypothesis] = {
                    "filename": filename, "categories": {c: omegas[c] for c in selected},
                    "missing": [c for c in used if c not in omegas],
                    "statistics": parse_scalar_report(filename)}

    # Format global category estimates, never interpreting categories as branches or alignment sites.
    def branch_omega_table(self, omegas, statistics):
        section = '<table class="backlit2"><thead><tr><th>Category</th><th>Mean</th><th>SD</th>'
        section += '<th>Median</th><th>95% credible interval</th></tr></thead><tbody>\n'
        for category, field in sorted(omegas.items()):
            median = statistics["median"].get(field, "NA")
            constant = field in statistics["constants"]
            mean = median if constant else statistics["mean"].get(field, "NA")
            sd = "—" if constant else statistics["stddev"].get(field, "NA")
            interval = "—" if constant else (
                f'({statistics["CI_low"].get(field, "NA")}, {statistics["CI_high"].get(field, "NA")})')
            section += '<tr>' + ''.join(f'<td>{html.escape(str(v))}</td>'
                                       for v in (category, mean, sd, median, interval)) + '</tr>\n'
        return section + '</tbody></table>\n'

    # Present model support separately from category estimates conditional on the sampled hypothesis.
    def section_branch_models(self):
        section = '<section id="branch-models" class="report-section"><h2>Branch models</h2>\n'
        section += ('<p>Category ω values describe model categories, not individual branches or '
                    'variation among sites. Category properties in the alignment viewer have the same meaning.</p>\n')
        for report in self.branch_model_reports:
            smodel = report["smodel"]
            partitions = [str(i + 1) for i, s in enumerate(self.get_smodel_indices()) if s == smodel - 1]
            section += f'<h3><a href="#S{smodel}">S{smodel}</a> — {report["kind"]}</h3>\n'
            section += '<p>Partitions: ' + ', '.join(f'<a href="#P{p}">P{p}</a>' for p in partitions) + '</p>\n'
            if report["kind"] == "BranchModel":
                if report["omegas"]:
                    section += self.branch_omega_table(report["omegas"], vars(self))
                else:
                    section += '<p>Category ω values were not logged; see the ordinary scalar summaries.</p>\n'
                continue
            section += '<h4>Hypothesis support</h4>\n'
            section += ('<p>Probabilities are Rao–Blackwellized estimates. Log odds compare each hypothesis '
                        'with all other hypotheses combined; they are not Bayes factors.</p>\n')
            section += ('<table class="backlit2"><thead><tr><th>Hypothesis</th><th>Posterior probability</th>'
                        '<th>Posterior log odds</th><th>Used categories</th></tr></thead><tbody>\n')
            for h, support in sorted(report["support"].items()):
                categories = ', '.join(map(str, report["usage"][h])) if report["usage"] is not None else 'Unknown'
                cells = (h, support.get("posterior_probability", "NA"), support.get("log_odds", "NA"), categories)
                section += '<tr>' + ''.join(f'<td>{html.escape(str(c))}</td>' for c in cells) + '</tr>\n'
            section += '</tbody></table>\n'
            section += '<p>Unconditional category summaries include samples where that category is unused.</p>\n'
            if report["usage"] is None:
                section += '<p>Category usage is unavailable; some entries may be unused.</p>\n'
            if not report["conditional"]:
                if not report["omegas"]:
                    section += '<p>Category ω values were not logged; conditional estimates are unavailable.</p>\n'
                if report["selector"] is None:
                    section += '<p>The hypothesis field was not logged; conditional estimates are unavailable.</p>\n'
            for h, conditional in sorted(report["conditional"].items()):
                statistics = conditional["statistics"]
                counts = statistics["matching_samples"]
                section += f'<h4>ω conditional on hypothesis {h}</h4>\n'
                if conditional["missing"]:
                    missing = ', '.join(map(str, conditional["missing"]))
                    section += f'<p>ω values were not logged for used categories {missing}; their estimates are unavailable.</p>\n'
                chain_counts = ', '.join(f'chain {c}: {n}' for c, n in counts.items() if c != "pooled")
                section += f'<p>Matching samples: {counts.get("pooled", 0)} ({chain_counts}).</p>\n'
                if counts.get("pooled", 0):
                    section += self.branch_omega_table(conditional["categories"], statistics)
                else:
                    section += '<p>No retained samples.</p>\n'
                name = html.escape(conditional["filename"].name)
                section += f'<p><a href="{name}">Complete conditional report</a></p>\n'
        return section + '</section>\n'

    # Recover positive-selection statistics' full model paths and the partitions using each model.
    # The column-name mapping is required because statreport intentionally shows Tracer-friendly names.
    def positive_selection_model_support(self):
        column_name_map = self.get_column_name_map()
        if column_name_map is None:
            return []
        short_to_long = {short_name: long_name for long_name, short_name in column_name_map.items()}
        if len(short_to_long) != len(column_name_map):
            raise Exception("TSV column names are not unique.")

        support = []
        for short_name, statistics in getattr(self, "model_support_statistics", {}).items():
            long_name = short_to_long.get(short_name, short_name)
            match = re.fullmatch(r'S([1-9][0-9]*)/(.+:LogOddsPosSelection)', long_name)
            if match is None:
                continue
            smodel = int(match.group(1))
            partitions = [
                partition + 1
                for partition, model_index in enumerate(self.get_smodel_indices())
                if model_index == smodel - 1
            ]
            support.append({
                "statistic": match.group(2),
                "smodel": smodel,
                "partitions": partitions,
                **statistics,
            })
        return support

    # Identify substitution models with an explicit positive-selection selector from their
    # complete logged-field names, independently of which state the selector happened to visit.
    def positive_selection_selector_models(self):
        column_name_map = self.get_column_name_map()
        if column_name_map is None:
            return None
        selector_models = set()
        for long_name in column_name_map:
            match = re.fullmatch(r'S([1-9][0-9]*)/.+:LogOddsPosSelection', long_name)
            if match is not None:
                selector_models.add(int(match.group(1)))
        return selector_models

    def has_positive_selection(self):
        return bool(getattr(self, "positive_selection_reports", [])) or bool(self.positive_selection_model_support())


    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.")

        topology_count = self.get_value_from_file(self.outdir / "consensus", "n_topologies =")
        if topology_count is None:
            raise Exception("The topology count is missing from the trees-consensus report.")
        return int(topology_count)

        # This routine computes mean branch lengths WITH NODE CIRCLES (*.ltree)
    def compute_mean_branch_lengths(self):
        print("\nComputing mean branch lengths: ",end='')
        for level in self.tree_consensus_levels:
            value = int(level*100)
            prefix = self.outdir / f"c{value}"
            tree = prefix.with_suffix(".tree")
            outfile = prefix.with_suffix(".ltree")
            if more_recent_than_all_of(outfile, self.get_trees_files()):
                continue

            cmd = ['tree-mean-lengths','--tree',tree,'--safe','--show-node-lengths']
            cmd += self.get_trees_files()

            if self.prune is not None:
                cmd.append("--ignore={}".format(self.prune))
            if self.subsample is not None and self.subsample != 1:
                cmd.append("--subsample={}".format(self.subsample))
            if self.burnin is not None:
                cmd.append("--skip={}".format(self.burnin))
            if self.until is not None:
                cmd.append("--until={}".format(self.until))

            self.exec_show(cmd,outfile=outfile)
        print(" done.")


    def get_alignment_info(self,filename,alphabet):
        output = self.exec_show(['alignment-info',filename,'--alphabet',alphabet])
        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(2)

            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'])

            treefile = tree + (".ltree" if self.node_bubbles else ".tree")
            if not more_recent_than(self.outdir / f"{tree}-tree.pdf", filename1):
                cmd = ['draw-tree', treefile, '--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', treefile, '--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):
        if not self.has_topology_variation():
            return

        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 = []
        if not self.has_topology_variation():
            return

        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:
            script = self.get_libexec_script("plot-SRQ.R")
            self.Rexec(script, [
                self.outdir / f"{srq}.SRQ",
                self.outdir / f"{srq}.SRQ.png",
                f"Scaled Regeneration Quantile (SRQ) plot: {srq}",
                srq,
            ])
        script = self.get_libexec_script("plot-support-levels.R")
        args = [self.outdir / "c-levels.plot", self.outdir / "c-levels.svg"]
        if self.subpartitions:
            args.append(self.outdir / "extended-c-levels.plot")
        self.Rexec(script, args)
    def compute_tree_MDS(self):
        if not self.has_topology_variation() or 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")
        x3dom_dir = Path("x3dom-1.8.3")
        self.copy_report_assets([
            x3dom_dir / "x3dom.js",
            x3dom_dir / "x3dom.css",
            x3dom_dir / "LICENSE.md",
            x3dom_dir / "README.md",
        ])
        with open(filename,"w+",encoding='utf-8') as x3d_file:
            print("""\
<!DOCTYPE html>
<html lang="en">
 <head>
   <meta charset="utf-8">
   <meta name="viewport" content="width=device-width, initial-scale=1">
   <title>MDS 3D Plot</title>
   <script src="bpy-summarize-assets/x3dom-1.8.3/x3dom.js"></script>
   <link rel="stylesheet" href="bpy-summarize-assets/x3dom-1.8.3/x3dom.css">
    <style>
      html, body { height:100%; margin:0; }
      x3d { border:2px solid darkorange; box-sizing:border-box; height:100%; width:100%; }
    </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)


    # Pool per-character property means only when every chain logged a partition.
    def summarize_character_properties(self):
        self.character_property_summaries = [None] * self.n_partitions()
        partition_files = [
            self.get_character_properties_for_partition(partition)
            for partition in range(self.n_partitions())
        ]
        if not any(any(filename is not None for filename in files)
                   for files in partition_files):
            return

        print("Summarizing character properties: ", end="", flush=True)
        for partition, filenames in enumerate(partition_files):
            present = sum(filename is not None for filename in filenames)
            if present == 0:
                continue
            if present != len(filenames):
                print(
                    f"\nWarning: character properties for partition {partition + 1} "
                    "were not logged by all chains; skipping this partition."
                )
                continue

            output = self.outdir / f"P{partition + 1}.site-property-summary.json"
            self.character_property_summaries[partition] = output
            if more_recent_than_all_of(output, filenames):
                continue

            command = ["character-properties", "summarize", *filenames]
            if self.burnin is not None:
                command.append(f"--skip={self.burnin}")
            if self.until is not None:
                command.append(f"--until={self.until}")
            if self.subsample is not None and self.subsample != 1:
                command.append(f"--subsample={self.subsample}")
            self.exec_show(command, outfile=output)
        print("done.")

    # Regenerate ranked positive-selection TSV tables and retain their report metadata and rows.
    # character-properties remains responsible for selecting, projecting, and sorting sites.
    def summarize_positive_selection(self):
        self.positive_selection_reports = []
        threshold = 0.5
        for partition, summary in enumerate(self.character_property_summaries):
            if summary is None:
                continue
            summary_data = get_json_from_file(summary)
            properties = sorted(
                name for name in summary_data["properties"]
                if name.endswith("posSelection")
            )
            conditioned = summary_data.get("conditioned", {}).get("positiveSelectionInModel")
            for property_name in properties:
                branch_models = {r["smodel"] for r in getattr(self, "branch_model_reports", [])}
                if (branch_models and re.fullmatch(r'branch[0-9]+-posSelection', property_name)
                        and self.get_smodel_indices()[partition] is not None
                        and self.get_smodel_indices()[partition] + 1 in branch_models):
                    continue
                use_condition = conditioned is not None and (
                    conditioned["retained_samples"] == 0 or property_name in conditioned["properties"]
                )
                report_name = property_name[:-len("posSelection")] + "positive-selection"
                report_name = re.sub(r"[^A-Za-z0-9_.-]", "-", report_name)
                output = self.outdir / f"P{partition + 1}.{report_name}.tsv"
                alignment = self.outdir / f"P{partition + 1}.initial.fasta"
                if self.get_imodel_for_partition(partition) is not None:
                    alignment = self.outdir / f"P{partition + 1}.consensus.pd-wsum.fasta"
                command = [
                    "character-properties", "positive-selection", summary, alignment, property_name,
                    "--highest=100%", "--sort=decreasing", "--format=tsv",
                    "--alphabet", self.get_alphabets()[partition],
                ]
                if not use_condition:
                    command.append("--unconditional")
                self.exec_show(command, outfile=output)

                with output.open(newline="", encoding="utf-8") as report_file:
                    rows = list(csv.DictReader(report_file, delimiter="\t"))
                model_averaged_count = sum(
                    float(row["model-averaged-probability"]) > threshold for row in rows
                )
                conditioned_count = sum(
                    bool(row["conditioned-probability"])
                    and float(row["conditioned-probability"]) > threshold
                    for row in rows
                )
                self.positive_selection_reports.append({
                    "partition": partition + 1,
                    "property": property_name,
                    "filename": output,
                    "alignment": alignment,
                    "rows": rows,
                    "conditioned": use_condition,
                    "threshold": threshold,
                    "model_averaged_count": model_averaged_count,
                    "conditioned_count": conditioned_count if use_condition else None,
                    "retained_samples": (
                        conditioned["retained_samples"] if use_condition else summary_data["retained_samples"]
                    ),
                    "total_samples": summary_data["retained_samples"],
                })


    def compute_wpd_alignments(self):
        print("\nComputing consensus alignments via posterior decoding: ", 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)

            for score,scorename in scores.items():
                name = f"P{i+1}.consensus.pd-{score}"

                self.alignments.append((name, self.get_alphabets()[i], f'Consensus (posterior-decoding: {scorename})'))

                unordered = self.outdir / f"Work/{name}-unordered.fasta"
                if not more_recent_than_all_of(unordered, 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}")

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

                    chop_cmd=['alignment-chop-internal','--tree',self.outdir / 'MAP.tree']
                    max_cmd=['alignment-max',f'--score={score}', '--alphabet', self.get_alphabets()[i]]
                    self.exec_pipeline([cut_cmd, chop_cmd, max_cmd], unordered)
                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,contrast=True):
        if alphabet == "Amino-Acids":
            return "AA" + ("+contrast" if contrast else "")
        elif re.match(r'Numeric\(\d+\)',alphabet):
            return "plain" + ("+contrast" if contrast else "")
        else:
            return "DNA" + ("+contrast" if contrast else "")

    # Map tip-only partition alignments to their shared character-property summary.
    def character_properties_for_alignment(self, alignment):
        match = re.match(r"^P([0-9]+)\.", alignment)
        excluded = re.match(r"^P[0-9]+\.ancestors", alignment) or re.search(
            r"(?:^|[.-])diff(?:[.-]|$)", alignment
        )
        if not match or excluded:
            return None
        partition = int(match.group(1)) - 1
        if partition < 0 or partition >= len(self.character_property_summaries):
            return None
        return self.character_property_summaries[partition]

    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]

        properties = kwargs.get("properties")
        alphabet = kwargs.get("alphabet")
        if properties is not None and alphabet is None:
            raise ValueError("an alphabet is required with alignment properties")
        if properties is not None:
            cmd += ['--properties', properties]
        if alphabet is not None:
            cmd += ['--alphabet', alphabet]

        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)
            properties = self.character_properties_for_alignment(alignment)
            self.draw_alignment(
                filename,
                color_scheme=color_scheme,
                ruler=True,
                properties=properties,
                alphabet=alphabet,
            )
            print('*',end='',flush=True)

        # The difference view deliberately uses text/nucleotide coordinates in both commands.
        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}.consensus.pd-wsum.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]'],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}.consensus.pd-wsum.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',
                   '--alphabet',self.get_alphabets()[i]]
            if self.burnin is not None:
                cmd.append(f'--skip={self.burnin}')
            if self.until is not None:
                cmd.append(f'--until={self.until}')
            if self.subsample is not None and self.subsample != 1:
                cmd.append(f'--thin={self.subsample}')
#           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 afile, tfile in zip(afiles, self.get_trees_files()):
                cmd += ['-A',afile,'-T',tfile]
            output = self.outdir / f"P{i+1}.ancestors.fasta"
            if not more_recent_than_all_of(output,[template,tree,*afiles,*self.get_trees_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"
                map_tree = self.outdir / 'MAP.tree'
                au_inputs = [*afiles, afile, map_tree]
                if not more_recent_than_all_of(AUfile, au_inputs):
                    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}")
                    if self.subsample is not None and self.subsample != 1:
                        cut_cmd.append(f"--subsample={self.subsample}")
                    chop_cmd=['alignment-chop-internal','--tree',map_tree]
                    gild_cmd = [
                        'alignment-gild',
                        afile,
                        map_tree,
                        '--max-alignments=500',
                        '--alphabet',
                        alphabet,
                    ]
                    self.exec_pipeline([cut_cmd, chop_cmd, gild_cmd], AUfile)

                html_file = self.outdir / f"{alignment}-AU.html"
                properties = self.character_properties_for_alignment(alignment)
                html_inputs = [afile, AUfile]
                if properties is not None:
                    html_inputs.append(properties)
                if not more_recent_than_all_of(html_file, html_inputs):
                    color_scheme=self.color_scheme_for_alphabet(alphabet, contrast=False)
                    color_scheme += "+fade+fade+fade+fade"
                    self.draw_alignment(
                        afile,
                        ruler=True,
                        AU=AUfile,
                        color_scheme=color_scheme,
                        outfile=html_file,
                        properties=properties,
                        alphabet=alphabet,
                    )
                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 += '<main id="report" 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 getattr(self, "branch_model_reports", []):
                    output += self.section_branch_models()
            if self.has_positive_selection():
                self.copy_report_assets([Path("table-sort.js"), Path("positive-selection.js")])
                output += '<script defer src="bpy-summarize-assets/table-sort.js"></script>\n'
                output += '<script defer src="bpy-summarize-assets/positive-selection.js"></script>\n'
                output += self.section_positive_selection()
            if self.has_trees():
                output += self.section_phylogeny_distribution()
            if self.has_alignments():
                output += self.section_alignment_distribution()
            output += self.section_ancestral_sequences() # FIXME!
            if self.has_mixing_results():
                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_glossary()
            output += self.section_end()
            print(output, file=index)

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

    def html_header(self,title):
        section = """\
<!DOCTYPE html>
<html lang="en">
  <head>
"""
        section += f"<meta charset=\"utf-8\">\n"
        section += '<meta name="viewport" content="width=device-width, initial-scale=1">\n'
        section += f"<title>BAli-Phy: {html.escape(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;}
      tbody th[scope="row"] {font-weight: normal;}
      #parameters tbody th[scope="row"], #alignment tbody th[scope="row"] {text-align:left;}
      .backlit2 tbody th {padding-left: 0.5em; padding-right: 1em;}

#topbar {
        align-items: center;
        background-color: rgb(201,217,233);
        border-bottom: 1px solid rgb(155,177,199);
        box-shadow: 0 1px 4px rgb(0 0 0 / 18%);
        box-sizing: border-box;
        display: flex;
        gap: 0.75em;
        margin: 0;
        padding: 0.45em 0.75em;
        position: sticky;
        top: 0;
        width: 100%;
        z-index: 100;
}

#topbar #menu {
        display: flex;
        flex: 1 1 auto;
        gap: 0.35em;
        min-width: 0;
        overflow-x: auto;
        padding: 0.1em 0;
        scrollbar-width: thin;
}

#topbar #menu a {
        background-color: rgb(235,242,248);
        border: 1px solid rgb(139,166,192);
        border-radius: 0.35em;
        color: rgb(25,63,97);
        flex: 0 0 auto;
        padding: 0.3em 0.6em;
        text-decoration: none;
        white-space: nowrap;
}

#topbar #menu a:first-child {margin-left: auto;}
#topbar #menu a:hover {background-color: white; border-color: rgb(73,113,150);}
#topbar #menu a:focus-visible {outline: 2px solid rgb(25,79,128); outline-offset: 2px;}

#topbar #path {
        font-weight: bold;
        flex: 0 0 auto;
        white-space: nowrap;
}

      .skip-link {
        background: white;
        left: 0.75em;
        padding: 0.5em;
        position: absolute;
        top: 0;
        transform: translateY(-120%);
        z-index: 200;
      }
      .skip-link:focus {transform: translateY(0);}
      .visually-hidden {
        border:0;
        clip:rect(0 0 0 0);
        height:1px;
        margin:-1px;
        overflow:hidden;
        padding:0;
        position:absolute;
        white-space:nowrap;
        width:1px;
      }

      html {scroll-padding-top: 4em;}
      .content {margin:1em;}
      .backlit td {background: rgb(220,220,220);}

      .report-section {
        border-top: 1px solid rgb(155,177,199);
        clear: both;
        display: flow-root;
        margin-top: 2.5em;
        padding-top: 1.25em;
      }

      h1 {font-size: 150%;}
      h2 {font-size: 130%; margin-top:2.5em; margin-bottom: 1em}
      .report-section > h2 {
        border-left: 0.3rem solid rgb(73,113,150);
        break-after: avoid;
        margin-top: 0;
        padding: 0.2em 0 0.2em 0.7em;
      }
      h3 {font-size: 110%; margin-bottom: 0.2em;}

      ul {margin-top: 0.4em;}

      .center :is(th, td) {text-align: center;}

      *[title] {cursor:help;}
      a[title] {vertical-align:top;color:#00f;font-size:70%; border-bottom:1px dotted;}
      .phylogeny-grid {
        align-items:start;
        display:grid;
        gap:1.5em;
        grid-template-columns:repeat(auto-fit, minmax(min(100%, 20em), 1fr));
      }
      .mixing-overview, .tree-mixing-grid {
        align-items:start;
        display:grid;
        gap:1.5em;
        grid-template-columns:repeat(auto-fit, minmax(min(100%, 26em), 1fr));
      }
      .mixing-overview + .tree-mixing-grid {margin-block-start:2em;}
      .srq-plots {
        align-items:start;
        display:grid;
        gap:1em;
        grid-template-columns:repeat(auto-fit, minmax(min(100%, 15em), 1fr));
      }
      .plot-stack {display:grid; gap:2em;}
      .plot-panel {min-width:0;}
      .plot-panel .svg-image {display:block; height:auto; max-width:100%; width:100%;}
      .srq-plots img {display:block; height:auto; max-width:100%; width:300px;}
      figure.plot-panel {margin:0;}
      figure.plot-panel figcaption {font-weight:bold; margin-bottom:0.5em; text-align:center;}
      .table-scroll {max-width:100%; overflow-x:auto;}
      .table-scroll:focus-visible {outline:2px solid rgb(25,79,128); outline-offset:2px;}
      .table-scroll > table {min-width:max-content;}
      table.backlit2 {border-collapse: collapse;}
      table.backlit2 thead th {
        border-bottom: 1px solid #C0C0C0;
        padding-bottom: 0.35em;
        padding-top: 0.35em;
      }
      table caption {font-weight:bold; margin-bottom:0.5em; text-align:left;}
      table.backlit2 tr:nth-child(odd of :has(> td)) > :is(th, td) {background-color: #F0F0F0; margin:0}
      table.convergence-statistics tbody th {font-weight:bold; text-align:left;}
      table.scalar-variables td {text-align:right;}
      table.scalar-variables td.scalar-na {text-align:center;}
      .scalar-mean, .scalar-sd, .scalar-median, .scalar-bci {font-variant-numeric:tabular-nums;}
      .backlit2 td.scalar-mean {padding-right:0.25em;}
      table.scalar-variables td.scalar-pm {padding-left:0; padding-right:0; text-align:center;}
      table.model td.model-variable {padding-left:0.5em; padding-right:0.5em; text-align:left;}
      .backlit2 td.scalar-sd {padding-left:0.25em;}
      .backlit2 td.scalar-median {padding-right:0.5em;}
      .backlit2 td.scalar-bci {padding-left:0;}
      table.positive-selection-support td:nth-child(n+4) {text-align:right; font-variant-numeric:tabular-nums;}
      table.positive-selection-sites :is(th, td).site-column,
      table.positive-selection-sites td.site-probability,
      table.positive-selection-sites td.site-dnds-mean,
      table.positive-selection-sites td.site-dnds-sd {text-align:right; font-variant-numeric:tabular-nums;}
      table.positive-selection-sites th.posterior-group {text-align:center;}
      .positive-selection-sites button {
        padding:0; border:0; background:transparent; color:inherit;
        font:inherit; font-weight:bold; cursor:pointer;
      }
      .positive-selection-sites button:focus-visible {
        outline:3px solid rgba(37, 99, 235, 0.48); outline-offset:2px;
      }
      table.positive-selection-sites td.site-dnds-mean {padding-right:0.25em;}
      table.positive-selection-sites td.site-dnds-pm {padding-left:0; padding-right:0; text-align:center;}
      table.positive-selection-sites td.site-dnds-sd {padding-left:0.25em;}
      table.tree-files {border-collapse:collapse;}
      table.alignment-summary th.alignment-name {padding-right:3em;}
      table.alignment-summary .alignment-stat {padding-left:0.5em; padding-right:0.5em; text-align:center;}
      tr.scalar-unestimated {font-style:italic;}
      .diagnostic-legend {font-size:90%; margin-bottom:0.5em;}
      .diagnostic-marker {display:inline-block; font-weight:bold; margin-right:0.25em; min-width:0.75em; text-align:center;}
      .diagnostic-bad {color:rgb(176,0,32); font-weight:600;}
      .diagnostic-caution {color:rgb(128,81,0); font-weight:600;}
      .glossary-link {color:inherit; text-decoration:underline dotted; text-underline-offset:0.15em;}
      dl.glossary {max-width:70em;}
      dl.glossary dt {font-weight:bold;}
      dl.glossary dd {margin-bottom:0.75em;}

      .modelname {font-weight: bold; color: blue}
      code.model-expression, code.model-identifier {
        font-family:ui-monospace,SFMono-Regular,Consolas,"Liberation Mono",monospace;
      }
      table.model {margin-left: 2em}
      table.model :is(th, td) {vertical-align: top}
      table.models :is(th, td) {vertical-align: top}

      @page {margin:1.5cm;}
      @media print {
        #topbar, .skip-link {display:none;}
        .content {margin:0;}
        .phylogeny-grid {grid-template-columns:repeat(3, minmax(0, 1fr));}
        .mixing-overview, .tree-mixing-grid {grid-template-columns:minmax(0, 1fr) minmax(0, 2fr);}
        .srq-plots {grid-template-columns:repeat(3, minmax(0, 1fr));}
        .table-scroll {overflow:visible;}
        .table-scroll > table {min-width:0;}
        figure, tr, img, object {break-inside:avoid;}
        thead {display:table-header-group;}
      }

    </style>
</head>
<body>
  <a class="skip-link" href="#report">Skip to report content</a>
"""
        return section

    def topbar(self):
        sections = [
            ("data", "Overview", self.has_model()),
            ("parameters", "Parameters", self.has_parameters()),
            ("positive-selection", "Positive selection", self.has_positive_selection()),
            ("branch-models", "Branch models", bool(getattr(self, "branch_model_reports", []))),
            ("topology", "Trees", self.has_trees()),
            ("alignment", "Alignments", self.has_alignments()),
            ("mixing", "Convergence &amp; mixing", self.has_mixing_results()),
            ("analysis", "Run details", True),
            ("models", "Model &amp; priors", self.has_model()),
            ("code", "Code", self.has_code()),
        ]
        links = "\n".join(f'      <a href="#{name}">{label}</a>' for name, label, include in sections if include)
        return f"""\
  <nav id="topbar" aria-label="Report sections">
    <span id="path">Sections:</span>
    <div id="menu">
{links}
    </div>
  </nav>
"""

    def section_data_and_model(self):
        section = """\
<section id="data" class="report-section">
<h2>Data and model</h2>
<div class="table-scroll" role="region" aria-label="Data and model table" tabindex="0">
<table class="backlit2">
  <thead>
    <tr><th scope="col">Partition</th><th scope="col">Sequences</th><th scope="col">Lengths</th><th scope="col">Alphabet</th><th scope="col">Substitution&nbsp;Model</th><th scope="col">Indel&nbsp;Model</th><th scope="col">Scale&nbsp;Model</th></tr>
  </thead>
  <tbody>
"""
        for i in range(self.n_partitions()):
            section += f'<tr id="P{i+1}">\n'
            section += f'  <th scope="row">{i+1}</th>\n'
            section += f'  <td>{html.escape(str(self.get_input_files()[i]))}</td>\n'

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

            alphabet = self.get_alphabets()[i]
            section += f'<td>{html.escape(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} {print_model_string(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} {print_model_string(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>'
                section += f'<td>{link} {print_model_string(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 += '</tbody>\n'
        section += '</table>\n'
        section += '</div>\n'
        section += '</section>\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):
        with inpath.open(encoding='utf-8') as infile:
            code = html.escape(infile.read())
        title = html.escape(inpath.name)

        with outpath.open('w',encoding='utf-8') as outfile:
            print(f"""\
<!DOCTYPE html>
<html lang="en">
  <head>
    <meta charset="utf-8">
    <meta name="viewport" content="width=device-width, initial-scale=1">
    <title>{title}</title>
    <style>
      body {{font-family:system-ui,sans-serif; margin:0;}}
      main {{margin:1em;}}
      h1 {{font-size:1.35em;}}
      pre {{background:#fafafa; border:1px solid #ddd; overflow:auto; padding:1em; tab-size:4;}}
      code {{font-family:ui-monospace,SFMono-Regular,Consolas,"Liberation Mono",monospace;}}
    </style>
    <link rel="stylesheet" href="bpy-summarize-assets/atom-one-light.min.css">
    <script defer src="bpy-summarize-assets/highlight.min.js"></script>
    <script defer src="bpy-summarize-assets/haskell.min.js"></script>
    <script>document.addEventListener("DOMContentLoaded", () => hljs.highlightAll());</script>
  </head>
  <body>
    <main>
      <h1>{title}</h1>
      <pre><code class="language-haskell">{code}</code></pre>
    </main>
  </body>
</html>
""", file=outfile)

    # Copy selected installed web assets while retaining their versioned directory structure.
    def copy_report_assets(self, relative_paths):
        source_assets = self.libexecdir / "bpy-summarize-assets"
        report_assets = self.outdir / "bpy-summarize-assets"
        for relative_path in relative_paths:
            destination = report_assets / relative_path
            destination.parent.mkdir(parents=True, exist_ok=True)
            shutil.copyfile(source_assets / relative_path, destination)

    def copy_code(self):
        self.copy_report_assets([
            Path("highlight.min.js"),
            Path("haskell.min.js"),
            Path("atom-one-light.min.css"),
            Path("LICENSE"),
        ])

        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 = """\
<section id="code" class="report-section">
<h2>Model Code</h2>
    <ul>
"""
        for file_path in self.get_code():
            url = file_path.with_suffix(".html").name
            href = html.escape(url, quote=True)
            label = html.escape(file_path.name)
            section += f'        <li><a href="{href}">{label}</a></li>\n'
        section += "    </ul>\n"
        section += '</section>\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 = """\
<section id="parameters" class="report-section">
<h2>Scalar variables</h2>
  <p class="diagnostic-legend"><strong>Diagnostic flags:</strong>
    <span class="diagnostic-bad"><span class="diagnostic-marker" aria-hidden="true">!</span>concerning</span>;
    <span class="diagnostic-caution"><span class="diagnostic-marker" aria-hidden="true">&#9651;</span>caution</span>.
  </p>
  <div class="table-scroll" role="region" aria-label="Scalar variables table" tabindex="0">
    <table class="backlit2 scalar-variables">
      <thead>
        <tr><th scope="col">Statistic</th><th scope="colgroup" colspan="3">Mean &plusmn; SD</th><th scope="colgroup" colspan="2">Median (95% BCI)</th><th scope="col"><a class="glossary-link" href="#glossary-act"><abbr title="Autocorrelation time">ACT</abbr></a></th><th scope="col"><a class="glossary-link" href="#glossary-ess"><abbr title="Effective sample size">ESS</abbr></a></th><th scope="col">burnin</th><th scope="col"><a class="glossary-link" href="#glossary-psrf-ci80"><abbr title="Potential scale reduction factor based on the 80% credible interval">PSRF-80%CI</abbr></a></th><th scope="col"><a class="glossary-link" href="#glossary-psrf-rcf"><abbr title="Potential scale reduction factor for discrete distributions">PSRF-RCF</abbr></a></th></tr>
      </thead>
      <tbody>
"""
        bad_marker = (
            '<span class="diagnostic-marker" aria-hidden="true">!</span>'
            '<span class="visually-hidden">Concerning: </span>'
        )
        caution_marker = (
            '<span class="diagnostic-marker" aria-hidden="true">&#9651;</span>'
            '<span class="visually-hidden">Caution: </span>'
        )
        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_class = "" if var in self.CI_low else ' class="scalar-unestimated"'
            section += f'<tr{row_class}>\n'
            section += f'  <th scope="row">{html.escape(var)}</th>\n'
            if var in self.constants:
                section += f'  <td class="scalar-mean">{html.escape(self.median[var])}</td>\n'
                section += '  <td class="scalar-pm"></td>\n'
                section += '  <td class="scalar-sd scalar-na">&mdash;</td>\n'
                section += f'  <td class="scalar-median">{html.escape(self.median[var])}</td>\n'
                section += '  <td class="scalar-bci scalar-na">&mdash;</td>\n'
                section += '  <td class="scalar-na">&mdash;</td>\n'
                section += '  <td class="scalar-na">&mdash;</td>\n'
                section += '  <td class="scalar-na">&mdash;</td>\n'
                section += '  <td class="scalar-na">&mdash;</td>\n'
                section += '  <td class="scalar-na">&mdash;</td>\n'
                section += '</tr>\n'
                continue

            if var in self.mean:
                section += f'  <td class="scalar-mean">{html.escape(self.mean[var])}</td>\n'
                section += '  <td class="scalar-pm">&plusmn;</td>\n'
                section += f'  <td class="scalar-sd">{html.escape(self.stddev[var])}</td>\n'
            else:
                section += '  <td class="scalar-mean"></td>\n'
                section += '  <td class="scalar-pm"></td>\n'
                section += '  <td class="scalar-sd"></td>\n'

            if var in self.CI_low:
                section += f'  <td class="scalar-median">{html.escape(self.median[var])}</td>\n'
                section += (
                    f'  <td class="scalar-bci">({html.escape(self.CI_low[var])}, '
                    f'{html.escape(self.CI_high[var])})</td>\n'
                )

                section += f'  <td>{self.ACT[var]}</td>\n'

                if self.ESS[var] < 100:
                    ESS_class = ' class="diagnostic-bad"'
                    ESS_marker = bad_marker
                elif self.ESS[var] < 300:
                    ESS_class = ' class="diagnostic-caution"'
                    ESS_marker = caution_marker
                else:
                    ESS_class = ''
                    ESS_marker = ''
                section += f'  <td{ESS_class}>{ESS_marker}{self.ESS[var]}</td>\n'

                if self.Burnin[var] == "Not Converged!":
                    burnin_class = ' class="diagnostic-bad"'
                    burnin_value = (
                        '<span class="diagnostic-marker" aria-hidden="true">!</span>'
                        '<span class="visually-hidden">Not converged</span>'
                    )
                else:
                    burnin_class = ''
                    burnin_value = html.escape(str(self.Burnin[var]))
                section += f'  <td{burnin_class}>{burnin_value}</td>\n'

                if var in self.PSRF_CI80:
                    if self.PSRF_CI80[var] >= 1.2:
                        CI80_class = ' class="diagnostic-bad"'
                        CI80_marker = bad_marker
                    elif self.PSRF_CI80[var] >= 1.05:
                        CI80_class = ' class="diagnostic-caution"'
                        CI80_marker = caution_marker
                    else:
                        CI80_class = ''
                        CI80_marker = ''
                    section += f'  <td{CI80_class}>{CI80_marker}{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_class = ' class="diagnostic-bad"'
                        RCF_marker = bad_marker
                    elif self.PSRF_RCF[var] >= 1.05:
                        RCF_class = ' class="diagnostic-caution"'
                        RCF_marker = caution_marker
                    else:
                        RCF_class = ''
                        RCF_marker = ''
                    section += f'  <td{RCF_class}>{RCF_marker}{self.PSRF_RCF[var]}</td>\n'
                else:
                    section += '  <td>NA</td>\n'

            else:
                section += f'  <td class="scalar-median">{html.escape(self.median[var])}</td>\n'
                section += '  <td class="scalar-bci"></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 += '      </tbody>\n'
        section += '</table>\n'
        section += '</div>\n'
        section += '</section>\n'

        return section

    # Present model support and the highest-ranked sites without duplicating site selection logic.
    # Every displayed row comes from the complete TSV table generated by character-properties.
    def section_positive_selection(self):
        section = """\
<section id="positive-selection" class="report-section">
<h2>Positive selection</h2>
"""
        reports = getattr(self, "positive_selection_reports", [])
        if any(report["conditioned"] and 0 < report["retained_samples"] < report["total_samples"]
               for report in reports):
            section += """\
<p>The overall posterior allows positive selection to be absent; the posterior with selection
assumes it is included in the model.</p>
"""
        model_support = self.positive_selection_model_support()
        selector_models = self.positive_selection_selector_models()
        if model_support:
            section += """\
<h3>Support for positive selection in the model</h3>
<div class="table-scroll" role="region" aria-label="Positive-selection model support table" tabindex="0">
<table class="backlit2 positive-selection-support">
  <thead>
    <tr><th scope="col">Statistic</th><th scope="col">Substitution model</th><th scope="col">Partitions</th><th scope="col">Posterior probability</th><th scope="col">Log odds</th></tr>
  </thead>
  <tbody>
"""
            for support in model_support:
                statistic = html.escape(support["statistic"])
                smodel = support["smodel"]
                partition_links = ", ".join(
                    f'<a href="#positive-selection-P{partition}">P{partition}</a>'
                    for partition in support["partitions"]
                )
                section += f"""\
  <tr>
    <th scope="row">{statistic}</th>
    <td><a href="#S{smodel}">S{smodel}</a></td>
    <td>{partition_links}</td>
    <td>{html.escape(support["posterior_probability"])}</td>
    <td>{html.escape(support["log_odds"])}</td>
  </tr>
"""
            section += "  </tbody>\n</table>\n</div>\n"

        anchored_partitions = set()
        for report in reports:
            partition = report["partition"]
            property_name = html.escape(report["property"])
            heading_id = ""
            if partition not in anchored_partitions:
                heading_id = f' id="positive-selection-P{partition}"'
                anchored_partitions.add(partition)
            smodel_index = self.get_smodel_indices()[partition - 1]
            smodel = ""
            if smodel_index is not None:
                smodel = f' — <a href="#S{smodel_index + 1}">S{smodel_index + 1}</a>'
            section += f"<h3{heading_id}>Partition {partition}{smodel} — {property_name}</h3>\n"
            alignment_name = html.escape(report["alignment"].name)
            section += (
                f'<p>Display alignment: <a href="{alignment_name}">{alignment_name}</a>. '
                'Column numbers refer to this alignment.</p>\n'
            )

            rows = report["rows"]
            threshold = report["threshold"]
            # Keep every row in the download; filter only the HTML table, using either posterior.
            rows = [
                row for row in rows
                if float(row["model-averaged-probability"]) > threshold
                or (row["conditioned-probability"]
                    and float(row["conditioned-probability"]) > threshold)
            ]
            threshold_text = f"{threshold:g}"
            model_averaged_count = report["model_averaged_count"]
            conditioned_count = report["conditioned_count"]
            retained_samples = report["retained_samples"]
            total_samples = report["total_samples"]
            has_distinct_selection_posterior = (
                report["conditioned"] and 0 < retained_samples < total_samples
            )
            fixed_selection = None
            if (selector_models is not None and smodel_index is not None
                    and smodel_index + 1 not in selector_models and report["conditioned"]
                    and total_samples > 0):
                if retained_samples == total_samples:
                    fixed_selection = True
                elif retained_samples == 0:
                    fixed_selection = False

            section += "<p>"
            if fixed_selection is True:
                section += "Positive selection is included in this model.<br>\n"
            elif fixed_selection is False:
                section += "Positive selection is not included in this model.<br>\n"

            if report["conditioned"] and total_samples > 0 and retained_samples == total_samples:
                count_text = str(model_averaged_count)
            else:
                count_text = f"{model_averaged_count} overall"
                if has_distinct_selection_posterior:
                    count_text += f", {conditioned_count} with selection"
            section += (
                f"Columns with Pr(dN/dS &gt; 1) &gt; {threshold_text}: {count_text}</p>\n"
            )

            if rows:
                section += '<div class="positive-selection-table">\n'
                if report["conditioned"] and retained_samples > 0:
                    section += (
                        '<p>Each column uses the letter with the highest probability of positive selection, '
                        'conditional on selection being present in the model.'
                    )
                    if has_distinct_selection_posterior:
                        section += ' Both summaries refer to that same letter.'
                    section += '</p>\n'
                else:
                    section += (
                        '<p>Each column uses the letter with the highest overall probability '
                        'of positive selection.</p>\n'
                    )
                if len(rows) > 15:
                    section += (
                        f'<p class="site-count">Showing 15 of {len(rows)} selected columns.</p>\n'
                        '<p class="site-expand" hidden><label><input type="checkbox"> '
                        'Show all selected columns</label></p>\n'
                    )

                has_model_averaged_dnds = any(row["model-averaged-dNdS-mean"] for row in rows)
                has_conditioned = has_distinct_selection_posterior
                has_conditioned_dnds = has_conditioned and any(row["conditioned-dNdS-mean"] for row in rows)
                model_averaged_span = 1 + 3 * has_model_averaged_dnds
                conditioned_span = 1 + 3 * has_conditioned_dnds
                section += (
                    f'<div class="table-scroll" role="region" '
                    f'aria-label="Partition {partition} positive-selection sites table" tabindex="0">\n'
                )
                section += '<table class="backlit2 positive-selection-sites">\n  <thead>\n  <tr>\n'
                section += (
                    '    <th scope="col" rowspan="2" data-sort-column="0" '
                    'data-default-direction="1">Column</th>'
                    '<th scope="col" rowspan="2">Codon</th><th scope="col" rowspan="2">AA</th>'
                )
                section += (
                    f'<th scope="colgroup" class="posterior-group" colspan="{model_averaged_span}">'
                    "Overall posterior</th>"
                )
                if has_conditioned:
                    section += (
                        f'<th scope="colgroup" class="posterior-group" colspan="{conditioned_span}">'
                        "Posterior with selection</th>"
                    )
                section += '<th scope="col" rowspan="2">Source letter</th>\n  </tr>\n  <tr>\n'
                section += (
                    '    <th scope="col" data-sort-column="3" data-default-direction="-1">'
                    'Pr(dN/dS&gt;1)</th>'
                )
                if has_model_averaged_dnds:
                    section += '<th scope="colgroup" colspan="3">Posterior dN/dS<br>mean &plusmn; SD</th>'
                if has_conditioned:
                    section += (
                        f'<th scope="col" data-sort-column="{3 + model_averaged_span}" '
                        'data-default-direction="-1">Pr(dN/dS&gt;1)</th>'
                    )
                    if has_conditioned_dnds:
                        section += (
                            '<th scope="colgroup" colspan="3">Posterior dN/dS<br>'
                            'mean &plusmn; SD</th>'
                        )
                section += "\n  </tr>\n  </thead>\n  <tbody>\n"
                for rank, row in enumerate(rows):
                    model_averaged_probability = f'{float(row["model-averaged-probability"]):.3f}'
                    source = f'{row["sequence"]}:{row["sequence-character"]}'
                    hidden = ' hidden' if rank >= 15 else ''
                    section += f'  <tr data-report-row="{rank}"{hidden}>\n'
                    section += (
                        f'    <th scope="row" class="site-column" '
                        f'data-sort-value="{html.escape(row["column"], quote=True)}">'
                        f'{html.escape(row["column"])}</th>\n'
                    )
                    section += f'    <td>{html.escape(row["symbol"])}</td>\n'
                    section += f'    <td>{html.escape(row["translation"] or "-")}</td>\n'
                    section += (
                        f'    <td class="site-probability" '
                        f'data-sort-value="{float(row["model-averaged-probability"])}">'
                        f'{model_averaged_probability}</td>\n'
                    )
                    if has_model_averaged_dnds:
                        model_averaged_dnds_mean = f'{float(row["model-averaged-dNdS-mean"]):.3f}'
                        model_averaged_dnds_sd = f'{float(row["model-averaged-dNdS-sd"]):.3f}'
                        section += f'    <td class="site-dnds-mean">{model_averaged_dnds_mean}</td>\n'
                        section += '    <td class="site-dnds-pm">&plusmn;</td>\n'
                        section += f'    <td class="site-dnds-sd">{model_averaged_dnds_sd}</td>\n'
                    if has_conditioned:
                        conditioned_probability = f'{float(row["conditioned-probability"]):.3f}'
                        section += (
                            f'    <td class="site-probability" '
                            f'data-sort-value="{float(row["conditioned-probability"])}">'
                            f'{conditioned_probability}</td>\n'
                        )
                        if has_conditioned_dnds:
                            conditioned_dnds_mean = f'{float(row["conditioned-dNdS-mean"]):.3f}'
                            conditioned_dnds_sd = f'{float(row["conditioned-dNdS-sd"]):.3f}'
                            section += f'    <td class="site-dnds-mean">{conditioned_dnds_mean}</td>\n'
                            section += '    <td class="site-dnds-pm">&plusmn;</td>\n'
                            section += f'    <td class="site-dnds-sd">{conditioned_dnds_sd}</td>\n'
                    section += f'    <td>{html.escape(source)}</td>\n'
                    section += "  </tr>\n"
                section += "  </tbody>\n</table>\n</div>\n</div>\n"

            links = [f'<a href="{html.escape(report["filename"].name)}">All rows (TSV)</a>']
            alignment_viewer = report["alignment"].with_suffix(".html")
            if alignment_viewer.exists():
                links.append(f'<a href="{alignment_viewer.name}">Alignment viewer</a>')
            section += "<p>" + " · ".join(links) + "</p>\n"
        section += '</section>\n'
        return section

    def section_phylogeny_distribution(self):
        section = '<section id="topology" class="report-section">\n'
        section += '<h2>Phylogeny Distribution</h2>\n'
        section += '<div class="phylogeny-grid">\n'
        section += '<figure class="plot-panel">\n'
        section += '<figcaption>50% consensus tree</figcaption>\n'
        if path.exists(self.outdir / 'c50-tree.svg'):
            section += self.html_svg('c50-tree.svg', '50% consensus tree')
        else:
            section += '<p>Tree image not generated.</p>\n'
        section += '</figure>\n'
        section += '<div class="phylogeny-files">\n'
        if self.has_topology_variation():
            section += (
                '<p><strong>Partition support:</strong> <a href="consensus">Summary</a> · '
                '<a href="partitions.bs">Across chains</a></p>\n'
            )
        section += (
            '<div class="table-scroll" role="region" aria-label="Tree files table" tabindex="0">\n'
            '<table class="backlit2 tree-files">\n'
            '<caption>Tree files</caption>\n'
            '<thead>\n'
            '  <tr><th scope="col" rowspan="2">Tree</th>'
            '<th scope="colgroup" colspan="4">Consensus tree</th>'
        )
        if self.subpartitions:
            section += '<th scope="colgroup" colspan="3">Markov-chain tree</th>'
        section += '</tr>\n  <tr><th scope="col">Newick</th><th scope="col">+PP</th>'
        section += '<th scope="col">PDF</th><th scope="col">SVG</th>'
        if self.subpartitions:
            section += '<th scope="col">Newick</th><th scope="col">PDF</th><th scope="col">SVG</th>'
        section += '</tr>\n</thead>\n<tbody>\n'
        for (tree,name) in self.trees:
            tree_url = html.escape(str(tree), quote=True)
            tree_label = html.escape(str(name))
            download_label = html.escape(str(name), quote=True)
            section += f"""\
  <tr>
    <th scope="row">{tree_label}</th>
    <td><a href="{tree_url}.tree" aria-label="Download {download_label} tree in Newick format">Newick</a></td>
    <td><a href="{tree_url}.PP.tree" aria-label="Download {download_label} Newick tree with posterior probabilities">+PP</a></td>
"""
            if path.exists(self.outdir / f"{tree}-tree.pdf"):
                section += (
                    f'    <td><a href="{tree_url}-tree.pdf" aria-label="Download {download_label} '
                    'tree as PDF">PDF</a></td>\n'
                )
            else:
                section += '    <td>&mdash;</td>\n'
            if path.exists(self.outdir / f"{tree}-tree.svg"):
                section += (
                    f'    <td><a href="{tree_url}-tree.svg" aria-label="Download {download_label} '
                    'tree as SVG">SVG</a></td>\n'
                )
            else:
                section += '    <td>&mdash;</td>\n'
            if self.subpartitions:
                if path.exists(self.outdir / f"{tree}.mtree"):
                    section += (
                        f'<td><a href="{tree_url}.mtree" aria-label="Download {download_label} '
                        'Markov-chain tree in Newick format">Newick</a></td>\n'
                    )
                else:
                    section += '<td></td>\n'

                if path.exists(self.outdir / f"{tree}-mctree.pdf"):
                    section += (
                        f'<td><a href="{tree_url}-mctree.pdf" aria-label="Download {download_label} '
                        'Markov-chain tree as PDF">PDF</a></td>\n'
                    )
                else:
                    section += '<td></td>\n'

                if path.exists(self.outdir / f"{tree}-mctree.svg"):
                    section += (
                        f'<td><a href="{tree_url}-mctree.svg" aria-label="Download {download_label} '
                        'Markov-chain tree as SVG">SVG</a></td>\n'
                    )
                else:
                    section += '<td></td>\n'

            section += '  </tr>\n'
        section += '</tbody>\n</table>\n</div>\n'
        section += '</div>\n'
        if self.has_topology_variation():
            section += '<figure class="plot-panel">\n'
            section += '<figcaption>Consensus-tree support levels</figcaption>\n'
            if path.exists(self.outdir / 'c-levels.svg'):
                section += self.html_svg('c-levels.svg', 'Consensus-tree support levels')
            else:
                section += '<p>Support-level plot not generated.</p>\n'
            section += '</figure>\n'
        section += '</div>\n'
        section += '</section>\n'
        return section

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

        section = '<section id="alignment" class="report-section">\n'
        section += '<h2>Alignment Distribution</h2>\n'
        for i in range(self.n_partitions()):
            section += f'<h3>Partition {i+1}</h3>\n'
            section += f"""\
  <div class="table-scroll" role="region" aria-label="Partition {i+1} alignment table" tabindex="0">
  <table class="alignment-summary">
    <thead>
    <tr>
       <th scope="col" class="alignment-name">Alignment</th>
       <th scope="col">FASTA</th>
       <th scope="col">HTML</th>
       <th scope="col"><span aria-hidden="true">Diff</span><span class="visually-hidden">Difference from the WPD alignment</span></th>
       <th scope="col">AU</th>
       <th scope="col" class="alignment-stat"><span aria-hidden="true">Min. %identity</span><span class="visually-hidden">Minimum pairwise sequence identity</span></th>
       <th scope="col" class="alignment-stat"><span aria-hidden="true"># Sites</span><span class="visually-hidden">Number of alignment columns</span></th>
       <th scope="col" class="alignment-stat"><span aria-hidden="true">Constant</span><span class="visually-hidden">Number of invariant columns</span></th>
       <th scope="col" class="alignment-stat"><span aria-hidden="true">Informative</span><span class="visually-hidden">Number of parsimony-informative columns</span></th>
    </tr>
    </thead>
    <tbody>
"""
            for (alignment,alphabet,name) in self.alignments:
                if not alignment.startswith(f'P{i+1}.'):
                    continue
                section +=  '    <tr>\n'
                section += (
                    f'      <th scope="row" class="alignment-name">{html.escape(name)}</th>\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", alphabet)
                section += f'<td class="alignment-stat">{features["min_p_identity"]}%</td>\n'
                section += f'<td class="alignment-stat">{features["length"]}</td>\n'
                section += (
                    f'<td class="alignment-stat">{features["n_const"]} ({features["p_const"]}%)</td>\n'
                )
                section += (
                    f'<td class="alignment-stat">{features["n_inform"]} '
                    f'({features["p_inform"]}%)</td>\n'
                )

                section += '    </tr>\n'
            section += '    </tbody>\n'
            section += '  </table>\n'
            section += '  </div>\n'
        section += '</section>\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

    # Embed an SVG with an accessible name and useful fallback link.
    # Callers may add layout classes without sharing mutable state between plots.
    def html_svg(self,url,description,width=None,height=None,classes=None,extra=None):
        classes = list(classes) if classes else []
        classes.append('svg-image')
        class_ = ' '.join(classes)
        escaped_url = html.escape(url, quote=True)
        escaped_description = html.escape(description, quote=True)
        svg = (
            f'<object data="{escaped_url}" type="image/svg+xml" class="{class_}" role="img" '
            f'aria-label="{escaped_description}"'
        )
        if width:
            svg += f' width="{html.escape(str(width), quote=True)}"'
        if height:
            svg += f' height="{html.escape(str(height), quote=True)}"'
        if extra:
            svg += " " + extra
        svg += f'><a href="{escaped_url}">View {escaped_description}</a></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"""\
<div class="tree-mixing-grid">
  <figure class="plot-panel">
    <figcaption>Projection of RF distances for {MDS_title} (<a href="https://doi.org/10.1080/10635150590946961">Hillis et al 2005</a>)</figcaption>
"""
        if path.exists( self.outdir / MDS_figure ):
            section += self.html_svg(MDS_figure, 'Projection of Robinson-Foulds tree distances')
            if path.exists(self.outdir / f"{MDS_figure_3d}.html"):
                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 += '</figure>\n'
        section += '<div class="plot-stack">\n'
        section += '<figure class="plot-panel">\n'
        section += (
            '<figcaption>Split posterior probabilities across chains '
            '(<a href="https://doi.org/10.1080/10635150600812544">Beiko et al 2006</a>)'
            '</figcaption>\n'
        )
        if path.exists(self.outdir / "convergence1-PP.svg"):
            section += self.html_svg(
                "convergence1-PP.svg", "Split posterior probabilities across chains"
            )
        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>"
        section += '</figure>\n'

        if path.exists(self.outdir / "convergence2-PP.svg"):
            section += '<figure class="plot-panel">\n'
            section += '<figcaption>Split posterior log odds across chains</figcaption>\n'
            section += self.html_svg(
                "convergence2-PP.svg",
                "Split posterior log odds across chains",
            )
            section += '</figure>\n'
        section += "</div>\n"
        section += "</div>\n"
        return section

    def section_mixing(self):
        section = """\
<section id="mixing" class="report-section">
<h2>Convergence and mixing</h2>
  <div class="mixing-overview">
    <div>
        <table class="backlit2 convergence-statistics">
          <caption>Statistics</caption>
          <tbody>
"""
        burnin_before = "NA"
        min_NE = "NA"
        min_NE_partition = None
        asdsf = None
        msdsf = None
        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><th scope="row">scalar burnin</th>'
                f'<td>{html.escape(str(burnin_before))}</td></tr>\n'
            )
            section += (
                '<tr><th scope="row">scalar <a class="glossary-link" href="#glossary-ess">'
                '<abbr title="Effective sample size">ESS</abbr></a></th>'
                f'<td>{html.escape(str(min_NE))}</td></tr>\n'
            )


        if self.has_trees() and self.has_topology_variation():
            min_NE_partition = self.get_value_from_file(self.outdir / 'partitions.bs','min Ne =')
            section += (
                '<tr><th scope="row">minimum topological '
                '<a class="glossary-link" href="#glossary-ess">'
                '<abbr title="Effective sample size">ESS</abbr></a></th>'
                f'<td>{html.escape(str(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:
                asdsf = asdsf.rstrip()
            section += (
                '<tr><th scope="row"><a class="glossary-link" href="#glossary-asdsf">'
                '<abbr title="Average standard deviation of split frequencies">ASDSF</abbr>'
                '</a></th>'
                f'<td>{html.escape(str(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 += (
                '<tr><th scope="row"><a class="glossary-link" href="#glossary-msdsf">'
                '<abbr title="Maximum standard deviation of split frequencies">MSDSF</abbr>'
                '</a></th>'
                f'<td>{html.escape(str(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 += (
                '<tr><th scope="row"><a class="glossary-link" href="#glossary-psrf-ci80">'
                '<abbr title="Potential scale reduction factor based on the 80% credible interval">'
                f'PSRF-80%CI</abbr></a></th><td>{html.escape(str(psrf_80))}</td></tr>\n'
            )
            section += (
                '<tr><th scope="row"><a class="glossary-link" href="#glossary-psrf-rcf">'
                '<abbr title="Potential scale reduction factor for discrete distributions">'
                f'PSRF-RCF</abbr></a></th><td>{html.escape(str(psrf_rcf))}</td></tr>\n'
            )

        section += '          </tbody>\n'
        section += '</table>\n'
        section += '</div>\n'

        if self.has_trees() and self.has_topology_variation():
            section += '    <div class="srq-plots">\n'
            srq_plots = [
                ("c50.SRQ.png", "50% consensus-tree split support", "support of the 50% consensus tree"),
                ("partitions.SRQ.png", "All sampled splits", "support of each sampled split"),
            ]
            for filename, caption, description in srq_plots:
                section += '      <figure class="plot-panel">\n'
                section += f'        <figcaption>{caption}</figcaption>\n'
                if path.exists(self.outdir / filename):
                    section += f'        <img src="{filename}" alt="SRQ plot for {description}.">\n'
                else:
                    section += '        <p>SRQ plot not generated.</p>\n'
                section += '      </figure>\n'
            section += '    </div>\n'
        section += '  </div>\n'

        if self.has_trees() and self.has_topology_variation():
            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() and self.has_topology_variation():
            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}")
        section += '</section>\n'
        return section

    # Define only the diagnostic abbreviations that can occur in this report.
    # Inline diagnostic links keep the main report terse while making details easy to find.
    def section_glossary(self):
        has_topology_diagnostics = self.has_trees() and self.has_topology_variation()
        if not (self.has_parameters() or has_topology_diagnostics):
            return ""

        section = """\
<section id="glossary" class="report-section">
<h2>Glossary</h2>
<dl class="glossary">
"""
        if self.has_parameters():
            section += """\
  <dt id="glossary-act"><abbr title="Autocorrelation time">ACT</abbr></dt>
  <dd>Autocorrelation time measures serial dependence. The effective sample size is the retained sample count divided by ACT.</dd>
"""
        if self.has_parameters() or has_topology_diagnostics:
            section += """\
  <dt id="glossary-ess"><abbr title="Effective sample size">ESS</abbr></dt>
  <dd>Effective sample size approximates the number of independent samples that would give the same precision as the correlated MCMC samples. Values below 100 are flagged as concerning, and values from 100 through 299 warrant caution; obtaining more samples should increase ESS.</dd>
"""
        if self.has_parameters():
            section += """\
  <dt id="glossary-psrf-ci80"><abbr title="Potential scale reduction factor based on the 80% credible interval">PSRF-80%CI</abbr></dt>
  <dd>The width of the pooled 80% credible interval divided by the average within-chain width. Values from 1.05 through 1.19 warrant caution, and values of 1.2 or greater are concerning. Values near 1 indicate similar widths across chains.</dd>
  <dt id="glossary-psrf-rcf"><abbr title="Potential scale reduction factor for discrete distributions">PSRF-RCF</abbr></dt>
  <dd>For each chain's 80% credible interval, this compares its probability in that chain with its probability in the pooled samples, then averages the ratios across chains. It uses the same warning thresholds as PSRF-80%CI.</dd>
"""
        if has_topology_diagnostics:
            section += """\
  <dt id="glossary-asdsf"><abbr title="Average standard deviation of split frequencies">ASDSF</abbr></dt>
  <dd>The average across splits of the standard deviation, across chains, of each split's posterior probability. Values below 0.01 are commonly considered acceptable; larger values suggest obtaining more samples or comparing additional chains.</dd>
  <dt id="glossary-msdsf"><abbr title="Maximum standard deviation of split frequencies">MSDSF</abbr></dt>
  <dd>The largest standard deviation across chains of any split's posterior probability. A large value identifies at least one split whose estimated support still differs substantially among chains.</dd>
"""
        section += "</dl>\n</section>\n"
        return section

    def section_analysis(self):
        section = '<section id="analysis" class="report-section">\n'

        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 += '<h2>Run details</h2>\n'
        if something_common:
            section += "<p>"
        if all_same(commands):
            section += f"<b>command line</b>: {html.escape(str(commands[0]))}<br>\n"
        if all_same(parent_dirs):
            section += f"<b>directory</b>: {html.escape(str(parent_dirs[0]))}<br>\n"
        if all_same(versions):
            section += f"<b>version</b>: {html.escape(str(versions[0]))}\n"
        if something_common:
            section += "</p>"
        #section += '<table style="width:100%">'."\n"

        #section += '<table class="backlit2 center" style="width:100%">'."\n"
        section += '<div class="table-scroll" role="region" aria-label="Run details table" tabindex="0">\n'
        section += '<table class="backlit2 center">\n'
        section += '<thead>\n<tr><th scope="col">chain #</th>'
        if not all_same(versions):
            section += '<th scope="col">version</th>'
        section += '<th scope="col">burnin</th><th scope="col">subsample</th><th scope="col">samples</th>'
        if not all_same(commands):
            section += '<th scope="col">command line</th>'
        section += '<th scope="col">subdirectory</th>'
        if not all_same(parent_dirs):
            section += '<th scope="col">directory</th>'
        section += '</tr>\n</thead>\n<tbody>\n'

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

            section += f'  <th scope="row">{i+1}</th>\n'
            if not all_same(versions):
                section += f"  <td>{html.escape(str(versions[i]))}</td>\n"

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

            n_iterations = self.run(i).n_iterations()
            last_iteration = min(n_iterations, self.until) if self.until is not None else n_iterations
            eligible_samples = max(0, last_iteration - self.burnin)
            # Selection first caps the input, then skips burn-in, and finally keeps every nth sample.
            retained_samples = (eligible_samples + self.subsample - 1) // self.subsample
            print(
                f"RUN {i+1}: directory = {self.run(i).get_dir()}     iterations = {n_iterations}     "
                f"burnin = {self.burnin}"
            )
            section += f"  <td>{retained_samples}</td>\n"

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

            section += "</tr>\n"

        section += "</tbody>\n</table>\n</div>\n"
        section += '</section>\n'
        return section

    def section_model_and_priors(self):
        section = '<section id="models" class="report-section">\n'
        section += '<h2>Model and priors</h2>\n'

        section += '<h3 id="tree">Tree (+priors)</h3>\n'
#    $section .= "<table class=\"backlit2\">\n";
        section += '<table class="tree-priors" role="presentation">\n<tbody>\n'
        if self.get_topology_prior():
            section += (
                '<tr><td class="modelname">topology</td>'
                f'<td>{print_model_string(self.get_topology_prior())}</td></tr>'
            )
        if self.get_branch_length_prior():
            section += (
                '<tr><td class="modelname">branch lengths</td>'
                f'<td>{print_model_string(self.get_branch_length_prior())}</td></tr>'
            )
        if self.get_tree_prior():
            section += (
                '<tr><td class="modelname">tree</td>'
                f'<td>{print_model_string(self.get_tree_prior())}</td></tr>'
            )
        if self.get_subst_rates():
            section += (
                '<tr><td class="modelname">subst rates</td>'
                f'<td>{print_model_string(self.get_subst_rates())}</td></tr>'
            )
        if self.get_indel_rates():
            section += (
                '<tr><td class="modelname">indel rates</td>'
                f'<td>{print_model_string(self.get_indel_rates())}</td></tr>'
            )
        section += '</tbody>\n</table>\n'

        section += '<h3 id="smodel">Substitution model (+priors)</h3>\n'
        section += print_models("S",self.get_smodels())

        section += '<h3 id="imodel">Indel model (+priors)</h3>\n'
        section += print_models("I",self.get_imodels())

        if (self.get_scale_models()):
            section += '<h3 id="scales">Scales (+priors)</h3>\n'
            section += print_models("scale",self.get_scale_models());

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

    def section_end(self):
        return """\
     </main>
   </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   bpy-summarize analysis-dir-1/ analysis-dir-2/\n   bpy-summarize trees1.trees trees2.trees\n   bpy-summarize log1.log log2.log",
                                     formatter_class=argparse.RawTextHelpFormatter)

    parser.add_argument("mcmc_outputs", default=['.'], help="Subdirectories with MCMC runs",nargs='*')
    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 NUMth stored sample")
    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("--sub-partitions",default=False,action='store_true')
    parser.add_argument("--node-bubbles",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)
