Repository: jjmccollum/teiphy

Issue:
Assignment of random starting rate parameter values, new prior distributions, and support for `--seed` option
For larger collations with increasingly varied substitution matrices at different sites/variation units, it becomes increasingly likely that at least one of the substitution matrices will be singular under the initial assignments of transcriptional rate parameters (which presently all default to 2.0). This ceases to be a problem after BEAST 2's first evaluation of the site likelihoods, because in subsequent likelihood calculations, it samples these rate parameter values randomly according to a prior distribution, and the resulting substitution matrices will almost never be singular. Since BEAST 2 XML inputs require that every `RealParameter` element have an initial `@value` attribute specified, we need to supply these values, but we should be able to avoid singular matrices if we set these values randomly in `to_beast`.

This can be done by using distributions available in `numpy`, which is already a dependency of `teiphy`. I suggest replacing the current offset log-normal prior distribution assumed for transcriptional rate parameters with a simple gamma distribution (perhaps with `alpha="5"` and `beta="2"` as default values). Then we can sample random initial values for the rate parameters from this distribution using something like

```python
rng = np.random.default_rng(seed)
sample = rng.gamma(5.0, 2.0)
```

To ensure replicability for generated outputs, we will need to add support for a `--seed` command-line argument.

Discussion:

Patch:
diff --git a/.github/workflows/beast.yml b/.github/workflows/beast.yml
index bfcf37c..118ece7 100644
--- a/.github/workflows/beast.yml
+++ b/.github/workflows/beast.yml
@@ -26,9 +26,9 @@ jobs:
     - run: poetry install
     - name: Testing
       run: |
-        poetry run teiphy -t reconstructed -t defective -t orthographic -t subreading -m overlap -m lac -s"*" -s T --fill-correctors --clock strict example/ubs_ephesians.xml beast_strict_example.xml
-        poetry run teiphy -t reconstructed -t defective -t orthographic -t subreading -m overlap -m lac -s"*" -s T --fill-correctors --clock uncorrelated example/ubs_ephesians.xml beast_uncorrelated_example.xml
-        poetry run teiphy -t reconstructed -t defective -t orthographic -t subreading -m overlap -m lac -s"*" -s T --fill-correctors --clock local example/ubs_ephesians.xml beast_local_example.xml
+        poetry run teiphy -t reconstructed -t defective -t orthographic -t subreading -m overlap -m lac -s"*" -s T --fill-correctors --clock strict --seed 1337 example/ubs_ephesians.xml beast_strict_example.xml
+        poetry run teiphy -t reconstructed -t defective -t orthographic -t subreading -m overlap -m lac -s"*" -s T --fill-correctors --clock uncorrelated --seed 1337 example/ubs_ephesians.xml beast_uncorrelated_example.xml
+        poetry run teiphy -t reconstructed -t defective -t orthographic -t subreading -m overlap -m lac -s"*" -s T --fill-correctors --clock local --seed 1337 example/ubs_ephesians.xml beast_local_example.xml
     - name: Install Phylogenetics Package
       run: |
         export JAVA_HOME=/opt/hostedtoolcache/Java_Zulu_jdk+fx/17.0.6-10/x64
diff --git a/teiphy/collation.py b/teiphy/collation.py
index 0dd9f83..a87f297 100644
--- a/teiphy/collation.py
+++ b/teiphy/collation.py
@@ -7,7 +7,7 @@ from datetime import datetime  # for calculating the current year (for dating an
 import time  # to time calculations for users
 import string  # for easy retrieval of character ranges
 from lxml import etree as et  # for reading TEI XML inputs
-import numpy as np  # for collation matrix outputs
+import numpy as np  # for random number sampling and collation matrix outputs
 import pandas as pd  # for writing to DataFrames, CSV, Excel, etc.
 from slugify import slugify  # for converting Unicode text from readings to ASCII for NEXUS
 from jinja2 import Environment, PackageLoader, select_autoescape  # for filling output XML templates
@@ -1291,7 +1291,11 @@ class Collation:
         return root_frequencies_string
 
     def to_beast(
-        self, file_addr: Union[Path, str], drop_constant: bool = False, clock_model: ClockModel = ClockModel.strict
+        self,
+        file_addr: Union[Path, str],
+        drop_constant: bool = False,
+        clock_model: ClockModel = ClockModel.strict,
+        seed: int = None,
     ):
         """Writes this Collation to a file in BEAST format with the given address.
 
@@ -1299,6 +1303,7 @@ class Collation:
             file_addr: A string representing the path to an output file.
             drop_constant (bool, optional): An optional flag indicating whether to ignore variation units with one substantive reading.
             clock_model: A ClockModel option indicating which clock model to use.
+            seed: A seed for random number generation (for setting initial values of unspecified transcriptional rates).
         """
         # Populate a list of sites that will correspond to columns of the sequence alignment:
         substantive_variation_unit_ids = self.variation_unit_ids
@@ -1453,15 +1458,16 @@ class Collation:
             intrinsic_category_object["odds"] = odds if odds is not None else 1.0
             intrinsic_category_object["estimate"] = "false" if odds is not None else "true"
             intrinsic_category_objects.append(intrinsic_category_object)
-        # The proceed to transcriptional rate categories:
+        # Then proceed to transcriptional rate categories:
+        rng = np.random.default_rng(seed)
         for transcriptional_category in self.transcriptional_categories:
             transcriptional_category_object = {}
             # Copy the ID of this transcriptional category:
             transcriptional_category_object["id"] = transcriptional_category
             # Then copy the rate of this transcriptional category,
-            # setting it to 2.0 if it is not specified and setting the estimate attribute accordingly:
+            # setting it to a random number sampled from a Gamma distribution if it is not specified and setting the estimate attribute accordingly:
             rate = self.transcriptional_rates_by_id[transcriptional_category]
-            transcriptional_category_object["rate"] = rate if rate is not None else 2.0
+            transcriptional_category_object["rate"] = rate if rate is not None else rng.gamma(5.0, 2.0)
             transcriptional_category_object["estimate"] = "false" if rate is not None else "true"
             transcriptional_category_objects.append(transcriptional_category_object)
         # Now render the output XML file using the Jinja template:
@@ -1893,6 +1899,7 @@ class Collation:
         calibrate_dates: bool = False,
         mrbayes: bool = False,
         clock_model: ClockModel = ClockModel.strict,
+        seed: int = None,
     ):
         """Writes this Collation to the file with the given address.
 
@@ -1933,6 +1940,7 @@ class Collation:
             clock_model: A ClockModel option indicating which type of clock model to use.
                 This option is intended for inputs to MrBayes and BEAST 2.
                 MrBayes does not presently support a local clock model, so it will default to a strict clock model if a local clock model is specified.
+            seed: A seed for random number generation (for setting initial values of unspecified transcriptional rates in BEAST 2 XML output).
         """
         file_addr = Path(file_addr)
         format = format or Format.infer(
@@ -1961,7 +1969,7 @@ class Collation:
             return self.to_fasta(file_addr, drop_constant=drop_constant)
 
         if format == format.BEAST:
-            return self.to_beast(file_addr, drop_constant=drop_constant, clock_model=clock_model)
+            return self.to_beast(file_addr, drop_constant=drop_constant, clock_model=clock_model, seed=seed)
 
         if format == Format.CSV:
             return self.to_csv(
diff --git a/teiphy/main.py b/teiphy/main.py
index 766b1f9..be29648 100644
--- a/teiphy/main.py
+++ b/teiphy/main.py
@@ -75,6 +75,10 @@ def to_file(
         False,
         help="Treat missing characters/variation units as having a contribution of 1 split over all states/readings; if False, then missing data is ignored (i.e., all states are 0). Not applicable for non-tabular formats.",
     ),
+    seed: int = typer.Option(
+        None,
+        help="Seed for random number generation (used for setting default initial values of transcriptional rate parameters for BEAST 2 XML output); if not specified, then the default seeding of the numpy.random.default_rng class will be used.",
+    ),
     verbose: bool = typer.Option(False, help="Enable verbose logging (mostly for debugging purposes)."),
     version: bool = typer.Option(
         False,
@@ -128,4 +132,5 @@ def to_file(
         clock_model=clock,
         long_table=long_table,
         split_missing=split_missing,
+        seed=seed,
     )
diff --git a/teiphy/templates/beast_template.xml b/teiphy/templates/beast_template.xml
index 59231f3..e0625bc 100644
--- a/teiphy/templates/beast_template.xml
+++ b/teiphy/templates/beast_template.xml
@@ -62,7 +62,7 @@
 
             <!-- Start transcriptional rate parameters -->
             {%- for tc in transcriptional_categories %}
-            <stateNode spec="parameter.RealParameter" id="{{ tc.id }}_rate" lower="1.0" upper="Infinity" value="{{ tc.rate }}" estimate="{{ tc.estimate }}"/>
+            <stateNode spec="parameter.RealParameter" id="{{ tc.id }}_rate" lower="0.0" upper="Infinity" value="{{ tc.rate }}" estimate="{{ tc.estimate }}"/>
             {%- endfor %}
             <!-- We include a "default" rate fixed at 1 that corresponds to (unlikely) transitions with no transcriptional explanation.
             All other rates in the substitution matrices will be estimated relative to this.
@@ -115,7 +115,7 @@
                 {%- for tc in transcriptional_categories %}
                 {%- if tc.estimate == "true" %}
                 <distribution spec="Prior" id="{{ tc.id }}_ratePrior" x="@{{ tc.id }}_rate">
-                    <distr spec="LogNormalDistributionModel" offset="1.0" M="1.0" S="1.0"/>
+                    <distr spec="Gamma" alpha="5.0" beta="2.0"/>
                 </distribution>
                 {%- endif %}
                 {%- endfor %}


Test patch:
diff --git a/tests/test_main.py b/tests/test_main.py
index 6b46abc..dbd89a0 100644
--- a/tests/test_main.py
+++ b/tests/test_main.py
@@ -696,6 +696,8 @@ def test_to_beast_variable_rates():
                 "-s*",
                 "-sT",
                 "--fill-correctors",
+                "--seed",
+                "1337",
                 str(input_example),
                 str(output),
             ],
@@ -703,13 +705,21 @@ def test_to_beast_variable_rates():
         assert result.exit_code == 0
         assert output.exists()
         beast_xml = et.parse(output, parser=parser)
-        for xml_transcriptional_category in xml_transcriptional_categories:
+        expected_rates = [
+            9.499649179019467,
+            8.750762619939191,
+            25.00430409157674,
+            19.873724085272503,
+            9.988999535304949,
+            21.828078913817116,
+        ]  # Gamma(5,2)-distributed numbers generated with seed 1337
+        for i, xml_transcriptional_category in enumerate(xml_transcriptional_categories):
             transcriptional_category = xml_transcriptional_category.get("{%s}id" % xml_ns)
             beast_xml_transcriptional_rate_categories = beast_xml.xpath(
                 "//stateNode[@id=\"%s_rate\"]" % transcriptional_category
             )
             assert len(beast_xml_transcriptional_rate_categories) == 1
-            assert float(beast_xml_transcriptional_rate_categories[0].get("value")) == 2.0
+            assert float(beast_xml_transcriptional_rate_categories[0].get("value")) == expected_rates[i]
             assert beast_xml_transcriptional_rate_categories[0].get("estimate") == "true"