Coverage for teiphy/collation.py: 100.00%
1677 statements
« prev ^ index » next coverage.py v7.9.2, created at 2026-08-27 23:22 +0000
1#!/usr/bin/env python3
3from enum import Enum
4from typing import List, Union
5import os
6from pathlib import Path
7from datetime import datetime # for calculating the current year (for dating and tree height purposes)
8import math # for special functions
9import time # to time calculations for users
10import string # for easy retrieval of character ranges
11from lxml import etree as et # for reading TEI XML inputs
12import numpy as np # for random number sampling and collation matrix outputs
13import pandas as pd # for writing to DataFrames, CSV, Excel, etc.
14from slugify import slugify # for converting Unicode text from readings to ASCII for NEXUS
15from jinja2 import Environment, PackageLoader, select_autoescape # for filling output XML templates
16from tqdm import tqdm # for progress bars
18from .common import xml_ns, tei_ns
19from .format import Format
20from .witness import Witness
21from .variation_unit import VariationUnit
24class ParsingException(Exception):
25 pass
28class WitnessDateException(Exception):
29 pass
32class IntrinsicRelationsException(Exception):
33 pass
36class ClockModel(str, Enum):
37 strict = "strict"
38 uncorrelated = "uncorrelated"
39 local = "local"
42class AncestralLogger(str, Enum):
43 state = "state"
44 sequence = "sequence"
45 none = "none"
48class TableType(str, Enum):
49 matrix = "matrix"
50 distance = "distance"
51 similarity = "similarity"
52 idf = "idf"
53 mean_idf = "mean-idf"
54 mi = "mi"
55 mean_mi = "mean-mi"
56 nexus = "nexus"
57 long = "long"
60class SplitMissingType(str, Enum):
61 uniform = "uniform"
62 proportional = "proportional"
65class TransformMatrixType(str, Enum):
66 stddev = "stddev"
67 mad = "mad"
70class Collation:
71 """Base class for storing TEI XML collation data internally.
73 This corresponds to the entire XML tree, rooted at the TEI element of the collation.
75 Attributes:
76 manuscript_suffixes: A list of suffixes used to distinguish manuscript subwitnesses like first hands, correctors, main texts, alternate texts, and multiple attestations from their base witnesses.
77 trivial_reading_types: A set of reading types (e.g., "reconstructed", "defective", "orthographic", "subreading") whose readings should be collapsed under the previous substantive reading.
78 missing_reading_types: A set of reading types (e.g., "lac", "overlap") whose readings should be treated as missing data.
79 fill_corrector_lacunae: A boolean flag indicating whether or not to fill "lacunae" in witnesses with type "corrector".
80 fragmentary_threshold: A float representing the proportion such that all witnesses extant at fewer than this proportion of variation units are filtered out of the collation.
81 fill_correctors_threshold: A float representing the proportion such that all correctors extant at fewer than this proportion of variation units are not filled in.
82 witnesses: A list of Witness instances contained in this Collation.
83 witness_index_by_id: A dictionary mapping base witness ID strings to their int indices in the witnesses list.
84 variation_units: A list of VariationUnit instances contained in this Collation.
85 readings_by_witness: A dictionary mapping base witness ID strings to lists of reading support coefficients for all units (with at least two substantive readings).
86 substantive_variation_unit_ids: A list of ID strings for variation units with two or more substantive readings.
87 substantive_variation_unit_reading_tuples: A list of (variation unit ID, reading ID) tuples for substantive readings.
88 verbose: A boolean flag indicating whether or not to print timing and debugging details for the user.
89 """
91 def __init__(
92 self,
93 xml: et.ElementTree,
94 manuscript_suffixes: List[str] = [],
95 trivial_reading_types: List[str] = [],
96 missing_reading_types: List[str] = [],
97 fill_corrector_lacunae: bool = False,
98 fragmentary_threshold: float = None,
99 fill_correctors_threshold: float = None,
100 dates_file: Union[Path, str] = None,
101 verbose: bool = False,
102 ):
103 """Constructs a new Collation instance with the given settings.
105 Args:
106 xml: An lxml.etree.ElementTree representing an XML tree rooted at a TEI element.
107 manuscript_suffixes: An optional list of suffixes used to distinguish manuscript subwitnesses like first hands, correctors, main texts, alternate texts, and multiple attestations from their base witnesses.
108 trivial_reading_types: An optional set of reading types (e.g., "reconstructed", "defective", "orthographic", "subreading") whose readings should be collapsed under the previous substantive reading.
109 missing_reading_types: An optional set of reading types (e.g., "lac", "overlap") whose readings should be treated as missing data.
110 fill_corrector_lacunae: An optional flag indicating whether or not to fill "lacunae" in witnesses with type "corrector".
111 fragmentary_threshold: An optional float representing the proportion such that all witnesses extant at fewer than this proportion of variation units are filtered out of the collation.
112 fill_correctors_threshold: An optional float representing the proportion such that all correctors extant at fewer than this proportion of variation units are not filled in.
113 dates_file: An optional path to a CSV file containing witness IDs, minimum dates, and maximum dates. If specified, then for all witnesses in the first column, any existing date ranges for them in the TEI XML collation will be ignored.
114 verbose: An optional flag indicating whether or not to print timing and debugging details for the user.
115 """
116 self.manuscript_suffixes = manuscript_suffixes
117 self.trivial_reading_types = set(trivial_reading_types)
118 self.missing_reading_types = set(missing_reading_types)
119 self.fill_corrector_lacunae = fill_corrector_lacunae
120 self.fragmentary_threshold = fragmentary_threshold
121 self.fill_correctors_threshold = fill_correctors_threshold
122 self.verbose = verbose
123 self.witnesses = []
124 self.witness_index_by_id = {}
125 self.variation_units = []
126 self.readings_by_witness = {}
127 self.variation_unit_ids = []
128 self.substantive_variation_unit_reading_tuples = []
129 self.substantive_readings_by_variation_unit_id = {}
130 self.weight_categories = []
131 self.weights_by_id = {}
132 self.intrinsic_categories = []
133 self.intrinsic_odds_by_id = {}
134 self.transcriptional_categories = []
135 self.transcriptional_rates_by_id = {}
136 self.origin_date_range = []
137 # Now parse the XML tree to populate these data structures:
138 if self.verbose:
139 print("Initializing collation...")
140 t0 = time.time()
141 self.parse_origin_date_range(xml)
142 self.parse_list_wit(xml)
143 self.validate_wits(xml)
144 # If a dates file was specified, then update the witness date ranges manually:
145 if dates_file is not None:
146 self.update_witness_date_ranges_from_dates_file(dates_file)
147 # If the upper bound on a work's date of origin is not defined, then attempt to assign it an upper bound based on the witness dates;
148 # otherwise, attempt to assign lower bounds to witness dates based on it:
149 if self.origin_date_range[1] is None:
150 self.update_origin_date_range_from_witness_date_ranges()
151 else:
152 self.update_witness_date_ranges_from_origin_date_range()
153 self.parse_weights(xml)
154 self.parse_intrinsic_odds(xml)
155 self.parse_transcriptional_rates(xml)
156 self.parse_apps(xml)
157 self.validate_intrinsic_relations()
158 self.parse_readings_by_witness()
159 # If a threshold of readings for fragmentary witnesses is specified, then filter the witness list using the dictionary mapping witness IDs to readings:
160 if self.fragmentary_threshold is not None:
161 self.filter_fragmentary_witnesses(xml)
162 t1 = time.time()
163 if self.verbose:
164 print("Total time to initialize collation: %0.4fs." % (t1 - t0))
166 def parse_origin_date_range(self, xml: et.ElementTree):
167 """Given an XML tree for a collation, populates this Collation's list of origin date bounds.
169 Args:
170 xml: An lxml.etree.ElementTree representing an XML tree rooted at a TEI element.
171 """
172 if self.verbose:
173 print("Parsing origin date range...")
174 t0 = time.time()
175 self.origin_date_range = [None, None]
176 for date in xml.xpath(
177 "//tei:sourceDesc//tei:bibl//tei:date|//tei:sourceDesc//tei:biblStruct//tei:date|//tei:sourceDesc//tei:biblFull//tei:date",
178 namespaces={"tei": tei_ns},
179 ):
180 # Try the @when attribute first; if it is set, then it accounts for both ends of the date range:
181 if date.get("when") is not None:
182 self.origin_date_range[0] = int(date.get("when").split("-")[0])
183 self.origin_date_range[1] = self.origin_date_range[0]
184 # Failing that, if it has @from and @to attributes (indicating the period over which the work was completed),
185 # then the completion date of the work accounts for both ends of the date range:
186 elif date.get("to") is not None:
187 self.origin_date_range[0] = int(date.get("to").split("-")[0])
188 self.origin_date_range[1] = self.origin_date_range[0]
189 # Failing that, set lower and upper bounds on the origin date using the the @notBefore and @notAfter attributes:
190 elif date.get("notBefore") is not None or date.get("notAfter") is not None:
191 if date.get("notBefore") is not None:
192 self.origin_date_range[0] = int(date.get("notBefore").split("-")[0])
193 if date.get("notAfter") is not None:
194 self.origin_date_range[1] = int(date.get("notAfter").split("-")[0])
195 return
197 def get_base_wit(self, wit: str):
198 """Given a witness siglum, strips of the specified manuscript suffixes until the siglum matches one in the witness list or until no more suffixes can be stripped.
200 Args:
201 wit: A string representing a witness siglum, potentially including suffixes to be stripped.
202 """
203 base_wit = wit
204 # If our starting siglum corresponds to a siglum in the witness list, then just return it:
205 if base_wit in self.witness_index_by_id:
206 return base_wit
207 # Otherwise, strip any suffixes we find until the siglum corresponds to a base witness in the list
208 # or no more suffixes can be stripped:
209 suffix_found = True
210 while suffix_found:
211 suffix_found = False
212 for suffix in self.manuscript_suffixes:
213 if base_wit.endswith(suffix):
214 suffix_found = True
215 base_wit = base_wit[: -len(suffix)]
216 break # stop looking for other suffixes
217 # If the siglum stripped of this suffix now corresponds to a siglum in the witness list, then return it:
218 if base_wit in self.witness_index_by_id:
219 return base_wit
220 # If we get here, then all possible manuscript suffixes have been stripped, and the resulting siglum does not correspond to a siglum in the witness list:
221 return base_wit
223 def parse_list_wit(self, xml: et.ElementTree):
224 """Given an XML tree for a collation, populates its list of witnesses from its listWit element.
225 If the XML tree does not contain a listWit element, then a ParsingException is thrown listing all distinct witness sigla encountered in the collation.
227 Args:
228 xml: An lxml.etree.ElementTree representing an XML tree rooted at a TEI element.
229 """
230 if self.verbose:
231 print("Parsing witness list...")
232 t0 = time.time()
233 self.witnesses = []
234 self.witness_index_by_id = {}
235 list_wits = xml.xpath("/tei:TEI//tei:listWit", namespaces={"tei": tei_ns})
236 if len(list_wits) == 0:
237 # There is no listWit element: collect all distinct witness sigla in the collation and raise a ParsingException listing them:
238 distinct_sigla = set()
239 sigla = []
240 # Proceed for each rdg, rdgGrp, or witDetail element:
241 for rdg in xml.xpath("//tei:rdg|//tei:rdgGrp|//tei:witDetail", namespaces={"tei": tei_ns}):
242 wit_str = rdg.get("wit") if rdg.get("wit") is not None else ""
243 wits = wit_str.split()
244 for wit in wits:
245 siglum = wit.strip("#") # remove the URI prefix, if there is one
246 if siglum not in distinct_sigla:
247 distinct_sigla.add(siglum)
248 sigla.append(siglum)
249 sigla.sort()
250 msg = ""
251 msg += "An explicit listWit element must be included in the TEI XML collation.\n"
252 msg += "The following sigla occur in the collation and should be included as the @xml:id or @n attributes of witness elements under the listWit element:\n"
253 msg += ", ".join(sigla)
254 raise ParsingException(msg)
255 # Otherwise, take the first listWit element as the list of all witnesses and process it:
256 list_wit = list_wits[0]
257 for witness in list_wit.xpath("./tei:witness", namespaces={"tei": tei_ns}):
258 wit = Witness(witness, self.verbose)
259 self.witness_index_by_id[wit.id] = len(self.witnesses)
260 self.witnesses.append(wit)
261 t1 = time.time()
262 if self.verbose:
263 print("Finished processing %d witnesses in %0.4fs." % (len(self.witnesses), t1 - t0))
264 return
266 def validate_wits(self, xml: et.ElementTree):
267 """Given an XML tree for a collation, checks if any witness sigla listed in a rdg, rdgGrp, or witDetail element,
268 once stripped of ignored suffixes, is not found in the witness list.
269 A warning will be issued for each distinct siglum like this.
270 This method also checks if the upper bound of any witness's date is earlier than the lower bound on the collated work's date of origin
271 and throws an exception if so.
273 Args:
274 xml: An lxml.etree.ElementTree representing an XML tree rooted at a TEI element.
275 """
276 if self.verbose:
277 print("Validating witness list against collation...")
278 t0 = time.time()
279 # There is no listWit element: collect all distinct witness sigla in the collation and raise an exception listing them:
280 distinct_extra_sigla = set()
281 extra_sigla = []
282 # Proceed for each rdg, rdgGrp, or witDetail element:
283 for rdg in xml.xpath("//tei:rdg|//tei:rdgGrp|//tei:witDetail", namespaces={"tei": tei_ns}):
284 wit_str = rdg.get("wit") if rdg.get("wit") is not None else ""
285 wits = wit_str.split()
286 for wit in wits:
287 siglum = wit.strip("#") # remove the URI prefix, if there is one
288 base_siglum = self.get_base_wit(siglum)
289 if base_siglum not in self.witness_index_by_id:
290 if base_siglum not in distinct_extra_sigla:
291 distinct_extra_sigla.add(base_siglum)
292 extra_sigla.append(base_siglum)
293 if len(extra_sigla) > 0:
294 extra_sigla.sort()
295 msg = ""
296 msg += "WARNING: The following sigla occur in the collation that do not have corresponding witness entries in the listWit:\n"
297 msg += ", ".join(extra_sigla)
298 print(msg)
299 # If the lower bound on the date of origin is defined, then check each witness against it:
300 if self.origin_date_range[0] is not None:
301 bad_date_witness_sigla = []
302 bad_date_upper_bounds_by_witness = {}
303 for i, wit in enumerate(self.witnesses):
304 if wit.date_range[1] is not None and wit.date_range[1] < self.origin_date_range[0]:
305 bad_date_witness_sigla.append(wit.id)
306 bad_date_upper_bounds_by_witness[wit.id] = wit.date_range[1]
307 if len(bad_date_witness_sigla) > 0:
308 msg = ""
309 msg += "The following witnesses have their latest possible dates before the earliest date of origin %d specified for the collated work:\n"
310 msg += ", ".join(
311 [
312 (siglum + "(" + str(bad_date_upper_bounds_by_witness[siglum]) + ")")
313 for siglum in bad_date_witness_sigla
314 ]
315 )
316 raise WitnessDateException(msg)
317 t1 = time.time()
318 if self.verbose:
319 print("Finished witness validation in %0.4fs." % (t1 - t0))
320 return
322 def update_witness_date_ranges_from_dates_file(self, dates_file: Union[Path, str]):
323 """Given a CSV-formatted dates file, update the date ranges of all witnesses whose IDs are in the first column of the dates file
324 (overwriting existing date ranges if necessary).
325 """
326 if self.verbose:
327 print("Updating witness dates from file %s..." % (str(dates_file)))
328 t0 = time.time()
329 dates_df = pd.read_csv(dates_file, index_col=0, names=["id", "min", "max"])
330 for witness in self.witnesses:
331 wit_id = witness.id
332 if wit_id in dates_df.index:
333 # For every witness in the list whose ID is specified in the dates file,
334 # update their date ranges (as long as the date ranges in the file are are well-formed):
335 min_date = int(dates_df.loc[wit_id]["min"]) if not np.isnan(dates_df.loc[wit_id]["min"]) else None
336 max_date = (
337 int(dates_df.loc[wit_id]["max"])
338 if not np.isnan(dates_df.loc[wit_id]["max"])
339 else datetime.now().year
340 )
341 if min_date is not None and max_date is not None and min_date > max_date:
342 raise ParsingException(
343 "In dates file %s, for witness ID %s, the minimum date %d is greater than the maximum date %d."
344 % (str(dates_file), wit_id, min_date, max_date)
345 )
346 witness.date_range = [min_date, max_date]
347 t1 = time.time()
348 if self.verbose:
349 print("Finished witness date range updates in %0.4fs." % (t1 - t0))
350 return
352 def update_origin_date_range_from_witness_date_ranges(self):
353 """Conditionally updates the upper bound on the date of origin of the work represented by this Collation
354 based on the bounds on the witnesses' dates.
355 If none of the witnesses have bounds on their dates, then nothing is done.
356 This method is only invoked if the work's date of origin does not already have its upper bound defined.
357 """
358 if self.verbose:
359 print("Updating upper bound on origin date using witness dates...")
360 t0 = time.time()
361 # Set the origin date to the earliest witness date:
362 witness_date_lower_bounds = [wit.date_range[0] for wit in self.witnesses if wit.date_range[0] is not None]
363 witness_date_upper_bounds = [wit.date_range[1] for wit in self.witnesses if wit.date_range[1] is not None]
364 min_witness_date = (
365 min(witness_date_lower_bounds + witness_date_upper_bounds)
366 if len(witness_date_lower_bounds + witness_date_upper_bounds) > 0
367 else None
368 )
369 if min_witness_date is not None:
370 self.origin_date_range[1] = (
371 min(self.origin_date_range[1], min_witness_date)
372 if self.origin_date_range[1] is not None
373 else min_witness_date
374 )
375 t1 = time.time()
376 if self.verbose:
377 print("Finished updating upper bound on origin date in %0.4fs." % (t1 - t0))
378 return
380 def update_witness_date_ranges_from_origin_date_range(self):
381 """Attempts to update the lower bounds on the witnesses' dates of origin of the work represented by this Collation
382 using the upper bound on the date of origin of the work represented by this Collation.
383 This method is only invoked if the upper bound on the work's date of origin was not already defined
384 (i.e., if update_origin_date_range_from_witness_date_ranges was not invoked).
385 """
386 if self.verbose:
387 print("Updating lower bounds on witness dates using origin date...")
388 t0 = time.time()
389 # Proceed for every witness:
390 for i, wit in enumerate(self.witnesses):
391 # Ensure that the lower bound on this witness's date is no earlier than the upper bound on the date of the work's origin:
392 wit.date_range[0] = (
393 max(wit.date_range[0], self.origin_date_range[1])
394 if wit.date_range[0] is not None
395 else self.origin_date_range[1]
396 )
397 # Then ensure that the upper bound on this witness's date is no earlier than its lower bound, in case we updated it:
398 wit.date_range[1] = max(wit.date_range[0], wit.date_range[1])
399 t1 = time.time()
400 if self.verbose:
401 print("Finished updating lower bounds on witness dates in %0.4fs." % (t1 - t0))
402 return
404 def parse_weights(self, xml: et.ElementTree):
405 """Given an XML tree for a collation, populates this Collation's list of variation unit weight categories
406 (associated with types of variation that may have different expected frequencies)
407 and its dictionary mapping these categories to integer weights.
409 Args:
410 xml: An lxml.etree.ElementTree representing an XML tree rooted at a TEI element.
411 """
412 if self.verbose:
413 print("Parsing variation unit weight categories...")
414 t0 = time.time()
415 self.weight_categories = []
416 self.weights_by_id = {}
417 for interp in xml.xpath("//tei:interpGrp[@type=\"weight\"]/tei:interp", namespaces={"tei": tei_ns}):
418 # These must be indexed by the xml:id attribute, so skip any that do not have one:
419 if interp.get("{%s}id" % xml_ns) is None:
420 continue
421 weight_category = interp.get("{%s}id" % xml_ns)
422 # Retrieve this element's text value and coerce it to an integer, defaulting to 1 if it has none:
423 weight = 1
424 for certainty in interp.xpath("./tei:certainty", namespaces={"tei": tei_ns}):
425 if certainty.get("degree") is not None:
426 weight = int(certainty.get("degree"))
427 break
428 self.weight_categories.append(weight_category)
429 self.weights_by_id[weight_category] = weight
430 t1 = time.time()
431 if self.verbose:
432 print(
433 "Finished processing %d variation unit weight categories in %0.4fs."
434 % (len(self.weight_categories), t1 - t0)
435 )
436 return
438 def parse_intrinsic_odds(self, xml: et.ElementTree):
439 """Given an XML tree for a collation, populates this Collation's list of intrinsic probability categories
440 (e.g., "absolutely more likely," "highly more likely," "more likely," "slightly more likely," "equally likely")
441 and its dictionary mapping these categories to numerical odds.
442 If a category does not contain a certainty element specifying its number, then it will be assumed to be a parameter to be estimated.
444 Args:
445 xml: An lxml.etree.ElementTree representing an XML tree rooted at a TEI element.
446 """
447 if self.verbose:
448 print("Parsing intrinsic odds categories...")
449 t0 = time.time()
450 self.intrinsic_categories = []
451 self.intrinsic_odds_by_id = {}
452 for interp in xml.xpath("//tei:interpGrp[@type=\"intrinsic\"]/tei:interp", namespaces={"tei": tei_ns}):
453 # These must be indexed by the xml:id attribute, so skip any that do not have one:
454 if interp.get("{%s}id" % xml_ns) is None:
455 continue
456 odds_category = interp.get("{%s}id" % xml_ns)
457 # If this element contains a certainty subelement with a fixed odds value for this category, then set it:
458 odds = None
459 for certainty in interp.xpath("./tei:certainty", namespaces={"tei": tei_ns}):
460 if certainty.get("degree") is not None:
461 odds = float(certainty.get("degree"))
462 break
463 self.intrinsic_categories.append(odds_category)
464 self.intrinsic_odds_by_id[odds_category] = odds
465 t1 = time.time()
466 if self.verbose:
467 print(
468 "Finished processing %d intrinsic odds categories in %0.4fs."
469 % (len(self.intrinsic_categories), t1 - t0)
470 )
471 return
473 def parse_transcriptional_rates(self, xml: et.ElementTree):
474 """Given an XML tree for a collation, populates this Collation's dictionary mapping transcriptional change categories
475 (e.g., "aural confusion," "visual error," "clarification") to numerical rates.
476 If a category does not contain a certainty element specifying its number, then it will be assumed to be a parameter to be estimated.
478 Args:
479 xml: An lxml.etree.ElementTree representing an XML tree rooted at a TEI element.
480 """
481 if self.verbose:
482 print("Parsing transcriptional change categories...")
483 t0 = time.time()
484 self.transcriptional_categories = []
485 self.transcriptional_rates_by_id = {}
486 for interp in xml.xpath("//tei:interpGrp[@type=\"transcriptional\"]/tei:interp", namespaces={"tei": tei_ns}):
487 # These must be indexed by the xml:id attribute, so skip any that do not have one:
488 if interp.get("{%s}id" % xml_ns) is None:
489 continue
490 transcriptional_category = interp.get("{%s}id" % xml_ns)
491 # If this element contains a certainty subelement with a fixed rate for this category, then set it:
492 rate = None
493 for certainty in interp.xpath("./tei:certainty", namespaces={"tei": tei_ns}):
494 if certainty.get("degree") is not None:
495 rate = float(certainty.get("degree"))
496 break
497 self.transcriptional_categories.append(transcriptional_category)
498 self.transcriptional_rates_by_id[transcriptional_category] = rate
499 t1 = time.time()
500 if self.verbose:
501 print(
502 "Finished processing %d transcriptional change categories in %0.4fs."
503 % (len(self.transcriptional_rates_by_id), t1 - t0)
504 )
505 return
507 def validate_intrinsic_relations(self):
508 """Checks if any VariationUnit's intrinsic_relations map is not a forest.
509 If any is not, then an IntrinsicRelationsException is thrown describing the VariationUnit at fault.
510 """
511 if self.verbose:
512 print("Validating intrinsic relation graphs for variation units...")
513 t0 = time.time()
514 for vu in self.variation_units:
515 # Skip any variation units with an empty intrinsic_relations map:
516 if len(vu.intrinsic_relations) == 0:
517 continue
518 # For all others, start by identifying all reading IDs that are not related to by some other reading ID:
519 in_degree_by_reading = {}
520 for edge in vu.intrinsic_relations:
521 s = edge[0]
522 t = edge[1]
523 if s not in in_degree_by_reading:
524 in_degree_by_reading[s] = 0
525 if t not in in_degree_by_reading:
526 in_degree_by_reading[t] = 0
527 in_degree_by_reading[t] += 1
528 # If any reading has more than one relation pointing to it, then the intrinsic relations graph is not a forest:
529 excessive_in_degree_readings = [
530 rdg_id for rdg_id in in_degree_by_reading if in_degree_by_reading[rdg_id] > 1
531 ]
532 if len(excessive_in_degree_readings) > 0:
533 msg = ""
534 msg += (
535 "In variation unit %s, the following readings have more than one intrinsic relation pointing to them: %s.\n"
536 % (vu.id, ", ".join(excessive_in_degree_readings))
537 )
538 msg += "Please ensure that at least one reading has no relations pointing to it and that every reading has no more than one relation pointing to it."
539 raise IntrinsicRelationsException(msg)
540 # If every reading has another reading pointing to it, then the intrinsic relations graph contains a cycle and is not a forest:
541 starting_nodes = [rdg_id for rdg_id in in_degree_by_reading if in_degree_by_reading[rdg_id] == 0]
542 if len(starting_nodes) == 0:
543 msg = ""
544 msg += "In variation unit %s, the intrinsic relations contain a cycle.\n" % vu.id
545 msg += "Please ensure that at least one reading has no relations pointing to it and that every reading has no more than one relation pointing to it."
546 raise IntrinsicRelationsException(msg)
547 t1 = time.time()
548 if self.verbose:
549 print("Finished intrinsic relations validation in %0.4fs." % (t1 - t0))
550 return
552 def parse_apps(self, xml: et.ElementTree):
553 """Given an XML tree for a collation, populates its list of variation units from its app elements.
555 Args:
556 xml: An lxml.etree.ElementTree representing an XML tree rooted at a TEI element.
557 """
558 if self.verbose:
559 print("Parsing variation units...")
560 t0 = time.time()
561 for a in xml.xpath('//tei:app', namespaces={'tei': tei_ns}):
562 vu = VariationUnit(a, self.verbose)
563 self.variation_units.append(vu)
564 t1 = time.time()
565 if self.verbose:
566 print("Finished processing %d variation units in %0.4fs." % (len(self.variation_units), t1 - t0))
567 return
569 def get_readings_by_witness_for_unit(self, vu: VariationUnit):
570 """Returns a dictionary mapping witness IDs to a list of their reading coefficients for a given variation unit.
572 Args:
573 vu: A VariationUnit to be processed.
575 Returns:
576 A dictionary mapping witness ID strings to a list of their coefficients for all substantive readings in this VariationUnit.
577 """
578 # In a first pass, populate lists of substantive (variation unit ID, reading ID) tuples and reading labels
579 # and a map from reading IDs to the indices of their parent substantive reading in this unit:
580 reading_id_to_index = {}
581 self.substantive_readings_by_variation_unit_id[vu.id] = []
582 for rdg in vu.readings:
583 # If this reading is missing (e.g., lacunose or inapplicable due to an overlapping variant) or targets another reading, then skip it:
584 if rdg.type in self.missing_reading_types or len(rdg.certainties) > 0:
585 continue
586 # If this reading is trivial, then map it to the last substantive index:
587 if rdg.type in self.trivial_reading_types:
588 reading_id_to_index[rdg.id] = len(self.substantive_readings_by_variation_unit_id[vu.id]) - 1
589 continue
590 # Otherwise, the reading is substantive: add it to the map and update the last substantive index:
591 self.substantive_readings_by_variation_unit_id[vu.id].append(rdg.id)
592 self.substantive_variation_unit_reading_tuples.append(tuple([vu.id, rdg.id]))
593 reading_id_to_index[rdg.id] = len(self.substantive_readings_by_variation_unit_id[vu.id]) - 1
594 if len(self.substantive_readings_by_variation_unit_id[vu.id]) == 0:
595 raise ValueError(f"Variation unit {vu.id} has no substantive readings.")
596 # If the list of substantive readings only contains one entry, then this variation unit is not informative;
597 # return an empty dictionary and add nothing to the list of substantive reading labels:
598 if self.verbose:
599 print(
600 "Variation unit %s has %d substantive readings."
601 % (vu.id, len(self.substantive_readings_by_variation_unit_id[vu.id]))
602 )
603 readings_by_witness_for_unit = {}
604 # Initialize the output dictionary with empty sets for all base witnesses:
605 for wit in self.witnesses:
606 readings_by_witness_for_unit[wit.id] = [0] * len(self.substantive_readings_by_variation_unit_id[vu.id])
607 # In a second pass, assign each base witness a set containing the readings it supports in this unit:
608 for rdg in vu.readings:
609 # Initialize the dictionary indicating support for this reading (or its disambiguations):
610 rdg_support = [0] * len(self.substantive_readings_by_variation_unit_id[vu.id])
611 # If this is a missing reading (e.g., a lacuna or an overlap), then we can skip it, as its corresponding set will be empty:
612 if rdg.type in self.missing_reading_types:
613 continue
614 # Otherwise, if this reading is trivial, then it will contain an entry for the index of its parent substantive reading:
615 elif rdg.type in self.trivial_reading_types:
616 rdg_support[reading_id_to_index[rdg.id]] = 1
617 # Otherwise, if this reading has one or more nonzero certainty degrees,
618 # then set the entries for these readings to their degrees:
619 elif sum(rdg.certainties.values()) > 0:
620 for t in rdg.certainties:
621 # Skip any reading whose ID is unrecognized in this unit:
622 if t in reading_id_to_index:
623 rdg_support[reading_id_to_index[t]] = rdg.certainties[t]
624 # Otherwise, if this reading has one or more targets (i.e., if it is an ambiguous reading),
625 # then set the entries for each of its targets to 1:
626 elif len(rdg.targets) > 0:
627 for t in rdg.targets:
628 # Skip any reading whose ID is unrecognized in this unit:
629 if t in reading_id_to_index:
630 rdg_support[reading_id_to_index[t]] = 1
631 # Otherwise, this reading is itself substantive; set the entry for the index of this reading to 1:
632 else:
633 rdg_support[reading_id_to_index[rdg.id]] = 1
634 # Proceed for each witness siglum in the support for this reading:
635 for wit in rdg.wits:
636 # Is this siglum a base siglum?
637 base_wit = self.get_base_wit(wit)
638 if base_wit not in self.witness_index_by_id:
639 # If it is not, then it is probably just because we've encountered a corrector or some other secondary witness not included in the witness list;
640 # report this if we're in verbose mode and move on:
641 if self.verbose:
642 print(
643 "Skipping unknown witness siglum %s (base siglum %s) in variation unit %s, reading %s..."
644 % (wit, base_wit, vu.id, rdg.id)
645 )
646 continue
647 # If we've found a base siglum, then add this reading's contribution to the base witness's reading set for this unit;
648 # normally the existing set will be empty, but if we reduce two suffixed sigla to the same base witness,
649 # then that witness may attest to multiple readings in the same unit:
650 readings_by_witness_for_unit[base_wit] = [
651 (min(readings_by_witness_for_unit[base_wit][i] + rdg_support[i], 1))
652 for i in range(len(rdg_support))
653 ]
654 return readings_by_witness_for_unit
656 def parse_readings_by_witness(self):
657 """Populates the internal dictionary mapping witness IDs to a list of their reading support sets for all variation units, and then fills the empty reading support sets for witnesses of type "corrector" with the entries of the previous witness."""
658 if self.verbose:
659 print("Populating internal dictionary of witness readings...")
660 t0 = time.time()
661 # Initialize the data structures to be populated here:
662 self.readings_by_witness = {}
663 self.variation_unit_ids = []
664 for wit in self.witnesses:
665 self.readings_by_witness[wit.id] = []
666 # Populate them for each variation unit:
667 for vu in self.variation_units:
668 readings_by_witness_for_unit = self.get_readings_by_witness_for_unit(vu)
669 if len(readings_by_witness_for_unit) > 0:
670 self.variation_unit_ids.append(vu.id)
671 for wit in readings_by_witness_for_unit:
672 self.readings_by_witness[wit].append(readings_by_witness_for_unit[wit])
673 # Optionally, fill the lacunae of the correctors:
674 if self.fill_corrector_lacunae:
675 filled_readings = []
676 for i, wit in enumerate(self.witnesses):
677 # If this is the first witness, then it shouldn't be a corrector (since there is no previous witness against which to compare it):
678 if i == 0:
679 filled_readings = list(self.readings_by_witness[wit.id])
680 continue
681 # Otherwise, if this witness is not a corrector, then skip it:
682 if wit.type != "corrector":
683 filled_readings = list(self.readings_by_witness[wit.id])
684 continue
685 # Otherwise, add this corrector's extant readings to the filled readings list:
686 for j in range(len(self.readings_by_witness[wit.id])):
687 if sum(self.readings_by_witness[wit.id][j]) != 0:
688 filled_readings[j] = list(self.readings_by_witness[wit.id][j])
689 # If a threshold of extant readings is specified, then check if this corrector meets it, and skip it if not:
690 if self.fill_correctors_threshold is not None:
691 # If there is a threshold, then first check the proportion of variation units at which this witness is not lacunose:
692 proportion_extant = sum(
693 [
694 1
695 for j in range(len(self.readings_by_witness[wit.id]))
696 if sum(self.readings_by_witness[wit.id][j]) != 0
697 ]
698 ) / len(self.readings_by_witness[wit.id])
699 # If this corrector does not exceed the threshold, then don't fill it,
700 # but do update the running set of readings for the previous corrector and first hand:
701 if proportion_extant < self.fill_correctors_threshold:
702 continue
703 # Otherwise, fill every lacuna in this corrector based on the filled readings list:
704 for j in range(len(self.readings_by_witness[wit.id])):
705 if sum(self.readings_by_witness[wit.id][j]) == 0:
706 self.readings_by_witness[wit.id][j] = list(filled_readings[j])
707 t1 = time.time()
708 if self.verbose:
709 print(
710 "Populated dictionary for %d witnesses over %d substantive variation units in %0.4fs."
711 % (len(self.witnesses), len(self.variation_unit_ids), t1 - t0)
712 )
713 return
715 def filter_fragmentary_witnesses(self, xml):
716 """Filters the original witness list and readings by witness dictionary to exclude witnesses whose proportions of extant passages fall below the fragmentary readings threshold."""
717 if self.verbose:
718 print(
719 "Filtering fragmentary witnesses (extant in < %f of all variation units) out of internal witness list and dictionary of witness readings..."
720 % self.fragmentary_threshold
721 )
722 t0 = time.time()
723 fragmentary_witness_set = set()
724 # Proceed for each witness in order:
725 for wit in self.witnesses:
726 wit_id = wit.id
727 # We count the number of variation units at which this witness has an extant (i.e., non-missing) reading:
728 extant_reading_count = 0
729 total_reading_count = len(self.readings_by_witness[wit.id])
730 # Proceed through all reading support lists:
731 for rdg_support in self.readings_by_witness[wit_id]:
732 # If the current reading support list is not all zeroes, then increment this witness's count of extant readings:
733 if sum(rdg_support) != 0:
734 extant_reading_count += 1
735 # If the proportion of extant readings falls below the threshold, then add this witness to the list of fragmentary witnesses:
736 if extant_reading_count / total_reading_count < self.fragmentary_threshold:
737 fragmentary_witness_set.add(wit_id)
738 # Then filter the witness list to exclude the fragmentary witnesses:
739 filtered_witnesses = [wit for wit in self.witnesses if wit.id not in fragmentary_witness_set]
740 self.witnesses = filtered_witnesses
741 # Then remove the entries for the fragmentary witnesses from the witnesses-to-readings dictionary:
742 for wit_id in fragmentary_witness_set:
743 del self.readings_by_witness[wit_id]
744 t1 = time.time()
745 if self.verbose:
746 print(
747 "Filtered out %d fragmentary witness(es) (%s) in %0.4fs."
748 % (len(fragmentary_witness_set), str(list(fragmentary_witness_set)), t1 - t0)
749 )
750 return
752 def get_nexus_symbols(self):
753 """Returns a list of one-character symbols needed to represent the states of all substantive readings in NEXUS.
755 The number of symbols equals the maximum number of substantive readings at any variation unit.
757 Returns:
758 A list of individual characters representing states in readings.
759 """
760 # NOTE: IQTREE does not appear to support symbols outside of 0-9 and a-z, and its base symbols must be case-insensitive.
761 # The official version of MrBayes is likewise limited to 32 symbols.
762 # But PAUP* allows up to 64 symbols, and Andrew Edmondson's fork of MrBayes does, as well.
763 # So this method will support symbols from 0-9, a-z, and A-Z (for a total of 62 states)
764 possible_symbols = list(string.digits) + list(string.ascii_lowercase) + list(string.ascii_uppercase)
765 # The number of symbols needed is equal to the length of the longest substantive reading vector:
766 nsymbols = 0
767 # If there are no witnesses, then no symbols are needed at all:
768 if len(self.witnesses) == 0:
769 return []
770 wit_id = self.witnesses[0].id
771 for i, vu_id in enumerate(self.variation_unit_ids):
772 rdg_support = self.readings_by_witness[wit_id][i]
773 if len(rdg_support) > len(possible_symbols):
774 raise ValueError(
775 f"ERROR: too many substantive readings at variation unit '{vu_id}' to represent in NEXUS format."
776 )
777 nsymbols = max(nsymbols, len(rdg_support))
778 nexus_symbols = possible_symbols[:nsymbols]
779 return nexus_symbols
781 def to_nexus(
782 self,
783 file_addr: Union[Path, str],
784 drop_constant: bool = False,
785 char_state_labels: bool = True,
786 frequency: bool = False,
787 ambiguous_as_missing: bool = False,
788 calibrate_dates: bool = False,
789 mrbayes: bool = False,
790 clock_model: ClockModel = ClockModel.strict,
791 ):
792 """Writes this Collation to a NEXUS file with the given address.
794 Args:
795 file_addr: A string representing the path to an output NEXUS file; the file type should be .nex, .nexus, or .nxs.
796 drop_constant: An optional flag indicating whether to ignore variation units with one substantive reading.
797 char_state_labels: An optional flag indicating whether or not to include the CharStateLabels block.
798 frequency: An optional flag indicating whether to use the StatesFormat=Frequency setting
799 instead of the StatesFormat=StatesPresent setting
800 (and thus represent all states with frequency vectors rather than symbols).
801 Note that this setting is necessary to make use of certainty degrees assigned to multiple ambiguous states in the collation.
802 ambiguous_as_missing: An optional flag indicating whether to treat all ambiguous states as missing data.
803 If this flag is set, then only base symbols will be generated for the NEXUS file.
804 It is only applied if the frequency option is False.
805 calibrate_dates: An optional flag indicating whether to add an Assumptions block that specifies date distributions for witnesses.
806 This option is intended for inputs to BEAST 2.
807 mrbayes: An optional flag indicating whether to add a MrBayes block that specifies model settings and age calibrations for witnesses.
808 This option is intended for inputs to MrBayes.
809 clock_model: A ClockModel option indicating which type of clock model to use.
810 This option is intended for inputs to MrBayes and BEAST 2.
811 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.
812 """
813 # Populate a list of sites that will correspond to columns of the sequence alignment:
814 substantive_variation_unit_ids = self.variation_unit_ids
815 if drop_constant:
816 substantive_variation_unit_ids = [
817 vu_id
818 for vu_id in self.variation_unit_ids
819 if len(self.substantive_readings_by_variation_unit_id[vu_id]) > 1
820 ]
821 substantive_variation_unit_ids_set = set(substantive_variation_unit_ids)
822 substantive_variation_unit_reading_tuples_set = set(self.substantive_variation_unit_reading_tuples)
823 # Start by calculating the values we will be using here:
824 ntax = len(self.witnesses)
825 nchar = len(substantive_variation_unit_ids)
826 taxlabels = [slugify(wit.id, lowercase=False, separator='_') for wit in self.witnesses]
827 max_taxlabel_length = max(
828 [len(taxlabel) for taxlabel in taxlabels]
829 ) # keep track of the longest taxon label for tabular alignment purposes
830 charlabels = [slugify(vu_id, lowercase=False, separator='_') for vu_id in substantive_variation_unit_ids]
831 missing_symbol = '?'
832 symbols = self.get_nexus_symbols()
833 # Generate all parent folders for this file that don't already exist:
834 Path(file_addr).parent.mkdir(parents=True, exist_ok=True)
835 # Then write the file:
836 pbar = tqdm()
837 with open(file_addr, "w", encoding="utf-8") as f:
838 # Start with the NEXUS header:
839 f.write("#NEXUS\n\n")
840 # Then begin the data block:
841 f.write("Begin DATA;\n")
842 # Write the collation matrix dimensions:
843 f.write("\tDimensions ntax=%d nchar=%d;\n" % (ntax, nchar))
844 # Write the format subblock:
845 f.write("\tFormat\n")
846 f.write("\t\tDataType=Standard\n")
847 f.write("\t\tMissing=%s\n" % missing_symbol)
848 if frequency:
849 f.write("\t\tStatesFormat=Frequency\n")
850 f.write("\t\tSymbols=\"%s\";\n" % (" ".join(symbols)))
851 # If the char_state_labels is set, then write the labels for character-state labels, with each on its own line:
852 if char_state_labels:
853 f.write("\tCharStateLabels")
854 vu_ind = 1
855 for vu in self.variation_units:
856 if vu.id not in substantive_variation_unit_ids_set:
857 continue
858 if vu_ind == 1:
859 f.write("\n\t\t%d %s /" % (vu_ind, slugify(vu.id, lowercase=False, separator='_')))
860 else:
861 f.write(",\n\t\t%d %s /" % (vu_ind, slugify(vu.id, lowercase=False, separator='_')))
862 rdg_ind = 0
863 for rdg in vu.readings:
864 key = tuple([vu.id, rdg.id])
865 if key not in substantive_variation_unit_reading_tuples_set:
866 continue
867 ascii_rdg_text = slugify(
868 rdg.text, lowercase=False, separator='_', replacements=[['η', 'h'], ['ω', 'w']]
869 )
870 if ascii_rdg_text == "":
871 ascii_rdg_text = "om."
872 f.write(" %s" % ascii_rdg_text)
873 rdg_ind += 1
874 if rdg_ind > 0:
875 vu_ind += 1
876 f.write(";\n")
877 # Write the matrix subblock:
878 f.write("\tMatrix")
879 with tqdm(total=len(self.witnesses)) as pbar:
880 for i, wit in enumerate(self.witnesses):
881 taxlabel = taxlabels[i]
882 if frequency:
883 sequence = "\n\t\t" + taxlabel
884 for j, vu_id in enumerate(self.variation_unit_ids):
885 if vu_id not in substantive_variation_unit_ids_set:
886 continue
887 rdg_support = self.readings_by_witness[wit.id][j]
888 sequence += "\n\t\t\t"
889 # If this reading is lacunose in this witness, then use the missing character:
890 if sum(rdg_support) == 0:
891 sequence += missing_symbol
892 continue
893 # Otherwise, print out its frequencies for different readings in parentheses:
894 sequence += "("
895 for k, w in enumerate(rdg_support):
896 sequence += "%s:%0.4f" % (symbols[k], w)
897 if k < len(rdg_support) - 1:
898 sequence += " "
899 sequence += ")"
900 else:
901 sequence = "\n\t\t" + taxlabel
902 # Add enough space after this label ensure that all sequences are nicely aligned:
903 sequence += " " * (max_taxlabel_length - len(taxlabel) + 1)
904 for j, vu_id in enumerate(self.variation_unit_ids):
905 if vu_id not in substantive_variation_unit_ids_set:
906 continue
907 rdg_support = self.readings_by_witness[wit.id][j]
908 # If this reading is lacunose in this witness, then use the missing character:
909 if sum(rdg_support) == 0:
910 sequence += missing_symbol
911 continue
912 rdg_inds = [
913 k for k, w in enumerate(rdg_support) if w > 0
914 ] # the index list consists of the indices of all readings with any degree of certainty assigned to them
915 # For singleton readings, just print the symbol:
916 if len(rdg_inds) == 1:
917 sequence += symbols[rdg_inds[0]]
918 continue
919 # For multiple readings, print the corresponding readings in braces or the missing symbol depending on input settings:
920 if ambiguous_as_missing:
921 sequence += missing_symbol
922 else:
923 sequence += "{%s}" % "".join([str(rdg_ind) for rdg_ind in rdg_inds])
924 f.write("%s" % (sequence))
925 pbar.update(1)
926 f.write(";\n")
927 # End the data block:
928 f.write("End;")
929 # If calibrate_dates is set, then add the assumptions block:
930 if calibrate_dates:
931 f.write("\n\n")
932 f.write("Begin ASSUMPTIONS;\n")
933 # Set the scale to years:
934 f.write("\tOPTIONS SCALE = years;\n\n")
935 # Then calibrate the witness ages:
936 calibrate_strings = []
937 for i, wit in enumerate(self.witnesses):
938 taxlabel = taxlabels[i]
939 date_range = wit.date_range
940 if date_range[0] is not None:
941 # If there is a lower bound on the witness's date, then use either a fixed or uniform distribution,
942 # depending on whether the upper and lower bounds match:
943 min_age = datetime.now().year - date_range[1]
944 max_age = datetime.now().year - date_range[0]
945 if min_age == max_age:
946 calibrate_string = "\tCALIBRATE %s = fixed(%d)" % (taxlabel, min_age)
947 calibrate_strings.append(calibrate_string)
948 else:
949 calibrate_string = "\tCALIBRATE %s = uniform(%d,%d)" % (taxlabel, min_age, max_age)
950 calibrate_strings.append(calibrate_string)
951 else:
952 # If there is no lower bound on the witness's date, then use an offset log-normal distribution:
953 min_age = datetime.now().year - date_range[1]
954 calibrate_string = "\tCALIBRATE %s = offsetlognormal(%d,0.0,1.0)" % (taxlabel, min_age)
955 calibrate_strings.append(calibrate_string)
956 # Then print the calibrate strings, separated by commas and line breaks and terminated by a semicolon:
957 f.write("%s;\n\n" % ",\n".join(calibrate_strings))
958 # End the assumptions block:
959 f.write("End;")
960 # If mrbayes is set, then add the mrbayes block:
961 if mrbayes:
962 f.write("\n\n")
963 f.write("Begin MRBAYES;\n")
964 # Turn on the autoclose feature by default:
965 f.write("\tset autoclose=yes;\n")
966 # Set the branch lengths to be governed by a birth-death clock model, and set up the parameters for this model:
967 f.write("\n")
968 f.write("\tprset brlenspr = clock:birthdeath;\n")
969 f.write("\tprset speciationpr = uniform(0.0,10.0);\n")
970 f.write("\tprset extinctionpr = beta(2.0,4.0);\n")
971 f.write("\tprset sampleprob = 0.01;\n")
972 # Use the specified clock model:
973 f.write("\n")
974 if clock_model == clock_model.uncorrelated:
975 f.write("\tprset clockvarpr=igr;\n")
976 f.write("\tprset clockratepr=lognormal(0.0,1.0);\n")
977 f.write("\tprset igrvarpr=exponential(1.0);\n")
978 else:
979 f.write("\tprset clockvarpr=strict;\n")
980 f.write("\tprset clockratepr=lognormal(0.0,1.0);\n")
981 # Set the priors on the tree age depending on the date range for the origin of the collated work:
982 f.write("\n")
983 if self.origin_date_range[0] is not None:
984 min_tree_age = (
985 datetime.now().year - self.origin_date_range[1]
986 if self.origin_date_range[1] is not None
987 else 0.0
988 )
989 max_tree_age = datetime.now().year - self.origin_date_range[0]
990 f.write("\tprset treeagepr = uniform(%d,%d);\n" % (min_tree_age, max_tree_age))
991 else:
992 min_tree_age = (
993 datetime.now().year - self.origin_date_range[1]
994 if self.origin_date_range[1] is not None
995 else 0.0
996 )
997 f.write("\tprset treeagepr = offsetgamma(%d,1.0,1.0);\n" % (min_tree_age))
998 # Then calibrate the witness ages:
999 f.write("\n")
1000 f.write("\tprset nodeagepr = calibrated;\n")
1001 for i, wit in enumerate(self.witnesses):
1002 taxlabel = taxlabels[i]
1003 date_range = wit.date_range
1004 if date_range[0] is not None:
1005 # If there is a lower bound on the witness's date, then use either a fixed or uniform distribution,
1006 # depending on whether the upper and lower bounds match:
1007 min_age = datetime.now().year - date_range[1]
1008 max_age = datetime.now().year - date_range[0]
1009 if min_age == max_age:
1010 f.write("\tcalibrate %s = fixed(%d);\n" % (taxlabel, min_age))
1011 else:
1012 f.write("\tcalibrate %s = uniform(%d,%d);\n" % (taxlabel, min_age, max_age))
1013 else:
1014 # If there is no lower bound on the witness's date, then use an offset gamma distribution:
1015 min_age = datetime.now().year - date_range[1]
1016 f.write("\tcalibrate %s = offsetgamma(%d,1.0,1.0);\n" % (taxlabel, min_age))
1017 f.write("\n")
1018 # Add default settings for MCMC estimation of posterior distribution:
1019 f.write("\tmcmcp ngen=100000;\n")
1020 # Write the command to run MrBayes:
1021 f.write("\tmcmc;\n")
1022 # End the assumptions block:
1023 f.write("End;")
1024 return
1026 def get_hennig86_symbols(self):
1027 """Returns a list of one-character symbols needed to represent the states of all substantive readings in Hennig86 format.
1029 The number of symbols equals the maximum number of substantive readings at any variation unit.
1031 Returns:
1032 A list of individual characters representing states in readings.
1033 """
1034 possible_symbols = (
1035 list(string.digits) + list(string.ascii_uppercase)[:22]
1036 ) # NOTE: the maximum number of symbols allowed in Hennig86 format is 32
1037 # The number of symbols needed is equal to the length of the longest substantive reading vector:
1038 nsymbols = 0
1039 # If there are no witnesses, then no symbols are needed at all:
1040 if len(self.witnesses) == 0:
1041 return []
1042 wit_id = self.witnesses[0].id
1043 for i, vu_id in enumerate(self.variation_unit_ids):
1044 rdg_support = self.readings_by_witness[wit_id][i]
1045 if len(rdg_support) > len(possible_symbols):
1046 raise ValueError(
1047 f"ERROR: too many substantive readings at variation unit '{vu_id}' to represent in Hennig86 format."
1048 )
1049 nsymbols = max(nsymbols, len(rdg_support))
1050 hennig86_symbols = possible_symbols[:nsymbols]
1051 return hennig86_symbols
1053 def to_hennig86(self, file_addr: Union[Path, str], drop_constant: bool = False):
1054 """Writes this Collation to a file in Hennig86 format with the given address.
1055 Note that because Hennig86 format does not support NEXUS-style ambiguities, such ambiguities will be treated as missing data.
1057 Args:
1058 file_addr: A string representing the path to an output file.
1059 drop_constant (bool, optional): An optional flag indicating whether to ignore variation units with one substantive reading.
1060 """
1061 # Populate a list of sites that will correspond to columns of the sequence alignment:
1062 substantive_variation_unit_ids = self.variation_unit_ids
1063 if drop_constant:
1064 substantive_variation_unit_ids = [
1065 vu_id
1066 for vu_id in self.variation_unit_ids
1067 if len(self.substantive_readings_by_variation_unit_id[vu_id]) > 1
1068 ]
1069 substantive_variation_unit_ids_set = set(substantive_variation_unit_ids)
1070 substantive_variation_unit_reading_tuples_set = set(self.substantive_variation_unit_reading_tuples)
1071 # Start by calculating the values we will be using here:
1072 ntax = len(self.witnesses)
1073 nchar = len(substantive_variation_unit_ids)
1074 taxlabels = []
1075 for wit in self.witnesses:
1076 taxlabel = wit.id
1077 # Hennig86 format requires taxon names to start with a letter, so if this is not the case, then append "WIT_" to the start of the name:
1078 if taxlabel[0] not in string.ascii_letters:
1079 taxlabel = "WIT_" + taxlabel
1080 # Then replace any disallowed characters in the string with an underscore:
1081 taxlabel = slugify(taxlabel, lowercase=False, separator='_')
1082 taxlabels.append(taxlabel)
1083 max_taxlabel_length = max(
1084 [len(taxlabel) for taxlabel in taxlabels]
1085 ) # keep track of the longest taxon label for tabular alignment purposes
1086 missing_symbol = '?'
1087 symbols = self.get_hennig86_symbols()
1088 # Generate all parent folders for this file that don't already exist:
1089 Path(file_addr).parent.mkdir(parents=True, exist_ok=True)
1090 with open(file_addr, "w", encoding="ascii") as f:
1091 # Start with the nstates header:
1092 f.write("nstates %d;\n" % len(symbols))
1093 # Then begin the xread block:
1094 f.write("xread\n")
1095 # Write the dimensions:
1096 f.write("%d %d\n" % (nchar, ntax))
1097 # Now write the matrix:
1098 with tqdm(total=len(self.witnesses)) as pbar:
1099 for i, wit in enumerate(self.witnesses):
1100 taxlabel = taxlabels[i]
1101 # Add enough space after this label ensure that all sequences are nicely aligned:
1102 sequence = taxlabel + (" " * (max_taxlabel_length - len(taxlabel) + 1))
1103 for j, vu_id in enumerate(self.variation_unit_ids):
1104 if vu_id not in substantive_variation_unit_ids_set:
1105 continue
1106 rdg_support = self.readings_by_witness[wit.id][j]
1107 # If this reading is lacunose in this witness, then use the missing character:
1108 if sum(rdg_support) == 0:
1109 sequence += missing_symbol
1110 continue
1111 rdg_inds = [
1112 k for k, w in enumerate(rdg_support) if w > 0
1113 ] # the index list consists of the indices of all readings with any degree of certainty assigned to them
1114 # For singleton readings, just print the symbol:
1115 if len(rdg_inds) == 1:
1116 sequence += symbols[rdg_inds[0]]
1117 continue
1118 # For multiple readings, print the missing symbol:
1119 sequence += missing_symbol
1120 f.write("%s\n" % (sequence))
1121 pbar.update(1)
1122 f.write(";")
1123 return
1125 def get_phylip_symbols(self):
1126 """Returns a list of one-character symbols needed to represent the states of all substantive readings in PHYLIP format.
1128 The number of symbols equals the maximum number of substantive readings at any variation unit.
1130 Returns:
1131 A list of individual characters representing states in readings.
1132 """
1133 possible_symbols = (
1134 list(string.digits) + list(string.ascii_lowercase)[:22]
1135 ) # NOTE: for RAxML, multistate characters with an alphabet sizes up to 32 are supported
1136 # The number of symbols needed is equal to the length of the longest substantive reading vector:
1137 nsymbols = 0
1138 # If there are no witnesses, then no symbols are needed at all:
1139 if len(self.witnesses) == 0:
1140 return []
1141 wit_id = self.witnesses[0].id
1142 for i, vu_id in enumerate(self.variation_unit_ids):
1143 rdg_support = self.readings_by_witness[wit_id][i]
1144 if len(rdg_support) > len(possible_symbols):
1145 raise ValueError(
1146 f"ERROR: too many substantive readings at variation unit '{vu_id}' to represent in PHYLIP format."
1147 )
1148 nsymbols = max(nsymbols, len(rdg_support))
1149 phylip_symbols = possible_symbols[:nsymbols]
1150 return phylip_symbols
1152 def to_phylip(self, file_addr: Union[Path, str], drop_constant: bool = False):
1153 """Writes this Collation to a file in PHYLIP format with the given address.
1154 Note that because PHYLIP format does not support NEXUS-style ambiguities, such ambiguities will be treated as missing data.
1156 Args:
1157 file_addr: A string representing the path to an output file.
1158 drop_constant (bool, optional): An optional flag indicating whether to ignore variation units with one substantive reading.
1159 """
1160 # Populate a list of sites that will correspond to columns of the sequence alignment:
1161 substantive_variation_unit_ids = self.variation_unit_ids
1162 if drop_constant:
1163 substantive_variation_unit_ids = [
1164 vu_id
1165 for vu_id in self.variation_unit_ids
1166 if len(self.substantive_readings_by_variation_unit_id[vu_id]) > 1
1167 ]
1168 substantive_variation_unit_ids_set = set(substantive_variation_unit_ids)
1169 substantive_variation_unit_reading_tuples_set = set(self.substantive_variation_unit_reading_tuples)
1170 # Start by calculating the values we will be using here:
1171 ntax = len(self.witnesses)
1172 nchar = len(substantive_variation_unit_ids)
1173 taxlabels = []
1174 for wit in self.witnesses:
1175 taxlabel = wit.id
1176 # Then replace any disallowed characters in the string with an underscore:
1177 taxlabel = slugify(taxlabel, lowercase=False, separator='_')
1178 taxlabels.append(taxlabel)
1179 max_taxlabel_length = max(
1180 [len(taxlabel) for taxlabel in taxlabels]
1181 ) # keep track of the longest taxon label for tabular alignment purposes
1182 missing_symbol = '?'
1183 symbols = self.get_phylip_symbols()
1184 # Generate all parent folders for this file that don't already exist:
1185 Path(file_addr).parent.mkdir(parents=True, exist_ok=True)
1186 with open(file_addr, "w", encoding="ascii") as f:
1187 # Write the dimensions:
1188 f.write("%d %d\n" % (ntax, nchar))
1189 # Now write the matrix:
1190 for i, wit in enumerate(self.witnesses):
1191 taxlabel = taxlabels[i]
1192 # Add enough space after this label ensure that all sequences are nicely aligned:
1193 sequence = taxlabel + (" " * (max_taxlabel_length - len(taxlabel))) + "\t"
1194 for j, vu_id in enumerate(self.variation_unit_ids):
1195 if vu_id not in substantive_variation_unit_ids_set:
1196 continue
1197 rdg_support = self.readings_by_witness[wit.id][j]
1198 # If this reading is lacunose in this witness, then use the missing character:
1199 if sum(rdg_support) == 0:
1200 sequence += missing_symbol
1201 continue
1202 rdg_inds = [
1203 k for k, w in enumerate(rdg_support) if w > 0
1204 ] # the index list consists of the indices of all readings with any degree of certainty assigned to them
1205 # For singleton readings, just print the symbol:
1206 if len(rdg_inds) == 1:
1207 sequence += symbols[rdg_inds[0]]
1208 continue
1209 # For multiple readings, print the missing symbol:
1210 sequence += missing_symbol
1211 f.write("%s\n" % (sequence))
1212 return
1214 def get_fasta_symbols(self):
1215 """Returns a list of one-character symbols needed to represent the states of all substantive readings in FASTA format.
1217 The number of symbols equals the maximum number of substantive readings at any variation unit.
1219 Returns:
1220 A list of individual characters representing states in readings.
1221 """
1222 possible_symbols = (
1223 list(string.digits) + list(string.ascii_lowercase)[:22]
1224 ) # NOTE: for RAxML, multistate characters with an alphabet sizes up to 32 are supported
1225 # The number of symbols needed is equal to the length of the longest substantive reading vector:
1226 nsymbols = 0
1227 # If there are no witnesses, then no symbols are needed at all:
1228 if len(self.witnesses) == 0:
1229 return []
1230 wit_id = self.witnesses[0].id
1231 for i, vu_id in enumerate(self.variation_unit_ids):
1232 rdg_support = self.readings_by_witness[wit_id][i]
1233 if len(rdg_support) > len(possible_symbols):
1234 raise ValueError(
1235 f"ERROR: too many substantive readings at variation unit '{vu_id}' to represent in FASTA format."
1236 )
1237 nsymbols = max(nsymbols, len(rdg_support))
1238 fasta_symbols = possible_symbols[:nsymbols]
1239 return fasta_symbols
1241 def to_fasta(self, file_addr: Union[Path, str], drop_constant: bool = False):
1242 """Writes this Collation to a file in FASTA format with the given address.
1243 Note that because FASTA format does not support NEXUS-style ambiguities, such ambiguities will be treated as missing data.
1245 Args:
1246 file_addr: A string representing the path to an output file.
1247 drop_constant (bool, optional): An optional flag indicating whether to ignore variation units with one substantive reading.
1248 """
1249 # Populate a list of sites that will correspond to columns of the sequence alignment:
1250 substantive_variation_unit_ids = self.variation_unit_ids
1251 if drop_constant:
1252 substantive_variation_unit_ids = [
1253 vu_id
1254 for vu_id in self.variation_unit_ids
1255 if len(self.substantive_readings_by_variation_unit_id[vu_id]) > 1
1256 ]
1257 substantive_variation_unit_ids_set = set(substantive_variation_unit_ids)
1258 substantive_variation_unit_reading_tuples_set = set(self.substantive_variation_unit_reading_tuples)
1259 # Start by calculating the values we will be using here:
1260 ntax = len(self.witnesses)
1261 nchar = len(substantive_variation_unit_ids)
1262 taxlabels = []
1263 for wit in self.witnesses:
1264 taxlabel = wit.id
1265 # Then replace any disallowed characters in the string with an underscore:
1266 taxlabel = slugify(taxlabel, lowercase=False, separator='_')
1267 taxlabels.append(taxlabel)
1268 max_taxlabel_length = max(
1269 [len(taxlabel) for taxlabel in taxlabels]
1270 ) # keep track of the longest taxon label for tabular alignment purposes
1271 missing_symbol = '?'
1272 symbols = self.get_fasta_symbols()
1273 # Generate all parent folders for this file that don't already exist:
1274 Path(file_addr).parent.mkdir(parents=True, exist_ok=True)
1275 with open(file_addr, "w", encoding="ascii") as f:
1276 # Now write the matrix:
1277 with tqdm(total=len(self.witnesses)) as pbar:
1278 for i, wit in enumerate(self.witnesses):
1279 taxlabel = taxlabels[i]
1280 # Add enough space after this label ensure that all sequences are nicely aligned:
1281 sequence = ">%s\n" % taxlabel
1282 for j, vu_id in enumerate(self.variation_unit_ids):
1283 if vu_id not in substantive_variation_unit_ids_set:
1284 continue
1285 rdg_support = self.readings_by_witness[wit.id][j]
1286 # If this reading is lacunose in this witness, then use the missing character:
1287 if sum(rdg_support) == 0:
1288 sequence += missing_symbol
1289 continue
1290 rdg_inds = [
1291 k for k, w in enumerate(rdg_support) if w > 0
1292 ] # the index list consists of the indices of all readings with any degree of certainty assigned to them
1293 # For singleton readings, just print the symbol:
1294 if len(rdg_inds) == 1:
1295 sequence += symbols[rdg_inds[0]]
1296 continue
1297 # For multiple readings, print the missing symbol:
1298 sequence += missing_symbol
1299 f.write("%s\n" % (sequence))
1300 pbar.update(1)
1301 return
1303 def get_beast_symbols(self):
1304 """Returns a list of one-character symbols needed to represent the states of all substantive readings in BEAST format.
1306 The number of symbols equals the maximum number of substantive readings at any variation unit.
1308 Returns:
1309 A list of individual characters representing states in readings.
1310 """
1311 possible_symbols = (
1312 list(string.digits) + list(string.ascii_lowercase) + list(string.ascii_uppercase)
1313 ) # NOTE: for BEAST, any number of states should theoretically be permissible, but since code maps are required for some reason, we will limit the number of symbols to 62 for now
1314 # The number of symbols needed is equal to the length of the longest substantive reading vector:
1315 nsymbols = 0
1316 # If there are no witnesses, then no symbols are needed at all:
1317 if len(self.witnesses) == 0:
1318 return []
1319 wit_id = self.witnesses[0].id
1320 for i, vu_id in enumerate(self.variation_unit_ids):
1321 rdg_support = self.readings_by_witness[wit_id][i]
1322 if len(rdg_support) > len(possible_symbols):
1323 raise ValueError(
1324 f"ERROR: too many substantive readings at variation unit '{vu_id}' to represent in BEAST format."
1325 )
1326 nsymbols = max(nsymbols, len(rdg_support))
1327 beast_symbols = possible_symbols[:nsymbols]
1328 return beast_symbols
1330 def get_tip_date_range(self):
1331 """Gets the minimum and maximum dates attested among the witnesses.
1332 Also checks if the witness with the latest possible date has a fixed date
1333 (i.e, if the lower and upper bounds for its date are the same)
1334 and issues a warning if not, as this will cause unusual behavior in BEAST 2.
1336 Returns:
1337 A tuple containing the earliest and latest possible tip dates.
1338 """
1339 earliest_date = None
1340 earliest_wit = None
1341 latest_date = None
1342 latest_wit = None
1343 for wit in self.witnesses:
1344 wit_id = wit.id
1345 date_range = wit.date_range
1346 if date_range[0] is not None:
1347 if earliest_date is not None:
1348 earliest_wit = wit if date_range[0] < earliest_date else earliest_wit
1349 earliest_date = min(date_range[0], earliest_date)
1350 else:
1351 earliest_wit = wit
1352 earliest_date = date_range[0]
1353 if date_range[1] is not None:
1354 if latest_date is not None:
1355 latest_wit = (
1356 wit
1357 if (date_range[1] > latest_date or (date_range[0] == date_range[1] == latest_date))
1358 else latest_wit
1359 ) # the second check ensures that a witness with a fixed date is preferred to a witness with a date range that ends with the same date
1360 latest_date = max(date_range[1], latest_date)
1361 else:
1362 latest_wit = wit
1363 latest_date = date_range[1]
1364 if latest_wit.date_range[0] is None or latest_wit.date_range[0] != latest_wit.date_range[1]:
1365 print(
1366 "WARNING: the latest witness, %s, has a variable date range; this will result in problems with time-dependent substitution models and misalignment of trees in BEAST DensiTree outputs! Please ensure that witness %s has a fixed date."
1367 % (latest_wit.id, latest_wit.id)
1368 )
1369 return (earliest_date, latest_date)
1371 def get_beast_origin_span(self, tip_date_range):
1372 """Returns a tuple containing the lower and upper bounds for the height of the origin of the Birth-Death Skyline model.
1373 The upper bound on the height of the tree is the difference between the latest tip date
1374 and the lower bound on the date of the original work, if both are defined;
1375 otherwise, it is left undefined.
1376 The lower bound on the height of the tree is the difference between the latest tip date
1377 and the upper bound on the date of the original work, if both are defined;
1378 otherwise, it is the difference between the earliest tip date and the latest, if both are defined.
1380 Args:
1381 tip_date_range: A tuple containing the earliest and latest possible tip dates.
1383 Returns:
1384 A tuple containing lower and upper bounds on the origin height for the Birth-Death Skyline model.
1385 """
1386 origin_span = [0, None]
1387 # If the upper bound on the date of the work's composition is defined, then set the lower bound on the height of the origin using it and the latest tip date
1388 # (note that if it had to be defined in terms of witness date lower bounds, then this would have happened already):
1389 if self.origin_date_range[1] is not None:
1390 origin_span[0] = tip_date_range[1] - self.origin_date_range[1]
1391 # If the lower bound on the date of the work's composition is defined, then set the upper bound on the height of the origin using it and the latest tip date:
1392 if self.origin_date_range[0] is not None:
1393 origin_span[1] = tip_date_range[1] - self.origin_date_range[0]
1394 return tuple(origin_span)
1396 def get_beast_date_map(self, taxlabels):
1397 """Returns a string representing witness-to-date mappings in BEAST format.
1399 Since this format requires single dates as opposed to date ranges,
1400 witnesses with closed date ranges will be mapped to the average of their lower and upper bounds,
1401 and witnesses with open date ranges will not be mapped.
1403 Args:
1404 taxlabels: A list of slugified taxon labels.
1406 Returns:
1407 A string containing comma-separated date calibrations of the form witness_id=date.
1408 """
1409 calibrate_strings = []
1410 for i, wit in enumerate(self.witnesses):
1411 taxlabel = taxlabels[i]
1412 date_range = wit.date_range
1413 # If either end of this witness's date range is empty, then do not include it:
1414 if date_range[0] is None or date_range[1] is None:
1415 continue
1416 # Otherwise, take the midpoint of its date range as its date:
1417 date = int((date_range[0] + date_range[1]) / 2)
1418 calibrate_string = "%s=%d" % (taxlabel, date)
1419 calibrate_strings.append(calibrate_string)
1420 # Then output the full date map string:
1421 date_map = ",".join(calibrate_strings)
1422 return date_map
1424 def get_beast_code_map_for_unit(self, symbols, missing_symbol, vu_ind):
1425 """Returns a string containing state/reading code mappings in BEAST format using the given single-state and missing state symbols for the character/variation unit at the given index.
1426 If the variation unit at the given index is a singleton unit (i.e., if it has only one substantive reading), then a code for a dummy state will be included.
1428 Args:
1429 vu_ind: An integer index for the desired unit.
1431 Returns:
1432 A string containing comma-separated code mappings.
1433 """
1434 vu = self.variation_units[vu_ind]
1435 vu_id = vu.id
1436 code_map = {}
1437 for k in range(len(self.substantive_readings_by_variation_unit_id[vu.id])):
1438 code_map[symbols[k]] = str(k)
1439 # If this site is a singleton site, then add a code mapping for the dummy state:
1440 if len(self.substantive_readings_by_variation_unit_id[vu.id]) == 1:
1441 code_map[symbols[1]] = str(1)
1442 # Then add a mapping for the missing state, including a dummy state if this is a singleton site:
1443 code_map[missing_symbol] = " ".join(
1444 str(k) for k in range(len(self.substantive_readings_by_variation_unit_id[vu.id]))
1445 )
1446 # If this site is a singleton site, then add the dummy state to the missing state mapping:
1447 if len(self.substantive_readings_by_variation_unit_id[vu.id]) == 1:
1448 code_map[missing_symbol] = code_map[missing_symbol] + " " + str(1)
1449 # Then combine all of the mappings into a single string:
1450 code_map_string = ", ".join([code + "=" + code_map[code] for code in code_map])
1451 return code_map_string
1453 def get_beast_equilibrium_frequencies_for_unit(self, vu_ind):
1454 """Returns a string containing state/reading equilibrium frequencies in BEAST format for the character/variation unit at the given index.
1455 Since the equilibrium frequencies are not used with the substitution models, the equilibrium frequencies simply correspond to a uniform distribution over the states.
1456 If the variation unit at the given index is a singleton unit (i.e., if it has only one substantive reading), then an equilibrium frequency of 0 will be added for a dummy state.
1458 Args:
1459 vu_ind: An integer index for the desired unit.
1461 Returns:
1462 A string containing space-separated equilibrium frequencies.
1463 """
1464 vu = self.variation_units[vu_ind]
1465 vu_id = vu.id
1466 # If this unit is a singleton, then return the string "0.5 0.5":
1467 if len(self.substantive_readings_by_variation_unit_id[vu_id]) == 1:
1468 return "0.5 0.5"
1469 # Otherwise, set the equilibrium frequencies according to a uniform distribution:
1470 equilibrium_frequencies = [1.0 / len(self.substantive_readings_by_variation_unit_id[vu_id])] * len(
1471 self.substantive_readings_by_variation_unit_id[vu_id]
1472 )
1473 equilibrium_frequencies_string = " ".join([str(w) for w in equilibrium_frequencies])
1474 return equilibrium_frequencies_string
1476 def get_beast_root_frequencies_for_unit(self, vu_ind):
1477 """Returns a string containing state/reading root frequencies in BEAST format for the character/variation unit at the given index.
1478 The root frequencies are calculated from the intrinsic odds at this unit.
1479 If the variation unit at the given index is a singleton unit (i.e., if it has only one substantive reading), then a root frequency of 0 will be added for a dummy state.
1480 If no intrinsic odds are specified, then a uniform distribution over all states is assumed.
1482 Args:
1483 vu_ind: An integer index for the desired unit.
1485 Returns:
1486 A string containing space-separated root frequencies.
1487 """
1488 vu = self.variation_units[vu_ind]
1489 vu_id = vu.id
1490 intrinsic_relations = vu.intrinsic_relations
1491 intrinsic_odds_by_id = self.intrinsic_odds_by_id
1492 # If this unit is a singleton, then return the string "1 0":
1493 if len(self.substantive_readings_by_variation_unit_id[vu_id]) == 1:
1494 return "1 0"
1495 # If this unit has no intrinsic odds, then assume a uniform distribution over all readings:
1496 if len(intrinsic_relations) == 0:
1497 root_frequencies = [1.0 / len(self.substantive_readings_by_variation_unit_id[vu_id])] * len(
1498 self.substantive_readings_by_variation_unit_id[vu_id]
1499 )
1500 root_frequencies_string = " ".join([str(w) for w in root_frequencies])
1501 return root_frequencies_string
1502 # We will populate the root frequencies based on the intrinsic odds of the readings:
1503 root_frequencies_by_id = {}
1504 for rdg_id in self.substantive_readings_by_variation_unit_id[vu_id]:
1505 root_frequencies_by_id[rdg_id] = 0
1506 # First, construct an adjacency list for efficient edge iteration:
1507 neighbors_by_source = {}
1508 for edge in intrinsic_relations:
1509 s = edge[0]
1510 t = edge[1]
1511 if s not in neighbors_by_source:
1512 neighbors_by_source[s] = []
1513 if t not in neighbors_by_source:
1514 neighbors_by_source[t] = []
1515 neighbors_by_source[s].append(t)
1516 # Next, identify all readings that are not targeted by any intrinsic odds relation:
1517 in_degree_by_reading = {}
1518 for edge in intrinsic_relations:
1519 s = edge[0]
1520 t = edge[1]
1521 if s not in in_degree_by_reading:
1522 in_degree_by_reading[s] = 0
1523 if t not in in_degree_by_reading:
1524 in_degree_by_reading[t] = 0
1525 in_degree_by_reading[t] += 1
1526 starting_nodes = [t for t in in_degree_by_reading if in_degree_by_reading[t] == 0]
1527 # Set the root frequencies for these readings to 1 (they will be normalized later):
1528 for starting_node in starting_nodes:
1529 root_frequencies_by_id[starting_node] = 1.0
1530 # Next, set the frequencies for the remaining readings recursively using the adjacency list:
1531 def update_root_frequencies(s):
1532 for t in neighbors_by_source[s]:
1533 intrinsic_category = intrinsic_relations[(s, t)]
1534 odds = (
1535 intrinsic_odds_by_id[intrinsic_category]
1536 if intrinsic_odds_by_id[intrinsic_category] is not None
1537 else 1.0
1538 ) # TODO: This needs to be handled using parameters once we have it implemented in BEAST
1539 root_frequencies_by_id[t] = root_frequencies_by_id[s] / odds
1540 update_root_frequencies(t)
1541 return
1543 for starting_node in starting_nodes:
1544 update_root_frequencies(starting_node)
1545 # Then produce a normalized vector of root frequencies that corresponds to a probability distribution:
1546 root_frequencies = [
1547 root_frequencies_by_id[rdg_id] for rdg_id in self.substantive_readings_by_variation_unit_id[vu_id]
1548 ]
1549 total_frequencies = sum(root_frequencies)
1550 for k in range(len(root_frequencies)):
1551 root_frequencies[k] = root_frequencies[k] / total_frequencies
1552 root_frequencies_string = " ".join([str(w) for w in root_frequencies])
1553 return root_frequencies_string
1555 def to_beast(
1556 self,
1557 file_addr: Union[Path, str],
1558 drop_constant: bool = False,
1559 clock_model: ClockModel = ClockModel.strict,
1560 ancestral_logger: AncestralLogger = AncestralLogger.state,
1561 seed: int = None,
1562 ):
1563 """Writes this Collation to a file in BEAST format with the given address.
1565 Args:
1566 file_addr: A string representing the path to an output file.
1567 drop_constant (bool, optional): An optional flag indicating whether to ignore variation units with one substantive reading.
1568 clock_model: A ClockModel option indicating which clock model to use.
1569 ancestral_logger: An AncestralLogger option indicating which class of logger (if any) to use for ancestral states.
1570 seed: A seed for random number generation (for setting initial values of unspecified transcriptional rates).
1571 """
1572 # Populate a list of sites that will correspond to columns of the sequence alignment:
1573 substantive_variation_unit_ids = self.variation_unit_ids
1574 if drop_constant:
1575 substantive_variation_unit_ids = [
1576 vu_id
1577 for vu_id in self.variation_unit_ids
1578 if len(self.substantive_readings_by_variation_unit_id[vu_id]) > 1
1579 ]
1580 # Populate sets of substantive variation unit IDs and substantive variant reading tuples:
1581 substantive_variation_unit_ids_set = set(substantive_variation_unit_ids)
1582 substantive_variation_unit_reading_tuples_set = set(self.substantive_variation_unit_reading_tuples)
1583 # First, calculate the values we will be using for the main template:
1584 taxlabels = [slugify(wit.id, lowercase=False, separator='_') for wit in self.witnesses]
1585 missing_symbol = '?'
1586 symbols = self.get_beast_symbols()
1587 tip_date_range = self.get_tip_date_range()
1588 origin_span = self.get_beast_origin_span(tip_date_range)
1589 date_map = self.get_beast_date_map(taxlabels)
1590 # Then populate the necessary objects for the BEAST XML Jinja template:
1591 witness_objects = []
1592 variation_unit_objects = []
1593 intrinsic_category_objects = []
1594 transcriptional_category_objects = []
1595 with tqdm(total=len(self.witnesses)) as pbar:
1596 # Start with witnesses:
1597 for i, wit in enumerate(self.witnesses):
1598 witness_object = {}
1599 # Copy the ID for this witness:
1600 witness_object["id"] = wit.id
1601 # Copy its date bounds:
1602 witness_object["min_date"] = wit.date_range[0]
1603 witness_object["max_date"] = wit.date_range[1]
1604 # Populate its sequence from its entries in the witness's readings dictionary:
1605 sequence = ""
1606 for j, rdg_support in enumerate(self.readings_by_witness[wit.id]):
1607 vu_id = self.variation_unit_ids[j]
1608 # Skip any variation units deemed non-substantive:
1609 if vu_id not in substantive_variation_unit_ids:
1610 continue
1611 # If this witness has a certainty of 0 for all readings, then it is a gap; assign a likelihood of 1 to each reading:
1612 if sum(rdg_support) == 0:
1613 for k, w in enumerate(rdg_support):
1614 sequence += "1"
1615 if k < len(rdg_support) - 1:
1616 sequence += ", "
1617 else:
1618 if len(rdg_support) > 1:
1619 sequence += "; "
1620 else:
1621 # If this site is a singleton site, then add a dummy state:
1622 sequence += ", 0; "
1623 # Otherwise, read the probabilities as they are given:
1624 else:
1625 for k, w in enumerate(rdg_support):
1626 sequence += str(w)
1627 if k < len(rdg_support) - 1:
1628 sequence += ", "
1629 else:
1630 if len(rdg_support) > 1:
1631 sequence += "; "
1632 else:
1633 # If this site is a singleton site, then add a dummy state:
1634 sequence += ", 0; "
1635 # Strip the final semicolon and space from the sequence:
1636 sequence = sequence.strip("; ")
1637 # Then set the witness object's sequence attribute to this string:
1638 witness_object["sequence"] = sequence
1639 witness_objects.append(witness_object)
1640 pbar.update(1)
1641 # Then proceed to variation units:
1642 for j, vu in enumerate(self.variation_units):
1643 if vu.id not in substantive_variation_unit_ids_set:
1644 continue
1645 variation_unit_object = {}
1646 # Copy the one-based index of this variation unit:
1647 variation_unit_object["index"] = j + 1
1648 # Copy the ID of this variation unit:
1649 variation_unit_object["id"] = vu.id
1650 # Set a flag indicating if this variation unit is constant:
1651 variation_unit_object["is_constant"] = (
1652 True if len(self.substantive_readings_by_variation_unit_id[vu.id]) == 1 else False
1653 )
1654 # Copy this variation unit's number of substantive readings,
1655 # setting it to 2 if it is a singleton unit:
1656 variation_unit_object["nstates"] = (
1657 len(self.substantive_readings_by_variation_unit_id[vu.id])
1658 if len(self.substantive_readings_by_variation_unit_id[vu.id]) > 1
1659 else 2
1660 )
1661 # Then construct the code map for this unit:
1662 variation_unit_object["code_map"] = self.get_beast_code_map_for_unit(symbols, missing_symbol, j)
1663 # Then populate a comma-separated string of reading labels for this unit:
1664 rdg_texts = []
1665 vu_label = vu.id
1666 for rdg in vu.readings:
1667 key = tuple([vu.id, rdg.id])
1668 if key not in substantive_variation_unit_reading_tuples_set:
1669 continue
1670 rdg_text = slugify(rdg.text, lowercase=False, allow_unicode=True, separator='_')
1671 # Replace any empty reading text with an omission marker:
1672 if rdg_text == "":
1673 rdg_text = "om."
1674 rdg_texts.append(rdg_text)
1675 # If this site is a singleton site, then add a dummy reading for the dummy state:
1676 if len(self.substantive_readings_by_variation_unit_id[vu.id]) == 1:
1677 rdg_texts.append("DUMMY")
1678 rdg_texts_string = ", ".join(rdg_texts)
1679 variation_unit_object["rdg_texts"] = rdg_texts_string
1680 # Then populate this unit's equilibrium frequency string and its root frequency string:
1681 variation_unit_object["equilibrium_frequencies"] = self.get_beast_equilibrium_frequencies_for_unit(j)
1682 variation_unit_object["root_frequencies"] = self.get_beast_root_frequencies_for_unit(j)
1683 # Then populate a dictionary mapping epoch height ranges to lists of off-diagonal entries for substitution models:
1684 rate_objects_by_epoch_height_range = {}
1685 epoch_height_ranges = []
1686 # Then proceed based on whether the transcriptional relations for this variation unit have been defined:
1687 if len(vu.transcriptional_relations_by_date_range) == 0:
1688 # If there are no transcriptional relations, then map the epoch range of (None, None) to their list of off-diagonal entries:
1689 epoch_height_ranges.append((None, None))
1690 rate_objects_by_epoch_height_range[(None, None)] = []
1691 rate_objects = rate_objects_by_epoch_height_range[(None, None)]
1692 if len(self.substantive_readings_by_variation_unit_id[vu.id]) == 1:
1693 # If this is a singleton site, then use an arbitrary 2x2 rate matrix:
1694 rate_objects.append({"transcriptional_categories": ["default"], "expression": None})
1695 rate_objects.append({"transcriptional_categories": ["default"], "expression": None})
1696 else:
1697 # If this is a site with multiple substantive readings, but no transcriptional relations list,
1698 # then use a Lewis Mk substitution matrix with the appropriate number of states:
1699 for k_1, rdg_id_1 in enumerate(self.substantive_readings_by_variation_unit_id[vu.id]):
1700 for k_2, rdg_id_2 in enumerate(self.substantive_readings_by_variation_unit_id[vu.id]):
1701 # Skip diagonal elements:
1702 if k_1 == k_2:
1703 continue
1704 rate_objects.append({"transcriptional_categories": ["default"], "expression": None})
1705 else:
1706 # Otherwise, proceed for every date range:
1707 for date_range in vu.transcriptional_relations_by_date_range:
1708 # Get the map of transcriptional relations for reference later:
1709 transcriptional_relations = vu.transcriptional_relations_by_date_range[date_range]
1710 # Now get the epoch height range corresponding to this date range, and initialize its list in the dictionary:
1711 epoch_height_range = [None, None]
1712 epoch_height_range[0] = tip_date_range[1] - date_range[1] if date_range[1] is not None else None
1713 epoch_height_range[1] = tip_date_range[1] - date_range[0] if date_range[0] is not None else None
1714 epoch_height_range = tuple(epoch_height_range)
1715 epoch_height_ranges.append(epoch_height_range)
1716 rate_objects_by_epoch_height_range[epoch_height_range] = []
1717 rate_objects = rate_objects_by_epoch_height_range[epoch_height_range]
1718 # Then proceed for every pair of readings in this unit:
1719 for k_1, rdg_id_1 in enumerate(self.substantive_readings_by_variation_unit_id[vu.id]):
1720 for k_2, rdg_id_2 in enumerate(self.substantive_readings_by_variation_unit_id[vu.id]):
1721 # Skip diagonal elements:
1722 if k_1 == k_2:
1723 continue
1724 # If the first reading has no transcriptional relation to the second in this unit, then use the default rate:
1725 if (rdg_id_1, rdg_id_2) not in transcriptional_relations:
1726 rate_objects.append({"transcriptional_categories": ["default"], "expression": None})
1727 continue
1728 # Otherwise, if only one category of transcriptional relations holds between the first and second readings,
1729 # then use its rate:
1730 if len(transcriptional_relations[(rdg_id_1, rdg_id_2)]) == 1:
1731 # If there is only one such category, then add its rate as a standalone var element:
1732 transcriptional_category = list(transcriptional_relations[(rdg_id_1, rdg_id_2)])[0]
1733 rate_objects.append(
1734 {"transcriptional_categories": [transcriptional_category], "expression": None}
1735 )
1736 continue
1737 # If there is more than one, then add a var element that is a sum of the individual categories' rates:
1738 transcriptional_categories = list(transcriptional_relations[(rdg_id_1, rdg_id_2)])
1739 args = []
1740 for transcriptional_category in transcriptional_categories:
1741 args.append("%s_rate" % transcriptional_category)
1742 args_string = " ".join(args)
1743 ops = ["+"] * (len(args) - 1)
1744 ops_string = " ".join(ops)
1745 expression_string = " ".join([args_string, ops_string])
1746 rate_objects.append(
1747 {
1748 "transcriptional_categories": transcriptional_categories,
1749 "expression": expression_string,
1750 }
1751 )
1752 # Now reorder the list of epoch height ranges, and get a list of non-null epoch dates in ascending order from the dictionary:
1753 epoch_height_ranges.reverse()
1754 epoch_heights = [
1755 epoch_height_range[0] for epoch_height_range in epoch_height_ranges if epoch_height_range[0] is not None
1756 ]
1757 # Then add all of these data structures to the variation unit object:
1758 variation_unit_object["epoch_heights"] = epoch_heights
1759 variation_unit_object["epoch_heights_string"] = " ".join(
1760 [str(epoch_height) for epoch_height in epoch_heights]
1761 )
1762 variation_unit_object["epoch_height_ranges"] = epoch_height_ranges
1763 variation_unit_object["epoch_rates"] = [
1764 rate_objects_by_epoch_height_range[epoch_height_range] for epoch_height_range in epoch_height_ranges
1765 ]
1766 variation_unit_objects.append(variation_unit_object)
1767 # Then proceed to intrinsic odds categories:
1768 for intrinsic_category in self.intrinsic_categories:
1769 intrinsic_category_object = {}
1770 # Copy the ID of this intrinsic category:
1771 intrinsic_category_object["id"] = intrinsic_category
1772 # Then copy the odds factors associated with this intrinsic category,
1773 # setting it to 1.0 if it is not specified and setting the estimate attribute accordingly:
1774 odds = self.intrinsic_odds_by_id[intrinsic_category]
1775 intrinsic_category_object["odds"] = odds if odds is not None else 1.0
1776 intrinsic_category_object["estimate"] = "false" if odds is not None else "true"
1777 intrinsic_category_objects.append(intrinsic_category_object)
1778 # Then proceed to transcriptional rate categories:
1779 rng = np.random.default_rng(seed)
1780 for transcriptional_category in self.transcriptional_categories:
1781 transcriptional_category_object = {}
1782 # Copy the ID of this transcriptional category:
1783 transcriptional_category_object["id"] = transcriptional_category
1784 # Then copy the rate of this transcriptional category,
1785 # setting it to a random number sampled from a Gamma distribution if it is not specified and setting the estimate attribute accordingly:
1786 rate = self.transcriptional_rates_by_id[transcriptional_category]
1787 transcriptional_category_object["rate"] = rate if rate is not None else rng.gamma(5.0, 2.0)
1788 transcriptional_category_object["estimate"] = "false" if rate is not None else "true"
1789 transcriptional_category_objects.append(transcriptional_category_object)
1790 # Now render the output XML file using the Jinja template:
1791 env = Environment(loader=PackageLoader("teiphy", "templates"), autoescape=select_autoescape())
1792 template = env.get_template("beast_template.xml")
1793 rendered = template.render(
1794 nsymbols=len(symbols),
1795 date_map=date_map,
1796 origin_span=origin_span,
1797 clock_model=clock_model.value,
1798 clock_rate_categories=2 * len(self.witnesses) - 1,
1799 ancestral_logger=ancestral_logger.value,
1800 witnesses=witness_objects,
1801 variation_units=variation_unit_objects,
1802 non_constant_variation_units=[
1803 variation_unit_object
1804 for variation_unit_object in variation_unit_objects
1805 if not variation_unit_object["is_constant"]
1806 ],
1807 constant_variation_unit_filter=",".join(
1808 [
1809 str(variation_unit_object["index"])
1810 for variation_unit_object in variation_unit_objects
1811 if variation_unit_object["is_constant"]
1812 ]
1813 ),
1814 intrinsic_categories=intrinsic_category_objects,
1815 transcriptional_categories=transcriptional_category_objects,
1816 )
1817 # Generate all parent folders for this file that don't already exist:
1818 Path(file_addr).parent.mkdir(parents=True, exist_ok=True)
1819 with open(file_addr, "w", encoding="utf-8") as f:
1820 f.write(rendered)
1821 return
1823 def to_numpy(self, drop_constant: bool = False, split_missing: SplitMissingType = None):
1824 """Returns this Collation in the form of a NumPy array, along with arrays of its row and column labels.
1826 Args:
1827 drop_constant (bool, optional): An optional flag indicating whether to ignore variation units with one substantive reading.
1828 split_missing (SplitMissingType, optional): An option indicating whether or not to treat missing characters/variation units as having a contribution of 1 split over all states/readings.
1829 If not specified, then missing data is ignored (i.e., all states are 0).
1830 If "uniform", then the contribution of 1 is divided evenly over all substantive readings.
1831 If "proportional", then the contribution of 1 is divided between the readings in proportion to their support among the witnesses that are not missing.
1833 Returns:
1834 A NumPy array with a row for each substantive reading and a column for each witness.
1835 A list of substantive reading ID strings.
1836 A list of witness ID strings.
1837 """
1838 # Populate a list of sites that will correspond to columns of the sequence alignment:
1839 substantive_variation_unit_ids = self.variation_unit_ids
1840 if drop_constant:
1841 substantive_variation_unit_ids = [
1842 vu_id
1843 for vu_id in self.variation_unit_ids
1844 if len(self.substantive_readings_by_variation_unit_id[vu_id]) > 1
1845 ]
1846 substantive_variation_unit_ids_set = set(substantive_variation_unit_ids)
1847 substantive_variation_unit_reading_tuples_set = set(self.substantive_variation_unit_reading_tuples)
1848 # Initialize the output array with the appropriate dimensions:
1849 reading_labels = []
1850 for vu in self.variation_units:
1851 if vu.id not in substantive_variation_unit_ids_set:
1852 continue
1853 for rdg in vu.readings:
1854 key = tuple([vu.id, rdg.id])
1855 if key in substantive_variation_unit_reading_tuples_set:
1856 reading_labels.append(vu.id + ", " + rdg.text)
1857 witness_labels = [wit.id for wit in self.witnesses]
1858 matrix = np.zeros((len(reading_labels), len(witness_labels)), dtype=float)
1859 # For each variation unit, keep a record of the proportion of non-missing witnesses supporting the substantive variant readings:
1860 support_proportions_by_unit = {}
1861 for j, vu_id in enumerate(self.variation_unit_ids):
1862 if vu_id not in substantive_variation_unit_ids_set:
1863 continue
1864 support_proportions = [0.0] * len(self.substantive_readings_by_variation_unit_id[vu_id])
1865 for i, wit in enumerate(self.witnesses):
1866 rdg_support = self.readings_by_witness[wit.id][j]
1867 for l, w in enumerate(rdg_support):
1868 support_proportions[l] += w
1869 norm = (
1870 sum(support_proportions) if sum(support_proportions) > 0 else 1.0
1871 ) # if this variation unit has no extant witnesses (e.g., if its only witnesses are fragmentary and we have excluded them), then assume a norm of 1 to avoid division by zero
1872 for l in range(len(support_proportions)):
1873 support_proportions[l] = support_proportions[l] / norm
1874 support_proportions_by_unit[vu_id] = support_proportions
1875 # Then populate it with the appropriate values:
1876 col_ind = 0
1877 with tqdm(total=len(self.witnesses)) as pbar:
1878 for i, wit in enumerate(self.witnesses):
1879 row_ind = 0
1880 for j, vu_id in enumerate(self.variation_unit_ids):
1881 if vu_id not in substantive_variation_unit_ids_set:
1882 continue
1883 rdg_support = self.readings_by_witness[wit.id][j]
1884 # If this reading support vector sums to 0, then this is missing data; handle it as specified:
1885 if sum(rdg_support) == 0:
1886 if split_missing == SplitMissingType.uniform:
1887 for l in range(len(rdg_support)):
1888 matrix[row_ind, col_ind] = 1 / len(rdg_support)
1889 row_ind += 1
1890 elif split_missing == SplitMissingType.proportional:
1891 for l in range(len(rdg_support)):
1892 matrix[row_ind, col_ind] = support_proportions_by_unit[vu_id][l]
1893 row_ind += 1
1894 else:
1895 row_ind += len(rdg_support)
1896 # Otherwise, add its coefficients normally:
1897 else:
1898 for l in range(len(rdg_support)):
1899 matrix[row_ind, col_ind] = rdg_support[l]
1900 row_ind += 1
1901 col_ind += 1
1902 pbar.update(1)
1903 return matrix, reading_labels, witness_labels
1905 def get_ext_matrix(self, drop_constant: bool = False, split_missing: SplitMissingType = None):
1906 """Returns a NumPy matrix containing a row and column for each witness and the number of variation units shared by the row and column witnesses in each cell.
1907 Note that if the split_missing option is specified, all variation units are counted.
1909 Args:
1910 drop_constant (bool, optional): An optional flag indicating whether to ignore variation units with one substantive reading.
1911 Default value is False.
1912 split_missing (SplitMissingType, optional): An option indicating whether or not to treat missing characters/variation units as having a contribution of 1 split over all states/readings.
1913 If not specified, then missing data is ignored (i.e., all states are 0).
1914 If "uniform", then the contribution of 1 is divided evenly over all substantive readings.
1915 If "proportional", then the contribution of 1 is divided between the readings in proportion to their support among the witnesses that are not missing.
1917 Returns:
1918 A NumPy matrix with a row and column for each witness and the number of variation units shared by the row and column witnesses in each cell.
1919 """
1920 # Populate a list of sites that will correspond to columns of the sequence alignment:
1921 substantive_variation_unit_ids = self.variation_unit_ids
1922 if drop_constant:
1923 substantive_variation_unit_ids = [
1924 vu_id
1925 for vu_id in self.variation_unit_ids
1926 if len(self.substantive_readings_by_variation_unit_id[vu_id]) > 1
1927 ]
1928 substantive_variation_unit_ids_set = set(substantive_variation_unit_ids)
1929 # Then initialize the output matrix:
1930 witness_labels = [wit.id for wit in self.witnesses]
1931 ext_matrix = ext_matrix = np.full((len(witness_labels), len(witness_labels)), 0, dtype=int)
1932 # If the split_missing option has been specified, then all entries in the matrix will be the number of substantive variation units,
1933 # so we can just fill the matrix with this value and return it:
1934 if split_missing is not None:
1935 ext_matrix = ext_matrix = np.full(
1936 (len(witness_labels), len(witness_labels)), len(substantive_variation_unit_ids), dtype=int
1937 )
1938 return ext_matrix
1939 # Otherwise, populate the matrix for all pairs of witnesses:
1940 with tqdm(total=len(self.witnesses) ** 2) as pbar:
1941 # Then calculate the mutual information contribution for each pair of witnesses:
1942 for i, wit_1 in enumerate(witness_labels):
1943 for j, wit_2 in enumerate(witness_labels):
1944 shared_ext_units = 0
1945 # The contribution to the entry for these witnesses will be identical regardless of the order in which they are specified,
1946 # so we only have to calculate it once:
1947 if i > j:
1948 pbar.update(1)
1949 continue
1950 # Otherwise, calculate the number of substantive variation units at which both of these witnesses are extant:
1951 for k, vu_id in enumerate(self.variation_unit_ids):
1952 if vu_id not in substantive_variation_unit_ids_set:
1953 continue
1954 wit_1_rdg_support = self.readings_by_witness[wit_1][k]
1955 wit_2_rdg_support = self.readings_by_witness[wit_2][k]
1956 if sum(wit_1_rdg_support) == 0.0 or sum(wit_2_rdg_support) == 0.0:
1957 continue
1958 shared_ext_units += 1
1959 ext_matrix[i][j] = shared_ext_units
1960 ext_matrix[j][i] = shared_ext_units
1961 pbar.update(1)
1962 return ext_matrix
1964 def transform_matrix(self, matrix: np.ndarray, transform_matrix: TransformMatrixType = None):
1965 """Transforms a given matrix's columns based on the specified transform_matrix option.
1967 Args:
1968 matrix (ndarray): The matrix whose columns are to be transformed.
1969 transform_matrix (TransformMatrixType, optional): A TransformMatrixType option indicating how the columns of a witness-to-witness matrix output should be transformed.
1970 Only applicable for tabular outputs in which the rows and columns correspond to the witnesses in the collation.
1972 Returns:
1973 A Pandas DataFrame corresponding to a collation matrix with reading frequencies or a long table with discrete reading states.
1974 """
1975 # If no transform_matrix option was supplied, then return the matrix as-is:
1976 if transform_matrix is None:
1977 return matrix
1978 # Otherwise, apply the specified column transformation:
1979 if transform_matrix == TransformMatrixType.stddev:
1980 means = np.mean(matrix, 0) # get the means of the columns (axis 0)
1981 stddevs = np.std(matrix, 0) # get the standard deviations of the columns (axis 0)
1982 transformed_matrix = np.full(
1983 matrix.shape, 0.0, dtype=float
1984 ) # the standard deviation can only be 0 if all values equal the mean, so a default value of 0 is acceptable
1985 np.divide(matrix - means, stddevs, out=transformed_matrix, where=(stddevs != 0))
1986 return transformed_matrix
1987 if transform_matrix == TransformMatrixType.mad:
1988 medians = np.median(matrix, 0) # get the medians of the columns (axis 0)
1989 mads = np.median(np.abs(matrix - medians), 0) # get the median absolute deviations of the columns (axis 0)
1990 transformed_matrix = np.full(
1991 matrix.shape, np.nan, dtype=float
1992 ) # the MAD can be 0 even if all values are not equal to the median, so NaN should be used if a column's MAD is 0
1993 np.divide(matrix - medians, mads, out=transformed_matrix, where=(mads != 0))
1994 return transformed_matrix
1996 def to_distance_matrix(
1997 self,
1998 drop_constant: bool = False,
1999 proportion: bool = False,
2000 show_ext: bool = False,
2001 transform_matrix: TransformMatrixType = None,
2002 ):
2003 """Transforms this Collation into a NumPy distance matrix between witnesses, along with an array of its labels for the witnesses.
2004 Distances can be computed either as counts of disagreements (the default setting), or as proportions of disagreements over all variation units where both witnesses have singleton readings.
2005 Optionally, the count of units where both witnesses have singleton readings can be included after the count/proportion of disagreements.
2007 Args:
2008 drop_constant (bool, optional): An optional flag indicating whether to ignore variation units with one substantive reading.
2009 Default value is False.
2010 proportion (bool, optional): An optional flag indicating whether or not to calculate distances as proportions over extant, unambiguous variation units.
2011 Default value is False.
2012 show_ext: An optional flag indicating whether each cell in the matrix
2013 should include the number of their extant, unambiguous variation units after the number of their disagreements.
2014 Default value is False.
2015 transform_matrix (TransformMatrixType, optional): A TransformMatrixType option indicating how the columns of the matrix should be transformed.
2017 Returns:
2018 A NumPy distance matrix with a row and column for each witness.
2019 A list of witness ID strings.
2020 """
2021 # Populate a list of sites that will correspond to columns of the sequence alignment:
2022 substantive_variation_unit_ids = self.variation_unit_ids
2023 if drop_constant:
2024 substantive_variation_unit_ids = [
2025 vu_id
2026 for vu_id in self.variation_unit_ids
2027 if len(self.substantive_readings_by_variation_unit_id[vu_id]) > 1
2028 ]
2029 substantive_variation_unit_ids_set = set(substantive_variation_unit_ids)
2030 substantive_variation_unit_reading_tuples_set = set(self.substantive_variation_unit_reading_tuples)
2031 # Initialize the output array with the appropriate dimensions:
2032 witness_labels = [wit.id for wit in self.witnesses]
2033 matrix = np.full((len(witness_labels), len(witness_labels)), 0, dtype=int) # ints of the form disagreements
2034 with tqdm(total=len(self.witnesses) ** 2) as pbar:
2035 for i, wit_1 in enumerate(witness_labels):
2036 for j, wit_2 in enumerate(witness_labels):
2037 disagreements = 0
2038 # The contribution to the entry for these witnesses will be identical regardless of the order in which they are specified,
2039 # so we only have to calculate it once:
2040 if i > j:
2041 pbar.update(1)
2042 continue
2043 # Otherwise, calculate the number of units where both witnesses disagree:
2044 for k, vu_id in enumerate(self.variation_unit_ids):
2045 if vu_id not in substantive_variation_unit_ids_set:
2046 continue
2047 wit_1_rdg_support = self.readings_by_witness[wit_1][k]
2048 wit_2_rdg_support = self.readings_by_witness[wit_2][k]
2049 # If either witness is lacunose, then move on:
2050 if sum(wit_1_rdg_support) == 0.0 or sum(wit_2_rdg_support) == 0.0:
2051 continue
2052 # Otherwise, if the (potential) readings of the two witnesses do not overlap, then count them as disagreeing:
2053 if (
2054 sum([wit_1_rdg_support[l] * wit_2_rdg_support[l] for l in range(len(wit_1_rdg_support))])
2055 == 0.0
2056 ):
2057 disagreements += 1
2058 matrix[i, j] = disagreements
2059 matrix[j, i] = disagreements
2060 pbar.update(1)
2061 # Initialize a matrix for shared extant variation units for witnesses, and populate it if the proportion or show_ext option is specified:
2062 ext_matrix = None
2063 if proportion or show_ext:
2064 ext_matrix = self.get_ext_matrix(drop_constant=drop_constant)
2065 # If the proportion option is set, then divide every value in the matrix by the corresponding entry in the matrix of shared extant variation units:
2066 if proportion:
2067 proportion_matrix = np.full((len(witness_labels), len(witness_labels)), 0.0, dtype=float)
2068 np.divide(
2069 matrix, ext_matrix, out=proportion_matrix, where=(ext_matrix != 0)
2070 ) # division by 0 can occur if two witnesses have no overlapping units; leave their proportion as 0.0
2071 matrix = proportion_matrix
2072 # Then transform the columns of the main matrix as specified:
2073 matrix = self.transform_matrix(matrix, transform_matrix)
2074 # If the show_ext option is set, then append the number of shared extant variation units after the matrix's values:
2075 if show_ext:
2076 serialized_values = []
2077 for i, wit_1 in enumerate(witness_labels):
2078 serialized_values.append([])
2079 for j, wit_2 in enumerate(witness_labels):
2080 serialized_values[-1].append("/".join([str(matrix[i][j]), str(ext_matrix[i][j])]))
2081 matrix = np.array(serialized_values)
2082 return matrix, witness_labels
2084 def to_similarity_matrix(
2085 self,
2086 drop_constant: bool = False,
2087 proportion: bool = False,
2088 show_ext: bool = False,
2089 transform_matrix: TransformMatrixType = None,
2090 ):
2091 """Transforms this Collation into a NumPy similarity matrix between witnesses, along with an array of its labels for the witnesses.
2092 Similarities can be computed either as counts of agreements (the default setting), or as proportions of agreements over all variation units where both witnesses have singleton readings.
2093 Optionally, the count of units where both witnesses have singleton readings can be included after the count/proportion of agreements.
2095 Args:
2096 drop_constant (bool, optional): An optional flag indicating whether to ignore variation units with one substantive reading.
2097 Default value is False.
2098 proportion (bool, optional): An optional flag indicating whether or not to calculate similarities as proportions over extant, unambiguous variation units.
2099 Default value is False.
2100 show_ext: An optional flag indicating whether each cell in the matrix
2101 should include the number of their extant, unambiguous variation units after the number of agreements.
2102 Default value is False.
2103 transform_matrix (TransformMatrixType, optional): A TransformMatrixType option indicating how the columns of the matrix should be transformed.
2105 Returns:
2106 A NumPy agreement matrix with a row and column for each witness.
2107 A list of witness ID strings.
2108 """
2109 # Populate a list of sites that will correspond to columns of the sequence alignment:
2110 substantive_variation_unit_ids = self.variation_unit_ids
2111 if drop_constant:
2112 substantive_variation_unit_ids = [
2113 vu_id
2114 for vu_id in self.variation_unit_ids
2115 if len(self.substantive_readings_by_variation_unit_id[vu_id]) > 1
2116 ]
2117 substantive_variation_unit_ids_set = set(substantive_variation_unit_ids)
2118 substantive_variation_unit_reading_tuples_set = set(self.substantive_variation_unit_reading_tuples)
2119 # Initialize the output array with the appropriate dimensions:
2120 witness_labels = [wit.id for wit in self.witnesses]
2121 matrix = np.full((len(witness_labels), len(witness_labels)), 0, dtype=int) # ints of the form agreements
2122 with tqdm(total=len(self.witnesses) ** 2) as pbar:
2123 for i, wit_1 in enumerate(witness_labels):
2124 for j, wit_2 in enumerate(witness_labels):
2125 agreements = 0
2126 # The contribution to the entry for these witnesses will be identical regardless of the order in which they are specified,
2127 # so we only have to calculate it once:
2128 if i > j:
2129 pbar.update(1)
2130 continue
2131 # Otherwise, calculate the number of units where both witnesses unambiguously agree:
2132 for k, vu_id in enumerate(self.variation_unit_ids):
2133 if vu_id not in substantive_variation_unit_ids_set:
2134 continue
2135 wit_1_rdg_support = self.readings_by_witness[wit_1][k]
2136 wit_2_rdg_support = self.readings_by_witness[wit_2][k]
2137 wit_1_rdg_inds = [l for l, w in enumerate(wit_1_rdg_support) if w > 0]
2138 wit_2_rdg_inds = [l for l, w in enumerate(wit_2_rdg_support) if w > 0]
2139 if len(wit_1_rdg_inds) != 1 or len(wit_2_rdg_inds) != 1:
2140 continue
2141 if wit_1_rdg_inds[0] == wit_2_rdg_inds[0]:
2142 agreements += 1
2143 matrix[i, j] = agreements
2144 matrix[j, i] = agreements
2145 pbar.update(1)
2146 # Initialize a matrix for shared extant variation units for witnesses, and populate it if the proportion or show_ext option is specified:
2147 ext_matrix = None
2148 if proportion or show_ext:
2149 ext_matrix = self.get_ext_matrix(drop_constant=drop_constant)
2150 # If the proportion option is set, then divide every value in the matrix by the corresponding entry in the matrix of shared extant variation units:
2151 if proportion:
2152 proportion_matrix = np.full((len(witness_labels), len(witness_labels)), 0.0, dtype=float)
2153 np.divide(
2154 matrix, ext_matrix, out=proportion_matrix, where=(ext_matrix != 0)
2155 ) # division by 0 can occur if two witnesses have no overlapping units; leave their proportion as 0.0
2156 matrix = proportion_matrix
2157 # Then transform the columns of the main matrix as specified:
2158 matrix = self.transform_matrix(matrix, transform_matrix)
2159 # If the show_ext option is set, then append the number of shared extant variation units after the matrix's values:
2160 if show_ext:
2161 serialized_values = []
2162 for i, wit_1 in enumerate(witness_labels):
2163 serialized_values.append([])
2164 for j, wit_2 in enumerate(witness_labels):
2165 serialized_values[-1].append("/".join([str(matrix[i][j]), str(ext_matrix[i][j])]))
2166 matrix = np.array(serialized_values)
2167 return matrix, witness_labels
2169 def to_idf_matrix(
2170 self,
2171 drop_constant: bool = False,
2172 split_missing: SplitMissingType = None,
2173 proportion: bool = False,
2174 show_ext: bool = False,
2175 transform_matrix: TransformMatrixType = None,
2176 ):
2177 """Transforms this Collation into a NumPy matrix of agreements between witnesses weighted by inverse document frequency (IDF), along with an array of its labels for the witnesses.
2178 The IDF weight of an agreement on a given reading is the information content -log(Pr(R)) of the event R of randomly sampling a witness with that reading.
2179 The IDF-weighted agreement score for two witnesses is the sum of the IDF weights for the readings at which they agree.
2180 Where any witness is ambiguous, it contributes to its potential readings' sampling probabilities in proportion to its degrees of support for those readings.
2181 Similarly, where one or both target witnesses are ambiguous, the expected information content of their agreement is calculated based on the probabilities of their having the same reading.
2182 If a split_missing argument is supplied, then lacunae are handled in the same way.
2183 If the proportion option is set to True, then each cell will contain the mean IDF weight for its corresponding pair of witnesses over all variation units at which both are extant.
2184 (If the split_missing argument is also specified, then the mean is taken over all substantive variation units.)
2186 Args:
2187 drop_constant (bool, optional): An optional flag indicating whether to ignore variation units with one substantive reading.
2188 Default value is False.
2189 split_missing (SplitMissingType, optional): An option indicating whether or not to treat missing characters/variation units as having a contribution of 1 split over all states/readings.
2190 If not specified, then missing data is ignored (i.e., all states are 0).
2191 If "uniform", then the contribution of 1 is divided evenly over all substantive readings.
2192 If "proportional", then the contribution of 1 is divided between the readings in proportion to their support among the witnesses that are not missing.
2193 proportion (bool, optional): An optional flag indicating whether or not to calculate the mean IDF weight of each pair of witnesses' agreements as opposed to the total IDF weight of their agreements.
2194 Default value is False.
2195 show_ext: An optional flag indicating whether each cell in the matrix
2196 should include the number of their extant, unambiguous variation units after the number of agreements.
2197 Default value is False.
2198 transform_matrix (TransformMatrixType, optional): A TransformMatrixType option indicating how the columns of the matrix should be transformed.
2200 Returns:
2201 A NumPy IDF-weighted agreement matrix with a row and column for each witness.
2202 A list of witness ID strings.
2203 """
2204 # Populate a list of sites that will correspond to columns of the sequence alignment:
2205 substantive_variation_unit_ids = self.variation_unit_ids
2206 if drop_constant:
2207 substantive_variation_unit_ids = [
2208 vu_id
2209 for vu_id in self.variation_unit_ids
2210 if len(self.substantive_readings_by_variation_unit_id[vu_id]) > 1
2211 ]
2212 substantive_variation_unit_ids_set = set(substantive_variation_unit_ids)
2213 substantive_variation_unit_reading_tuples_set = set(self.substantive_variation_unit_reading_tuples)
2214 # Initialize the output array with the appropriate dimensions:
2215 witness_labels = [wit.id for wit in self.witnesses]
2216 # If the split_missing option is "proportional", then for each variation unit, keep a record of the proportion of non-missing witnesses supporting the substantive variant readings:
2217 support_proportions_by_unit = {}
2218 if split_missing == SplitMissingType.proportional:
2219 for k, vu_id in enumerate(self.variation_unit_ids):
2220 # Skip this variation unit if it is a dropped constant site:
2221 if vu_id not in substantive_variation_unit_ids_set:
2222 continue
2223 support_proportions = [0.0] * len(self.substantive_readings_by_variation_unit_id[vu_id])
2224 for i, wit in enumerate(witness_labels):
2225 rdg_support = self.readings_by_witness[wit][k]
2226 for l, w in enumerate(rdg_support):
2227 support_proportions[l] += w
2228 norm = (
2229 sum(support_proportions) if sum(support_proportions) > 0 else 1.0
2230 ) # if this variation unit has no extant witnesses (e.g., if its only witnesses are fragmentary and we have excluded them), then assume a norm of 1 to avoid division by zero
2231 for l in range(len(support_proportions)):
2232 support_proportions[l] = support_proportions[l] / norm
2233 support_proportions_by_unit[vu_id] = support_proportions
2234 # Then populate data structures mapping each variation unit's ID to normalized reading support dictionaries
2235 # and vectors of sampling probabilities for its substantive readings:
2236 normalized_reading_support_dicts_by_vu_id = {}
2237 sampling_probabilities_by_vu_id = {}
2238 for k, vu_id in enumerate(self.variation_unit_ids):
2239 # Skip this variation unit if it is a dropped constant site:
2240 if vu_id not in substantive_variation_unit_ids_set:
2241 continue
2242 # Otherwise, populate normalized reading support vector dictionaries and sampling probability vectors in this unit:
2243 normalized_reading_support_by_wit = {}
2244 sampling_probabilities = [0.0] * len(self.substantive_readings_by_variation_unit_id[vu_id])
2245 for i, wit in enumerate(witness_labels):
2246 rdg_support = self.readings_by_witness[wit][k]
2247 # Check if this reading support vector represents missing data:
2248 norm = sum(rdg_support)
2249 if norm == 0:
2250 # If this reading support vector sums to 0, then this is missing data; handle it as specified:
2251 if split_missing == SplitMissingType.uniform:
2252 rdg_support = [1 / len(rdg_support) for l in range(len(rdg_support))]
2253 elif split_missing == SplitMissingType.proportional:
2254 rdg_support = [support_proportions_by_unit[vu_id][l] for l in range(len(rdg_support))]
2255 else:
2256 # Otherwise, the data is present, though it may be ambiguous; normalize the reading probabilities to sum to 1:
2257 rdg_support = [w / norm for l, w in enumerate(rdg_support)]
2258 normalized_reading_support_by_wit[wit] = rdg_support
2259 # Then add this witness's contributions to the readings' sampling probabilities:
2260 for l, w in enumerate(normalized_reading_support_by_wit[wit]):
2261 sampling_probabilities[l] += w
2262 norm = (
2263 sum(sampling_probabilities) if sum(sampling_probabilities) > 0 else 1.0
2264 ) # if this variation unit has no extant witnesses (e.g., if its only witnesses are fragmentary and we have excluded them), then assume a norm of 1 to avoid division by zero
2265 # Otherwise, normalize the sampling probabilities so they sum to 1:
2266 sampling_probabilities = [w / norm for w in sampling_probabilities]
2267 normalized_reading_support_dicts_by_vu_id[vu_id] = normalized_reading_support_by_wit
2268 sampling_probabilities_by_vu_id[vu_id] = sampling_probabilities
2269 # Then populate the matrix with the total expected information content for agreements between each pair of witnesses:
2270 matrix = np.full((len(witness_labels), len(witness_labels)), 0, dtype=float)
2271 with tqdm(total=len(self.witnesses) ** 2) as pbar:
2272 for i, wit_1 in enumerate(witness_labels):
2273 for j, wit_2 in enumerate(witness_labels):
2274 total_information_content = 0.0
2275 # The contribution to the entry for these witnesses will be identical regardless of the order in which they are specified,
2276 # so we only have to calculate it once:
2277 if i > j:
2278 pbar.update(1)
2279 continue
2280 # Otherwise, calculate the expected information content of agreements between these witnesses in each substantive variation unit
2281 # based on the sampling probabilities of the substantive readings in the unit:
2282 for k, vu_id in enumerate(self.variation_unit_ids):
2283 if vu_id not in substantive_variation_unit_ids_set:
2284 continue
2285 wit_1_rdg_support = normalized_reading_support_dicts_by_vu_id[vu_id][wit_1]
2286 wit_2_rdg_support = normalized_reading_support_dicts_by_vu_id[vu_id][wit_2]
2287 sampling_probabilities = sampling_probabilities_by_vu_id[vu_id]
2288 # First, calculate the probability that these two witnesses agree:
2289 probability_of_agreement = sum(
2290 [wit_1_rdg_support[l] * wit_2_rdg_support[l] for l in range(len(sampling_probabilities))]
2291 )
2292 # If these witnesses do not agree at this variation unit, then this unit contributes nothing to their total score:
2293 if probability_of_agreement == 0.0:
2294 continue
2295 # Otherwise, calculate the expected information content (in bits) of their agreement given their agreement on that reading
2296 # (skipping readings with a sampling probability of 0):
2297 expected_information_content = sum(
2298 [
2299 -math.log2(sampling_probabilities[l])
2300 * (wit_1_rdg_support[l] * wit_2_rdg_support[l] / probability_of_agreement)
2301 for l in range(len(sampling_probabilities))
2302 if sampling_probabilities[l] > 0.0
2303 ]
2304 )
2305 # Then add this contribution to the total score for these two witnesses:
2306 total_information_content += expected_information_content
2307 matrix[i, j] = total_information_content
2308 matrix[j, i] = total_information_content
2309 pbar.update(1)
2310 # Initialize a matrix for shared extant variation units for witnesses, and populate it if the proportion or show_ext option is specified:
2311 ext_matrix = None
2312 if proportion or show_ext:
2313 ext_matrix = self.get_ext_matrix(drop_constant=drop_constant, split_missing=split_missing)
2314 # If the proportion option is set, then divide every value in the matrix by the corresponding entry in the matrix of shared extant variation units:
2315 if proportion:
2316 proportion_matrix = np.full((len(witness_labels), len(witness_labels)), 0.0, dtype=float)
2317 np.divide(
2318 matrix, ext_matrix, out=proportion_matrix, where=(ext_matrix != 0)
2319 ) # division by 0 can occur if two witnesses have no overlapping units; leave their proportion as 0.0
2320 matrix = proportion_matrix
2321 # Then transform the columns of the main matrix as specified:
2322 matrix = self.transform_matrix(matrix, transform_matrix)
2323 # If the show_ext option is set, then append the number of shared extant variation units after the matrix's values:
2324 if show_ext:
2325 serialized_values = []
2326 for i, wit_1 in enumerate(witness_labels):
2327 serialized_values.append([])
2328 for j, wit_2 in enumerate(witness_labels):
2329 serialized_values[-1].append("/".join([str(matrix[i][j]), str(ext_matrix[i][j])]))
2330 matrix = np.array(serialized_values)
2331 return matrix, witness_labels
2333 def to_mi_matrix(
2334 self,
2335 drop_constant: bool = False,
2336 split_missing: SplitMissingType = None,
2337 proportion: bool = False,
2338 show_ext: bool = False,
2339 transform_matrix: TransformMatrixType = None,
2340 ):
2341 """Transforms this Collation into a NumPy matrix of the total mutual information (MI), in bits, between witnesses over all variation units, along with an array of its labels for the witnesses.
2342 This is equivalent to the total Kullback-Leibler divergence of the joint distribution of the witnesses' observed readings
2343 from the joint distribution of their expected readings under the assumption that the witnesses are independent, taken over all variation units.
2344 The value of 0 if and only if the witnesses are completely independent.
2346 Args:
2347 drop_constant (bool, optional): An optional flag indicating whether to ignore variation units with one substantive reading.
2348 Default value is False.
2349 split_missing (SplitMissingType, optional): An option indicating whether or not to treat missing characters/variation units as having a contribution of 1 split over all states/readings.
2350 If not specified, then missing data is ignored (i.e., all states are 0).
2351 If "uniform", then the contribution of 1 is divided evenly over all substantive readings.
2352 If "proportional", then the contribution of 1 is divided between the readings in proportion to their support among the witnesses that are not missing.
2353 proportion (bool, optional): An optional flag indicating whether or not to calculate the mean IDF weight of each pair of witnesses' agreements as opposed to the total IDF weight of their agreements.
2354 Default value is False.
2355 show_ext: An optional flag indicating whether each cell in the matrix
2356 should include the number of their extant, unambiguous variation units after the number of agreements.
2357 Default value is False.
2358 transform_matrix (TransformMatrixType, optional): A TransformMatrixType option indicating how the columns of the matrix should be transformed.
2360 Returns:
2361 A NumPy MI matrix with a row and column for each witness.
2362 A list of witness ID strings.
2363 """
2364 # Populate a list of sites that will correspond to columns of the sequence alignment:
2365 substantive_variation_unit_ids = self.variation_unit_ids
2366 if drop_constant:
2367 substantive_variation_unit_ids = [
2368 vu_id
2369 for vu_id in self.variation_unit_ids
2370 if len(self.substantive_readings_by_variation_unit_id[vu_id]) > 1
2371 ]
2372 substantive_variation_unit_ids_set = set(substantive_variation_unit_ids)
2373 substantive_variation_unit_reading_tuples_set = set(self.substantive_variation_unit_reading_tuples)
2374 # Initialize the output array with the appropriate dimensions:
2375 witness_labels = [wit.id for wit in self.witnesses]
2376 # If the split_missing option is "proportional", then for each variation unit, keep a record of the proportion of non-missing witnesses supporting the substantive variant readings:
2377 support_proportions_by_unit = {}
2378 if split_missing == SplitMissingType.proportional:
2379 for k, vu_id in enumerate(self.variation_unit_ids):
2380 # Skip this variation unit if it is a dropped constant site:
2381 if vu_id not in substantive_variation_unit_ids_set:
2382 continue
2383 support_proportions = [0.0] * len(self.substantive_readings_by_variation_unit_id[vu_id])
2384 for i, wit in enumerate(witness_labels):
2385 rdg_support = self.readings_by_witness[wit][k]
2386 for l, w in enumerate(rdg_support):
2387 support_proportions[l] += w
2388 norm = (
2389 sum(support_proportions) if sum(support_proportions) > 0 else 1.0
2390 ) # if this variation unit has no extant witnesses (e.g., if its only witnesses are fragmentary and we have excluded them), then assume a norm of 1 to avoid division by zero
2391 for l in range(len(support_proportions)):
2392 support_proportions[l] = support_proportions[l] / norm
2393 support_proportions_by_unit[vu_id] = support_proportions
2394 # Then populate data structures mapping each variation unit's ID to normalized reading support dictionaries,
2395 # vectors of sampling probabilities for its substantive readings, and expected joint probability matrices:
2396 normalized_reading_support_dicts_by_vu_id = {}
2397 sampling_probabilities_by_vu_id = {}
2398 expected_joint_probabilities_by_vu_id = {}
2399 for k, vu_id in enumerate(self.variation_unit_ids):
2400 # Skip this variation unit if it is a dropped constant site:
2401 if vu_id not in substantive_variation_unit_ids_set:
2402 continue
2403 # Otherwise, populate normalized reading support vector dictionaries and sampling probability vectors in this unit:
2404 normalized_reading_support_by_wit = {}
2405 sampling_probabilities = [0.0] * len(self.substantive_readings_by_variation_unit_id[vu_id])
2406 for i, wit in enumerate(witness_labels):
2407 rdg_support = self.readings_by_witness[wit][k]
2408 # Check if this reading support vector represents missing data:
2409 norm = sum(rdg_support)
2410 if norm == 0:
2411 # If this reading support vector sums to 0, then this is missing data; handle it as specified:
2412 if split_missing == SplitMissingType.uniform:
2413 rdg_support = [1 / len(rdg_support) for l in range(len(rdg_support))]
2414 elif split_missing == SplitMissingType.proportional:
2415 rdg_support = [support_proportions_by_unit[vu_id][l] for l in range(len(rdg_support))]
2416 else:
2417 # Otherwise, the data is present, though it may be ambiguous; normalize the reading probabilities to sum to 1:
2418 rdg_support = [w / norm for l, w in enumerate(rdg_support)]
2419 normalized_reading_support_by_wit[wit] = rdg_support
2420 # Then add this witness's contributions to the readings' sampling probabilities:
2421 for l, w in enumerate(normalized_reading_support_by_wit[wit]):
2422 sampling_probabilities[l] += w
2423 norm = (
2424 sum(sampling_probabilities) if sum(sampling_probabilities) > 0 else 1.0
2425 ) # if this variation unit has no extant witnesses (e.g., if its only witnesses are fragmentary and we have excluded them), then assume a norm of 1 to avoid division by zero
2426 # Otherwise, normalize the sampling probabilities so they sum to 1:
2427 sampling_probabilities = [w / norm for w in sampling_probabilities]
2428 normalized_reading_support_dicts_by_vu_id[vu_id] = normalized_reading_support_by_wit
2429 sampling_probabilities_by_vu_id[vu_id] = sampling_probabilities
2430 # Then populate a contingency table for the expected probabilities of joint support in this unit:
2431 expected_joint_probabilities = np.full(
2432 (len(sampling_probabilities), len(sampling_probabilities)), 0, dtype=float
2433 )
2434 for l1 in range(len(sampling_probabilities)):
2435 for l2 in range(len(sampling_probabilities)):
2436 expected_joint_probabilities[l1, l2] = sampling_probabilities[l1] * sampling_probabilities[l2]
2437 expected_joint_probabilities_by_vu_id[vu_id] = expected_joint_probabilities
2438 # Then populate the matrix one variation unit at a time:
2439 matrix = np.full((len(witness_labels), len(witness_labels)), 0, dtype=float)
2440 with tqdm(total=len(self.witnesses) ** 2) as pbar:
2441 # Then calculate the mutual information contribution for each pair of witnesses:
2442 for i, wit_1 in enumerate(witness_labels):
2443 for j, wit_2 in enumerate(witness_labels):
2444 total_mutual_information = 0.0
2445 # The contribution to the entry for these witnesses will be identical regardless of the order in which they are specified,
2446 # so we only have to calculate it once:
2447 if i > j:
2448 pbar.update(1)
2449 continue
2450 # Otherwise, calculate the mutual information between these witnesses in each substantive variation unit:
2451 for k, vu_id in enumerate(self.variation_unit_ids):
2452 if vu_id not in substantive_variation_unit_ids_set:
2453 continue
2454 wit_1_rdg_support = normalized_reading_support_dicts_by_vu_id[vu_id][wit_1]
2455 wit_2_rdg_support = normalized_reading_support_dicts_by_vu_id[vu_id][wit_2]
2456 sampling_probabilities = sampling_probabilities_by_vu_id[vu_id]
2457 expected_joint_probabilities = expected_joint_probabilities_by_vu_id[vu_id]
2458 # If either witness has an all-zeroes vector (because it is lacunose in this unit), then we can skip these witnesses here:
2459 if sum(wit_1_rdg_support) == 0 or sum(wit_2_rdg_support) == 0:
2460 continue
2461 # Otherwise, populate a contingency table for the observed probabilities of joint support in this unit:
2462 observed_joint_probabilities = np.full(
2463 (len(sampling_probabilities), len(sampling_probabilities)), 0, dtype=float
2464 )
2465 for l1, w1 in enumerate(wit_1_rdg_support):
2466 for l2, w2 in enumerate(wit_2_rdg_support):
2467 observed_joint_probabilities[l1, l2] = w1 * w2
2468 # Then calculate the mutual information using the expected and observed distribution matrices:
2469 mutual_information = 0.0
2470 for l1 in range(len(sampling_probabilities)):
2471 for l2 in range(len(sampling_probabilities)):
2472 observed = observed_joint_probabilities[l1, l2]
2473 expected = expected_joint_probabilities[l1, l2]
2474 if observed == 0:
2475 continue
2476 mutual_information += observed * math.log2(observed / expected)
2477 # Then add this mutual information to the total for these two witnesses:
2478 total_mutual_information += mutual_information
2479 matrix[i, j] += total_mutual_information
2480 matrix[j, i] += total_mutual_information
2481 pbar.update(1)
2482 # Initialize a matrix for shared extant variation units for witnesses, and populate it if the proportion or show_ext option is specified:
2483 ext_matrix = None
2484 if proportion or show_ext:
2485 ext_matrix = self.get_ext_matrix(drop_constant=drop_constant, split_missing=split_missing)
2486 # If the proportion option is set, then divide every value in the matrix by the corresponding entry in the matrix of shared extant variation units:
2487 if proportion:
2488 proportion_matrix = np.full((len(witness_labels), len(witness_labels)), 0.0, dtype=float)
2489 np.divide(
2490 matrix, ext_matrix, out=proportion_matrix, where=(ext_matrix != 0)
2491 ) # division by 0 can occur if two witnesses have no overlapping units; leave their proportion as 0.0
2492 matrix = proportion_matrix
2493 # Then transform the columns of the main matrix as specified:
2494 matrix = self.transform_matrix(matrix, transform_matrix)
2495 # If the show_ext option is set, then append the number of shared extant variation units after the matrix's values:
2496 if show_ext:
2497 serialized_values = []
2498 for i, wit_1 in enumerate(witness_labels):
2499 serialized_values.append([])
2500 for j, wit_2 in enumerate(witness_labels):
2501 serialized_values[-1].append("/".join([str(matrix[i][j]), str(ext_matrix[i][j])]))
2502 matrix = np.array(serialized_values)
2503 return matrix, witness_labels
2505 def to_nexus_table(self, drop_constant: bool = False, ambiguous_as_missing: bool = False):
2506 """Returns this Collation in the form of a table with rows for taxa, columns for characters, and reading IDs in cells.
2508 Args:
2509 drop_constant (bool, optional): An optional flag indicating whether to ignore variation units with one substantive reading.
2510 Default value is False.
2511 ambiguous_as_missing (bool, optional): An optional flag indicating whether to treat all ambiguous states as missing data.
2512 Default value is False.
2514 Returns:
2515 A NumPy array with rows for taxa, columns for characters, and reading IDs in cells.
2516 A list of substantive reading ID strings.
2517 A list of witness ID strings.
2518 """
2519 # Populate a list of sites that will correspond to columns of the sequence alignment:
2520 substantive_variation_unit_ids = self.variation_unit_ids
2521 if drop_constant:
2522 substantive_variation_unit_ids = [
2523 vu_id
2524 for vu_id in self.variation_unit_ids
2525 if len(self.substantive_readings_by_variation_unit_id[vu_id]) > 1
2526 ]
2527 substantive_variation_unit_ids_set = set(substantive_variation_unit_ids)
2528 substantive_variation_unit_reading_tuples_set = set(self.substantive_variation_unit_reading_tuples)
2529 # In a first pass, populate a dictionary mapping (variation unit index, reading index) tuples from the readings_by_witness dictionary
2530 # to the readings' IDs:
2531 reading_ids_by_indices = {}
2532 for j, vu in enumerate(self.variation_units):
2533 if vu.id not in substantive_variation_unit_ids_set:
2534 continue
2535 k = 0
2536 for rdg in vu.readings:
2537 key = tuple([vu.id, rdg.id])
2538 if key not in substantive_variation_unit_reading_tuples_set:
2539 continue
2540 indices = tuple([j, k])
2541 reading_ids_by_indices[indices] = rdg.id
2542 k += 1
2543 # Initialize the output array with the appropriate dimensions:
2544 missing_symbol = '?'
2545 witness_labels = [wit.id for wit in self.witnesses]
2546 matrix = np.full(
2547 (len(witness_labels), len(substantive_variation_unit_ids)), missing_symbol, dtype=object
2548 ) # use dtype=object because the maximum string length is not known up front
2549 # Then populate it with the appropriate values:
2550 with tqdm(total=len(self.witnesses)) as pbar:
2551 row_ind = 0
2552 for i, wit in enumerate(self.witnesses):
2553 col_ind = 0
2554 for j, vu in enumerate(self.variation_units):
2555 if vu.id not in substantive_variation_unit_ids_set:
2556 continue
2557 rdg_support = self.readings_by_witness[wit.id][j]
2558 # If this reading support vector sums to 0, then this is missing data; handle it as specified:
2559 if sum(rdg_support) == 0:
2560 matrix[row_ind, col_ind] = missing_symbol
2561 # Otherwise, add its coefficients normally:
2562 else:
2563 rdg_inds = [
2564 k for k, w in enumerate(rdg_support) if w > 0
2565 ] # the index list consists of the indices of all readings with any degree of certainty assigned to them
2566 # For singleton readings, just print the reading ID:
2567 if len(rdg_inds) == 1:
2568 k = rdg_inds[0]
2569 matrix[row_ind, col_ind] = reading_ids_by_indices[(j, k)]
2570 # For multiple readings, print the corresponding reading IDs in braces or the missing symbol depending on input settings:
2571 else:
2572 if ambiguous_as_missing:
2573 matrix[row_ind, col_ind] = missing_symbol
2574 else:
2575 matrix[row_ind, col_ind] = "{%s}" % " ".join(
2576 [reading_ids_by_indices[(j, k)] for k in rdg_inds]
2577 )
2578 col_ind += 1
2579 row_ind += 1
2580 pbar.update(1)
2581 return matrix, witness_labels, substantive_variation_unit_ids
2583 def to_long_table(self, drop_constant: bool = False):
2584 """Returns this Collation in the form of a long table with columns for taxa, characters, reading indices, and reading values.
2585 Note that this method treats ambiguous readings as missing data.
2587 Args:
2588 drop_constant (bool, optional): An optional flag indicating whether to ignore variation units with one substantive reading.
2589 Default value is False.
2591 Returns:
2592 A NumPy array with columns for taxa, characters, reading indices, and reading values, and rows for each combination of these values in the matrix.
2593 A list of column label strings.
2594 """
2595 # Populate a list of sites that will correspond to columns of the sequence alignment:
2596 substantive_variation_unit_ids = self.variation_unit_ids
2597 if drop_constant:
2598 substantive_variation_unit_ids = [
2599 vu_id
2600 for vu_id in self.variation_unit_ids
2601 if len(self.substantive_readings_by_variation_unit_id[vu_id]) > 1
2602 ]
2603 substantive_variation_unit_ids_set = set(substantive_variation_unit_ids)
2604 substantive_variation_unit_reading_tuples_set = set(self.substantive_variation_unit_reading_tuples)
2605 # Initialize the outputs:
2606 column_labels = ["taxon", "character", "state", "value"]
2607 long_table_list = []
2608 # Populate a dictionary mapping (variation unit index, reading index) tuples to reading texts:
2609 reading_texts_by_indices = {}
2610 for j, vu in enumerate(self.variation_units):
2611 if vu.id not in substantive_variation_unit_ids_set:
2612 continue
2613 k = 0
2614 for rdg in vu.readings:
2615 key = tuple([vu.id, rdg.id])
2616 if key not in substantive_variation_unit_reading_tuples_set:
2617 continue
2618 indices = tuple([j, k])
2619 reading_texts_by_indices[indices] = rdg.text
2620 k += 1
2621 # Then populate the output list with the appropriate values:
2622 witness_labels = [wit.id for wit in self.witnesses]
2623 missing_symbol = '?'
2624 with tqdm(total=len(self.witnesses)) as pbar:
2625 for i, wit in enumerate(self.witnesses):
2626 row_ind = 0
2627 for j, vu_id in enumerate(self.variation_unit_ids):
2628 if vu_id not in substantive_variation_unit_ids_set:
2629 continue
2630 rdg_support = self.readings_by_witness[wit.id][j]
2631 # Populate a list of nonzero coefficients for this reading support vector:
2632 rdg_inds = [k for k, w in enumerate(rdg_support) if w > 0]
2633 # If this list does not consist of exactly one reading, then treat it as missing data:
2634 if len(rdg_inds) != 1:
2635 long_table_list.append([wit.id, vu_id, missing_symbol, missing_symbol])
2636 continue
2637 k = rdg_inds[0]
2638 rdg_text = reading_texts_by_indices[(j, k)]
2639 # Replace empty reading texts with the omission placeholder:
2640 if rdg_text == "":
2641 rdg_text = "om."
2642 long_table_list.append([wit.id, vu_id, k, rdg_text])
2643 pbar.update(1)
2644 # Then convert the long table entries list to a NumPy array:
2645 long_table = np.array(long_table_list)
2646 return long_table, column_labels
2648 def to_dataframe(
2649 self,
2650 drop_constant: bool = False,
2651 ambiguous_as_missing: bool = False,
2652 proportion: bool = False,
2653 table_type: TableType = TableType.matrix,
2654 split_missing: SplitMissingType = None,
2655 transform_matrix: TransformMatrixType = None,
2656 show_ext: bool = False,
2657 ):
2658 """Returns this Collation in the form of a Pandas DataFrame array, including the appropriate row and column labels.
2660 Args:
2661 drop_constant (bool, optional): An optional flag indicating whether to ignore variation units with one substantive reading.
2662 Default value is False.
2663 ambiguous_as_missing (bool, optional): An optional flag indicating whether to treat all ambiguous states as missing data.
2664 Default value is False.
2665 proportion (bool, optional): An optional flag indicating whether or not to calculate distances as proportions over extant, unambiguous variation units.
2666 Default value is False.
2667 table_type (TableType, optional): A TableType option indicating which type of tabular output to generate.
2668 Only applicable for tabular outputs.
2669 Default value is "matrix".
2670 split_missing (SplitMissingType, optional): An option indicating whether or not to treat missing characters/variation units as having a contribution of 1 split over all states/readings.
2671 If not specified, then missing data is ignored (i.e., all states are 0).
2672 If "uniform", then the contribution of 1 is divided evenly over all substantive readings.
2673 If "proportional", then the contribution of 1 is divided between the readings in proportion to their support among the witnesses that are not missing.
2674 Only applicable for table types "matrix" and "idf".
2675 transform_matrix (TransformMatrixType, optional): A TransformMatrixType option indicating how the columns of a witness-to-witness matrix output should be transformed.
2676 Only applicable for tabular outputs in which the rows and columns correspond to the witnesses in the collation.
2677 show_ext: An optional flag indicating whether each cell in the matrix
2678 should include the number of their extant, unambiguous variation units after the number of their disagreements/agreements.
2679 Only applicable for tabular output formats of type \"distance\" or \"similarity\".
2680 Default value is False.
2682 Returns:
2683 A Pandas DataFrame corresponding to a collation matrix with reading frequencies or a long table with discrete reading states.
2684 """
2685 df = None
2686 # Proceed based on the table type:
2687 if table_type == TableType.matrix:
2688 # Convert the collation to a NumPy array and get its row and column labels first:
2689 matrix, reading_labels, witness_labels = self.to_numpy(
2690 drop_constant=drop_constant, split_missing=split_missing
2691 )
2692 df = pd.DataFrame(matrix, index=reading_labels, columns=witness_labels)
2693 elif table_type == TableType.distance:
2694 # Convert the collation to a NumPy array and get its row and column labels first:
2695 matrix, witness_labels = self.to_distance_matrix(
2696 drop_constant=drop_constant, proportion=proportion, transform_matrix=transform_matrix, show_ext=show_ext
2697 )
2698 df = pd.DataFrame(matrix, index=witness_labels, columns=witness_labels)
2699 elif table_type == TableType.similarity:
2700 # Convert the collation to a NumPy array and get its row and column labels first:
2701 matrix, witness_labels = self.to_similarity_matrix(
2702 drop_constant=drop_constant, proportion=proportion, transform_matrix=transform_matrix, show_ext=show_ext
2703 )
2704 df = pd.DataFrame(matrix, index=witness_labels, columns=witness_labels)
2705 elif table_type == TableType.idf:
2706 # Convert the collation to a NumPy array and get its row and column labels first:
2707 matrix, witness_labels = self.to_idf_matrix(
2708 drop_constant=drop_constant,
2709 split_missing=split_missing,
2710 proportion=proportion,
2711 transform_matrix=transform_matrix,
2712 show_ext=show_ext,
2713 )
2714 df = pd.DataFrame(matrix, index=witness_labels, columns=witness_labels)
2715 elif table_type == TableType.mi:
2716 # Convert the collation to a NumPy array and get its row and column labels first:
2717 matrix, witness_labels = self.to_mi_matrix(
2718 drop_constant=drop_constant,
2719 split_missing=split_missing,
2720 proportion=proportion,
2721 transform_matrix=transform_matrix,
2722 show_ext=show_ext,
2723 )
2724 df = pd.DataFrame(matrix, index=witness_labels, columns=witness_labels)
2725 elif table_type == TableType.nexus:
2726 # Convert the collation to a NumPy array and get its row and column labels first:
2727 matrix, witness_labels, vu_labels = self.to_nexus_table(
2728 drop_constant=drop_constant, ambiguous_as_missing=ambiguous_as_missing
2729 )
2730 df = pd.DataFrame(matrix, index=witness_labels, columns=vu_labels)
2731 elif table_type == TableType.long:
2732 # Convert the collation to a long table and get its column labels first:
2733 long_table, column_labels = self.to_long_table(drop_constant=drop_constant)
2734 df = pd.DataFrame(long_table, columns=column_labels)
2735 return df
2737 def to_csv(
2738 self,
2739 file_addr: Union[Path, str],
2740 drop_constant: bool = False,
2741 ambiguous_as_missing: bool = False,
2742 proportion: bool = False,
2743 table_type: TableType = TableType.matrix,
2744 split_missing: SplitMissingType = None,
2745 transform_matrix: TransformMatrixType = None,
2746 show_ext: bool = False,
2747 **kwargs
2748 ):
2749 """Writes this Collation to a comma-separated value (CSV) file with the given address.
2751 If your witness IDs are numeric (e.g., Gregory-Aland numbers), then they will be written in full to the CSV file, but Excel will likely interpret them as numbers and truncate any leading zeroes!
2753 Args:
2754 file_addr: A string representing the path to an output CSV file; the file type should be .csv.
2755 drop_constant: An optional flag indicating whether to ignore variation units with one substantive reading.
2756 Default value is False.
2757 ambiguous_as_missing: An optional flag indicating whether to treat all ambiguous states as missing data.
2758 Default value is False.
2759 proportion: An optional flag indicating whether or not to calculate distances as proportions over extant, unambiguous variation units.
2760 Default value is False.
2761 table_type: A TableType option indicating which type of tabular output to generate.
2762 Only applicable for tabular outputs.
2763 Default value is "matrix".
2764 split_missing: An option indicating whether or not to treat missing characters/variation units as having a contribution of 1 split over all states/readings.
2765 If not specified, then missing data is ignored (i.e., all states are 0).
2766 If "uniform", then the contribution of 1 is divided evenly over all substantive readings.
2767 If "proportional", then the contribution of 1 is divided between the readings in proportion to their support among the witnesses that are not missing.
2768 Only applicable for table types "matrix" and "idf".
2769 transform_matrix: A TransformMatrixType option indicating how the columns of a witness-to-witness matrix output should be transformed.
2770 Only applicable for tabular outputs in which the rows and columns correspond to the witnesses in the collation.
2771 show_ext: An optional flag indicating whether each cell in the matrix
2772 should include the number of their extant, unambiguous variation units after the number of their disagreements/agreements.
2773 Only applicable for tabular output formats of type \"distance\" or \"similarity\".
2774 Default value is False.
2775 **kwargs: Keyword arguments for pandas.DataFrame.to_csv.
2776 """
2777 # Convert the collation to a Pandas DataFrame first:
2778 df = self.to_dataframe(
2779 drop_constant=drop_constant,
2780 ambiguous_as_missing=ambiguous_as_missing,
2781 proportion=proportion,
2782 table_type=table_type,
2783 split_missing=split_missing,
2784 show_ext=show_ext,
2785 transform_matrix=transform_matrix,
2786 )
2787 # Generate all parent folders for this file that don't already exist:
2788 Path(file_addr).parent.mkdir(parents=True, exist_ok=True)
2789 # Proceed based on the table type:
2790 if table_type == TableType.long:
2791 return df.to_csv(
2792 file_addr, encoding="utf-8-sig", index=False, **kwargs
2793 ) # add BOM to start of file so that Excel will know to read it as Unicode
2794 return df.to_csv(
2795 file_addr, encoding="utf-8-sig", **kwargs
2796 ) # add BOM to start of file so that Excel will know to read it as Unicode
2798 def to_excel(
2799 self,
2800 file_addr: Union[Path, str],
2801 drop_constant: bool = False,
2802 ambiguous_as_missing: bool = False,
2803 proportion: bool = False,
2804 table_type: TableType = TableType.matrix,
2805 split_missing: SplitMissingType = None,
2806 transform_matrix: TransformMatrixType = None,
2807 show_ext: bool = False,
2808 ):
2809 """Writes this Collation to an Excel (.xlsx) file with the given address.
2811 Since Pandas is deprecating its support for xlwt, specifying an output in old Excel (.xls) output is not recommended.
2813 Args:
2814 file_addr: A string representing the path to an output Excel file; the file type should be .xlsx.
2815 drop_constant: An optional flag indicating whether to ignore variation units with one substantive reading.
2816 Default value is False.
2817 ambiguous_as_missing: An optional flag indicating whether to treat all ambiguous states as missing data.
2818 Default value is False.
2819 proportion: An optional flag indicating whether or not to calculate distances as proportions over extant, unambiguous variation units.
2820 Default value is False.
2821 table_type: A TableType option indicating which type of tabular output to generate.
2822 Only applicable for tabular outputs.
2823 Default value is "matrix".
2824 split_missing: An option indicating whether or not to treat missing characters/variation units as having a contribution of 1 split over all states/readings.
2825 If not specified, then missing data is ignored (i.e., all states are 0).
2826 If "uniform", then the contribution of 1 is divided evenly over all substantive readings.
2827 If "proportional", then the contribution of 1 is divided between the readings in proportion to their support among the witnesses that are not missing.
2828 Only applicable for table types "matrix" and "idf".
2829 transform_matrix: A TransformMatrixType option indicating how the columns of a witness-to-witness matrix output should be transformed.
2830 Only applicable for tabular outputs in which the rows and columns correspond to the witnesses in the collation.
2831 show_ext: An optional flag indicating whether each cell in the matrix
2832 should include the number of their extant, unambiguous variation units after the number of their disagreements/agreements.
2833 Only applicable for tabular output formats of type \"distance\" or \"similarity\".
2834 Default value is False.
2835 """
2836 # Convert the collation to a Pandas DataFrame first:
2837 df = self.to_dataframe(
2838 drop_constant=drop_constant,
2839 ambiguous_as_missing=ambiguous_as_missing,
2840 proportion=proportion,
2841 table_type=table_type,
2842 split_missing=split_missing,
2843 show_ext=show_ext,
2844 transform_matrix=transform_matrix,
2845 )
2846 # Generate all parent folders for this file that don't already exist:
2847 Path(file_addr).parent.mkdir(parents=True, exist_ok=True)
2848 # Proceed based on the table type:
2849 if table_type == TableType.long:
2850 return df.to_excel(file_addr, index=False)
2851 return df.to_excel(file_addr)
2853 def to_phylip_matrix(
2854 self,
2855 file_addr: Union[Path, str],
2856 drop_constant: bool = False,
2857 proportion: bool = False,
2858 table_type: TableType = TableType.distance,
2859 show_ext: bool = False,
2860 ):
2861 """Writes this Collation as a PHYLIP-formatted distance/similarity matrix to the file with the given address.
2863 Args:
2864 file_addr: A string representing the path to an output PHYLIP file; the file type should be .ph or .phy.
2865 drop_constant: An optional flag indicating whether to ignore variation units with one substantive reading.
2866 Default value is False.
2867 proportion: An optional flag indicating whether or not to calculate distances as proportions over extant, unambiguous variation units.
2868 Default value is False.
2869 table_type: A TableType option indicating which type of tabular output to generate.
2870 For PHYLIP-formatted outputs, distance and similarity matrices are the only supported table types.
2871 Default value is "distance".
2872 show_ext: An optional flag indicating whether each cell in the matrix
2873 should include the number of their extant, unambiguous variation units after the number of their disagreements/agreements.
2874 Only applicable for tabular output formats of type \"distance\" or \"similarity\".
2875 Default value is False.
2876 """
2877 # Convert the collation to a Pandas DataFrame first:
2878 matrix = None
2879 witness_labels = []
2880 # Proceed based on the table type:
2881 if table_type == TableType.distance:
2882 # Convert the collation to a NumPy array and get its row and column labels first:
2883 matrix, witness_labels = self.to_distance_matrix(
2884 drop_constant=drop_constant, proportion=proportion, show_ext=show_ext
2885 )
2886 elif table_type == TableType.similarity:
2887 # Convert the collation to a NumPy array and get its row and column labels first:
2888 matrix, witness_labels = self.to_similarity_matrix(
2889 drop_constant=drop_constant, proportion=proportion, show_ext=show_ext
2890 )
2891 # Generate all parent folders for this file that don't already exist:
2892 Path(file_addr).parent.mkdir(parents=True, exist_ok=True)
2893 with open(file_addr, "w", encoding="utf-8") as f:
2894 # The first line contains the number of taxa:
2895 f.write("%d\n" % len(witness_labels))
2896 # Every subsequent line contains a witness label, followed by the values in its row of the matrix:
2897 for i, wit_id in enumerate(witness_labels):
2898 wit_label = slugify(wit_id, lowercase=False, allow_unicode=True, separator='_')
2899 f.write("%s %s\n" % (wit_label, " ".join([str(v) for v in matrix[i]])))
2900 return
2902 def get_stemma_symbols(self):
2903 """Returns a list of one-character symbols needed to represent the states of all substantive readings in stemma format.
2905 The number of symbols equals the maximum number of substantive readings at any variation unit.
2907 Returns:
2908 A list of individual characters representing states in readings.
2909 """
2910 possible_symbols = (
2911 list(string.digits) + list(string.ascii_lowercase) + list(string.ascii_uppercase)
2912 ) # NOTE: the maximum number of symbols allowed in stemma format (other than "?" and "-") is 62
2913 # The number of symbols needed is equal to the length of the longest substantive reading vector:
2914 nsymbols = 0
2915 # If there are no witnesses, then no symbols are needed at all:
2916 if len(self.witnesses) == 0:
2917 return []
2918 wit_id = self.witnesses[0].id
2919 for i, vu_id in enumerate(self.variation_unit_ids):
2920 rdg_support = self.readings_by_witness[wit_id][i]
2921 if len(rdg_support) > len(possible_symbols):
2922 raise ValueError(
2923 f"ERROR: too many substantive readings at variation unit '{vu_id}' to represent in stemma format."
2924 )
2925 nsymbols = max(nsymbols, len(rdg_support))
2926 stemma_symbols = possible_symbols[:nsymbols]
2927 return stemma_symbols
2929 def to_stemma(self, file_addr: Union[Path, str]):
2930 """Writes this Collation to a stemma file without an extension and a Chron file (containing low, middle, and high dates for all witnesses) without an extension.
2932 Since this format does not support ambiguous states, all reading vectors with anything other than one nonzero entry will be interpreted as lacunose.
2933 If an interpGrp for weights is specified in the TEI XML collation, then the weights for the interp elements will be used as weights
2934 for the variation units that specify them in their ana attribute.
2936 Args:
2937 file_addr: A string representing the path to an output stemma prep file; the file should have no extension.
2938 The accompanying chron file will match this file name, except that it will have "_chron" appended to the end.
2939 drop_constant: An optional flag indicating whether to ignore variation units with one substantive reading.
2940 """
2941 # Populate a list of sites that will correspond to columns of the sequence alignment
2942 # (by default, constant sites are dropped):
2943 substantive_variation_unit_ids = [
2944 vu_id for vu_id in self.variation_unit_ids if len(self.substantive_readings_by_variation_unit_id[vu_id]) > 1
2945 ]
2946 substantive_variation_unit_ids_set = set(substantive_variation_unit_ids)
2947 substantive_variation_unit_reading_tuples_set = set(self.substantive_variation_unit_reading_tuples)
2948 # In a first pass, populate a dictionary mapping (variation unit index, reading index) tuples from the readings_by_witness dictionary
2949 # to the readings' texts:
2950 reading_texts_by_indices = {}
2951 for j, vu in enumerate(self.variation_units):
2952 if vu.id not in substantive_variation_unit_ids_set:
2953 continue
2954 k = 0
2955 for rdg in vu.readings:
2956 key = tuple([vu.id, rdg.id])
2957 if key not in substantive_variation_unit_reading_tuples_set:
2958 continue
2959 indices = tuple([j, k])
2960 reading_texts_by_indices[indices] = rdg.text
2961 k += 1
2962 # In a second pass, populate another dictionary mapping (variation unit index, reading index) tuples from the readings_by_witness dictionary
2963 # to the witnesses exclusively supporting those readings:
2964 reading_wits_by_indices = {}
2965 for indices in reading_texts_by_indices:
2966 reading_wits_by_indices[indices] = []
2967 for i, wit in enumerate(self.witnesses):
2968 for j, vu_id in enumerate(self.variation_unit_ids):
2969 if vu_id not in substantive_variation_unit_ids_set:
2970 continue
2971 rdg_support = self.readings_by_witness[wit.id][j]
2972 # If this witness does not exclusively support exactly one reading at this unit, then treat it as lacunose:
2973 if len([k for k, w in enumerate(rdg_support) if w > 0]) != 1:
2974 continue
2975 k = rdg_support.index(1)
2976 indices = tuple([j, k])
2977 reading_wits_by_indices[indices].append(wit.id)
2978 # In a third pass, write to the stemma file:
2979 symbols = self.get_stemma_symbols()
2980 Path(file_addr).parent.mkdir(
2981 parents=True, exist_ok=True
2982 ) # generate all parent folders for this file that don't already exist
2983 chron_file_addr = str(file_addr) + "_chron"
2984 with open(file_addr, "w", encoding="utf-8") as f:
2985 # Start with the witness list:
2986 f.write(
2987 "* %s ;\n\n"
2988 % " ".join(
2989 [slugify(wit.id, lowercase=False, allow_unicode=True, separator='_') for wit in self.witnesses]
2990 )
2991 )
2992 # f.write("^ %s\n\n" % chron_file_addr) #write the relative path to the chron file
2993 f.write(
2994 "^ %s\n\n" % ("." + os.sep + Path(chron_file_addr).name)
2995 ) # write the relative path to the chron file
2996 # Then add a line indicating that all witnesses are lacunose unless they are specified explicitly:
2997 f.write("= $? $* ;\n\n")
2998 with tqdm(total=len(self.variation_unit_ids)) as pbar:
2999 # Then proceed for each variation unit:
3000 for j, vu_id in enumerate(self.variation_unit_ids):
3001 if vu_id not in substantive_variation_unit_ids_set:
3002 pbar.update(1)
3003 continue
3004 # Print the variation unit ID first:
3005 f.write("@ %s\n" % vu_id)
3006 # In a first pass, print the texts of all readings enclosed in brackets:
3007 f.write("[ ")
3008 k = 0
3009 while True:
3010 indices = tuple([j, k])
3011 if indices not in reading_texts_by_indices:
3012 break
3013 text = slugify(
3014 reading_texts_by_indices[indices], lowercase=False, allow_unicode=True, separator='.'
3015 )
3016 # Denote omissions by en-dashes:
3017 if text == "":
3018 text = "\u2013"
3019 # The first reading should not be preceded by anything:
3020 if k == 0:
3021 f.write(text)
3022 f.write(" |")
3023 # Add the weight of this variation unit after the pipe by comparing its analysis categories to their weights:
3024 weight = 1
3025 vu = self.variation_units[j]
3026 if len(vu.analysis_categories) > 0:
3027 weight = int(
3028 sum(
3029 [
3030 self.weights_by_id[ana] if ana in self.weights_by_id else 1
3031 for ana in vu.analysis_categories
3032 ]
3033 )
3034 / len(vu.analysis_categories)
3035 )
3036 f.write("*%d" % weight)
3037 # Every subsequent reading should be preceded by a space:
3038 elif k > 0:
3039 f.write(" %s" % text)
3040 k += 1
3041 f.write(" ]\n")
3042 # In a second pass, print the indices and witnesses for all readings enclosed in angle brackets:
3043 k = 0
3044 f.write("\t< ")
3045 while True:
3046 indices = tuple([j, k])
3047 if indices not in reading_wits_by_indices:
3048 break
3050 rdg_symbol = symbols[k] # get the one-character alphanumeric code for this state
3051 wits = " ".join(reading_wits_by_indices[indices])
3052 # Open the variant reading support block with an angle bracket:
3053 if k == 0:
3054 f.write("%s %s" % (rdg_symbol, wits))
3055 # Open all subsequent variant reading support blocks with pipes on the next line:
3056 else:
3057 f.write("\n\t| %s %s" % (rdg_symbol, wits))
3058 k += 1
3059 f.write(" >\n")
3060 pbar.update(1)
3061 # In a fourth pass, write to the chron file:
3062 max_id_length = max(
3063 [len(slugify(wit.id, lowercase=False, allow_unicode=True, separator='_')) for wit in self.witnesses]
3064 )
3065 max_date_length = 0
3066 for wit in self.witnesses:
3067 if wit.date_range[0] is not None:
3068 max_date_length = max(max_date_length, len(str(wit.date_range[0])))
3069 if wit.date_range[1] is not None:
3070 max_date_length = max(max_date_length, len(str(wit.date_range[1])))
3071 # Attempt to get the minimum and maximum dates for witnesses; if we can't do this, then don't write a chron file:
3072 min_date = None
3073 max_date = None
3074 try:
3075 min_date = min([wit.date_range[0] for wit in self.witnesses if wit.date_range[0] is not None])
3076 max_date = max([wit.date_range[1] for wit in self.witnesses if wit.date_range[1] is not None])
3077 except Exception as e:
3078 print("WARNING: no witnesses have date ranges; no chron file will be written!")
3079 return
3080 with open(chron_file_addr, "w", encoding="utf-8") as f:
3081 for wit in self.witnesses:
3082 wit_label = slugify(wit.id, lowercase=False, allow_unicode=True, separator='_')
3083 f.write(wit_label)
3084 f.write(" " * (max_id_length - len(wit.id) + 1))
3085 # If either the lower bound on this witness's date is empty, then use the min and max dates over all witnesses as defaults:
3086 date_range = wit.date_range
3087 if date_range[0] is None:
3088 date_range = tuple([min_date, date_range[1]])
3089 # Then write the date range minimum, average, and maximum to the chron file:
3090 low_date = str(date_range[0])
3091 f.write(" " * (max_date_length - len(low_date) + 2))
3092 f.write(low_date)
3093 avg_date = str(int(((date_range[0] + date_range[1]) / 2)))
3094 f.write(" " * (max_date_length - len(str(avg_date)) + 2))
3095 f.write(avg_date)
3096 high_date = str(date_range[1])
3097 f.write(" " * (max_date_length - len(high_date) + 2))
3098 f.write(high_date)
3099 f.write("\n")
3100 return
3102 def to_file(
3103 self,
3104 file_addr: Union[Path, str],
3105 format: Format = None,
3106 drop_constant: bool = False,
3107 split_missing: SplitMissingType = None,
3108 char_state_labels: bool = True,
3109 frequency: bool = False,
3110 ambiguous_as_missing: bool = False,
3111 proportion: bool = False,
3112 calibrate_dates: bool = False,
3113 mrbayes: bool = False,
3114 clock_model: ClockModel = ClockModel.strict,
3115 ancestral_logger: AncestralLogger = AncestralLogger.state,
3116 table_type: TableType = TableType.matrix,
3117 transform_matrix: TransformMatrixType = None,
3118 show_ext: bool = False,
3119 seed: int = None,
3120 ):
3121 """Writes this Collation to the file with the given address.
3123 Args:
3124 file_addr (Union[Path, str]): The path to the output file.
3125 format (Format, optional): The desired output format.
3126 If None then it is infered from the file suffix.
3127 Defaults to None.
3128 drop_constant (bool, optional): An optional flag indicating whether to ignore variation units with one substantive reading.
3129 Default value is False.
3130 split_missing (SplitMissingType, optional): An option indicating whether or not to treat missing characters/variation units as having a contribution of 1 split over all states/readings.
3131 If not specified, then missing data is ignored (i.e., all states are 0).
3132 If "uniform", then the contribution of 1 is divided evenly over all substantive readings.
3133 If "proportional", then the contribution of 1 is divided between the readings in proportion to their support among the witnesses that are not missing.
3134 Only applicable for tabular outputs of type "matrix", "idf", "mean-idf", "mi", and "mean-mi".
3135 char_state_labels (bool, optional): An optional flag indicating whether to print
3136 the CharStateLabels block in NEXUS output.
3137 Default value is True.
3138 frequency (bool, optional): An optional flag indicating whether to use the StatesFormat=Frequency setting
3139 instead of the StatesFormat=StatesPresent setting
3140 (and thus represent all states with frequency vectors rather than symbols)
3141 in NEXUS output.
3142 Note that this setting is necessary to make use of certainty degrees assigned to multiple ambiguous states in the collation.
3143 Default value is False.
3144 ambiguous_as_missing (bool, optional): An optional flag indicating whether to treat all ambiguous states as missing data.
3145 If this flag is set, then only base symbols will be generated for the NEXUS file.
3146 It is only applied if the frequency option is False.
3147 Default value is False.
3148 proportion (bool, optional): An optional flag indicating whether to populate a distance matrix's cells
3149 with a proportion of disagreements to variation units where both witnesses are extant.
3150 It is only applied if the table_type option is "distance".
3151 Default value is False.
3152 calibrate_dates (bool, optional): An optional flag indicating whether to add an Assumptions block that specifies date distributions for witnesses
3153 in NEXUS output.
3154 This option is intended for inputs to BEAST 2.
3155 Default value is False.
3156 mrbayes (bool, optional): An optional flag indicating whether to add a MrBayes block that specifies model settings and age calibrations for witnesses
3157 in NEXUS output.
3158 This option is intended for inputs to MrBayes.
3159 Default value is False.
3160 clock_model (ClockModel, optional): A ClockModel option indicating which type of clock model to use.
3161 This option is intended for inputs to MrBayes and BEAST 2.
3162 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.
3163 Default value is "strict".
3164 ancestral_logger (AncestralLogger, optional): An AncestralLogger option indicating which class of logger (if any) to use for ancestral states.
3165 This option is intended for inputs to BEAST 2.
3166 table_type (TableType, optional): A TableType option indicating which type of tabular output to generate.
3167 Only applicable for tabular outputs and PHYLIP outputs.
3168 If the output is a PHYLIP file, then the type of tabular output must be "distance" or "similarity"; otherwise, it will be ignored.
3169 Default value is "matrix".
3170 transform_matrix (TransformMatrixType, optional): A TransformMatrixType option indicating how the columns of a witness-to-witness matrix output should be transformed.
3171 Only applicable for tabular outputs in which the rows and columns correspond to the witnesses in the collation.
3172 show_ext (bool, optional): An optional flag indicating whether each cell in the matrix
3173 should include the number of variation units where both witnesses are extant after the number of their disagreements/agreements.
3174 Only applicable for tabular output formats of type "distance" or "similarity".
3175 Default value is False.
3176 seed (optional, int): A seed for random number generation (for setting initial values of unspecified transcriptional rates in BEAST 2 XML output).
3177 """
3178 file_addr = Path(file_addr)
3179 format = format or Format.infer(
3180 file_addr.suffix
3181 ) # an exception will be raised here if the format or suffix is invalid
3183 if format == Format.NEXUS:
3184 return self.to_nexus(
3185 file_addr,
3186 drop_constant=drop_constant,
3187 char_state_labels=char_state_labels,
3188 frequency=frequency,
3189 ambiguous_as_missing=ambiguous_as_missing,
3190 calibrate_dates=calibrate_dates,
3191 mrbayes=mrbayes,
3192 clock_model=clock_model,
3193 )
3195 if format == format.HENNIG86:
3196 return self.to_hennig86(file_addr, drop_constant=drop_constant)
3198 if format == format.PHYLIP:
3199 if table_type in [TableType.distance, TableType.similarity]:
3200 return self.to_phylip_matrix(
3201 file_addr,
3202 drop_constant=drop_constant,
3203 proportion=proportion,
3204 table_type=table_type,
3205 show_ext=show_ext,
3206 )
3207 return self.to_phylip(file_addr, drop_constant=drop_constant)
3209 if format == format.FASTA:
3210 return self.to_fasta(file_addr, drop_constant=drop_constant)
3212 if format == format.BEAST:
3213 return self.to_beast(
3214 file_addr,
3215 drop_constant=drop_constant,
3216 clock_model=clock_model,
3217 ancestral_logger=ancestral_logger,
3218 seed=seed,
3219 )
3221 if format == Format.CSV:
3222 return self.to_csv(
3223 file_addr,
3224 drop_constant=drop_constant,
3225 ambiguous_as_missing=ambiguous_as_missing,
3226 proportion=proportion,
3227 table_type=table_type,
3228 split_missing=split_missing,
3229 transform_matrix=transform_matrix,
3230 show_ext=show_ext,
3231 )
3233 if format == Format.TSV:
3234 return self.to_csv(
3235 file_addr,
3236 drop_constant=drop_constant,
3237 ambiguous_as_missing=ambiguous_as_missing,
3238 proportion=proportion,
3239 table_type=table_type,
3240 split_missing=split_missing,
3241 transform_matrix=transform_matrix,
3242 show_ext=show_ext,
3243 sep="\t",
3244 )
3246 if format == Format.EXCEL:
3247 return self.to_excel(
3248 file_addr,
3249 drop_constant=drop_constant,
3250 ambiguous_as_missing=ambiguous_as_missing,
3251 proportion=proportion,
3252 table_type=table_type,
3253 split_missing=split_missing,
3254 transform_matrix=transform_matrix,
3255 show_ext=show_ext,
3256 )
3258 if format == Format.STEMMA:
3259 return self.to_stemma(file_addr)