Fix bugs, update README for beta, clean up build artifacts

- Fix perm_to_dist[5] summing to 4/3 instead of 1 in discretizer.py
- Fix emission_prob normalization operator precedence bug in DTO and KMeans
- Fix psi validation (np.all -> np.any) catching any negative entry
- Fix limited_int isinstance check that always failed
- Remove unused --num-classes-min/max args and dead comments in phlag.py
- Update README: new CLI flags, expanded toy example, hyperparameter table
- Add tensorflow-probability to requirements.txt
- Delete build.sh and setup.py (qqs.pyx replaced by pure Python qqs.py)
- Add TODO.md with remaining issues

Restructure as installable Python package for PyPI

- Move modules into phlag/ package (hmm, qqs, fxy, utils, main)
- Convert imports to relative
- Add pyproject.toml with hatchling build, console script entry point
- Add phlag/__init__.py with version 0.0.2
- Clean discretizer.py dead code, vectorize DTO, rename to fxy.py
- Fix concentration typos in hmm.py, add explicit imports
- Implement validate_parameters in main.py
- Add test/check_output.py for output validation
- Update README: pip install, phlag command, expanded docs
- Remove sbag.py

Author: bo1929

README.md

@@ -4,66 +4,91 @@
 
 Given a species tree and a sequence of gene trees, Phlag aims to detect strong deviations from the multi-species coalescent (MSC) process using quadripartition quartet scores (QQS) and a hidden Markov model (HMM).
 
-The input is a sequence of gene trees, ideally sampled uniformly across a chromosome. The expected output is a set of "flagged" subsequences corresponding to model violations (e.g., admixture, recombination suppression, or dramatic changes in effective population size) and an updated species tree with new branch length (in coalescent units), estimated by ignoring the flagged gene trees.
+The input is a sequence of gene trees, ideally sampled uniformly across a chromosome.
+The expected output is a set of "flagged" subsequences corresponding to model violations (e.g., introgression, recombination suppression, or dramatic changes in effective population size) and an updated species tree with branch lengths in coalescent units (CU), estimated by ignoring the flagged gene trees.
 
 ## Installation
 
-You will need a Python 3.9 environment. To install, follow the steps below:
-
+Phlag requires Python == 3.9.
 ```shell
 micromamba create -n phlag python=3.9 -y
 micromamba activate phlag
 git clone https://github.com/bo1929/phlag.git
 cd phlag
-pip install -r requirements.txt
-./build.sh
+pip install .
 ```
 
+<!-- You can simply use pip. -->
+<!-- ```shell -->
+<!-- pip install phlag -->
+<!-- ``` -->
+
+<!-- Alternatively, install Phlag from the source. -->
+<!-- ```shell -->
+<!-- git clone https://github.com/bo1929/phlag.git -->
+<!-- cd phlag -->
+<!-- pip install . -->
+<!-- ``` -->
+
 ## Quickstart with a toy example
 
-Suppose we have a species tree (`test/neoaves.nwk`) and an ordered set of gene trees (`test/emission.gtrees`).
-This gene tree sequence is a mixture of two coalescent with recombination processes, with one key difference: the effective population size for certain *Charadriiformes* increased 10-fold.
+The `test/` directory contains a simulated dataset based on a Neoaves species tree with 191 taxa.
+We use this to demonstrate how Phlag detects a genomic region where the MSC model is violated.
 
-The mixture happens in a way that only a specific loci are affected by this change.
-In this particular example, 150 out of 1500 consecutive gene trees are generated by the alternative process.
-Our goal is to detect this region (gene trees in indices [913, 1062]) using Phlag.
+Suppose, we have a species tree (`test/neoaves.nwk`) and a sequence of 1500 gene trees (`test/emission.gtrees`), ordered along a chromosome.
+This gene tree sequence is a mixture of two coalescent with recombination (CwR) processes:
+- **Background (default)**: Gene trees simulated under the standard MSC with the original species tree parameters.
+- **Anomalous**: Gene trees simulated with a 10-fold increase in the effective population size for the branch leading to the *Charadriiformes* clade (labeled `N159`).
 
+Out of 1500 consecutive gene trees, 150 (10%) come from the anomalous process.
+These occupy a single contiguous block at indices [913, 1062] (0-indexed, inclusive).
+Our goal is to detect this region using Phlag.
 
 ### Input
-* **Species tree**: A Newick species tree with labeled internal nodes.
-* **Gene trees**: One Newick tree per line, ordered along the genome (not necessarily with labeled internal nodes).
-* **Focal branch(es)**: The label of the clade(s) to be targeted. In this example, the clade under suspicion is labeled `N159`.
+* **Species tree** (`-s`): A Newick species tree with labeled internal nodes.
+* **Gene trees** (`-g`): One Newick tree per line, ordered along the genome (not necessarily with labeled internal nodes).
+* **Focal edge(s)** (`-e`): The label of the internal node(s) defining the edge(s) to target. In this example, the clade under suspicion is labeled `N159`.
 
 ### Usage
 ```shell
-python phlag.py \
-  --species-tree test/neoaves.nwk \
-  --gene-trees  test/emission.gtrees \
-  -c N159 \
-  --num-iters 10 \
-  --read-qqs-freqs test/qqs.tsv \
-  --output-file results-neoaves-N159.txt
+phlag \
+  -s test/neoaves.nwk \
+  -g test/emission.gtrees \
+  -e N159 \
+  -L 10 \
+  --read-qqs-path test/qqs.tsv \
+  -o results-neoaves-N159.txt
 ```
-The quadripartition quartet scores are already computed and stored in `test/qqs.tsv`.
-You can directly read those to avoid extra computation for this toy example.
-If you do not give `--read-qqs-freqs`, QQS values need to be computed, which adds to the running time.
-Note that you can also save QQS values to read later for different runs (i.e., changing the focal branch or the hyperparameters) using `write-qqs-freqs`.
 
-For other options, run `python phlag.py --help`.
+The quadripartition quartet scores for this dataset are precomputed and stored in `test/qqs.tsv`.
+Reading them with `--read-qqs-path` avoids recomputing QQS values from scratch.
+If you omit this flag, QQS values will be computed on the fly, which adds to the running time.
+You can also precompute and save QQS values for later runs (e.g., with different focal edges or hyperparameters) using `--write-qqs-path`.
+
+For the full list of options, run `phlag --help`.
+
+### Key hyperparameters
+
+| Flag | Default | Description |
+|------|---------|-------------|
+| `--rho` | 0.9 | Controls sensitivity; reduce to flag more regions |
+| `--beta` | 5 | Controls contiguity of flagged intervals; reduce to merge nearby flags |
+| `--emission-lambda` | 1.0 | Controls expected deviation of anomalies from MSC |
+| `--eta` | 0.5 | Occupancy bias penalty on the marginal log-likelihood |
 
 **Note**: If your data lacks strong deviations from the MSC, you may see few or no flagged trees.
-Conversely, an incorrect/unreliable species tree may result in too many flagged gene trees.
-In such cases, you can attempt to adjust the hyperparameters to detect anomalies at the resolution you desire.
+Conversely, an incorrect or unreliable species tree may result in too many flagged gene trees.
+In such cases, you can adjust the hyperparameters to detect anomalies at the resolution you desire.
 
 A high-level and generic suggestion is running Phlag with varying hyperparameters (in a sensible range tested in the paper), and focusing on regions that are consistently detected. You may want to increase `--expected-num-anomalies` if you are looking for more granular results or if large contigous chunks are flagged. Similarly, `--expected-anamoly-proportion` can be varied between 0 and 0.5 to change to prior on the stationary distribution. 
 
 ### Output
-The output is rather simple to interpret:
-- **Header (Lines 1-5)**: Auxiliary information prefixed with `#`:
-    * Simply the invoked command and the version information.
+The output file contains:
+- **Header lines** (prefixed with `#`):
+    * The invoked command.
     * The initial species tree with branch lengths in CU, estimated using the entire gene tree sequence.
-    * The same species tree topology, but branch lengths are estimated only based on the gene trees that are not flagged.
-    * The labels of the clades for the focal branches that are used, and their final branch lengths.
-    * Distances between the emission distributions of the two states, separately for all the focal branches following the order given in the above line.
-- **State label**: `1` corresponds to flagged gene trees; `0` corresponds to MSC-compliant gene trees.
-- **Posterior probabilities**: The decoding probability for each gene tree.
+    * The labels and final branch lengths of the focal edges.
+    * The final species tree with branch lengths re-estimated after excluding the flagged gene trees.
+    * The final focal edge lengths, and distances between the emission distributions of the two HMM states for each focal edge.
+- **State labels**: A comma-separated sequence where `1` = flagged (anomalous) and `0` = MSC-compliant. Gene trees excluded due to missing data are labeled `nan`.
+- **Posterior probabilities**: The smoothed posterior probability of being in the anomalous state for each gene tree.

TODO.md

@@ -0,0 +1,11 @@
+# TODO
+
+* Make sure that both labels around the root gives you the same QQS values.
+
+## Code Issues
+
+### Design Improvements
+
+- **`phlag.py` (heavy `__init__`)**: Constructing a `Phlag` object triggers tree parsing, QQS computation, edge length estimation, discretization, HMM initialization, and output header writing. This makes it impossible to test or use any step in isolation. A more conventional pattern would be a lightweight `__init__` followed by explicit method calls.
+- **`hmm.py`**: Wildcard import `from dynamax.hidden_markov_model.inference import *` still present. Key symbols (`partial`, `jit`, `lax`, `Callable`) now have explicit imports, but `HMMPosterior`, `HMMPosteriorFiltered`, `get_trans_mat`, `compute_transition_probs` still come from the wildcard.
+- **`qqs.py` (performance)**: `MSC.compute_qqs` is O(n_edges * n_gt * n_taxa) in pure Python. The inner loop iterates every gene tree for every edge, calling `QQS.compute()` which traverses the tree. This was previously Cython. For large datasets this will be a major bottleneck. Consider Cython/Numba/JAX acceleration for the inner loop.