aarmey · Sep 30, 2026 · Sep 30, 2026 · Sep 30, 2026
diff --git a/.gitignore b/.gitignore
@@ -25,3 +25,6 @@ lineage/data/taxol
 output/
 docs/_build/
 
+
+# KP-Tracer analysis intermediates (large; regenerate with analysis/kptracer scripts)
+analysis/kptracer/results/*.pkl
diff --git a/analysis/kptracer/README.md b/analysis/kptracer/README.md
@@ -0,0 +1,34 @@
+# Latent-state tHMM on KP-Tracer phylogenies (issue #1014)
+
+The findings are in [`REPORT.md`](REPORT.md). This directory holds the scripts that produce them.
+
+## Data
+
+Download the processed KP-Tracer release (Yang et al. 2022, Zenodo record 5847462, CC-BY, 1.3 GB):
+
+```sh
+mkdir -p ~/data/kptracer && cd ~/data/kptracer
+curl -L -o KPTracer-Data.tar.gz https://zenodo.org/api/records/5847462/files/KPTracer-Data.tar.gz/content
+tar xzf KPTracer-Data.tar.gz
+```
+
+The scripts look for `~/data/kptracer/KPTracer-Data` or `$KPTRACER_DATA`.
+
+## Pipeline
+
+Run from the repository root. `anndata` is only needed to read the h5ad, and `matplotlib` only for figures.
+
+```sh
+uv run python analysis/kptracer/simulation_study.py                 # results/simulation.csv
+uv run python analysis/kptracer/scalability.py                      # results/scalability.csv (run from this dir)
+uv run --with anndata python analysis/kptracer/prepare_data.py      # results/kptracer_inputs.pkl
+cd analysis/kptracer
+uv run python kptracer_analysis.py select        # K and branch-length selection  (~45 min on 46 cores)
+uv run python kptracer_analysis.py fit 12        # anchored (fixed emissions), unsupervised K=12, sensitivity fits
+uv run python kptracer_analysis.py bootstrap anchored
+uv run python kptracer_analysis.py compare
+uv run --with matplotlib python make_figures.py
+uv run python bdmm_comparison.py                 # tHMM side of the BDMM-Prime comparison (bdmm/ holds the BEAST XML)
+```
+
+`results/*.pkl` are large intermediates and are git-ignored. The CSV summaries and `figures/` are committed.
diff --git a/analysis/kptracer/REPORT.md b/analysis/kptracer/REPORT.md
diff --git a/analysis/kptracer/bdmm/bdmm_fixedtree.xml b/analysis/kptracer/bdmm/bdmm_fixedtree.xml
@@ -0,0 +1,281 @@
+<?xml version="1.0" encoding="UTF-8" standalone="no"?>
+<beast version="2.7" namespace="beast.base.inference:beast.base.inference.parameter:beast.base.inference.distribution:beast.base.inference.operator:beast.base.evolution.tree:beast.base.evolution.alignment:beast.base.evolution.operator:beast.base.inference.util:bdmmprime.distribution:bdmmprime.parameterization:bdmmprime.util">
+
+  <taxonset id="taxa" spec="TaxonSet">
+    <taxon id="n123" spec="Taxon"/>
+    <taxon id="n185" spec="Taxon"/>
+    <taxon id="n186" spec="Taxon"/>
+    <taxon id="n130" spec="Taxon"/>
+    <taxon id="n66" spec="Taxon"/>
+    <taxon id="n175" spec="Taxon"/>
+    <taxon id="n176" spec="Taxon"/>
+    <taxon id="n159" spec="Taxon"/>
+    <taxon id="n160" spec="Taxon"/>
+    <taxon id="n137" spec="Taxon"/>
+    <taxon id="n33" spec="Taxon"/>
+    <taxon id="n34" spec="Taxon"/>
+    <taxon id="n149" spec="Taxon"/>
+    <taxon id="n150" spec="Taxon"/>
+    <taxon id="n47" spec="Taxon"/>
+    <taxon id="n96" spec="Taxon"/>
+    <taxon id="n100" spec="Taxon"/>
+    <taxon id="n99" spec="Taxon"/>
+    <taxon id="n197" spec="Taxon"/>
+    <taxon id="n198" spec="Taxon"/>
+    <taxon id="n177" spec="Taxon"/>
+    <taxon id="n178" spec="Taxon"/>
+    <taxon id="n191" spec="Taxon"/>
+    <taxon id="n192" spec="Taxon"/>
+    <taxon id="n139" spec="Taxon"/>
+    <taxon id="n181" spec="Taxon"/>
+    <taxon id="n182" spec="Taxon"/>
+    <taxon id="n131" spec="Taxon"/>
+    <taxon id="n106" spec="Taxon"/>
+    <taxon id="n127" spec="Taxon"/>
+    <taxon id="n121" spec="Taxon"/>
+    <taxon id="n122" spec="Taxon"/>
+    <taxon id="n135" spec="Taxon"/>
+    <taxon id="n136" spec="Taxon"/>
+    <taxon id="n143" spec="Taxon"/>
+    <taxon id="n144" spec="Taxon"/>
+    <taxon id="n195" spec="Taxon"/>
+    <taxon id="n196" spec="Taxon"/>
+    <taxon id="n147" spec="Taxon"/>
+    <taxon id="n148" spec="Taxon"/>
+    <taxon id="n45" spec="Taxon"/>
+    <taxon id="n46" spec="Taxon"/>
+    <taxon id="n37" spec="Taxon"/>
+    <taxon id="n165" spec="Taxon"/>
+    <taxon id="n110" spec="Taxon"/>
+    <taxon id="n171" spec="Taxon"/>
+    <taxon id="n172" spec="Taxon"/>
+    <taxon id="n69" spec="Taxon"/>
+    <taxon id="n71" spec="Taxon"/>
+    <taxon id="n193" spec="Taxon"/>
+    <taxon id="n194" spec="Taxon"/>
+    <taxon id="n187" spec="Taxon"/>
+    <taxon id="n188" spec="Taxon"/>
+    <taxon id="n102" spec="Taxon"/>
+    <taxon id="n103" spec="Taxon"/>
+    <taxon id="n104" spec="Taxon"/>
+    <taxon id="n53" spec="Taxon"/>
+    <taxon id="n85" spec="Taxon"/>
+    <taxon id="n86" spec="Taxon"/>
+    <taxon id="n88" spec="Taxon"/>
+    <taxon id="n125" spec="Taxon"/>
+    <taxon id="n73" spec="Taxon"/>
+    <taxon id="n151" spec="Taxon"/>
+    <taxon id="n152" spec="Taxon"/>
+    <taxon id="n77" spec="Taxon"/>
+    <taxon id="n78" spec="Taxon"/>
+    <taxon id="n120" spec="Taxon"/>
+    <taxon id="n133" spec="Taxon"/>
+    <taxon id="n134" spec="Taxon"/>
+    <taxon id="n141" spec="Taxon"/>
+    <taxon id="n142" spec="Taxon"/>
+    <taxon id="n183" spec="Taxon"/>
+    <taxon id="n184" spec="Taxon"/>
+    <taxon id="n117" spec="Taxon"/>
+    <taxon id="n118" spec="Taxon"/>
+    <taxon id="n60" spec="Taxon"/>
+    <taxon id="n91" spec="Taxon"/>
+    <taxon id="n92" spec="Taxon"/>
+    <taxon id="n161" spec="Taxon"/>
+    <taxon id="n157" spec="Taxon"/>
+    <taxon id="n158" spec="Taxon"/>
+    <taxon id="n167" spec="Taxon"/>
+    <taxon id="n168" spec="Taxon"/>
+    <taxon id="n63" spec="Taxon"/>
+    <taxon id="n113" spec="Taxon"/>
+    <taxon id="n179" spec="Taxon"/>
+    <taxon id="n180" spec="Taxon"/>
+    <taxon id="n108" spec="Taxon"/>
+    <taxon id="n163" spec="Taxon"/>
+    <taxon id="n164" spec="Taxon"/>
+    <taxon id="n84" spec="Taxon"/>
+    <taxon id="n89" spec="Taxon"/>
+    <taxon id="n90" spec="Taxon"/>
+    <taxon id="n153" spec="Taxon"/>
+    <taxon id="n154" spec="Taxon"/>
+    <taxon id="n155" spec="Taxon"/>
+    <taxon id="n156" spec="Taxon"/>
+    <taxon id="n170" spec="Taxon"/>
+    <taxon id="n189" spec="Taxon"/>
+    <taxon id="n190" spec="Taxon"/>
+  </taxonset>
+
+  <!-- FIXED tree: not in the state, no tree operators -->
+  <tree id="tree" spec="TreeParser" IsLabelledNewick="true" adjustTipHeights="true" taxonset="@taxa"
+        newick="((((n185:0.105435,n186:0.105435)n75:0.758506,((n175:0.169748,n176:0.169748)n145:0.112940,(n159:0.234216,n160:0.234216)n146:0.048472)n76:0.581254)n43:0.449058,(n123:0.450563,((n137:0.324692,(n149:0.272999,n150:0.272999)n138:0.051694)n129:0.099779,n130:0.424471)n124:0.026092)n44:0.862437)n1:1.968832,(((((n47:1.159471,(n181:0.117813,n182:0.117813)n48:1.041658)n23:1.098771,((n131:0.409898,(n165:0.217494,(n193:0.047186,n194:0.047186)n166:0.170309)n132:0.192403)n55:0.683481,(((n187:0.104629,n188:0.104629)n109:0.501520,n110:0.606149)n105:0.017891,n106:0.624040)n56:0.469339)n24:1.164863)n15:0.438788,(((n127:0.431499,(n171:0.173877,n172:0.173877)n128:0.257622)n95:0.332616,n96:0.764115)n79:0.093977,(n100:0.734049,n99:0.734049)n80:0.124044)n16:1.838938)n5:0.300383,((n33:1.783769,n34:1.783769)n21:0.587406,(((n121:0.490003,n122:0.490003)n111:0.102243,(n135:0.369630,n136:0.369630)n112:0.222616)n97:0.169138,(n197:0.012586,n198:0.012586)n98:0.748797)n22:1.609791)n6:0.626239)n3:0.015266,((((n177:0.135231,n178:0.135231)n173:0.037077,(n191:0.052557,n192:0.052557)n174:0.119751)n65:0.745767,n66:0.918074)n7:2.050193,(((n139:0.323601,(n143:0.288174,n144:0.288174)n140:0.035427)n41:1.038820,((n195:0.015571,n196:0.015571)n115:0.514336,(n147:0.281240,n148:0.281240)n116:0.248668)n42:0.832513)n10:1.591565,((((n69:0.909870,((n125:0.436570,(n183:0.108761,n184:0.108761)n126:0.327809)n101:0.212951,n102:0.649521)n70:0.260349)n19:1.594987,((n103:0.645152,n104:0.645152)n49:0.464789,(n53:1.095520,(n73:0.887923,((n157:0.235615,n158:0.235615)n81:0.607617,(n167:0.214877,n168:0.214877)n82:0.628355)n74:0.044692)n54:0.207596)n50:0.014421)n20:1.394916)n17:0.101622,((n71:0.903043,(n85:0.829271,n86:0.829271)n72:0.073772)n25:1.128496,(((n151:0.267928,n152:0.267928)n31:1.569555,(n77:0.859500,n78:0.859500)n32:0.977983)n27:0.069964,(((n117:0.510409,n118:0.510409)n35:1.142729,((n63:0.927585,((n153:0.264503,n154:0.264503)n107:0.346663,n108:0.611166)n64:0.316418)n59:0.114692,n60:1.042276)n36:0.610861)n29:0.219068,((n91:0.777811,n92:0.777811)n57:0.302458,((n113:0.531419,(n163:0.218871,n164:0.218871)n114:0.312548)n61:0.413352,((((n155:0.236850,n156:0.236850)n93:0.527531,((n189:0.086787,n190:0.086787)n169:0.097964,n170:0.184751)n94:0.579630)n83:0.066969,n84:0.831350)n67:0.078867,(n89:0.779814,n90:0.779814)n68:0.130402)n62:0.034554)n58:0.135498)n30:0.791937)n28:0.035241)n26:0.124093)n18:0.574939)n11:0.316059,((n45:1.179748,n46:1.179748)n13:1.605255,(n37:1.583666,((((n161:0.221276,(n179:0.121742,n180:0.121742)n162:0.099534)n119:0.287553,n120:0.508829)n51:0.592224,(n133:0.382097,n134:0.382097)n52:0.718956)n39:0.284872,((n141:0.313129,n142:0.313129)n87:0.509341,n88:0.822470)n40:0.563455)n38:0.197741)n14:1.201337)n12:0.137534)n9:0.031449)n8:0.014281)n4:0.044412)n2:0.269154)n0:0.000000;"/>
+
+  <typeTraitSet id="typeTraitSet" spec="TraitSet" traitname="type" taxa="@taxa">
+n123=1,
+n185=2,
+n186=1,
+n130=0,
+n66=1,
+n175=0,
+n176=0,
+n159=0,
+n160=0,
+n137=0,
+n33=0,
+n34=1,
+n149=0,
+n150=0,
+n47=0,
+n96=2,
+n100=2,
+n99=2,
+n197=2,
+n198=2,
+n177=0,
+n178=1,
+n191=1,
+n192=1,
+n139=2,
+n181=0,
+n182=0,
+n131=2,
+n106=1,
+n127=2,
+n121=2,
+n122=2,
+n135=2,
+n136=2,
+n143=2,
+n144=2,
+n195=1,
+n196=1,
+n147=1,
+n148=1,
+n45=2,
+n46=1,
+n37=2,
+n165=2,
+n110=1,
+n171=3,
+n172=3,
+n69=1,
+n71=1,
+n193=2,
+n194=2,
+n187=1,
+n188=1,
+n102=0,
+n103=1,
+n104=1,
+n53=2,
+n85=1,
+n86=1,
+n88=0,
+n125=0,
+n73=0,
+n151=3,
+n152=3,
+n77=2,
+n78=1,
+n120=1,
+n133=0,
+n134=1,
+n141=1,
+n142=1,
+n183=0,
+n184=0,
+n117=1,
+n118=1,
+n60=1,
+n91=1,
+n92=0,
+n161=1,
+n157=0,
+n158=0,
+n167=0,
+n168=0,
+n63=2,
+n113=0,
+n179=1,
+n180=1,
+n108=1,
+n163=0,
+n164=0,
+n84=0,
+n89=1,
+n90=1,
+n153=1,
+n154=1,
+n155=0,
+n156=1,
+n170=0,
+n189=1,
+n190=1
+  </typeTraitSet>
+
+  <typeSet id="typeSet" spec="TypeSet" typeTraitSet="@typeTraitSet"/>
+  <processLength id="processLength" spec="ProcessLength" tree="@tree"/>
+
+  <run id="mcmc" spec="MCMC" chainLength="2500000">
+    <state id="state" storeEvery="5000">
+      <stateNode id="birthRate" spec="RealParameter" lower="0.0" dimension="1">1.0</stateNode>
+      <stateNode id="migrationRate" spec="RealParameter" lower="0.0" dimension="12">0.1</stateNode>
+      <stateNode id="rootFreqs" spec="RealParameter" lower="0.0" upper="1.0" dimension="4">0.25</stateNode>
+    </state>
+
+    <distribution id="posterior" spec="CompoundDistribution">
+      <distribution id="prior" spec="CompoundDistribution">
+        <distribution id="birthRatePrior" spec="Prior" x="@birthRate">
+          <distr spec="LogNormalDistributionModel" meanInRealSpace="false" M="0.0" S="1.0"/>
+        </distribution>
+        <distribution id="migrationRatePrior" spec="Prior" x="@migrationRate">
+          <distr spec="Exponential" mean="0.3"/>
+        </distribution>
+        <distribution id="rootFreqsPrior" spec="Prior" x="@rootFreqs">
+          <distr spec="Dirichlet"><parameter spec="RealParameter" name="alpha" value="1.0 1.0 1.0 1.0"/></distr>
+        </distribution>
+      </distribution>
+      <distribution id="likelihood" spec="CompoundDistribution">
+        <distribution id="treePrior" spec="BirthDeathMigrationDistribution" tree="@tree"
+                      typeTraitSet="@typeTraitSet" startTypePriorProbs="@rootFreqs"
+                      conditionOnRoot="true" parallelize="true">
+          <parameterization id="param" spec="CanonicalParameterization" typeSet="@typeSet" processLength="@processLength">
+            <birthRate spec="SkylineVectorParameter" typeSet="@typeSet" processLength="@processLength" skylineValues="@birthRate"/>
+            <deathRate spec="SkylineVectorParameter" typeSet="@typeSet" processLength="@processLength">
+              <skylineValues spec="RealParameter" value="0.0"/>
+            </deathRate>
+            <samplingRate spec="SkylineVectorParameter" typeSet="@typeSet" processLength="@processLength">
+              <skylineValues spec="RealParameter" value="0.0"/>
+            </samplingRate>
+            <removalProb spec="SkylineVectorParameter" typeSet="@typeSet" processLength="@processLength">
+              <skylineValues spec="RealParameter" value="1.0"/>
+            </removalProb>
+            <rhoSampling spec="TimedParameter" typeSet="@typeSet" processLength="@processLength" timesAreAges="true">
+              <times spec="RealParameter" value="0.0"/>
+              <values spec="RealParameter" value="1.0"/>
+            </rhoSampling>
+            <migrationRate spec="SkylineMatrixParameter" typeSet="@typeSet" processLength="@processLength" skylineValues="@migrationRate"/>
+          </parameterization>
+        </distribution>
+      </distribution>
+    </distribution>
+
+    <operator id="birthRateScale" spec="ScaleOperator" parameter="@birthRate" scaleFactor="0.5" weight="3.0"/>
+    <operator id="migRateScale" spec="ScaleOperator" parameter="@migrationRate" scaleFactor="0.5" weight="12.0"/>
+    <operator id="migRateScaleAll" spec="ScaleOperator" parameter="@migrationRate" scaleAll="true" scaleFactor="0.75" weight="2.0"/>
+    <operator id="rootFreqsExchange" spec="DeltaExchangeOperator" parameter="@rootFreqs" delta="0.2" weight="2.0"/>
+
+    <logger id="tracelog" spec="Logger" fileName="bdmm_fixedtree_chain1.log" logEvery="500" sanitiseHeaders="true">
+      <log idref="posterior"/>
+      <log idref="likelihood"/>
+      <log idref="prior"/>
+      <log idref="treePrior"/>
+      <log idref="birthRate"/>
+      <log idref="migrationRate"/>
+      <log idref="rootFreqs"/>
+    </logger>
+    <logger id="screenlog" spec="Logger" logEvery="25000">
+      <log idref="posterior"/>
+      <log idref="likelihood"/>
+      <log spec="ESS" arg="@posterior"/>
+    </logger>
+  </run>
+</beast>