"""
Appendix code: computational comparison of selected APIs with APIs present in
FDA-approved amorphous solid dispersion (ASD) products.

The script performs the following analyses using one consistent molecular-
structure standardisation procedure:

1. Conversion of input SMILES strings to standardised neutral parent structures.
2. Explicit demonstration of the conversion of naproxen sodium to naproxen.
3. Calculation of nine molecular descriptors using RDKit.
4. Principal component analysis (PCA).
5. Tanimoto-similarity heatmap for the selected API parent structures.
6. Hierarchical clustering based on Tanimoto distance.
7. Molecular-weight versus calculated-LogP plot.
8. Histograms comparing Lipinski and additional molecular descriptors.

Required input workbook
-----------------------
File: "FDA ASD smiles.xlsx"
Sheet 1 (index 0): APIs present in FDA-approved ASD products
Sheet 2 (index 1): selected APIs

Each sheet must contain a name column ("Name" or "Drug name") and a SMILES
column ("SMILES" or "Smiles").

Required Python packages
------------------------
pandas, numpy, matplotlib, seaborn, scipy, scikit-learn, RDKit and openpyxl
"""

from pathlib import Path

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns

from matplotlib.lines import Line2D
from scipy.cluster.hierarchy import dendrogram, linkage
from scipy.spatial.distance import squareform
from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler

from rdkit import Chem, DataStructs
from rdkit.Chem import AllChem, Crippen, Descriptors, Lipinski, rdMolDescriptors
from rdkit.Chem.MolStandardize import rdMolStandardize


# =============================================================================
# 1. FILES, ANALYSIS SETTINGS AND FIGURE STYLE
# =============================================================================

INPUT_FILE = Path("FDA ASD smiles.xlsx")
FDA_SHEET = 0
SELECTED_SHEET = 1

OUTPUT_DIRECTORY = Path("computational_analysis_outputs")
OUTPUT_DIRECTORY.mkdir(parents=True, exist_ok=True)

MORGAN_RADIUS = 2
MORGAN_BITS = 2048

DESCRIPTOR_COLUMNS = [
    "MW",
    "cLogP",
    "HBD",
    "HBA",
    "TPSA",
    "Rotatable_Bonds",
    "Aromatic_Rings",
    "Fraction_Csp3",
    "Ring_Count",
]

# Colour-blind-friendly palette used consistently in every figure.
SELECTED_COLOURS = {
    "Mefenamic acid": "#0072B2",
    "Allopurinol": "#E69F00",
    "Naproxen acid/naproxen sodium": "#CC79A7",
}

FDA_COLOUR = "#B8B8B8"
FDA_EDGE_COLOUR = "#4D4D4D"

sns.set_context("talk")
sns.set_style("white")

plt.rcParams.update(
    {
        "font.family": "Arial",
        "font.size": 16,
        "axes.titlesize": 20,
        "axes.labelsize": 18,
        "xtick.labelsize": 15,
        "ytick.labelsize": 15,
        "legend.fontsize": 14,
        "pdf.fonttype": 42,
        "ps.fonttype": 42,
        "savefig.dpi": 600,
    }
)


def save_figure(fig, filename_stem):
    """Save each figure as a 600-dpi PNG and a vector PDF."""

    fig.savefig(
        OUTPUT_DIRECTORY / f"{filename_stem}.png",
        dpi=600,
        bbox_inches="tight",
    )
    fig.savefig(
        OUTPUT_DIRECTORY / f"{filename_stem}.pdf",
        format="pdf",
        bbox_inches="tight",
    )


# =============================================================================
# 2. LOAD AND CLEAN THE SMILES WORKBOOK
# =============================================================================

def identify_column(dataframe, accepted_names):
    """Return a column matching one of the accepted labels, ignoring case."""

    normalised_columns = {
        str(column).strip().casefold(): column for column in dataframe.columns
    }

    for accepted_name in accepted_names:
        match = normalised_columns.get(accepted_name.casefold())
        if match is not None:
            return match

    raise KeyError(
        f"Could not find any of these columns: {accepted_names}. "
        f"Available columns: {list(dataframe.columns)}"
    )


def prepare_input_table(dataframe, dataset_label):
    """Create a consistent Name/SMILES/Dataset input table."""

    dataframe = dataframe.copy()
    dataframe.columns = dataframe.columns.astype(str).str.strip()

    name_column = identify_column(
        dataframe,
        ["Name", "Drug name", "API name", "API"],
    )
    smiles_column = identify_column(
        dataframe,
        ["SMILES", "Smiles"],
    )

    cleaned = dataframe[[name_column, smiles_column]].copy()
    cleaned.columns = ["Name", "SMILES"]
    cleaned["Name"] = cleaned["Name"].astype(str).str.strip()
    cleaned["SMILES"] = cleaned["SMILES"].astype(str).str.strip()
    cleaned["Dataset"] = dataset_label

    return cleaned


def standardise_selected_name(name):
    """Apply consistent display names to the four experimental API materials."""

    normalised_name = str(name).strip().casefold()

    if normalised_name == "mefenamic acid":
        return "Mefenamic acid"
    if normalised_name == "allopurinol":
        return "Allopurinol"
    if normalised_name in {"naproxen", "naproxen acid"}:
        return "Naproxen acid"
    if normalised_name == "naproxen sodium":
        return "Naproxen sodium"

    return str(name).strip()


fda_raw = prepare_input_table(
    pd.read_excel(INPUT_FILE, sheet_name=FDA_SHEET),
    "FDA-approved ASD API",
)

selected_raw = prepare_input_table(
    pd.read_excel(INPUT_FILE, sheet_name=SELECTED_SHEET),
    "Selected API",
)

selected_raw["Name"] = selected_raw["Name"].apply(standardise_selected_name)

print("Input record counts:\n")
print(
    pd.concat([fda_raw, selected_raw], ignore_index=True)["Dataset"]
    .value_counts()
    .to_string()
)


# =============================================================================
# 3. CONVERT SMILES TO STANDARDISED NEUTRAL PARENT STRUCTURES
# =============================================================================
#
# The same procedure is applied to every API:
#   1. Parse the input SMILES string.
#   2. Retain the largest molecular fragment, removing counterions.
#   3. Neutralise formal charges where chemically possible.
#
# Naproxen acid is already neutral and is therefore unchanged. For naproxen
# sodium, the Na+ counterion is removed and the carboxylate is protonated,
# producing the same neutral parent structure as naproxen acid.

fragment_chooser = rdMolStandardize.LargestFragmentChooser()
uncharger = rdMolStandardize.Uncharger()


def standardise_structure(smiles):
    """Convert a SMILES string into its standardised neutral parent molecule."""

    molecule = Chem.MolFromSmiles(smiles)

    if molecule is None:
        return None

    largest_fragment = fragment_chooser.choose(molecule)
    neutral_parent = uncharger.uncharge(largest_fragment)

    return neutral_parent


def add_standardised_structures(dataframe):
    """Add RDKit molecule objects and canonical parent SMILES strings."""

    dataframe = dataframe.copy()
    dataframe["Mol_standardised"] = dataframe["SMILES"].apply(
        standardise_structure
    )

    invalid = dataframe["Mol_standardised"].isna()

    if invalid.any():
        invalid_records = dataframe.loc[invalid, ["Name", "SMILES"]]
        raise ValueError(
            "The following SMILES strings could not be interpreted:\n"
            + invalid_records.to_string(index=False)
        )

    dataframe["Canonical_SMILES"] = dataframe["Mol_standardised"].apply(
        lambda molecule: Chem.MolToSmiles(
            molecule,
            canonical=True,
            isomericSmiles=True,
        )
    )

    return dataframe


fda = add_standardised_structures(fda_raw)
selected = add_standardised_structures(selected_raw)

print("\nAll molecular structures were successfully standardised.\n")
print("Standardised selected API structures:\n")
print(selected[["Name", "Canonical_SMILES"]].to_string(index=False))


def show_naproxen_sodium_conversion(selected_dataframe):
    """Print each step in the naproxen-sodium parent conversion."""

    sodium_rows = selected_dataframe[
        selected_dataframe["Name"].str.casefold() == "naproxen sodium"
    ]
    acid_rows = selected_dataframe[
        selected_dataframe["Name"].str.casefold() == "naproxen acid"
    ]

    if sodium_rows.empty:
        raise ValueError("Naproxen sodium was not found in the selected-API sheet.")
    if acid_rows.empty:
        raise ValueError("Naproxen acid was not found in the selected-API sheet.")

    sodium_smiles = sodium_rows.iloc[0]["SMILES"]
    original_molecule = Chem.MolFromSmiles(sodium_smiles)
    largest_fragment = fragment_chooser.choose(original_molecule)
    neutral_parent = uncharger.uncharge(largest_fragment)

    original_canonical = Chem.MolToSmiles(
        original_molecule,
        canonical=True,
        isomericSmiles=True,
    )
    fragment_canonical = Chem.MolToSmiles(
        largest_fragment,
        canonical=True,
        isomericSmiles=True,
    )
    parent_canonical = Chem.MolToSmiles(
        neutral_parent,
        canonical=True,
        isomericSmiles=True,
    )
    acid_parent = acid_rows.iloc[0]["Canonical_SMILES"]

    print("\nNaproxen sodium standardisation steps:\n")
    print(f"Input structure:              {original_canonical}")
    print(f"Largest fragment:             {fragment_canonical}")
    print(f"Neutral parent structure:     {parent_canonical}")
    print(f"Naproxen acid parent:         {acid_parent}")
    print(f"Parent structures identical:  {parent_canonical == acid_parent}")


show_naproxen_sodium_conversion(selected)


# Retain one record for each unique parent structure in the reference set.
fda_unique = (
    fda.drop_duplicates(subset="Canonical_SMILES", keep="first")
    .reset_index(drop=True)
    .copy()
)

# Naproxen acid and naproxen sodium become one computational parent entry.
naproxen_parent_smiles = selected.loc[
    selected["Name"].str.contains("naproxen", case=False, na=False),
    "Canonical_SMILES",
].unique()

if len(naproxen_parent_smiles) != 1:
    raise ValueError(
        "Naproxen acid and naproxen sodium did not produce one identical "
        "standardised parent structure. Check the input SMILES strings."
    )

selected_unique = (
    selected.drop_duplicates(subset="Canonical_SMILES", keep="first")
    .reset_index(drop=True)
    .copy()
)

selected_unique["Display_name"] = selected_unique["Name"]
selected_unique.loc[
    selected_unique["Canonical_SMILES"].isin(naproxen_parent_smiles),
    "Display_name",
] = "Naproxen acid/naproxen sodium"

print(f"\nNumber of FDA API records: {len(fda)}")
print(f"Number of unique FDA parent structures: {len(fda_unique)}")
print(f"Number of experimental selected API forms: {len(selected)}")
print(f"Number of unique selected parent structures: {len(selected_unique)}")


# =============================================================================
# 4. CALCULATE MOLECULAR DESCRIPTORS
# =============================================================================

def calculate_descriptors(molecule):
    """Calculate the nine descriptors used for PCA and descriptor plots."""

    return {
        "MW": Descriptors.MolWt(molecule),
        "cLogP": Crippen.MolLogP(molecule),
        "HBD": Lipinski.NumHDonors(molecule),
        "HBA": Lipinski.NumHAcceptors(molecule),
        "TPSA": rdMolDescriptors.CalcTPSA(molecule),
        "Rotatable_Bonds": Lipinski.NumRotatableBonds(molecule),
        "Aromatic_Rings": rdMolDescriptors.CalcNumAromaticRings(molecule),
        "Fraction_Csp3": rdMolDescriptors.CalcFractionCSP3(molecule),
        "Ring_Count": rdMolDescriptors.CalcNumRings(molecule),
    }


def add_descriptors(dataframe):
    """Append the molecular-descriptor columns to a structure table."""

    descriptor_table = pd.DataFrame(
        dataframe["Mol_standardised"].apply(calculate_descriptors).tolist(),
        index=dataframe.index,
    )

    return pd.concat([dataframe, descriptor_table], axis=1)


fda_unique = add_descriptors(fda_unique)
selected_unique = add_descriptors(selected_unique)

print("\nSelected API descriptor values:\n")
print(
    selected_unique[["Display_name"] + DESCRIPTOR_COLUMNS]
    .round(3)
    .to_string(index=False)
)


# =============================================================================
# 5. PRINCIPAL COMPONENT ANALYSIS
# =============================================================================
#
# The FDA-reference APIs define the descriptor scaling and PCA axes. The
# selected API parent structures are transformed using the fitted scaler and
# projected into the resulting PCA space. They are not used to fit the PCA.

scaler = StandardScaler()
fda_scaled = scaler.fit_transform(fda_unique[DESCRIPTOR_COLUMNS])
selected_scaled = scaler.transform(selected_unique[DESCRIPTOR_COLUMNS])

pca = PCA(n_components=2, svd_solver="full")
fda_scores_unoriented = pca.fit_transform(fda_scaled)
selected_scores_unoriented = pca.transform(selected_scaled)

# The mathematical signs of PCA axes are arbitrary. Orient PC1 so its MW
# coefficient is positive and PC2 so its Fraction Csp3 coefficient is positive.
orientation = np.ones(2)

mw_index = DESCRIPTOR_COLUMNS.index("MW")
fraction_csp3_index = DESCRIPTOR_COLUMNS.index("Fraction_Csp3")

if pca.components_[0, mw_index] < 0:
    orientation[0] = -1

if pca.components_[1, fraction_csp3_index] < 0:
    orientation[1] = -1

fda_scores = fda_scores_unoriented * orientation
selected_scores = selected_scores_unoriented * orientation
oriented_components = pca.components_ * orientation[:, np.newaxis]

explained_variance = pca.explained_variance_ratio_ * 100

print("\nFDA-reference explained variance:")
print(f"PC1: {explained_variance[0]:.1f}%")
print(f"PC2: {explained_variance[1]:.1f}%")
print(f"Cumulative variance: {explained_variance.sum():.1f}%")

pca_coefficients = pd.DataFrame(
    oriented_components.T,
    index=DESCRIPTOR_COLUMNS,
    columns=["PC1_coefficient", "PC2_coefficient"],
)

# For standardised variables, percentage contribution can be calculated from
# the squared component coefficients. The values sum to 100% for each PC.
pca_coefficients["PC1_contribution_percent"] = (
    pca_coefficients["PC1_coefficient"] ** 2 * 100
)
pca_coefficients["PC2_contribution_percent"] = (
    pca_coefficients["PC2_coefficient"] ** 2 * 100
)
pca_coefficients["PC1_rank"] = (
    pca_coefficients["PC1_contribution_percent"]
    .rank(ascending=False, method="min")
    .astype(int)
)
pca_coefficients["PC2_rank"] = (
    pca_coefficients["PC2_contribution_percent"]
    .rank(ascending=False, method="min")
    .astype(int)
)

fda_score_table = pd.DataFrame(
    {
        "Name": fda_unique["Name"].to_numpy(),
        "Canonical_SMILES": fda_unique["Canonical_SMILES"].to_numpy(),
        "PC1": fda_scores[:, 0],
        "PC2": fda_scores[:, 1],
    }
)

selected_score_table = pd.DataFrame(
    {
        "Name": selected_unique["Display_name"].to_numpy(),
        "Canonical_SMILES": selected_unique["Canonical_SMILES"].to_numpy(),
        "PC1": selected_scores[:, 0],
        "PC2": selected_scores[:, 1],
    }
)

explained_variance_table = pd.DataFrame(
    {
        "Principal_component": ["PC1", "PC2"],
        "Explained_variance_percent": explained_variance,
    }
)

print("\nPCA coefficients and contributions:\n")
print(pca_coefficients.round(4).to_string())

print("\nSelected API PCA coordinates:\n")
print(selected_score_table[["Name", "PC1", "PC2"]].round(3).to_string(index=False))


def plot_pca():
    """Plot the FDA-reference PCA scores and projected selected APIs."""

    fig, ax = plt.subplots(figsize=(12, 8))

    ax.scatter(
        fda_scores[:, 0],
        fda_scores[:, 1],
        s=70,
        color=FDA_COLOUR,
        edgecolor="none",
        alpha=0.95,
        label="FDA-approved ASD APIs",
    )

    for row_number, row in selected_unique.iterrows():
        display_name = row["Display_name"]
        ax.scatter(
            selected_scores[row_number, 0],
            selected_scores[row_number, 1],
            s=190,
            marker="X",
            color=SELECTED_COLOURS[display_name],
            edgecolor="black",
            linewidth=1.1,
            label=display_name,
            zorder=5,
        )

    ax.set_xlabel(f"PC1 ({explained_variance[0]:.1f}% variance)")
    ax.set_ylabel(f"PC2 ({explained_variance[1]:.1f}% variance)")
    ax.set_title(
        "Principal component analysis of FDA-approved ASD APIs and selected APIs",
        pad=18,
    )

    ax.legend(
        loc="upper center",
        bbox_to_anchor=(0.5, -0.14),
        ncol=2,
        frameon=False,
    )

    sns.despine(ax=ax)
    fig.subplots_adjust(bottom=0.23)
    save_figure(fig, "PCA_FDA_ASD_and_selected_APIs")
    plt.show()


plot_pca()


# =============================================================================
# 6. MORGAN FINGERPRINTS AND TANIMOTO SIMILARITY
# =============================================================================

def calculate_morgan_fingerprint(molecule):
    """Generate a 2048-bit radius-2 Morgan fingerprint."""

    return AllChem.GetMorganFingerprintAsBitVect(
        molecule,
        radius=MORGAN_RADIUS,
        nBits=MORGAN_BITS,
    )


def calculate_tanimoto_matrix(fingerprints):
    """Calculate a symmetric pairwise Tanimoto-similarity matrix."""

    number_of_structures = len(fingerprints)
    similarity_matrix = np.zeros(
        (number_of_structures, number_of_structures),
        dtype=float,
    )

    for row in range(number_of_structures):
        for column in range(number_of_structures):
            similarity_matrix[row, column] = DataStructs.TanimotoSimilarity(
                fingerprints[row],
                fingerprints[column],
            )

    return similarity_matrix


selected_unique["Morgan_FP"] = selected_unique["Mol_standardised"].apply(
    calculate_morgan_fingerprint
)

selected_tanimoto_matrix = calculate_tanimoto_matrix(
    selected_unique["Morgan_FP"].tolist()
)

selected_tanimoto_table = pd.DataFrame(
    selected_tanimoto_matrix,
    index=selected_unique["Display_name"],
    columns=selected_unique["Display_name"],
)

print("\nSelected parent-structure Tanimoto matrix:\n")
print(selected_tanimoto_table.round(3))


def plot_selected_tanimoto_heatmap():
    """Plot the 3 x 3 Tanimoto heatmap of unique selected parent structures."""

    fig, ax = plt.subplots(figsize=(10, 8))

    sns.heatmap(
        selected_tanimoto_table,
        annot=True,
        fmt=".2f",
        cmap="viridis",
        vmin=0,
        vmax=1,
        linewidths=0.75,
        linecolor="white",
        square=True,
        annot_kws={"size": 15, "family": "Arial"},
        cbar_kws={"label": "Tanimoto similarity", "shrink": 0.90},
        ax=ax,
    )

    ax.set_title(
        "Tanimoto similarity of selected API parent structures",
        pad=18,
    )
    ax.set_xlabel("")
    ax.set_ylabel("")

    plt.setp(
        ax.get_xticklabels(),
        rotation=45,
        ha="right",
        rotation_mode="anchor",
    )
    plt.setp(ax.get_yticklabels(), rotation=0)

    colour_bar = ax.collections[0].colorbar
    colour_bar.ax.tick_params(labelsize=15)
    colour_bar.set_label(
        "Tanimoto similarity",
        fontsize=17,
        fontname="Arial",
        labelpad=12,
    )

    fig.tight_layout()
    save_figure(fig, "Tanimoto_heatmap_selected_API_parent_structures")
    plt.show()


plot_selected_tanimoto_heatmap()


# =============================================================================
# 7. HIERARCHICAL CLUSTERING USING TANIMOTO DISTANCE
# =============================================================================

fda_for_clustering = fda_unique.copy()
fda_for_clustering["Plot_name"] = fda_for_clustering["Name"]

selected_for_clustering = selected_unique.copy()
selected_for_clustering["Plot_name"] = selected_for_clustering["Display_name"]

clustering_data = pd.concat(
    [fda_for_clustering, selected_for_clustering],
    ignore_index=True,
)

clustering_data["Morgan_FP"] = clustering_data["Mol_standardised"].apply(
    calculate_morgan_fingerprint
)

clustering_similarity = calculate_tanimoto_matrix(
    clustering_data["Morgan_FP"].tolist()
)

clustering_distance = 1.0 - clustering_similarity
np.fill_diagonal(clustering_distance, 0.0)

condensed_distance = squareform(clustering_distance, checks=False)

linkage_matrix = linkage(
    condensed_distance,
    method="average",
    optimal_ordering=True,
)


def plot_tanimoto_dendrogram():
    """Plot average-linkage clustering based on Tanimoto distance."""

    fig, ax = plt.subplots(figsize=(13, 16))

    dendrogram(
        linkage_matrix,
        labels=clustering_data["Plot_name"].tolist(),
        orientation="right",
        leaf_font_size=14,
        color_threshold=0,
        above_threshold_color="#777777",
        link_color_func=lambda _: "#777777",
        ax=ax,
    )

    for collection in ax.collections:
        collection.set_linewidth(1.0)

    for label in ax.get_yticklabels():
        label.set_fontname("Arial")
        label.set_fontsize(14)

        if label.get_text() in SELECTED_COLOURS:
            label.set_color(SELECTED_COLOURS[label.get_text()])
            label.set_fontweight("bold")

    ax.set_xlabel("Tanimoto distance")
    ax.set_ylabel("")
    ax.set_title(
        "Hierarchical clustering of FDA-approved ASD APIs and selected APIs",
        pad=18,
    )
    ax.set_xlim(left=0)

    sns.despine(ax=ax, left=True)
    fig.tight_layout()
    save_figure(fig, "Tanimoto_dendrogram_FDA_ASD_and_selected_APIs")
    plt.show()


plot_tanimoto_dendrogram()


# =============================================================================
# 8. MOLECULAR WEIGHT VERSUS CALCULATED LOGP
# =============================================================================

def plot_mw_vs_clogp():
    """Plot MW against RDKit-calculated cLogP and show Lipinski cut-offs."""

    fig, ax = plt.subplots(figsize=(13, 8.5))

    ax.scatter(
        fda_unique["MW"],
        fda_unique["cLogP"],
        s=75,
        facecolor="#D9D9D9",
        edgecolor=FDA_EDGE_COLOUR,
        linewidth=1.3,
        label="FDA-approved ASD APIs",
    )

    for _, row in selected_unique.iterrows():
        display_name = row["Display_name"]
        ax.scatter(
            row["MW"],
            row["cLogP"],
            s=220,
            marker="X",
            color=SELECTED_COLOURS[display_name],
            edgecolor="black",
            linewidth=1.1,
            label=display_name,
            zorder=5,
        )

    ax.axvline(
        500,
        color="black",
        linestyle=":",
        linewidth=2.3,
        label=r"Lipinski MW cut-off (500 g mol$^{-1}$)",
    )
    ax.axhline(
        5,
        color="black",
        linestyle="--",
        linewidth=2.3,
        label="Lipinski cLogP cut-off (5)",
    )

    ax.set_xlim(0, 1200)
    ax.set_ylim(-2, 10)
    ax.set_xlabel(r"Molecular weight (g mol$^{-1}$)")
    ax.set_ylabel("Calculated LogP (cLogP)")
    ax.set_title(
        "Molecular weight and calculated LogP of FDA-approved ASD APIs and selected APIs",
        pad=18,
    )

    # Reorder the legend so the two reference lines appear first.
    handles, labels = ax.get_legend_handles_labels()
    desired_order = [4, 5, 0, 1, 2, 3]
    ax.legend(
        [handles[index] for index in desired_order],
        [labels[index] for index in desired_order],
        loc="upper center",
        bbox_to_anchor=(0.5, -0.14),
        ncol=3,
        frameon=False,
    )

    sns.despine(ax=ax)
    fig.subplots_adjust(bottom=0.25)
    save_figure(fig, "MW_vs_cLogP_FDA_ASD_and_selected_APIs")
    plt.show()


plot_mw_vs_clogp()


# =============================================================================
# 9. HISTOGRAM HELPER FUNCTIONS
# =============================================================================

def add_selected_descriptor_lines(ax, descriptor_name):
    """
    Add selected-API values as dashed vertical lines.

    When values coincide, progressively narrower lines are overlaid at the
    same coordinate. This displays each colour without shifting the true value.
    """

    entries_by_value = {}

    for _, row in selected_unique.iterrows():
        value = float(row[descriptor_name])
        entry = (row["Display_name"], SELECTED_COLOURS[row["Display_name"]])
        entries_by_value.setdefault(value, []).append(entry)

    for value, entries in entries_by_value.items():
        if len(entries) == 1:
            line_widths = [2.8]
        else:
            line_widths = np.linspace(6.0, 2.5, len(entries))

        for (display_name, colour), line_width in zip(entries, line_widths):
            ax.axvline(
                value,
                color=colour,
                linestyle="--",
                linewidth=float(line_width),
                zorder=5,
            )


def add_ro5_cutoff(ax, cutoff):
    """Add a black dotted Lipinski rule-of-five cut-off line."""

    ax.axvline(
        cutoff,
        color="black",
        linestyle=":",
        linewidth=2.3,
        zorder=4,
    )

    x_minimum, x_maximum = ax.get_xlim()
    offset = 0.012 * (x_maximum - x_minimum)

    if cutoff > x_minimum + 0.8 * (x_maximum - x_minimum):
        text_x = cutoff - offset
        horizontal_alignment = "right"
    else:
        text_x = cutoff + offset
        horizontal_alignment = "left"

    y_minimum, y_maximum = ax.get_ylim()

    ax.text(
        text_x,
        y_maximum * 0.97,
        "Ro5 cut-off",
        rotation=90,
        ha=horizontal_alignment,
        va="top",
        fontsize=12,
        color="black",
    )


def selected_line_legend_handles():
    """Return consistent legend handles for selected API descriptor lines."""

    return [
        Line2D(
            [0],
            [0],
            color=colour,
            linestyle="--",
            linewidth=2.8,
            label=display_name,
        )
        for display_name, colour in SELECTED_COLOURS.items()
    ]


def draw_histogram(ax, descriptor_name, bins, title, x_label, x_limits=None):
    """Draw one FDA-reference descriptor histogram and selected-API lines."""

    ax.hist(
        fda_unique[descriptor_name],
        bins=bins,
        color="#D9D9D9",
        edgecolor=FDA_EDGE_COLOUR,
        linewidth=1.0,
    )

    add_selected_descriptor_lines(ax, descriptor_name)

    if x_limits is not None:
        ax.set_xlim(*x_limits)

    ax.set_title(title)
    ax.set_xlabel(x_label)
    ax.set_ylabel("Count")
    sns.despine(ax=ax)


# =============================================================================
# 10. LIPINSKI-DESCRIPTOR HISTOGRAMS
# =============================================================================

def plot_lipinski_histograms():
    """Compare MW, cLogP, HBD and HBA distributions."""

    fig, axes = plt.subplots(2, 2, figsize=(16, 11.5))

    draw_histogram(
        axes[0, 0],
        "MW",
        bins=20,
        title="Molecular weight (MW)",
        x_label=r"Molecular weight (g mol$^{-1}$)",
        x_limits=(0, 1200),
    )
    add_ro5_cutoff(axes[0, 0], 500)

    draw_histogram(
        axes[0, 1],
        "cLogP",
        bins=20,
        title="Calculated LogP (cLogP)",
        x_label="Calculated LogP (cLogP)",
        x_limits=(-2, 10),
    )
    add_ro5_cutoff(axes[0, 1], 5)

    draw_histogram(
        axes[1, 0],
        "HBD",
        bins=np.arange(-0.5, 6.5, 1),
        title="Hydrogen-bond donors (HBD)",
        x_label="HBD",
        x_limits=(-0.7, 5.3),
    )
    add_ro5_cutoff(axes[1, 0], 5)

    maximum_hba = int(max(fda_unique["HBA"].max(), 18))
    draw_histogram(
        axes[1, 1],
        "HBA",
        bins=np.arange(-0.5, maximum_hba + 1.5, 1),
        title="Hydrogen-bond acceptors (HBA)",
        x_label="HBA",
        x_limits=(-0.7, 18.5),
    )
    add_ro5_cutoff(axes[1, 1], 10)

    fig.suptitle(
        "Comparison of Lipinski descriptors for FDA-approved ASD APIs and selected APIs",
        fontsize=22,
        y=0.995,
    )

    fig.legend(
        handles=selected_line_legend_handles(),
        loc="lower center",
        bbox_to_anchor=(0.5, 0.015),
        ncol=3,
        frameon=False,
    )

    fig.tight_layout(rect=[0, 0.08, 1, 0.96])
    save_figure(fig, "Lipinski_descriptor_histograms")
    plt.show()


plot_lipinski_histograms()


# =============================================================================
# 11. ADDITIONAL MOLECULAR-DESCRIPTOR HISTOGRAMS
# =============================================================================

def plot_additional_descriptor_histograms():
    """Compare TPSA, rotatable bonds, fraction Csp3 and aromatic rings."""

    fig, axes = plt.subplots(2, 2, figsize=(16, 11.5))

    draw_histogram(
        axes[0, 0],
        "TPSA",
        bins=20,
        title="Topological polar surface area (TPSA)",
        x_label=r"TPSA ($\mathrm{\AA^2}$)",
        x_limits=(0, 280),
    )

    maximum_rotatable = int(max(fda_unique["Rotatable_Bonds"].max(), 18))
    draw_histogram(
        axes[0, 1],
        "Rotatable_Bonds",
        bins=np.arange(-0.5, maximum_rotatable + 1.5, 1),
        title="Rotatable bonds",
        x_label="Number of rotatable bonds",
        x_limits=(-0.7, 18.5),
    )

    draw_histogram(
        axes[1, 0],
        "Fraction_Csp3",
        bins=np.linspace(0, 0.8, 17),
        title="Fraction Csp3",
        x_label="Fraction Csp3",
        x_limits=(-0.04, 0.82),
    )

    maximum_aromatic = int(max(fda_unique["Aromatic_Rings"].max(), 7))
    draw_histogram(
        axes[1, 1],
        "Aromatic_Rings",
        bins=np.arange(-0.5, maximum_aromatic + 1.5, 1),
        title="Aromatic-ring count",
        x_label="Number of aromatic rings",
        x_limits=(-0.5, 7),
    )

    fig.suptitle(
        "Comparison of additional molecular descriptors for FDA-approved ASD APIs and selected APIs",
        fontsize=22,
        y=0.995,
    )

    fig.legend(
        handles=selected_line_legend_handles(),
        loc="lower center",
        bbox_to_anchor=(0.5, 0.015),
        ncol=3,
        frameon=False,
    )

    fig.tight_layout(rect=[0, 0.08, 1, 0.96])
    save_figure(fig, "Additional_molecular_descriptor_histograms")
    plt.show()


plot_additional_descriptor_histograms()


# =============================================================================
# 12. EXPORT NUMERICAL RESULTS
# =============================================================================

fda_descriptor_export = fda_unique[
    ["Name", "Canonical_SMILES"] + DESCRIPTOR_COLUMNS
].copy()

selected_descriptor_export = selected_unique[
    ["Display_name", "Canonical_SMILES"] + DESCRIPTOR_COLUMNS
].copy()
selected_descriptor_export = selected_descriptor_export.rename(
    columns={"Display_name": "Name"}
)

with pd.ExcelWriter(
    OUTPUT_DIRECTORY / "Computational_analysis_results.xlsx"
) as writer:
    fda_descriptor_export.to_excel(
        writer,
        sheet_name="FDA descriptors",
        index=False,
    )
    selected_descriptor_export.to_excel(
        writer,
        sheet_name="Selected descriptors",
        index=False,
    )
    fda_score_table.to_excel(
        writer,
        sheet_name="FDA PCA scores",
        index=False,
    )
    selected_score_table.to_excel(
        writer,
        sheet_name="Selected PCA scores",
        index=False,
    )
    pca_coefficients.reset_index(names="Descriptor").to_excel(
        writer,
        sheet_name="PCA coefficients",
        index=False,
    )
    explained_variance_table.to_excel(
        writer,
        sheet_name="Explained variance",
        index=False,
    )
    selected_tanimoto_table.to_excel(
        writer,
        sheet_name="Selected Tanimoto",
    )

print(
    "\nAnalysis complete. Figures and numerical tables were saved in:\n"
    f"{OUTPUT_DIRECTORY.resolve()}"
)