#!/usr/bin/env python3
"""
Tool: phylo_upgma_distance_matrix.py
Category: Phylogenetics & Evolutionary Biology
Description: Computes pairwise Jukes-Cantor and Hamming genetic distance matrices from
             Multiple Sequence Alignments (MSA) and executes UPGMA hierarchical clustering
             to construct phylogenetic trees in Newick format.
Usage:
    python phylo_upgma_distance_matrix.py -i alignment.fasta
"""

import argparse
import math
import sys

def hamming_distance(seq1: str, seq2: str) -> float:
    """Calculates p-distance (proportion of nucleotide mismatches)."""
    mismatches = 0
    valid_len = 0
    for b1, b2 in zip(seq1.upper(), seq2.upper()):
        if b1 not in ('-', 'N') and b2 not in ('-', 'N'):
            valid_len += 1
            if b1 != b2:
                mismatches += 1
    return mismatches / valid_len if valid_len > 0 else 0.0

def jukes_cantor_distance(p_dist: float) -> float:
    """Calculates evolutionary distance under the Jukes-Cantor (1969) substitution model."""
    if p_dist >= 0.75:
        return 3.0  # Infinite saturation cap
    return -0.75 * math.log(1 - (4.0 / 3.0) * p_dist)

def run_upgma(names, dist_matrix):
    """Executes UPGMA hierarchical clustering to construct a Newick tree string."""
    clusters = {i: names[i] for i in range(len(names))}
    active = list(range(len(names)))
    dists = { (i, j): dist_matrix[i][j] for i in range(len(names)) for j in range(len(names)) }

    while len(active) > 1:
        # Find minimum distance pair
        min_d = float('inf')
        pair = None
        for i in range(len(active)):
            for j in range(i + 1, len(active)):
                u, v = active[i], active[j]
                if dists.get((u, v), float('inf')) < min_d:
                    min_d = dists[(u, v)]
                    pair = (u, v)

        u, v = pair
        new_node = max(clusters.keys()) + 1
        clusters[new_node] = f"({clusters[u]}:{min_d/2:.3f},{clusters[v]}:{min_d/2:.3f})"
        
        # Calculate new cluster distances
        active.remove(u)
        active.remove(v)
        for w in active:
            avg_d = (dists.get((min(u, w), max(u, w)), 0) + dists.get((min(v, w), max(v, w)), 0)) / 2.0
            dists[(min(w, new_node), max(w, new_node))] = avg_d
            
        active.append(new_node)
        
    return clusters[active[0]] + ";"

def main():
    print("=== Phylogenetics: Jukes-Cantor Distance & UPGMA Tree ===")
    
    taxa = ["Human", "Chimp", "Gorilla", "Orangutan"]
    msa = [
        "ATGCGATCGATCGATCGATCGATCGATC",
        "ATGCGATCGATCGATCGATCGATCGATC", # 0 diffs with human
        "ATGCGATCGATCGATAGATCGATCGATC", # 1 diff
        "ATGTGATCGATCGTTAGATCGATCGATC"  # 3 diffs
    ]
    
    n = len(taxa)
    dist_matrix = [[0.0] * n for _ in range(n)]
    
    for i in range(n):
        for j in range(i + 1, n):
            p = hamming_distance(msa[i], msa[j])
            d = jukes_cantor_distance(p)
            dist_matrix[i][j] = d
            dist_matrix[j][i] = d

    print("\nPairwise Jukes-Cantor Distance Matrix:")
    header = "          " + "".join(f"{t:>12}" for t in taxa)
    print(header)
    for i, t in enumerate(taxa):
        row = f"{t:<10}" + "".join(f"{dist_matrix[i][j]:12.4f}" for j in range(n))
        print(row)

    tree = run_upgma(taxa, dist_matrix)
    print(f"\nConstructed UPGMA Newick Tree:")
    print(f"  -> {tree}")

if __name__ == "__main__":
    main()