Coverage for teiphy/collation.py: 100.00%
1657 statements
« prev ^ index » next coverage.py v7.9.2, created at 2026-07-20 02:19 +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 the list of substantive readings only contains one entry, then this variation unit is not informative;
595 # return an empty dictionary and add nothing to the list of substantive reading labels:
596 if self.verbose:
597 print(
598 "Variation unit %s has %d substantive readings."
599 % (vu.id, len(self.substantive_readings_by_variation_unit_id[vu.id]))
600 )
601 readings_by_witness_for_unit = {}
602 # Initialize the output dictionary with empty sets for all base witnesses:
603 for wit in self.witnesses:
604 readings_by_witness_for_unit[wit.id] = [0] * len(self.substantive_readings_by_variation_unit_id[vu.id])
605 # In a second pass, assign each base witness a set containing the readings it supports in this unit:
606 for rdg in vu.readings:
607 # Initialize the dictionary indicating support for this reading (or its disambiguations):
608 rdg_support = [0] * len(self.substantive_readings_by_variation_unit_id[vu.id])
609 # 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:
610 if rdg.type in self.missing_reading_types:
611 continue
612 # Otherwise, if this reading is trivial, then it will contain an entry for the index of its parent substantive reading:
613 elif rdg.type in self.trivial_reading_types:
614 rdg_support[reading_id_to_index[rdg.id]] = 1
615 # Otherwise, if this reading has one or more nonzero certainty degrees,
616 # then set the entries for these readings to their degrees:
617 elif sum(rdg.certainties.values()) > 0:
618 for t in rdg.certainties:
619 # Skip any reading whose ID is unrecognized in this unit:
620 if t in reading_id_to_index:
621 rdg_support[reading_id_to_index[t]] = rdg.certainties[t]
622 # Otherwise, if this reading has one or more targets (i.e., if it is an ambiguous reading),
623 # then set the entries for each of its targets to 1:
624 elif len(rdg.targets) > 0:
625 for t in rdg.targets:
626 # Skip any reading whose ID is unrecognized in this unit:
627 if t in reading_id_to_index:
628 rdg_support[reading_id_to_index[t]] = 1
629 # Otherwise, this reading is itself substantive; set the entry for the index of this reading to 1:
630 else:
631 rdg_support[reading_id_to_index[rdg.id]] = 1
632 # Proceed for each witness siglum in the support for this reading:
633 for wit in rdg.wits:
634 # Is this siglum a base siglum?
635 base_wit = self.get_base_wit(wit)
636 if base_wit not in self.witness_index_by_id:
637 # 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;
638 # report this if we're in verbose mode and move on:
639 if self.verbose:
640 print(
641 "Skipping unknown witness siglum %s (base siglum %s) in variation unit %s, reading %s..."
642 % (wit, base_wit, vu.id, rdg.id)
643 )
644 continue
645 # If we've found a base siglum, then add this reading's contribution to the base witness's reading set for this unit;
646 # normally the existing set will be empty, but if we reduce two suffixed sigla to the same base witness,
647 # then that witness may attest to multiple readings in the same unit:
648 readings_by_witness_for_unit[base_wit] = [
649 (min(readings_by_witness_for_unit[base_wit][i] + rdg_support[i], 1))
650 for i in range(len(rdg_support))
651 ]
652 return readings_by_witness_for_unit
654 def parse_readings_by_witness(self):
655 """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."""
656 if self.verbose:
657 print("Populating internal dictionary of witness readings...")
658 t0 = time.time()
659 # Initialize the data structures to be populated here:
660 self.readings_by_witness = {}
661 self.variation_unit_ids = []
662 for wit in self.witnesses:
663 self.readings_by_witness[wit.id] = []
664 # Populate them for each variation unit:
665 for vu in self.variation_units:
666 readings_by_witness_for_unit = self.get_readings_by_witness_for_unit(vu)
667 if len(readings_by_witness_for_unit) > 0:
668 self.variation_unit_ids.append(vu.id)
669 for wit in readings_by_witness_for_unit:
670 self.readings_by_witness[wit].append(readings_by_witness_for_unit[wit])
671 # Optionally, fill the lacunae of the correctors:
672 if self.fill_corrector_lacunae:
673 filled_readings = []
674 for i, wit in enumerate(self.witnesses):
675 # If this is the first witness, then it shouldn't be a corrector (since there is no previous witness against which to compare it):
676 if i == 0:
677 filled_readings = list(self.readings_by_witness[wit.id])
678 continue
679 # Otherwise, if this witness is not a corrector, then skip it:
680 if wit.type != "corrector":
681 filled_readings = list(self.readings_by_witness[wit.id])
682 continue
683 # Otherwise, add this corrector's extant readings to the filled readings list:
684 for j in range(len(self.readings_by_witness[wit.id])):
685 if sum(self.readings_by_witness[wit.id][j]) != 0:
686 filled_readings[j] = list(self.readings_by_witness[wit.id][j])
687 # If a threshold of extant readings is specified, then check if this corrector meets it, and skip it if not:
688 if self.fill_correctors_threshold is not None:
689 # If there is a threshold, then first check the proportion of variation units at which this witness is not lacunose:
690 proportion_extant = sum(
691 [
692 1
693 for j in range(len(self.readings_by_witness[wit.id]))
694 if sum(self.readings_by_witness[wit.id][j]) != 0
695 ]
696 ) / len(self.readings_by_witness[wit.id])
697 # If this corrector does not exceed the threshold, then don't fill it,
698 # but do update the running set of readings for the previous corrector and first hand:
699 if proportion_extant < self.fill_correctors_threshold:
700 continue
701 # Otherwise, fill every lacuna in this corrector based on the filled readings list:
702 for j in range(len(self.readings_by_witness[wit.id])):
703 if sum(self.readings_by_witness[wit.id][j]) == 0:
704 self.readings_by_witness[wit.id][j] = list(filled_readings[j])
705 t1 = time.time()
706 if self.verbose:
707 print(
708 "Populated dictionary for %d witnesses over %d substantive variation units in %0.4fs."
709 % (len(self.witnesses), len(self.variation_unit_ids), t1 - t0)
710 )
711 return
713 def filter_fragmentary_witnesses(self, xml):
714 """Filters the original witness list and readings by witness dictionary to exclude witnesses whose proportions of extant passages fall below the fragmentary readings threshold."""
715 if self.verbose:
716 print(
717 "Filtering fragmentary witnesses (extant in < %f of all variation units) out of internal witness list and dictionary of witness readings..."
718 % self.fragmentary_threshold
719 )
720 t0 = time.time()
721 fragmentary_witness_set = set()
722 # Proceed for each witness in order:
723 for wit in self.witnesses:
724 wit_id = wit.id
725 # We count the number of variation units at which this witness has an extant (i.e., non-missing) reading:
726 extant_reading_count = 0
727 total_reading_count = len(self.readings_by_witness[wit.id])
728 # Proceed through all reading support lists:
729 for rdg_support in self.readings_by_witness[wit_id]:
730 # If the current reading support list is not all zeroes, then increment this witness's count of extant readings:
731 if sum(rdg_support) != 0:
732 extant_reading_count += 1
733 # If the proportion of extant readings falls below the threshold, then add this witness to the list of fragmentary witnesses:
734 if extant_reading_count / total_reading_count < self.fragmentary_threshold:
735 fragmentary_witness_set.add(wit_id)
736 # Then filter the witness list to exclude the fragmentary witnesses:
737 filtered_witnesses = [wit for wit in self.witnesses if wit.id not in fragmentary_witness_set]
738 self.witnesses = filtered_witnesses
739 # Then remove the entries for the fragmentary witnesses from the witnesses-to-readings dictionary:
740 for wit_id in fragmentary_witness_set:
741 del self.readings_by_witness[wit_id]
742 t1 = time.time()
743 if self.verbose:
744 print(
745 "Filtered out %d fragmentary witness(es) (%s) in %0.4fs."
746 % (len(fragmentary_witness_set), str(list(fragmentary_witness_set)), t1 - t0)
747 )
748 return
750 def get_nexus_symbols(self):
751 """Returns a list of one-character symbols needed to represent the states of all substantive readings in NEXUS.
753 The number of symbols equals the maximum number of substantive readings at any variation unit.
755 Returns:
756 A list of individual characters representing states in readings.
757 """
758 # NOTE: IQTREE does not appear to support symbols outside of 0-9 and a-z, and its base symbols must be case-insensitive.
759 # The official version of MrBayes is likewise limited to 32 symbols.
760 # But PAUP* allows up to 64 symbols, and Andrew Edmondson's fork of MrBayes does, as well.
761 # So this method will support symbols from 0-9, a-z, and A-Z (for a total of 62 states)
762 possible_symbols = list(string.digits) + list(string.ascii_lowercase) + list(string.ascii_uppercase)
763 # The number of symbols needed is equal to the length of the longest substantive reading vector:
764 nsymbols = 0
765 # If there are no witnesses, then no symbols are needed at all:
766 if len(self.witnesses) == 0:
767 return []
768 wit_id = self.witnesses[0].id
769 for rdg_support in self.readings_by_witness[wit_id]:
770 nsymbols = max(nsymbols, len(rdg_support))
771 nexus_symbols = possible_symbols[:nsymbols]
772 return nexus_symbols
774 def to_nexus(
775 self,
776 file_addr: Union[Path, str],
777 drop_constant: bool = False,
778 char_state_labels: bool = True,
779 frequency: bool = False,
780 ambiguous_as_missing: bool = False,
781 calibrate_dates: bool = False,
782 mrbayes: bool = False,
783 clock_model: ClockModel = ClockModel.strict,
784 ):
785 """Writes this Collation to a NEXUS file with the given address.
787 Args:
788 file_addr: A string representing the path to an output NEXUS file; the file type should be .nex, .nexus, or .nxs.
789 drop_constant: An optional flag indicating whether to ignore variation units with one substantive reading.
790 char_state_labels: An optional flag indicating whether or not to include the CharStateLabels block.
791 frequency: An optional flag indicating whether to use the StatesFormat=Frequency setting
792 instead of the StatesFormat=StatesPresent setting
793 (and thus represent all states with frequency vectors rather than symbols).
794 Note that this setting is necessary to make use of certainty degrees assigned to multiple ambiguous states in the collation.
795 ambiguous_as_missing: An optional flag indicating whether to treat all ambiguous states as missing data.
796 If this flag is set, then only base symbols will be generated for the NEXUS file.
797 It is only applied if the frequency option is False.
798 calibrate_dates: An optional flag indicating whether to add an Assumptions block that specifies date distributions for witnesses.
799 This option is intended for inputs to BEAST 2.
800 mrbayes: An optional flag indicating whether to add a MrBayes block that specifies model settings and age calibrations for witnesses.
801 This option is intended for inputs to MrBayes.
802 clock_model: A ClockModel option indicating which type of clock model to use.
803 This option is intended for inputs to MrBayes and BEAST 2.
804 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.
805 """
806 # Populate a list of sites that will correspond to columns of the sequence alignment:
807 substantive_variation_unit_ids = self.variation_unit_ids
808 if drop_constant:
809 substantive_variation_unit_ids = [
810 vu_id
811 for vu_id in self.variation_unit_ids
812 if len(self.substantive_readings_by_variation_unit_id[vu_id]) > 1
813 ]
814 substantive_variation_unit_ids_set = set(substantive_variation_unit_ids)
815 substantive_variation_unit_reading_tuples_set = set(self.substantive_variation_unit_reading_tuples)
816 # Start by calculating the values we will be using here:
817 ntax = len(self.witnesses)
818 nchar = len(substantive_variation_unit_ids)
819 taxlabels = [slugify(wit.id, lowercase=False, separator='_') for wit in self.witnesses]
820 max_taxlabel_length = max(
821 [len(taxlabel) for taxlabel in taxlabels]
822 ) # keep track of the longest taxon label for tabular alignment purposes
823 charlabels = [slugify(vu_id, lowercase=False, separator='_') for vu_id in substantive_variation_unit_ids]
824 missing_symbol = '?'
825 symbols = self.get_nexus_symbols()
826 # Generate all parent folders for this file that don't already exist:
827 Path(file_addr).parent.mkdir(parents=True, exist_ok=True)
828 # Then write the file:
829 pbar = tqdm()
830 with open(file_addr, "w", encoding="utf-8") as f:
831 # Start with the NEXUS header:
832 f.write("#NEXUS\n\n")
833 # Then begin the data block:
834 f.write("Begin DATA;\n")
835 # Write the collation matrix dimensions:
836 f.write("\tDimensions ntax=%d nchar=%d;\n" % (ntax, nchar))
837 # Write the format subblock:
838 f.write("\tFormat\n")
839 f.write("\t\tDataType=Standard\n")
840 f.write("\t\tMissing=%s\n" % missing_symbol)
841 if frequency:
842 f.write("\t\tStatesFormat=Frequency\n")
843 f.write("\t\tSymbols=\"%s\";\n" % (" ".join(symbols)))
844 # If the char_state_labels is set, then write the labels for character-state labels, with each on its own line:
845 if char_state_labels:
846 f.write("\tCharStateLabels")
847 vu_ind = 1
848 for vu in self.variation_units:
849 if vu.id not in substantive_variation_unit_ids_set:
850 continue
851 if vu_ind == 1:
852 f.write("\n\t\t%d %s /" % (vu_ind, slugify(vu.id, lowercase=False, separator='_')))
853 else:
854 f.write(",\n\t\t%d %s /" % (vu_ind, slugify(vu.id, lowercase=False, separator='_')))
855 rdg_ind = 0
856 for rdg in vu.readings:
857 key = tuple([vu.id, rdg.id])
858 if key not in substantive_variation_unit_reading_tuples_set:
859 continue
860 ascii_rdg_text = slugify(
861 rdg.text, lowercase=False, separator='_', replacements=[['η', 'h'], ['ω', 'w']]
862 )
863 if ascii_rdg_text == "":
864 ascii_rdg_text = "om."
865 f.write(" %s" % ascii_rdg_text)
866 rdg_ind += 1
867 if rdg_ind > 0:
868 vu_ind += 1
869 f.write(";\n")
870 # Write the matrix subblock:
871 f.write("\tMatrix")
872 with tqdm(total=len(self.witnesses)) as pbar:
873 for i, wit in enumerate(self.witnesses):
874 taxlabel = taxlabels[i]
875 if frequency:
876 sequence = "\n\t\t" + taxlabel
877 for j, vu_id in enumerate(self.variation_unit_ids):
878 if vu_id not in substantive_variation_unit_ids_set:
879 continue
880 rdg_support = self.readings_by_witness[wit.id][j]
881 sequence += "\n\t\t\t"
882 # If this reading is lacunose in this witness, then use the missing character:
883 if sum(rdg_support) == 0:
884 sequence += missing_symbol
885 continue
886 # Otherwise, print out its frequencies for different readings in parentheses:
887 sequence += "("
888 for k, w in enumerate(rdg_support):
889 sequence += "%s:%0.4f" % (symbols[k], w)
890 if k < len(rdg_support) - 1:
891 sequence += " "
892 sequence += ")"
893 else:
894 sequence = "\n\t\t" + taxlabel
895 # Add enough space after this label ensure that all sequences are nicely aligned:
896 sequence += " " * (max_taxlabel_length - len(taxlabel) + 1)
897 for j, vu_id in enumerate(self.variation_unit_ids):
898 if vu_id not in substantive_variation_unit_ids_set:
899 continue
900 rdg_support = self.readings_by_witness[wit.id][j]
901 # If this reading is lacunose in this witness, then use the missing character:
902 if sum(rdg_support) == 0:
903 sequence += missing_symbol
904 continue
905 rdg_inds = [
906 k for k, w in enumerate(rdg_support) if w > 0
907 ] # the index list consists of the indices of all readings with any degree of certainty assigned to them
908 # For singleton readings, just print the symbol:
909 if len(rdg_inds) == 1:
910 sequence += symbols[rdg_inds[0]]
911 continue
912 # For multiple readings, print the corresponding readings in braces or the missing symbol depending on input settings:
913 if ambiguous_as_missing:
914 sequence += missing_symbol
915 else:
916 sequence += "{%s}" % "".join([str(rdg_ind) for rdg_ind in rdg_inds])
917 f.write("%s" % (sequence))
918 pbar.update(1)
919 f.write(";\n")
920 # End the data block:
921 f.write("End;")
922 # If calibrate_dates is set, then add the assumptions block:
923 if calibrate_dates:
924 f.write("\n\n")
925 f.write("Begin ASSUMPTIONS;\n")
926 # Set the scale to years:
927 f.write("\tOPTIONS SCALE = years;\n\n")
928 # Then calibrate the witness ages:
929 calibrate_strings = []
930 for i, wit in enumerate(self.witnesses):
931 taxlabel = taxlabels[i]
932 date_range = wit.date_range
933 if date_range[0] is not None:
934 # If there is a lower bound on the witness's date, then use either a fixed or uniform distribution,
935 # depending on whether the upper and lower bounds match:
936 min_age = datetime.now().year - date_range[1]
937 max_age = datetime.now().year - date_range[0]
938 if min_age == max_age:
939 calibrate_string = "\tCALIBRATE %s = fixed(%d)" % (taxlabel, min_age)
940 calibrate_strings.append(calibrate_string)
941 else:
942 calibrate_string = "\tCALIBRATE %s = uniform(%d,%d)" % (taxlabel, min_age, max_age)
943 calibrate_strings.append(calibrate_string)
944 else:
945 # If there is no lower bound on the witness's date, then use an offset log-normal distribution:
946 min_age = datetime.now().year - date_range[1]
947 calibrate_string = "\tCALIBRATE %s = offsetlognormal(%d,0.0,1.0)" % (taxlabel, min_age)
948 calibrate_strings.append(calibrate_string)
949 # Then print the calibrate strings, separated by commas and line breaks and terminated by a semicolon:
950 f.write("%s;\n\n" % ",\n".join(calibrate_strings))
951 # End the assumptions block:
952 f.write("End;")
953 # If mrbayes is set, then add the mrbayes block:
954 if mrbayes:
955 f.write("\n\n")
956 f.write("Begin MRBAYES;\n")
957 # Turn on the autoclose feature by default:
958 f.write("\tset autoclose=yes;\n")
959 # Set the branch lengths to be governed by a birth-death clock model, and set up the parameters for this model:
960 f.write("\n")
961 f.write("\tprset brlenspr = clock:birthdeath;\n")
962 f.write("\tprset speciationpr = uniform(0.0,10.0);\n")
963 f.write("\tprset extinctionpr = beta(2.0,4.0);\n")
964 f.write("\tprset sampleprob = 0.01;\n")
965 # Use the specified clock model:
966 f.write("\n")
967 if clock_model == clock_model.uncorrelated:
968 f.write("\tprset clockvarpr=igr;\n")
969 f.write("\tprset clockratepr=lognormal(0.0,1.0);\n")
970 f.write("\tprset igrvarpr=exponential(1.0);\n")
971 else:
972 f.write("\tprset clockvarpr=strict;\n")
973 f.write("\tprset clockratepr=lognormal(0.0,1.0);\n")
974 # Set the priors on the tree age depending on the date range for the origin of the collated work:
975 f.write("\n")
976 if self.origin_date_range[0] is not None:
977 min_tree_age = (
978 datetime.now().year - self.origin_date_range[1]
979 if self.origin_date_range[1] is not None
980 else 0.0
981 )
982 max_tree_age = datetime.now().year - self.origin_date_range[0]
983 f.write("\tprset treeagepr = uniform(%d,%d);\n" % (min_tree_age, max_tree_age))
984 else:
985 min_tree_age = (
986 datetime.now().year - self.origin_date_range[1]
987 if self.origin_date_range[1] is not None
988 else 0.0
989 )
990 f.write("\tprset treeagepr = offsetgamma(%d,1.0,1.0);\n" % (min_tree_age))
991 # Then calibrate the witness ages:
992 f.write("\n")
993 f.write("\tprset nodeagepr = calibrated;\n")
994 for i, wit in enumerate(self.witnesses):
995 taxlabel = taxlabels[i]
996 date_range = wit.date_range
997 if date_range[0] is not None:
998 # If there is a lower bound on the witness's date, then use either a fixed or uniform distribution,
999 # depending on whether the upper and lower bounds match:
1000 min_age = datetime.now().year - date_range[1]
1001 max_age = datetime.now().year - date_range[0]
1002 if min_age == max_age:
1003 f.write("\tcalibrate %s = fixed(%d);\n" % (taxlabel, min_age))
1004 else:
1005 f.write("\tcalibrate %s = uniform(%d,%d);\n" % (taxlabel, min_age, max_age))
1006 else:
1007 # If there is no lower bound on the witness's date, then use an offset gamma distribution:
1008 min_age = datetime.now().year - date_range[1]
1009 f.write("\tcalibrate %s = offsetgamma(%d,1.0,1.0);\n" % (taxlabel, min_age))
1010 f.write("\n")
1011 # Add default settings for MCMC estimation of posterior distribution:
1012 f.write("\tmcmcp ngen=100000;\n")
1013 # Write the command to run MrBayes:
1014 f.write("\tmcmc;\n")
1015 # End the assumptions block:
1016 f.write("End;")
1017 return
1019 def get_hennig86_symbols(self):
1020 """Returns a list of one-character symbols needed to represent the states of all substantive readings in Hennig86 format.
1022 The number of symbols equals the maximum number of substantive readings at any variation unit.
1024 Returns:
1025 A list of individual characters representing states in readings.
1026 """
1027 possible_symbols = (
1028 list(string.digits) + list(string.ascii_uppercase)[:22]
1029 ) # NOTE: the maximum number of symbols allowed in Hennig86 format is 32
1030 # The number of symbols needed is equal to the length of the longest substantive reading vector:
1031 nsymbols = 0
1032 # If there are no witnesses, then no symbols are needed at all:
1033 if len(self.witnesses) == 0:
1034 return []
1035 wit_id = self.witnesses[0].id
1036 for rdg_support in self.readings_by_witness[wit_id]:
1037 nsymbols = max(nsymbols, len(rdg_support))
1038 hennig86_symbols = possible_symbols[:nsymbols]
1039 return hennig86_symbols
1041 def to_hennig86(self, file_addr: Union[Path, str], drop_constant: bool = False):
1042 """Writes this Collation to a file in Hennig86 format with the given address.
1043 Note that because Hennig86 format does not support NEXUS-style ambiguities, such ambiguities will be treated as missing data.
1045 Args:
1046 file_addr: A string representing the path to an output file.
1047 drop_constant (bool, optional): An optional flag indicating whether to ignore variation units with one substantive reading.
1048 """
1049 # Populate a list of sites that will correspond to columns of the sequence alignment:
1050 substantive_variation_unit_ids = self.variation_unit_ids
1051 if drop_constant:
1052 substantive_variation_unit_ids = [
1053 vu_id
1054 for vu_id in self.variation_unit_ids
1055 if len(self.substantive_readings_by_variation_unit_id[vu_id]) > 1
1056 ]
1057 substantive_variation_unit_ids_set = set(substantive_variation_unit_ids)
1058 substantive_variation_unit_reading_tuples_set = set(self.substantive_variation_unit_reading_tuples)
1059 # Start by calculating the values we will be using here:
1060 ntax = len(self.witnesses)
1061 nchar = len(substantive_variation_unit_ids)
1062 taxlabels = []
1063 for wit in self.witnesses:
1064 taxlabel = wit.id
1065 # 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:
1066 if taxlabel[0] not in string.ascii_letters:
1067 taxlabel = "WIT_" + taxlabel
1068 # Then replace any disallowed characters in the string with an underscore:
1069 taxlabel = slugify(taxlabel, lowercase=False, separator='_')
1070 taxlabels.append(taxlabel)
1071 max_taxlabel_length = max(
1072 [len(taxlabel) for taxlabel in taxlabels]
1073 ) # keep track of the longest taxon label for tabular alignment purposes
1074 missing_symbol = '?'
1075 symbols = self.get_hennig86_symbols()
1076 # Generate all parent folders for this file that don't already exist:
1077 Path(file_addr).parent.mkdir(parents=True, exist_ok=True)
1078 with open(file_addr, "w", encoding="ascii") as f:
1079 # Start with the nstates header:
1080 f.write("nstates %d;\n" % len(symbols))
1081 # Then begin the xread block:
1082 f.write("xread\n")
1083 # Write the dimensions:
1084 f.write("%d %d\n" % (nchar, ntax))
1085 # Now write the matrix:
1086 with tqdm(total=len(self.witnesses)) as pbar:
1087 for i, wit in enumerate(self.witnesses):
1088 taxlabel = taxlabels[i]
1089 # Add enough space after this label ensure that all sequences are nicely aligned:
1090 sequence = taxlabel + (" " * (max_taxlabel_length - len(taxlabel) + 1))
1091 for j, vu_id in enumerate(self.variation_unit_ids):
1092 if vu_id not in substantive_variation_unit_ids_set:
1093 continue
1094 rdg_support = self.readings_by_witness[wit.id][j]
1095 # If this reading is lacunose in this witness, then use the missing character:
1096 if sum(rdg_support) == 0:
1097 sequence += missing_symbol
1098 continue
1099 rdg_inds = [
1100 k for k, w in enumerate(rdg_support) if w > 0
1101 ] # the index list consists of the indices of all readings with any degree of certainty assigned to them
1102 # For singleton readings, just print the symbol:
1103 if len(rdg_inds) == 1:
1104 sequence += symbols[rdg_inds[0]]
1105 continue
1106 # For multiple readings, print the missing symbol:
1107 sequence += missing_symbol
1108 f.write("%s\n" % (sequence))
1109 pbar.update(1)
1110 f.write(";")
1111 return
1113 def get_phylip_symbols(self):
1114 """Returns a list of one-character symbols needed to represent the states of all substantive readings in PHYLIP format.
1116 The number of symbols equals the maximum number of substantive readings at any variation unit.
1118 Returns:
1119 A list of individual characters representing states in readings.
1120 """
1121 possible_symbols = (
1122 list(string.digits) + list(string.ascii_lowercase)[:22]
1123 ) # NOTE: for RAxML, multistate characters with an alphabet sizes up to 32 are supported
1124 # The number of symbols needed is equal to the length of the longest substantive reading vector:
1125 nsymbols = 0
1126 # If there are no witnesses, then no symbols are needed at all:
1127 if len(self.witnesses) == 0:
1128 return []
1129 wit_id = self.witnesses[0].id
1130 for rdg_support in self.readings_by_witness[wit_id]:
1131 nsymbols = max(nsymbols, len(rdg_support))
1132 phylip_symbols = possible_symbols[:nsymbols]
1133 return phylip_symbols
1135 def to_phylip(self, file_addr: Union[Path, str], drop_constant: bool = False):
1136 """Writes this Collation to a file in PHYLIP format with the given address.
1137 Note that because PHYLIP format does not support NEXUS-style ambiguities, such ambiguities will be treated as missing data.
1139 Args:
1140 file_addr: A string representing the path to an output file.
1141 drop_constant (bool, optional): An optional flag indicating whether to ignore variation units with one substantive reading.
1142 """
1143 # Populate a list of sites that will correspond to columns of the sequence alignment:
1144 substantive_variation_unit_ids = self.variation_unit_ids
1145 if drop_constant:
1146 substantive_variation_unit_ids = [
1147 vu_id
1148 for vu_id in self.variation_unit_ids
1149 if len(self.substantive_readings_by_variation_unit_id[vu_id]) > 1
1150 ]
1151 substantive_variation_unit_ids_set = set(substantive_variation_unit_ids)
1152 substantive_variation_unit_reading_tuples_set = set(self.substantive_variation_unit_reading_tuples)
1153 # Start by calculating the values we will be using here:
1154 ntax = len(self.witnesses)
1155 nchar = len(substantive_variation_unit_ids)
1156 taxlabels = []
1157 for wit in self.witnesses:
1158 taxlabel = wit.id
1159 # Then replace any disallowed characters in the string with an underscore:
1160 taxlabel = slugify(taxlabel, lowercase=False, separator='_')
1161 taxlabels.append(taxlabel)
1162 max_taxlabel_length = max(
1163 [len(taxlabel) for taxlabel in taxlabels]
1164 ) # keep track of the longest taxon label for tabular alignment purposes
1165 missing_symbol = '?'
1166 symbols = self.get_phylip_symbols()
1167 # Generate all parent folders for this file that don't already exist:
1168 Path(file_addr).parent.mkdir(parents=True, exist_ok=True)
1169 with open(file_addr, "w", encoding="ascii") as f:
1170 # Write the dimensions:
1171 f.write("%d %d\n" % (ntax, nchar))
1172 # Now write the matrix:
1173 for i, wit in enumerate(self.witnesses):
1174 taxlabel = taxlabels[i]
1175 # Add enough space after this label ensure that all sequences are nicely aligned:
1176 sequence = taxlabel + (" " * (max_taxlabel_length - len(taxlabel))) + "\t"
1177 for j, vu_id in enumerate(self.variation_unit_ids):
1178 if vu_id not in substantive_variation_unit_ids_set:
1179 continue
1180 rdg_support = self.readings_by_witness[wit.id][j]
1181 # If this reading is lacunose in this witness, then use the missing character:
1182 if sum(rdg_support) == 0:
1183 sequence += missing_symbol
1184 continue
1185 rdg_inds = [
1186 k for k, w in enumerate(rdg_support) if w > 0
1187 ] # the index list consists of the indices of all readings with any degree of certainty assigned to them
1188 # For singleton readings, just print the symbol:
1189 if len(rdg_inds) == 1:
1190 sequence += symbols[rdg_inds[0]]
1191 continue
1192 # For multiple readings, print the missing symbol:
1193 sequence += missing_symbol
1194 f.write("%s\n" % (sequence))
1195 return
1197 def get_fasta_symbols(self):
1198 """Returns a list of one-character symbols needed to represent the states of all substantive readings in FASTA format.
1200 The number of symbols equals the maximum number of substantive readings at any variation unit.
1202 Returns:
1203 A list of individual characters representing states in readings.
1204 """
1205 possible_symbols = (
1206 list(string.digits) + list(string.ascii_lowercase)[:22]
1207 ) # NOTE: for RAxML, multistate characters with an alphabet sizes up to 32 are supported
1208 # The number of symbols needed is equal to the length of the longest substantive reading vector:
1209 nsymbols = 0
1210 # If there are no witnesses, then no symbols are needed at all:
1211 if len(self.witnesses) == 0:
1212 return []
1213 wit_id = self.witnesses[0].id
1214 for rdg_support in self.readings_by_witness[wit_id]:
1215 nsymbols = max(nsymbols, len(rdg_support))
1216 fasta_symbols = possible_symbols[:nsymbols]
1217 return fasta_symbols
1219 def to_fasta(self, file_addr: Union[Path, str], drop_constant: bool = False):
1220 """Writes this Collation to a file in FASTA format with the given address.
1221 Note that because FASTA format does not support NEXUS-style ambiguities, such ambiguities will be treated as missing data.
1223 Args:
1224 file_addr: A string representing the path to an output file.
1225 drop_constant (bool, optional): An optional flag indicating whether to ignore variation units with one substantive reading.
1226 """
1227 # Populate a list of sites that will correspond to columns of the sequence alignment:
1228 substantive_variation_unit_ids = self.variation_unit_ids
1229 if drop_constant:
1230 substantive_variation_unit_ids = [
1231 vu_id
1232 for vu_id in self.variation_unit_ids
1233 if len(self.substantive_readings_by_variation_unit_id[vu_id]) > 1
1234 ]
1235 substantive_variation_unit_ids_set = set(substantive_variation_unit_ids)
1236 substantive_variation_unit_reading_tuples_set = set(self.substantive_variation_unit_reading_tuples)
1237 # Start by calculating the values we will be using here:
1238 ntax = len(self.witnesses)
1239 nchar = len(substantive_variation_unit_ids)
1240 taxlabels = []
1241 for wit in self.witnesses:
1242 taxlabel = wit.id
1243 # Then replace any disallowed characters in the string with an underscore:
1244 taxlabel = slugify(taxlabel, lowercase=False, separator='_')
1245 taxlabels.append(taxlabel)
1246 max_taxlabel_length = max(
1247 [len(taxlabel) for taxlabel in taxlabels]
1248 ) # keep track of the longest taxon label for tabular alignment purposes
1249 missing_symbol = '?'
1250 symbols = self.get_fasta_symbols()
1251 # Generate all parent folders for this file that don't already exist:
1252 Path(file_addr).parent.mkdir(parents=True, exist_ok=True)
1253 with open(file_addr, "w", encoding="ascii") as f:
1254 # Now write the matrix:
1255 with tqdm(total=len(self.witnesses)) as pbar:
1256 for i, wit in enumerate(self.witnesses):
1257 taxlabel = taxlabels[i]
1258 # Add enough space after this label ensure that all sequences are nicely aligned:
1259 sequence = ">%s\n" % taxlabel
1260 for j, vu_id in enumerate(self.variation_unit_ids):
1261 if vu_id not in substantive_variation_unit_ids_set:
1262 continue
1263 rdg_support = self.readings_by_witness[wit.id][j]
1264 # If this reading is lacunose in this witness, then use the missing character:
1265 if sum(rdg_support) == 0:
1266 sequence += missing_symbol
1267 continue
1268 rdg_inds = [
1269 k for k, w in enumerate(rdg_support) if w > 0
1270 ] # the index list consists of the indices of all readings with any degree of certainty assigned to them
1271 # For singleton readings, just print the symbol:
1272 if len(rdg_inds) == 1:
1273 sequence += symbols[rdg_inds[0]]
1274 continue
1275 # For multiple readings, print the missing symbol:
1276 sequence += missing_symbol
1277 f.write("%s\n" % (sequence))
1278 pbar.update(1)
1279 return
1281 def get_beast_symbols(self):
1282 """Returns a list of one-character symbols needed to represent the states of all substantive readings in BEAST format.
1284 The number of symbols equals the maximum number of substantive readings at any variation unit.
1286 Returns:
1287 A list of individual characters representing states in readings.
1288 """
1289 possible_symbols = (
1290 list(string.digits) + list(string.ascii_lowercase) + list(string.ascii_uppercase)
1291 ) # 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
1292 # The number of symbols needed is equal to the length of the longest substantive reading vector:
1293 nsymbols = 0
1294 # If there are no witnesses, then no symbols are needed at all:
1295 if len(self.witnesses) == 0:
1296 return []
1297 wit_id = self.witnesses[0].id
1298 for rdg_support in self.readings_by_witness[wit_id]:
1299 nsymbols = max(nsymbols, len(rdg_support))
1300 beast_symbols = possible_symbols[:nsymbols]
1301 return beast_symbols
1303 def get_tip_date_range(self):
1304 """Gets the minimum and maximum dates attested among the witnesses.
1305 Also checks if the witness with the latest possible date has a fixed date
1306 (i.e, if the lower and upper bounds for its date are the same)
1307 and issues a warning if not, as this will cause unusual behavior in BEAST 2.
1309 Returns:
1310 A tuple containing the earliest and latest possible tip dates.
1311 """
1312 earliest_date = None
1313 earliest_wit = None
1314 latest_date = None
1315 latest_wit = None
1316 for wit in self.witnesses:
1317 wit_id = wit.id
1318 date_range = wit.date_range
1319 if date_range[0] is not None:
1320 if earliest_date is not None:
1321 earliest_wit = wit if date_range[0] < earliest_date else earliest_wit
1322 earliest_date = min(date_range[0], earliest_date)
1323 else:
1324 earliest_wit = wit
1325 earliest_date = date_range[0]
1326 if date_range[1] is not None:
1327 if latest_date is not None:
1328 latest_wit = (
1329 wit
1330 if (date_range[1] > latest_date or (date_range[0] == date_range[1] == latest_date))
1331 else latest_wit
1332 ) # 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
1333 latest_date = max(date_range[1], latest_date)
1334 else:
1335 latest_wit = wit
1336 latest_date = date_range[1]
1337 if latest_wit.date_range[0] is None or latest_wit.date_range[0] != latest_wit.date_range[1]:
1338 print(
1339 "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."
1340 % (latest_wit.id, latest_wit.id)
1341 )
1342 return (earliest_date, latest_date)
1344 def get_beast_origin_span(self, tip_date_range):
1345 """Returns a tuple containing the lower and upper bounds for the height of the origin of the Birth-Death Skyline model.
1346 The upper bound on the height of the tree is the difference between the latest tip date
1347 and the lower bound on the date of the original work, if both are defined;
1348 otherwise, it is left undefined.
1349 The lower bound on the height of the tree is the difference between the latest tip date
1350 and the upper bound on the date of the original work, if both are defined;
1351 otherwise, it is the difference between the earliest tip date and the latest, if both are defined.
1353 Args:
1354 tip_date_range: A tuple containing the earliest and latest possible tip dates.
1356 Returns:
1357 A tuple containing lower and upper bounds on the origin height for the Birth-Death Skyline model.
1358 """
1359 origin_span = [0, None]
1360 # 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
1361 # (note that if it had to be defined in terms of witness date lower bounds, then this would have happened already):
1362 if self.origin_date_range[1] is not None:
1363 origin_span[0] = tip_date_range[1] - self.origin_date_range[1]
1364 # 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:
1365 if self.origin_date_range[0] is not None:
1366 origin_span[1] = tip_date_range[1] - self.origin_date_range[0]
1367 return tuple(origin_span)
1369 def get_beast_date_map(self, taxlabels):
1370 """Returns a string representing witness-to-date mappings in BEAST format.
1372 Since this format requires single dates as opposed to date ranges,
1373 witnesses with closed date ranges will be mapped to the average of their lower and upper bounds,
1374 and witnesses with open date ranges will not be mapped.
1376 Args:
1377 taxlabels: A list of slugified taxon labels.
1379 Returns:
1380 A string containing comma-separated date calibrations of the form witness_id=date.
1381 """
1382 calibrate_strings = []
1383 for i, wit in enumerate(self.witnesses):
1384 taxlabel = taxlabels[i]
1385 date_range = wit.date_range
1386 # If either end of this witness's date range is empty, then do not include it:
1387 if date_range[0] is None or date_range[1] is None:
1388 continue
1389 # Otherwise, take the midpoint of its date range as its date:
1390 date = int((date_range[0] + date_range[1]) / 2)
1391 calibrate_string = "%s=%d" % (taxlabel, date)
1392 calibrate_strings.append(calibrate_string)
1393 # Then output the full date map string:
1394 date_map = ",".join(calibrate_strings)
1395 return date_map
1397 def get_beast_code_map_for_unit(self, symbols, missing_symbol, vu_ind):
1398 """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.
1399 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.
1401 Args:
1402 vu_ind: An integer index for the desired unit.
1404 Returns:
1405 A string containing comma-separated code mappings.
1406 """
1407 vu = self.variation_units[vu_ind]
1408 vu_id = vu.id
1409 code_map = {}
1410 for k in range(len(self.substantive_readings_by_variation_unit_id[vu.id])):
1411 code_map[symbols[k]] = str(k)
1412 # If this site is a singleton site, then add a code mapping for the dummy state:
1413 if len(self.substantive_readings_by_variation_unit_id[vu.id]) == 1:
1414 code_map[symbols[1]] = str(1)
1415 # Then add a mapping for the missing state, including a dummy state if this is a singleton site:
1416 code_map[missing_symbol] = " ".join(
1417 str(k) for k in range(len(self.substantive_readings_by_variation_unit_id[vu.id]))
1418 )
1419 # If this site is a singleton site, then add the dummy state to the missing state mapping:
1420 if len(self.substantive_readings_by_variation_unit_id[vu.id]) == 1:
1421 code_map[missing_symbol] = code_map[missing_symbol] + " " + str(1)
1422 # Then combine all of the mappings into a single string:
1423 code_map_string = ", ".join([code + "=" + code_map[code] for code in code_map])
1424 return code_map_string
1426 def get_beast_equilibrium_frequencies_for_unit(self, vu_ind):
1427 """Returns a string containing state/reading equilibrium frequencies in BEAST format for the character/variation unit at the given index.
1428 Since the equilibrium frequencies are not used with the substitution models, the equilibrium frequencies simply correspond to a uniform distribution over the states.
1429 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.
1431 Args:
1432 vu_ind: An integer index for the desired unit.
1434 Returns:
1435 A string containing space-separated equilibrium frequencies.
1436 """
1437 vu = self.variation_units[vu_ind]
1438 vu_id = vu.id
1439 # If this unit is a singleton, then return the string "0.5 0.5":
1440 if len(self.substantive_readings_by_variation_unit_id[vu_id]) == 1:
1441 return "0.5 0.5"
1442 # Otherwise, set the equilibrium frequencies according to a uniform distribution:
1443 equilibrium_frequencies = [1.0 / len(self.substantive_readings_by_variation_unit_id[vu_id])] * len(
1444 self.substantive_readings_by_variation_unit_id[vu_id]
1445 )
1446 equilibrium_frequencies_string = " ".join([str(w) for w in equilibrium_frequencies])
1447 return equilibrium_frequencies_string
1449 def get_beast_root_frequencies_for_unit(self, vu_ind):
1450 """Returns a string containing state/reading root frequencies in BEAST format for the character/variation unit at the given index.
1451 The root frequencies are calculated from the intrinsic odds at this unit.
1452 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.
1453 If no intrinsic odds are specified, then a uniform distribution over all states is assumed.
1455 Args:
1456 vu_ind: An integer index for the desired unit.
1458 Returns:
1459 A string containing space-separated root frequencies.
1460 """
1461 vu = self.variation_units[vu_ind]
1462 vu_id = vu.id
1463 intrinsic_relations = vu.intrinsic_relations
1464 intrinsic_odds_by_id = self.intrinsic_odds_by_id
1465 # If this unit is a singleton, then return the string "1 0":
1466 if len(self.substantive_readings_by_variation_unit_id[vu_id]) == 1:
1467 return "1 0"
1468 # If this unit has no intrinsic odds, then assume a uniform distribution over all readings:
1469 if len(intrinsic_relations) == 0:
1470 root_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 root_frequencies_string = " ".join([str(w) for w in root_frequencies])
1474 return root_frequencies_string
1475 # We will populate the root frequencies based on the intrinsic odds of the readings:
1476 root_frequencies_by_id = {}
1477 for rdg_id in self.substantive_readings_by_variation_unit_id[vu_id]:
1478 root_frequencies_by_id[rdg_id] = 0
1479 # First, construct an adjacency list for efficient edge iteration:
1480 neighbors_by_source = {}
1481 for edge in intrinsic_relations:
1482 s = edge[0]
1483 t = edge[1]
1484 if s not in neighbors_by_source:
1485 neighbors_by_source[s] = []
1486 if t not in neighbors_by_source:
1487 neighbors_by_source[t] = []
1488 neighbors_by_source[s].append(t)
1489 # Next, identify all readings that are not targeted by any intrinsic odds relation:
1490 in_degree_by_reading = {}
1491 for edge in intrinsic_relations:
1492 s = edge[0]
1493 t = edge[1]
1494 if s not in in_degree_by_reading:
1495 in_degree_by_reading[s] = 0
1496 if t not in in_degree_by_reading:
1497 in_degree_by_reading[t] = 0
1498 in_degree_by_reading[t] += 1
1499 starting_nodes = [t for t in in_degree_by_reading if in_degree_by_reading[t] == 0]
1500 # Set the root frequencies for these readings to 1 (they will be normalized later):
1501 for starting_node in starting_nodes:
1502 root_frequencies_by_id[starting_node] = 1.0
1503 # Next, set the frequencies for the remaining readings recursively using the adjacency list:
1504 def update_root_frequencies(s):
1505 for t in neighbors_by_source[s]:
1506 intrinsic_category = intrinsic_relations[(s, t)]
1507 odds = (
1508 intrinsic_odds_by_id[intrinsic_category]
1509 if intrinsic_odds_by_id[intrinsic_category] is not None
1510 else 1.0
1511 ) # TODO: This needs to be handled using parameters once we have it implemented in BEAST
1512 root_frequencies_by_id[t] = root_frequencies_by_id[s] / odds
1513 update_root_frequencies(t)
1514 return
1516 for starting_node in starting_nodes:
1517 update_root_frequencies(starting_node)
1518 # Then produce a normalized vector of root frequencies that corresponds to a probability distribution:
1519 root_frequencies = [
1520 root_frequencies_by_id[rdg_id] for rdg_id in self.substantive_readings_by_variation_unit_id[vu_id]
1521 ]
1522 total_frequencies = sum(root_frequencies)
1523 for k in range(len(root_frequencies)):
1524 root_frequencies[k] = root_frequencies[k] / total_frequencies
1525 root_frequencies_string = " ".join([str(w) for w in root_frequencies])
1526 return root_frequencies_string
1528 def to_beast(
1529 self,
1530 file_addr: Union[Path, str],
1531 drop_constant: bool = False,
1532 clock_model: ClockModel = ClockModel.strict,
1533 ancestral_logger: AncestralLogger = AncestralLogger.state,
1534 seed: int = None,
1535 ):
1536 """Writes this Collation to a file in BEAST format with the given address.
1538 Args:
1539 file_addr: A string representing the path to an output file.
1540 drop_constant (bool, optional): An optional flag indicating whether to ignore variation units with one substantive reading.
1541 clock_model: A ClockModel option indicating which clock model to use.
1542 ancestral_logger: An AncestralLogger option indicating which class of logger (if any) to use for ancestral states.
1543 seed: A seed for random number generation (for setting initial values of unspecified transcriptional rates).
1544 """
1545 # Populate a list of sites that will correspond to columns of the sequence alignment:
1546 substantive_variation_unit_ids = self.variation_unit_ids
1547 if drop_constant:
1548 substantive_variation_unit_ids = [
1549 vu_id
1550 for vu_id in self.variation_unit_ids
1551 if len(self.substantive_readings_by_variation_unit_id[vu_id]) > 1
1552 ]
1553 # Populate sets of substantive variation unit IDs and substantive variant reading tuples:
1554 substantive_variation_unit_ids_set = set(substantive_variation_unit_ids)
1555 substantive_variation_unit_reading_tuples_set = set(self.substantive_variation_unit_reading_tuples)
1556 # First, calculate the values we will be using for the main template:
1557 taxlabels = [slugify(wit.id, lowercase=False, separator='_') for wit in self.witnesses]
1558 missing_symbol = '?'
1559 symbols = self.get_beast_symbols()
1560 tip_date_range = self.get_tip_date_range()
1561 origin_span = self.get_beast_origin_span(tip_date_range)
1562 date_map = self.get_beast_date_map(taxlabels)
1563 # Then populate the necessary objects for the BEAST XML Jinja template:
1564 witness_objects = []
1565 variation_unit_objects = []
1566 intrinsic_category_objects = []
1567 transcriptional_category_objects = []
1568 with tqdm(total=len(self.witnesses)) as pbar:
1569 # Start with witnesses:
1570 for i, wit in enumerate(self.witnesses):
1571 witness_object = {}
1572 # Copy the ID for this witness:
1573 witness_object["id"] = wit.id
1574 # Copy its date bounds:
1575 witness_object["min_date"] = wit.date_range[0]
1576 witness_object["max_date"] = wit.date_range[1]
1577 # Populate its sequence from its entries in the witness's readings dictionary:
1578 sequence = ""
1579 for j, rdg_support in enumerate(self.readings_by_witness[wit.id]):
1580 vu_id = self.variation_unit_ids[j]
1581 # Skip any variation units deemed non-substantive:
1582 if vu_id not in substantive_variation_unit_ids:
1583 continue
1584 # If this witness has a certainty of 0 for all readings, then it is a gap; assign a likelihood of 1 to each reading:
1585 if sum(rdg_support) == 0:
1586 for k, w in enumerate(rdg_support):
1587 sequence += "1"
1588 if k < len(rdg_support) - 1:
1589 sequence += ", "
1590 else:
1591 if len(rdg_support) > 1:
1592 sequence += "; "
1593 else:
1594 # If this site is a singleton site, then add a dummy state:
1595 sequence += ", 0; "
1596 # Otherwise, read the probabilities as they are given:
1597 else:
1598 for k, w in enumerate(rdg_support):
1599 sequence += str(w)
1600 if k < len(rdg_support) - 1:
1601 sequence += ", "
1602 else:
1603 if len(rdg_support) > 1:
1604 sequence += "; "
1605 else:
1606 # If this site is a singleton site, then add a dummy state:
1607 sequence += ", 0; "
1608 # Strip the final semicolon and space from the sequence:
1609 sequence = sequence.strip("; ")
1610 # Then set the witness object's sequence attribute to this string:
1611 witness_object["sequence"] = sequence
1612 witness_objects.append(witness_object)
1613 pbar.update(1)
1614 # Then proceed to variation units:
1615 for j, vu in enumerate(self.variation_units):
1616 if vu.id not in substantive_variation_unit_ids_set:
1617 continue
1618 variation_unit_object = {}
1619 # Copy the one-based index of this variation unit:
1620 variation_unit_object["index"] = j + 1
1621 # Copy the ID of this variation unit:
1622 variation_unit_object["id"] = vu.id
1623 # Set a flag indicating if this variation unit is constant:
1624 variation_unit_object["is_constant"] = (
1625 True if len(self.substantive_readings_by_variation_unit_id[vu.id]) == 1 else False
1626 )
1627 # Copy this variation unit's number of substantive readings,
1628 # setting it to 2 if it is a singleton unit:
1629 variation_unit_object["nstates"] = (
1630 len(self.substantive_readings_by_variation_unit_id[vu.id])
1631 if len(self.substantive_readings_by_variation_unit_id[vu.id]) > 1
1632 else 2
1633 )
1634 # Then construct the code map for this unit:
1635 variation_unit_object["code_map"] = self.get_beast_code_map_for_unit(symbols, missing_symbol, j)
1636 # Then populate a comma-separated string of reading labels for this unit:
1637 rdg_texts = []
1638 vu_label = vu.id
1639 for rdg in vu.readings:
1640 key = tuple([vu.id, rdg.id])
1641 if key not in substantive_variation_unit_reading_tuples_set:
1642 continue
1643 rdg_text = slugify(rdg.text, lowercase=False, allow_unicode=True, separator='_')
1644 # Replace any empty reading text with an omission marker:
1645 if rdg_text == "":
1646 rdg_text = "om."
1647 rdg_texts.append(rdg_text)
1648 # If this site is a singleton site, then add a dummy reading for the dummy state:
1649 if len(self.substantive_readings_by_variation_unit_id[vu.id]) == 1:
1650 rdg_texts.append("DUMMY")
1651 rdg_texts_string = ", ".join(rdg_texts)
1652 variation_unit_object["rdg_texts"] = rdg_texts_string
1653 # Then populate this unit's equilibrium frequency string and its root frequency string:
1654 variation_unit_object["equilibrium_frequencies"] = self.get_beast_equilibrium_frequencies_for_unit(j)
1655 variation_unit_object["root_frequencies"] = self.get_beast_root_frequencies_for_unit(j)
1656 # Then populate a dictionary mapping epoch height ranges to lists of off-diagonal entries for substitution models:
1657 rate_objects_by_epoch_height_range = {}
1658 epoch_height_ranges = []
1659 # Then proceed based on whether the transcriptional relations for this variation unit have been defined:
1660 if len(vu.transcriptional_relations_by_date_range) == 0:
1661 # If there are no transcriptional relations, then map the epoch range of (None, None) to their list of off-diagonal entries:
1662 epoch_height_ranges.append((None, None))
1663 rate_objects_by_epoch_height_range[(None, None)] = []
1664 rate_objects = rate_objects_by_epoch_height_range[(None, None)]
1665 if len(self.substantive_readings_by_variation_unit_id[vu.id]) == 1:
1666 # If this is a singleton site, then use an arbitrary 2x2 rate matrix:
1667 rate_objects.append({"transcriptional_categories": ["default"], "expression": None})
1668 rate_objects.append({"transcriptional_categories": ["default"], "expression": None})
1669 else:
1670 # If this is a site with multiple substantive readings, but no transcriptional relations list,
1671 # then use a Lewis Mk substitution matrix with the appropriate number of states:
1672 for k_1, rdg_id_1 in enumerate(self.substantive_readings_by_variation_unit_id[vu.id]):
1673 for k_2, rdg_id_2 in enumerate(self.substantive_readings_by_variation_unit_id[vu.id]):
1674 # Skip diagonal elements:
1675 if k_1 == k_2:
1676 continue
1677 rate_objects.append({"transcriptional_categories": ["default"], "expression": None})
1678 else:
1679 # Otherwise, proceed for every date range:
1680 for date_range in vu.transcriptional_relations_by_date_range:
1681 # Get the map of transcriptional relations for reference later:
1682 transcriptional_relations = vu.transcriptional_relations_by_date_range[date_range]
1683 # Now get the epoch height range corresponding to this date range, and initialize its list in the dictionary:
1684 epoch_height_range = [None, None]
1685 epoch_height_range[0] = tip_date_range[1] - date_range[1] if date_range[1] is not None else None
1686 epoch_height_range[1] = tip_date_range[1] - date_range[0] if date_range[0] is not None else None
1687 epoch_height_range = tuple(epoch_height_range)
1688 epoch_height_ranges.append(epoch_height_range)
1689 rate_objects_by_epoch_height_range[epoch_height_range] = []
1690 rate_objects = rate_objects_by_epoch_height_range[epoch_height_range]
1691 # Then proceed for every pair of readings in this unit:
1692 for k_1, rdg_id_1 in enumerate(self.substantive_readings_by_variation_unit_id[vu.id]):
1693 for k_2, rdg_id_2 in enumerate(self.substantive_readings_by_variation_unit_id[vu.id]):
1694 # Skip diagonal elements:
1695 if k_1 == k_2:
1696 continue
1697 # If the first reading has no transcriptional relation to the second in this unit, then use the default rate:
1698 if (rdg_id_1, rdg_id_2) not in transcriptional_relations:
1699 rate_objects.append({"transcriptional_categories": ["default"], "expression": None})
1700 continue
1701 # Otherwise, if only one category of transcriptional relations holds between the first and second readings,
1702 # then use its rate:
1703 if len(transcriptional_relations[(rdg_id_1, rdg_id_2)]) == 1:
1704 # If there is only one such category, then add its rate as a standalone var element:
1705 transcriptional_category = list(transcriptional_relations[(rdg_id_1, rdg_id_2)])[0]
1706 rate_objects.append(
1707 {"transcriptional_categories": [transcriptional_category], "expression": None}
1708 )
1709 continue
1710 # If there is more than one, then add a var element that is a sum of the individual categories' rates:
1711 transcriptional_categories = list(transcriptional_relations[(rdg_id_1, rdg_id_2)])
1712 args = []
1713 for transcriptional_category in transcriptional_categories:
1714 args.append("%s_rate" % transcriptional_category)
1715 args_string = " ".join(args)
1716 ops = ["+"] * (len(args) - 1)
1717 ops_string = " ".join(ops)
1718 expression_string = " ".join([args_string, ops_string])
1719 rate_objects.append(
1720 {
1721 "transcriptional_categories": transcriptional_categories,
1722 "expression": expression_string,
1723 }
1724 )
1725 # Now reorder the list of epoch height ranges, and get a list of non-null epoch dates in ascending order from the dictionary:
1726 epoch_height_ranges.reverse()
1727 epoch_heights = [
1728 epoch_height_range[0] for epoch_height_range in epoch_height_ranges if epoch_height_range[0] is not None
1729 ]
1730 # Then add all of these data structures to the variation unit object:
1731 variation_unit_object["epoch_heights"] = epoch_heights
1732 variation_unit_object["epoch_heights_string"] = " ".join(
1733 [str(epoch_height) for epoch_height in epoch_heights]
1734 )
1735 variation_unit_object["epoch_height_ranges"] = epoch_height_ranges
1736 variation_unit_object["epoch_rates"] = [
1737 rate_objects_by_epoch_height_range[epoch_height_range] for epoch_height_range in epoch_height_ranges
1738 ]
1739 variation_unit_objects.append(variation_unit_object)
1740 # Then proceed to intrinsic odds categories:
1741 for intrinsic_category in self.intrinsic_categories:
1742 intrinsic_category_object = {}
1743 # Copy the ID of this intrinsic category:
1744 intrinsic_category_object["id"] = intrinsic_category
1745 # Then copy the odds factors associated with this intrinsic category,
1746 # setting it to 1.0 if it is not specified and setting the estimate attribute accordingly:
1747 odds = self.intrinsic_odds_by_id[intrinsic_category]
1748 intrinsic_category_object["odds"] = odds if odds is not None else 1.0
1749 intrinsic_category_object["estimate"] = "false" if odds is not None else "true"
1750 intrinsic_category_objects.append(intrinsic_category_object)
1751 # Then proceed to transcriptional rate categories:
1752 rng = np.random.default_rng(seed)
1753 for transcriptional_category in self.transcriptional_categories:
1754 transcriptional_category_object = {}
1755 # Copy the ID of this transcriptional category:
1756 transcriptional_category_object["id"] = transcriptional_category
1757 # Then copy the rate of this transcriptional category,
1758 # setting it to a random number sampled from a Gamma distribution if it is not specified and setting the estimate attribute accordingly:
1759 rate = self.transcriptional_rates_by_id[transcriptional_category]
1760 transcriptional_category_object["rate"] = rate if rate is not None else rng.gamma(5.0, 2.0)
1761 transcriptional_category_object["estimate"] = "false" if rate is not None else "true"
1762 transcriptional_category_objects.append(transcriptional_category_object)
1763 # Now render the output XML file using the Jinja template:
1764 env = Environment(loader=PackageLoader("teiphy", "templates"), autoescape=select_autoescape())
1765 template = env.get_template("beast_template.xml")
1766 rendered = template.render(
1767 nsymbols=len(symbols),
1768 date_map=date_map,
1769 origin_span=origin_span,
1770 clock_model=clock_model.value,
1771 clock_rate_categories=2 * len(self.witnesses) - 1,
1772 ancestral_logger=ancestral_logger.value,
1773 witnesses=witness_objects,
1774 variation_units=variation_unit_objects,
1775 non_constant_variation_units=[
1776 variation_unit_object
1777 for variation_unit_object in variation_unit_objects
1778 if not variation_unit_object["is_constant"]
1779 ],
1780 constant_variation_unit_filter=",".join(
1781 [
1782 str(variation_unit_object["index"])
1783 for variation_unit_object in variation_unit_objects
1784 if variation_unit_object["is_constant"]
1785 ]
1786 ),
1787 intrinsic_categories=intrinsic_category_objects,
1788 transcriptional_categories=transcriptional_category_objects,
1789 )
1790 # Generate all parent folders for this file that don't already exist:
1791 Path(file_addr).parent.mkdir(parents=True, exist_ok=True)
1792 with open(file_addr, "w", encoding="utf-8") as f:
1793 f.write(rendered)
1794 return
1796 def to_numpy(self, drop_constant: bool = False, split_missing: SplitMissingType = None):
1797 """Returns this Collation in the form of a NumPy array, along with arrays of its row and column labels.
1799 Args:
1800 drop_constant (bool, optional): An optional flag indicating whether to ignore variation units with one substantive reading.
1801 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.
1802 If not specified, then missing data is ignored (i.e., all states are 0).
1803 If "uniform", then the contribution of 1 is divided evenly over all substantive readings.
1804 If "proportional", then the contribution of 1 is divided between the readings in proportion to their support among the witnesses that are not missing.
1806 Returns:
1807 A NumPy array with a row for each substantive reading and a column for each witness.
1808 A list of substantive reading ID strings.
1809 A list of witness ID strings.
1810 """
1811 # Populate a list of sites that will correspond to columns of the sequence alignment:
1812 substantive_variation_unit_ids = self.variation_unit_ids
1813 if drop_constant:
1814 substantive_variation_unit_ids = [
1815 vu_id
1816 for vu_id in self.variation_unit_ids
1817 if len(self.substantive_readings_by_variation_unit_id[vu_id]) > 1
1818 ]
1819 substantive_variation_unit_ids_set = set(substantive_variation_unit_ids)
1820 substantive_variation_unit_reading_tuples_set = set(self.substantive_variation_unit_reading_tuples)
1821 # Initialize the output array with the appropriate dimensions:
1822 reading_labels = []
1823 for vu in self.variation_units:
1824 if vu.id not in substantive_variation_unit_ids_set:
1825 continue
1826 for rdg in vu.readings:
1827 key = tuple([vu.id, rdg.id])
1828 if key in substantive_variation_unit_reading_tuples_set:
1829 reading_labels.append(vu.id + ", " + rdg.text)
1830 witness_labels = [wit.id for wit in self.witnesses]
1831 matrix = np.zeros((len(reading_labels), len(witness_labels)), dtype=float)
1832 # For each variation unit, keep a record of the proportion of non-missing witnesses supporting the substantive variant readings:
1833 support_proportions_by_unit = {}
1834 for j, vu_id in enumerate(self.variation_unit_ids):
1835 if vu_id not in substantive_variation_unit_ids_set:
1836 continue
1837 support_proportions = [0.0] * len(self.substantive_readings_by_variation_unit_id[vu_id])
1838 for i, wit in enumerate(self.witnesses):
1839 rdg_support = self.readings_by_witness[wit.id][j]
1840 for l, w in enumerate(rdg_support):
1841 support_proportions[l] += w
1842 norm = (
1843 sum(support_proportions) if sum(support_proportions) > 0 else 1.0
1844 ) # 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
1845 for l in range(len(support_proportions)):
1846 support_proportions[l] = support_proportions[l] / norm
1847 support_proportions_by_unit[vu_id] = support_proportions
1848 # Then populate it with the appropriate values:
1849 col_ind = 0
1850 with tqdm(total=len(self.witnesses)) as pbar:
1851 for i, wit in enumerate(self.witnesses):
1852 row_ind = 0
1853 for j, vu_id in enumerate(self.variation_unit_ids):
1854 if vu_id not in substantive_variation_unit_ids_set:
1855 continue
1856 rdg_support = self.readings_by_witness[wit.id][j]
1857 # If this reading support vector sums to 0, then this is missing data; handle it as specified:
1858 if sum(rdg_support) == 0:
1859 if split_missing == SplitMissingType.uniform:
1860 for l in range(len(rdg_support)):
1861 matrix[row_ind, col_ind] = 1 / len(rdg_support)
1862 row_ind += 1
1863 elif split_missing == SplitMissingType.proportional:
1864 for l in range(len(rdg_support)):
1865 matrix[row_ind, col_ind] = support_proportions_by_unit[vu_id][l]
1866 row_ind += 1
1867 else:
1868 row_ind += len(rdg_support)
1869 # Otherwise, add its coefficients normally:
1870 else:
1871 for l in range(len(rdg_support)):
1872 matrix[row_ind, col_ind] = rdg_support[l]
1873 row_ind += 1
1874 col_ind += 1
1875 pbar.update(1)
1876 return matrix, reading_labels, witness_labels
1878 def get_ext_matrix(self, drop_constant: bool = False, split_missing: SplitMissingType = None):
1879 """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.
1880 Note that if the split_missing option is specified, all variation units are counted.
1882 Args:
1883 drop_constant (bool, optional): An optional flag indicating whether to ignore variation units with one substantive reading.
1884 Default value is False.
1885 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.
1886 If not specified, then missing data is ignored (i.e., all states are 0).
1887 If "uniform", then the contribution of 1 is divided evenly over all substantive readings.
1888 If "proportional", then the contribution of 1 is divided between the readings in proportion to their support among the witnesses that are not missing.
1890 Returns:
1891 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.
1892 """
1893 # Populate a list of sites that will correspond to columns of the sequence alignment:
1894 substantive_variation_unit_ids = self.variation_unit_ids
1895 if drop_constant:
1896 substantive_variation_unit_ids = [
1897 vu_id
1898 for vu_id in self.variation_unit_ids
1899 if len(self.substantive_readings_by_variation_unit_id[vu_id]) > 1
1900 ]
1901 substantive_variation_unit_ids_set = set(substantive_variation_unit_ids)
1902 # Then initialize the output matrix:
1903 witness_labels = [wit.id for wit in self.witnesses]
1904 ext_matrix = ext_matrix = np.full((len(witness_labels), len(witness_labels)), 0, dtype=int)
1905 # If the split_missing option has been specified, then all entries in the matrix will be the number of substantive variation units,
1906 # so we can just fill the matrix with this value and return it:
1907 if split_missing is not None:
1908 ext_matrix = ext_matrix = np.full(
1909 (len(witness_labels), len(witness_labels)), len(substantive_variation_unit_ids), dtype=int
1910 )
1911 return ext_matrix
1912 # Otherwise, populate the matrix for all pairs of witnesses:
1913 with tqdm(total=len(self.witnesses) ** 2) as pbar:
1914 # Then calculate the mutual information contribution for each pair of witnesses:
1915 for i, wit_1 in enumerate(witness_labels):
1916 for j, wit_2 in enumerate(witness_labels):
1917 shared_ext_units = 0
1918 # The contribution to the entry for these witnesses will be identical regardless of the order in which they are specified,
1919 # so we only have to calculate it once:
1920 if i > j:
1921 pbar.update(1)
1922 continue
1923 # Otherwise, calculate the number of substantive variation units at which both of these witnesses are extant:
1924 for k, vu_id in enumerate(self.variation_unit_ids):
1925 if vu_id not in substantive_variation_unit_ids_set:
1926 continue
1927 wit_1_rdg_support = self.readings_by_witness[wit_1][k]
1928 wit_2_rdg_support = self.readings_by_witness[wit_2][k]
1929 if sum(wit_1_rdg_support) == 0.0 or sum(wit_2_rdg_support) == 0.0:
1930 continue
1931 shared_ext_units += 1
1932 ext_matrix[i][j] = shared_ext_units
1933 ext_matrix[j][i] = shared_ext_units
1934 pbar.update(1)
1935 return ext_matrix
1937 def transform_matrix(self, matrix: np.ndarray, transform_matrix: TransformMatrixType = None):
1938 """Transforms a given matrix's columns based on the specified transform_matrix option.
1940 Args:
1941 matrix (ndarray): The matrix whose columns are to be transformed.
1942 transform_matrix (TransformMatrixType, optional): A TransformMatrixType option indicating how the columns of a witness-to-witness matrix output should be transformed.
1943 Only applicable for tabular outputs in which the rows and columns correspond to the witnesses in the collation.
1945 Returns:
1946 A Pandas DataFrame corresponding to a collation matrix with reading frequencies or a long table with discrete reading states.
1947 """
1948 # If no transform_matrix option was supplied, then return the matrix as-is:
1949 if transform_matrix is None:
1950 return matrix
1951 # Otherwise, apply the specified column transformation:
1952 if transform_matrix == TransformMatrixType.stddev:
1953 means = np.mean(matrix, 0) # get the means of the columns (axis 0)
1954 stddevs = np.std(matrix, 0) # get the standard deviations of the columns (axis 0)
1955 transformed_matrix = np.full(
1956 matrix.shape, 0.0, dtype=float
1957 ) # the standard deviation can only be 0 if all values equal the mean, so a default value of 0 is acceptable
1958 np.divide(matrix - means, stddevs, out=transformed_matrix, where=(stddevs != 0))
1959 return transformed_matrix
1960 if transform_matrix == TransformMatrixType.mad:
1961 medians = np.median(matrix, 0) # get the medians of the columns (axis 0)
1962 mads = np.median(np.abs(matrix - medians), 0) # get the median absolute deviations of the columns (axis 0)
1963 transformed_matrix = np.full(
1964 matrix.shape, np.nan, dtype=float
1965 ) # 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
1966 np.divide(matrix - medians, mads, out=transformed_matrix, where=(mads != 0))
1967 return transformed_matrix
1969 def to_distance_matrix(
1970 self,
1971 drop_constant: bool = False,
1972 proportion: bool = False,
1973 show_ext: bool = False,
1974 transform_matrix: TransformMatrixType = None,
1975 ):
1976 """Transforms this Collation into a NumPy distance matrix between witnesses, along with an array of its labels for the witnesses.
1977 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.
1978 Optionally, the count of units where both witnesses have singleton readings can be included after the count/proportion of disagreements.
1980 Args:
1981 drop_constant (bool, optional): An optional flag indicating whether to ignore variation units with one substantive reading.
1982 Default value is False.
1983 proportion (bool, optional): An optional flag indicating whether or not to calculate distances as proportions over extant, unambiguous variation units.
1984 Default value is False.
1985 show_ext: An optional flag indicating whether each cell in the matrix
1986 should include the number of their extant, unambiguous variation units after the number of their disagreements.
1987 Default value is False.
1988 transform_matrix (TransformMatrixType, optional): A TransformMatrixType option indicating how the columns of the matrix should be transformed.
1990 Returns:
1991 A NumPy distance matrix with a row and column for each witness.
1992 A list of witness ID strings.
1993 """
1994 # Populate a list of sites that will correspond to columns of the sequence alignment:
1995 substantive_variation_unit_ids = self.variation_unit_ids
1996 if drop_constant:
1997 substantive_variation_unit_ids = [
1998 vu_id
1999 for vu_id in self.variation_unit_ids
2000 if len(self.substantive_readings_by_variation_unit_id[vu_id]) > 1
2001 ]
2002 substantive_variation_unit_ids_set = set(substantive_variation_unit_ids)
2003 substantive_variation_unit_reading_tuples_set = set(self.substantive_variation_unit_reading_tuples)
2004 # Initialize the output array with the appropriate dimensions:
2005 witness_labels = [wit.id for wit in self.witnesses]
2006 matrix = np.full((len(witness_labels), len(witness_labels)), 0, dtype=int) # ints of the form disagreements
2007 with tqdm(total=len(self.witnesses) ** 2) as pbar:
2008 for i, wit_1 in enumerate(witness_labels):
2009 for j, wit_2 in enumerate(witness_labels):
2010 disagreements = 0
2011 # The contribution to the entry for these witnesses will be identical regardless of the order in which they are specified,
2012 # so we only have to calculate it once:
2013 if i > j:
2014 pbar.update(1)
2015 continue
2016 # Otherwise, calculate the number of units where both witnesses disagree:
2017 for k, vu_id in enumerate(self.variation_unit_ids):
2018 if vu_id not in substantive_variation_unit_ids_set:
2019 continue
2020 wit_1_rdg_support = self.readings_by_witness[wit_1][k]
2021 wit_2_rdg_support = self.readings_by_witness[wit_2][k]
2022 # If either witness is lacunose, then move on:
2023 if sum(wit_1_rdg_support) == 0.0 or sum(wit_2_rdg_support) == 0.0:
2024 continue
2025 # Otherwise, if the (potential) readings of the two witnesses do not overlap, then count them as disagreeing:
2026 if (
2027 sum([wit_1_rdg_support[l] * wit_2_rdg_support[l] for l in range(len(wit_1_rdg_support))])
2028 == 0.0
2029 ):
2030 disagreements += 1
2031 matrix[i, j] = disagreements
2032 matrix[j, i] = disagreements
2033 pbar.update(1)
2034 # Initialize a matrix for shared extant variation units for witnesses, and populate it if the proportion or show_ext option is specified:
2035 ext_matrix = None
2036 if proportion or show_ext:
2037 ext_matrix = self.get_ext_matrix(drop_constant=drop_constant)
2038 # 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:
2039 if proportion:
2040 proportion_matrix = np.full((len(witness_labels), len(witness_labels)), 0.0, dtype=float)
2041 np.divide(
2042 matrix, ext_matrix, out=proportion_matrix, where=(ext_matrix != 0)
2043 ) # division by 0 can occur if two witnesses have no overlapping units; leave their proportion as 0.0
2044 matrix = proportion_matrix
2045 # Then transform the columns of the main matrix as specified:
2046 matrix = self.transform_matrix(matrix, transform_matrix)
2047 # If the show_ext option is set, then append the number of shared extant variation units after the matrix's values:
2048 if show_ext:
2049 serialized_values = []
2050 for i, wit_1 in enumerate(witness_labels):
2051 serialized_values.append([])
2052 for j, wit_2 in enumerate(witness_labels):
2053 serialized_values[-1].append("/".join([str(matrix[i][j]), str(ext_matrix[i][j])]))
2054 matrix = np.array(serialized_values)
2055 return matrix, witness_labels
2057 def to_similarity_matrix(
2058 self,
2059 drop_constant: bool = False,
2060 proportion: bool = False,
2061 show_ext: bool = False,
2062 transform_matrix: TransformMatrixType = None,
2063 ):
2064 """Transforms this Collation into a NumPy similarity matrix between witnesses, along with an array of its labels for the witnesses.
2065 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.
2066 Optionally, the count of units where both witnesses have singleton readings can be included after the count/proportion of agreements.
2068 Args:
2069 drop_constant (bool, optional): An optional flag indicating whether to ignore variation units with one substantive reading.
2070 Default value is False.
2071 proportion (bool, optional): An optional flag indicating whether or not to calculate similarities as proportions over extant, unambiguous variation units.
2072 Default value is False.
2073 show_ext: An optional flag indicating whether each cell in the matrix
2074 should include the number of their extant, unambiguous variation units after the number of agreements.
2075 Default value is False.
2076 transform_matrix (TransformMatrixType, optional): A TransformMatrixType option indicating how the columns of the matrix should be transformed.
2078 Returns:
2079 A NumPy agreement matrix with a row and column for each witness.
2080 A list of witness ID strings.
2081 """
2082 # Populate a list of sites that will correspond to columns of the sequence alignment:
2083 substantive_variation_unit_ids = self.variation_unit_ids
2084 if drop_constant:
2085 substantive_variation_unit_ids = [
2086 vu_id
2087 for vu_id in self.variation_unit_ids
2088 if len(self.substantive_readings_by_variation_unit_id[vu_id]) > 1
2089 ]
2090 substantive_variation_unit_ids_set = set(substantive_variation_unit_ids)
2091 substantive_variation_unit_reading_tuples_set = set(self.substantive_variation_unit_reading_tuples)
2092 # Initialize the output array with the appropriate dimensions:
2093 witness_labels = [wit.id for wit in self.witnesses]
2094 matrix = np.full((len(witness_labels), len(witness_labels)), 0, dtype=int) # ints of the form agreements
2095 with tqdm(total=len(self.witnesses) ** 2) as pbar:
2096 for i, wit_1 in enumerate(witness_labels):
2097 for j, wit_2 in enumerate(witness_labels):
2098 agreements = 0
2099 # The contribution to the entry for these witnesses will be identical regardless of the order in which they are specified,
2100 # so we only have to calculate it once:
2101 if i > j:
2102 pbar.update(1)
2103 continue
2104 # Otherwise, calculate the number of units where both witnesses unambiguously agree:
2105 for k, vu_id in enumerate(self.variation_unit_ids):
2106 if vu_id not in substantive_variation_unit_ids_set:
2107 continue
2108 wit_1_rdg_support = self.readings_by_witness[wit_1][k]
2109 wit_2_rdg_support = self.readings_by_witness[wit_2][k]
2110 wit_1_rdg_inds = [l for l, w in enumerate(wit_1_rdg_support) if w > 0]
2111 wit_2_rdg_inds = [l for l, w in enumerate(wit_2_rdg_support) if w > 0]
2112 if len(wit_1_rdg_inds) != 1 or len(wit_2_rdg_inds) != 1:
2113 continue
2114 if wit_1_rdg_inds[0] == wit_2_rdg_inds[0]:
2115 agreements += 1
2116 matrix[i, j] = agreements
2117 matrix[j, i] = agreements
2118 pbar.update(1)
2119 # Initialize a matrix for shared extant variation units for witnesses, and populate it if the proportion or show_ext option is specified:
2120 ext_matrix = None
2121 if proportion or show_ext:
2122 ext_matrix = self.get_ext_matrix(drop_constant=drop_constant)
2123 # 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:
2124 if proportion:
2125 proportion_matrix = np.full((len(witness_labels), len(witness_labels)), 0.0, dtype=float)
2126 np.divide(
2127 matrix, ext_matrix, out=proportion_matrix, where=(ext_matrix != 0)
2128 ) # division by 0 can occur if two witnesses have no overlapping units; leave their proportion as 0.0
2129 matrix = proportion_matrix
2130 # Then transform the columns of the main matrix as specified:
2131 matrix = self.transform_matrix(matrix, transform_matrix)
2132 # If the show_ext option is set, then append the number of shared extant variation units after the matrix's values:
2133 if show_ext:
2134 serialized_values = []
2135 for i, wit_1 in enumerate(witness_labels):
2136 serialized_values.append([])
2137 for j, wit_2 in enumerate(witness_labels):
2138 serialized_values[-1].append("/".join([str(matrix[i][j]), str(ext_matrix[i][j])]))
2139 matrix = np.array(serialized_values)
2140 return matrix, witness_labels
2142 def to_idf_matrix(
2143 self,
2144 drop_constant: bool = False,
2145 split_missing: SplitMissingType = None,
2146 proportion: bool = False,
2147 show_ext: bool = False,
2148 transform_matrix: TransformMatrixType = None,
2149 ):
2150 """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.
2151 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.
2152 The IDF-weighted agreement score for two witnesses is the sum of the IDF weights for the readings at which they agree.
2153 Where any witness is ambiguous, it contributes to its potential readings' sampling probabilities in proportion to its degrees of support for those readings.
2154 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.
2155 If a split_missing argument is supplied, then lacunae are handled in the same way.
2156 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.
2157 (If the split_missing argument is also specified, then the mean is taken over all substantive variation units.)
2159 Args:
2160 drop_constant (bool, optional): An optional flag indicating whether to ignore variation units with one substantive reading.
2161 Default value is False.
2162 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.
2163 If not specified, then missing data is ignored (i.e., all states are 0).
2164 If "uniform", then the contribution of 1 is divided evenly over all substantive readings.
2165 If "proportional", then the contribution of 1 is divided between the readings in proportion to their support among the witnesses that are not missing.
2166 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.
2167 Default value is False.
2168 show_ext: An optional flag indicating whether each cell in the matrix
2169 should include the number of their extant, unambiguous variation units after the number of agreements.
2170 Default value is False.
2171 transform_matrix (TransformMatrixType, optional): A TransformMatrixType option indicating how the columns of the matrix should be transformed.
2173 Returns:
2174 A NumPy IDF-weighted agreement matrix with a row and column for each witness.
2175 A list of witness ID strings.
2176 """
2177 # Populate a list of sites that will correspond to columns of the sequence alignment:
2178 substantive_variation_unit_ids = self.variation_unit_ids
2179 if drop_constant:
2180 substantive_variation_unit_ids = [
2181 vu_id
2182 for vu_id in self.variation_unit_ids
2183 if len(self.substantive_readings_by_variation_unit_id[vu_id]) > 1
2184 ]
2185 substantive_variation_unit_ids_set = set(substantive_variation_unit_ids)
2186 substantive_variation_unit_reading_tuples_set = set(self.substantive_variation_unit_reading_tuples)
2187 # Initialize the output array with the appropriate dimensions:
2188 witness_labels = [wit.id for wit in self.witnesses]
2189 # 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:
2190 support_proportions_by_unit = {}
2191 if split_missing == SplitMissingType.proportional:
2192 for k, vu_id in enumerate(self.variation_unit_ids):
2193 # Skip this variation unit if it is a dropped constant site:
2194 if vu_id not in substantive_variation_unit_ids_set:
2195 continue
2196 support_proportions = [0.0] * len(self.substantive_readings_by_variation_unit_id[vu_id])
2197 for i, wit in enumerate(witness_labels):
2198 rdg_support = self.readings_by_witness[wit][k]
2199 for l, w in enumerate(rdg_support):
2200 support_proportions[l] += w
2201 norm = (
2202 sum(support_proportions) if sum(support_proportions) > 0 else 1.0
2203 ) # 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
2204 for l in range(len(support_proportions)):
2205 support_proportions[l] = support_proportions[l] / norm
2206 support_proportions_by_unit[vu_id] = support_proportions
2207 # Then populate data structures mapping each variation unit's ID to normalized reading support dictionaries
2208 # and vectors of sampling probabilities for its substantive readings:
2209 normalized_reading_support_dicts_by_vu_id = {}
2210 sampling_probabilities_by_vu_id = {}
2211 for k, vu_id in enumerate(self.variation_unit_ids):
2212 # Skip this variation unit if it is a dropped constant site:
2213 if vu_id not in substantive_variation_unit_ids_set:
2214 continue
2215 # Otherwise, populate normalized reading support vector dictionaries and sampling probability vectors in this unit:
2216 normalized_reading_support_by_wit = {}
2217 sampling_probabilities = [0.0] * len(self.substantive_readings_by_variation_unit_id[vu_id])
2218 for i, wit in enumerate(witness_labels):
2219 rdg_support = self.readings_by_witness[wit][k]
2220 # Check if this reading support vector represents missing data:
2221 norm = sum(rdg_support)
2222 if norm == 0:
2223 # If this reading support vector sums to 0, then this is missing data; handle it as specified:
2224 if split_missing == SplitMissingType.uniform:
2225 rdg_support = [1 / len(rdg_support) for l in range(len(rdg_support))]
2226 elif split_missing == SplitMissingType.proportional:
2227 rdg_support = [support_proportions_by_unit[vu_id][l] for l in range(len(rdg_support))]
2228 else:
2229 # Otherwise, the data is present, though it may be ambiguous; normalize the reading probabilities to sum to 1:
2230 rdg_support = [w / norm for l, w in enumerate(rdg_support)]
2231 normalized_reading_support_by_wit[wit] = rdg_support
2232 # Then add this witness's contributions to the readings' sampling probabilities:
2233 for l, w in enumerate(normalized_reading_support_by_wit[wit]):
2234 sampling_probabilities[l] += w
2235 norm = (
2236 sum(sampling_probabilities) if sum(sampling_probabilities) > 0 else 1.0
2237 ) # 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
2238 # Otherwise, normalize the sampling probabilities so they sum to 1:
2239 sampling_probabilities = [w / norm for w in sampling_probabilities]
2240 normalized_reading_support_dicts_by_vu_id[vu_id] = normalized_reading_support_by_wit
2241 sampling_probabilities_by_vu_id[vu_id] = sampling_probabilities
2242 # Then populate the matrix with the total expected information content for agreements between each pair of witnesses:
2243 matrix = np.full((len(witness_labels), len(witness_labels)), 0, dtype=float)
2244 with tqdm(total=len(self.witnesses) ** 2) as pbar:
2245 for i, wit_1 in enumerate(witness_labels):
2246 for j, wit_2 in enumerate(witness_labels):
2247 total_information_content = 0.0
2248 # The contribution to the entry for these witnesses will be identical regardless of the order in which they are specified,
2249 # so we only have to calculate it once:
2250 if i > j:
2251 pbar.update(1)
2252 continue
2253 # Otherwise, calculate the expected information content of agreements between these witnesses in each substantive variation unit
2254 # based on the sampling probabilities of the substantive readings in the unit:
2255 for k, vu_id in enumerate(self.variation_unit_ids):
2256 if vu_id not in substantive_variation_unit_ids_set:
2257 continue
2258 wit_1_rdg_support = normalized_reading_support_dicts_by_vu_id[vu_id][wit_1]
2259 wit_2_rdg_support = normalized_reading_support_dicts_by_vu_id[vu_id][wit_2]
2260 sampling_probabilities = sampling_probabilities_by_vu_id[vu_id]
2261 # First, calculate the probability that these two witnesses agree:
2262 probability_of_agreement = sum(
2263 [wit_1_rdg_support[l] * wit_2_rdg_support[l] for l in range(len(sampling_probabilities))]
2264 )
2265 # If these witnesses do not agree at this variation unit, then this unit contributes nothing to their total score:
2266 if probability_of_agreement == 0.0:
2267 continue
2268 # Otherwise, calculate the expected information content (in bits) of their agreement given their agreement on that reading
2269 # (skipping readings with a sampling probability of 0):
2270 expected_information_content = sum(
2271 [
2272 -math.log2(sampling_probabilities[l])
2273 * (wit_1_rdg_support[l] * wit_2_rdg_support[l] / probability_of_agreement)
2274 for l in range(len(sampling_probabilities))
2275 if sampling_probabilities[l] > 0.0
2276 ]
2277 )
2278 # Then add this contribution to the total score for these two witnesses:
2279 total_information_content += expected_information_content
2280 matrix[i, j] = total_information_content
2281 matrix[j, i] = total_information_content
2282 pbar.update(1)
2283 # Initialize a matrix for shared extant variation units for witnesses, and populate it if the proportion or show_ext option is specified:
2284 ext_matrix = None
2285 if proportion or show_ext:
2286 ext_matrix = self.get_ext_matrix(drop_constant=drop_constant, split_missing=split_missing)
2287 # 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:
2288 if proportion:
2289 proportion_matrix = np.full((len(witness_labels), len(witness_labels)), 0.0, dtype=float)
2290 np.divide(
2291 matrix, ext_matrix, out=proportion_matrix, where=(ext_matrix != 0)
2292 ) # division by 0 can occur if two witnesses have no overlapping units; leave their proportion as 0.0
2293 matrix = proportion_matrix
2294 # Then transform the columns of the main matrix as specified:
2295 matrix = self.transform_matrix(matrix, transform_matrix)
2296 # If the show_ext option is set, then append the number of shared extant variation units after the matrix's values:
2297 if show_ext:
2298 serialized_values = []
2299 for i, wit_1 in enumerate(witness_labels):
2300 serialized_values.append([])
2301 for j, wit_2 in enumerate(witness_labels):
2302 serialized_values[-1].append("/".join([str(matrix[i][j]), str(ext_matrix[i][j])]))
2303 matrix = np.array(serialized_values)
2304 return matrix, witness_labels
2306 def to_mi_matrix(
2307 self,
2308 drop_constant: bool = False,
2309 split_missing: SplitMissingType = None,
2310 proportion: bool = False,
2311 show_ext: bool = False,
2312 transform_matrix: TransformMatrixType = None,
2313 ):
2314 """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.
2315 This is equivalent to the total Kullback-Leibler divergence of the joint distribution of the witnesses' observed readings
2316 from the joint distribution of their expected readings under the assumption that the witnesses are independent, taken over all variation units.
2317 The value of 0 if and only if the witnesses are completely independent.
2319 Args:
2320 drop_constant (bool, optional): An optional flag indicating whether to ignore variation units with one substantive reading.
2321 Default value is False.
2322 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.
2323 If not specified, then missing data is ignored (i.e., all states are 0).
2324 If "uniform", then the contribution of 1 is divided evenly over all substantive readings.
2325 If "proportional", then the contribution of 1 is divided between the readings in proportion to their support among the witnesses that are not missing.
2326 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.
2327 Default value is False.
2328 show_ext: An optional flag indicating whether each cell in the matrix
2329 should include the number of their extant, unambiguous variation units after the number of agreements.
2330 Default value is False.
2331 transform_matrix (TransformMatrixType, optional): A TransformMatrixType option indicating how the columns of the matrix should be transformed.
2333 Returns:
2334 A NumPy MI matrix with a row and column for each witness.
2335 A list of witness ID strings.
2336 """
2337 # Populate a list of sites that will correspond to columns of the sequence alignment:
2338 substantive_variation_unit_ids = self.variation_unit_ids
2339 if drop_constant:
2340 substantive_variation_unit_ids = [
2341 vu_id
2342 for vu_id in self.variation_unit_ids
2343 if len(self.substantive_readings_by_variation_unit_id[vu_id]) > 1
2344 ]
2345 substantive_variation_unit_ids_set = set(substantive_variation_unit_ids)
2346 substantive_variation_unit_reading_tuples_set = set(self.substantive_variation_unit_reading_tuples)
2347 # Initialize the output array with the appropriate dimensions:
2348 witness_labels = [wit.id for wit in self.witnesses]
2349 # 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:
2350 support_proportions_by_unit = {}
2351 if split_missing == SplitMissingType.proportional:
2352 for k, vu_id in enumerate(self.variation_unit_ids):
2353 # Skip this variation unit if it is a dropped constant site:
2354 if vu_id not in substantive_variation_unit_ids_set:
2355 continue
2356 support_proportions = [0.0] * len(self.substantive_readings_by_variation_unit_id[vu_id])
2357 for i, wit in enumerate(witness_labels):
2358 rdg_support = self.readings_by_witness[wit][k]
2359 for l, w in enumerate(rdg_support):
2360 support_proportions[l] += w
2361 norm = (
2362 sum(support_proportions) if sum(support_proportions) > 0 else 1.0
2363 ) # 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
2364 for l in range(len(support_proportions)):
2365 support_proportions[l] = support_proportions[l] / norm
2366 support_proportions_by_unit[vu_id] = support_proportions
2367 # Then populate data structures mapping each variation unit's ID to normalized reading support dictionaries,
2368 # vectors of sampling probabilities for its substantive readings, and expected joint probability matrices:
2369 normalized_reading_support_dicts_by_vu_id = {}
2370 sampling_probabilities_by_vu_id = {}
2371 expected_joint_probabilities_by_vu_id = {}
2372 for k, vu_id in enumerate(self.variation_unit_ids):
2373 # Skip this variation unit if it is a dropped constant site:
2374 if vu_id not in substantive_variation_unit_ids_set:
2375 continue
2376 # Otherwise, populate normalized reading support vector dictionaries and sampling probability vectors in this unit:
2377 normalized_reading_support_by_wit = {}
2378 sampling_probabilities = [0.0] * len(self.substantive_readings_by_variation_unit_id[vu_id])
2379 for i, wit in enumerate(witness_labels):
2380 rdg_support = self.readings_by_witness[wit][k]
2381 # Check if this reading support vector represents missing data:
2382 norm = sum(rdg_support)
2383 if norm == 0:
2384 # If this reading support vector sums to 0, then this is missing data; handle it as specified:
2385 if split_missing == SplitMissingType.uniform:
2386 rdg_support = [1 / len(rdg_support) for l in range(len(rdg_support))]
2387 elif split_missing == SplitMissingType.proportional:
2388 rdg_support = [support_proportions_by_unit[vu_id][l] for l in range(len(rdg_support))]
2389 else:
2390 # Otherwise, the data is present, though it may be ambiguous; normalize the reading probabilities to sum to 1:
2391 rdg_support = [w / norm for l, w in enumerate(rdg_support)]
2392 normalized_reading_support_by_wit[wit] = rdg_support
2393 # Then add this witness's contributions to the readings' sampling probabilities:
2394 for l, w in enumerate(normalized_reading_support_by_wit[wit]):
2395 sampling_probabilities[l] += w
2396 norm = (
2397 sum(sampling_probabilities) if sum(sampling_probabilities) > 0 else 1.0
2398 ) # 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
2399 # Otherwise, normalize the sampling probabilities so they sum to 1:
2400 sampling_probabilities = [w / norm for w in sampling_probabilities]
2401 normalized_reading_support_dicts_by_vu_id[vu_id] = normalized_reading_support_by_wit
2402 sampling_probabilities_by_vu_id[vu_id] = sampling_probabilities
2403 # Then populate a contingency table for the expected probabilities of joint support in this unit:
2404 expected_joint_probabilities = np.full(
2405 (len(sampling_probabilities), len(sampling_probabilities)), 0, dtype=float
2406 )
2407 for l1 in range(len(sampling_probabilities)):
2408 for l2 in range(len(sampling_probabilities)):
2409 expected_joint_probabilities[l1, l2] = sampling_probabilities[l1] * sampling_probabilities[l2]
2410 expected_joint_probabilities_by_vu_id[vu_id] = expected_joint_probabilities
2411 # Then populate the matrix one variation unit at a time:
2412 matrix = np.full((len(witness_labels), len(witness_labels)), 0, dtype=float)
2413 with tqdm(total=len(self.witnesses) ** 2) as pbar:
2414 # Then calculate the mutual information contribution for each pair of witnesses:
2415 for i, wit_1 in enumerate(witness_labels):
2416 for j, wit_2 in enumerate(witness_labels):
2417 total_mutual_information = 0.0
2418 # The contribution to the entry for these witnesses will be identical regardless of the order in which they are specified,
2419 # so we only have to calculate it once:
2420 if i > j:
2421 pbar.update(1)
2422 continue
2423 # Otherwise, calculate the mutual information between these witnesses in each substantive variation unit:
2424 for k, vu_id in enumerate(self.variation_unit_ids):
2425 if vu_id not in substantive_variation_unit_ids_set:
2426 continue
2427 wit_1_rdg_support = normalized_reading_support_dicts_by_vu_id[vu_id][wit_1]
2428 wit_2_rdg_support = normalized_reading_support_dicts_by_vu_id[vu_id][wit_2]
2429 sampling_probabilities = sampling_probabilities_by_vu_id[vu_id]
2430 expected_joint_probabilities = expected_joint_probabilities_by_vu_id[vu_id]
2431 # If either witness has an all-zeroes vector (because it is lacunose in this unit), then we can skip these witnesses here:
2432 if sum(wit_1_rdg_support) == 0 or sum(wit_2_rdg_support) == 0:
2433 continue
2434 # Otherwise, populate a contingency table for the observed probabilities of joint support in this unit:
2435 observed_joint_probabilities = np.full(
2436 (len(sampling_probabilities), len(sampling_probabilities)), 0, dtype=float
2437 )
2438 for l1, w1 in enumerate(wit_1_rdg_support):
2439 for l2, w2 in enumerate(wit_2_rdg_support):
2440 observed_joint_probabilities[l1, l2] = w1 * w2
2441 # Then calculate the mutual information using the expected and observed distribution matrices:
2442 mutual_information = 0.0
2443 for l1 in range(len(sampling_probabilities)):
2444 for l2 in range(len(sampling_probabilities)):
2445 observed = observed_joint_probabilities[l1, l2]
2446 expected = expected_joint_probabilities[l1, l2]
2447 if observed == 0:
2448 continue
2449 mutual_information += observed * math.log2(observed / expected)
2450 # Then add this mutual information to the total for these two witnesses:
2451 total_mutual_information += mutual_information
2452 matrix[i, j] += total_mutual_information
2453 matrix[j, i] += total_mutual_information
2454 pbar.update(1)
2455 # Initialize a matrix for shared extant variation units for witnesses, and populate it if the proportion or show_ext option is specified:
2456 ext_matrix = None
2457 if proportion or show_ext:
2458 ext_matrix = self.get_ext_matrix(drop_constant=drop_constant, split_missing=split_missing)
2459 # 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:
2460 if proportion:
2461 proportion_matrix = np.full((len(witness_labels), len(witness_labels)), 0.0, dtype=float)
2462 np.divide(
2463 matrix, ext_matrix, out=proportion_matrix, where=(ext_matrix != 0)
2464 ) # division by 0 can occur if two witnesses have no overlapping units; leave their proportion as 0.0
2465 matrix = proportion_matrix
2466 # Then transform the columns of the main matrix as specified:
2467 matrix = self.transform_matrix(matrix, transform_matrix)
2468 # If the show_ext option is set, then append the number of shared extant variation units after the matrix's values:
2469 if show_ext:
2470 serialized_values = []
2471 for i, wit_1 in enumerate(witness_labels):
2472 serialized_values.append([])
2473 for j, wit_2 in enumerate(witness_labels):
2474 serialized_values[-1].append("/".join([str(matrix[i][j]), str(ext_matrix[i][j])]))
2475 matrix = np.array(serialized_values)
2476 return matrix, witness_labels
2478 def to_nexus_table(self, drop_constant: bool = False, ambiguous_as_missing: bool = False):
2479 """Returns this Collation in the form of a table with rows for taxa, columns for characters, and reading IDs in cells.
2481 Args:
2482 drop_constant (bool, optional): An optional flag indicating whether to ignore variation units with one substantive reading.
2483 Default value is False.
2484 ambiguous_as_missing (bool, optional): An optional flag indicating whether to treat all ambiguous states as missing data.
2485 Default value is False.
2487 Returns:
2488 A NumPy array with rows for taxa, columns for characters, and reading IDs in cells.
2489 A list of substantive reading ID strings.
2490 A list of witness ID strings.
2491 """
2492 # Populate a list of sites that will correspond to columns of the sequence alignment:
2493 substantive_variation_unit_ids = self.variation_unit_ids
2494 if drop_constant:
2495 substantive_variation_unit_ids = [
2496 vu_id
2497 for vu_id in self.variation_unit_ids
2498 if len(self.substantive_readings_by_variation_unit_id[vu_id]) > 1
2499 ]
2500 substantive_variation_unit_ids_set = set(substantive_variation_unit_ids)
2501 substantive_variation_unit_reading_tuples_set = set(self.substantive_variation_unit_reading_tuples)
2502 # In a first pass, populate a dictionary mapping (variation unit index, reading index) tuples from the readings_by_witness dictionary
2503 # to the readings' IDs:
2504 reading_ids_by_indices = {}
2505 for j, vu in enumerate(self.variation_units):
2506 if vu.id not in substantive_variation_unit_ids_set:
2507 continue
2508 k = 0
2509 for rdg in vu.readings:
2510 key = tuple([vu.id, rdg.id])
2511 if key not in substantive_variation_unit_reading_tuples_set:
2512 continue
2513 indices = tuple([j, k])
2514 reading_ids_by_indices[indices] = rdg.id
2515 k += 1
2516 # Initialize the output array with the appropriate dimensions:
2517 missing_symbol = '?'
2518 witness_labels = [wit.id for wit in self.witnesses]
2519 matrix = np.full(
2520 (len(witness_labels), len(substantive_variation_unit_ids)), missing_symbol, dtype=object
2521 ) # use dtype=object because the maximum string length is not known up front
2522 # Then populate it with the appropriate values:
2523 with tqdm(total=len(self.witnesses)) as pbar:
2524 row_ind = 0
2525 for i, wit in enumerate(self.witnesses):
2526 col_ind = 0
2527 for j, vu in enumerate(self.variation_units):
2528 if vu.id not in substantive_variation_unit_ids_set:
2529 continue
2530 rdg_support = self.readings_by_witness[wit.id][j]
2531 # If this reading support vector sums to 0, then this is missing data; handle it as specified:
2532 if sum(rdg_support) == 0:
2533 matrix[row_ind, col_ind] = missing_symbol
2534 # Otherwise, add its coefficients normally:
2535 else:
2536 rdg_inds = [
2537 k for k, w in enumerate(rdg_support) if w > 0
2538 ] # the index list consists of the indices of all readings with any degree of certainty assigned to them
2539 # For singleton readings, just print the reading ID:
2540 if len(rdg_inds) == 1:
2541 k = rdg_inds[0]
2542 matrix[row_ind, col_ind] = reading_ids_by_indices[(j, k)]
2543 # For multiple readings, print the corresponding reading IDs in braces or the missing symbol depending on input settings:
2544 else:
2545 if ambiguous_as_missing:
2546 matrix[row_ind, col_ind] = missing_symbol
2547 else:
2548 matrix[row_ind, col_ind] = "{%s}" % " ".join(
2549 [reading_ids_by_indices[(j, k)] for k in rdg_inds]
2550 )
2551 col_ind += 1
2552 row_ind += 1
2553 pbar.update(1)
2554 return matrix, witness_labels, substantive_variation_unit_ids
2556 def to_long_table(self, drop_constant: bool = False):
2557 """Returns this Collation in the form of a long table with columns for taxa, characters, reading indices, and reading values.
2558 Note that this method treats ambiguous readings as missing data.
2560 Args:
2561 drop_constant (bool, optional): An optional flag indicating whether to ignore variation units with one substantive reading.
2562 Default value is False.
2564 Returns:
2565 A NumPy array with columns for taxa, characters, reading indices, and reading values, and rows for each combination of these values in the matrix.
2566 A list of column label strings.
2567 """
2568 # Populate a list of sites that will correspond to columns of the sequence alignment:
2569 substantive_variation_unit_ids = self.variation_unit_ids
2570 if drop_constant:
2571 substantive_variation_unit_ids = [
2572 vu_id
2573 for vu_id in self.variation_unit_ids
2574 if len(self.substantive_readings_by_variation_unit_id[vu_id]) > 1
2575 ]
2576 substantive_variation_unit_ids_set = set(substantive_variation_unit_ids)
2577 substantive_variation_unit_reading_tuples_set = set(self.substantive_variation_unit_reading_tuples)
2578 # Initialize the outputs:
2579 column_labels = ["taxon", "character", "state", "value"]
2580 long_table_list = []
2581 # Populate a dictionary mapping (variation unit index, reading index) tuples to reading texts:
2582 reading_texts_by_indices = {}
2583 for j, vu in enumerate(self.variation_units):
2584 if vu.id not in substantive_variation_unit_ids_set:
2585 continue
2586 k = 0
2587 for rdg in vu.readings:
2588 key = tuple([vu.id, rdg.id])
2589 if key not in substantive_variation_unit_reading_tuples_set:
2590 continue
2591 indices = tuple([j, k])
2592 reading_texts_by_indices[indices] = rdg.text
2593 k += 1
2594 # Then populate the output list with the appropriate values:
2595 witness_labels = [wit.id for wit in self.witnesses]
2596 missing_symbol = '?'
2597 with tqdm(total=len(self.witnesses)) as pbar:
2598 for i, wit in enumerate(self.witnesses):
2599 row_ind = 0
2600 for j, vu_id in enumerate(self.variation_unit_ids):
2601 if vu_id not in substantive_variation_unit_ids_set:
2602 continue
2603 rdg_support = self.readings_by_witness[wit.id][j]
2604 # Populate a list of nonzero coefficients for this reading support vector:
2605 rdg_inds = [k for k, w in enumerate(rdg_support) if w > 0]
2606 # If this list does not consist of exactly one reading, then treat it as missing data:
2607 if len(rdg_inds) != 1:
2608 long_table_list.append([wit.id, vu_id, missing_symbol, missing_symbol])
2609 continue
2610 k = rdg_inds[0]
2611 rdg_text = reading_texts_by_indices[(j, k)]
2612 # Replace empty reading texts with the omission placeholder:
2613 if rdg_text == "":
2614 rdg_text = "om."
2615 long_table_list.append([wit.id, vu_id, k, rdg_text])
2616 pbar.update(1)
2617 # Then convert the long table entries list to a NumPy array:
2618 long_table = np.array(long_table_list)
2619 return long_table, column_labels
2621 def to_dataframe(
2622 self,
2623 drop_constant: bool = False,
2624 ambiguous_as_missing: bool = False,
2625 proportion: bool = False,
2626 table_type: TableType = TableType.matrix,
2627 split_missing: SplitMissingType = None,
2628 transform_matrix: TransformMatrixType = None,
2629 show_ext: bool = False,
2630 ):
2631 """Returns this Collation in the form of a Pandas DataFrame array, including the appropriate row and column labels.
2633 Args:
2634 drop_constant (bool, optional): An optional flag indicating whether to ignore variation units with one substantive reading.
2635 Default value is False.
2636 ambiguous_as_missing (bool, optional): An optional flag indicating whether to treat all ambiguous states as missing data.
2637 Default value is False.
2638 proportion (bool, optional): An optional flag indicating whether or not to calculate distances as proportions over extant, unambiguous variation units.
2639 Default value is False.
2640 table_type (TableType, optional): A TableType option indicating which type of tabular output to generate.
2641 Only applicable for tabular outputs.
2642 Default value is "matrix".
2643 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.
2644 If not specified, then missing data is ignored (i.e., all states are 0).
2645 If "uniform", then the contribution of 1 is divided evenly over all substantive readings.
2646 If "proportional", then the contribution of 1 is divided between the readings in proportion to their support among the witnesses that are not missing.
2647 Only applicable for table types "matrix" and "idf".
2648 transform_matrix (TransformMatrixType, optional): A TransformMatrixType option indicating how the columns of a witness-to-witness matrix output should be transformed.
2649 Only applicable for tabular outputs in which the rows and columns correspond to the witnesses in the collation.
2650 show_ext: An optional flag indicating whether each cell in the matrix
2651 should include the number of their extant, unambiguous variation units after the number of their disagreements/agreements.
2652 Only applicable for tabular output formats of type \"distance\" or \"similarity\".
2653 Default value is False.
2655 Returns:
2656 A Pandas DataFrame corresponding to a collation matrix with reading frequencies or a long table with discrete reading states.
2657 """
2658 df = None
2659 # Proceed based on the table type:
2660 if table_type == TableType.matrix:
2661 # Convert the collation to a NumPy array and get its row and column labels first:
2662 matrix, reading_labels, witness_labels = self.to_numpy(
2663 drop_constant=drop_constant, split_missing=split_missing
2664 )
2665 df = pd.DataFrame(matrix, index=reading_labels, columns=witness_labels)
2666 elif table_type == TableType.distance:
2667 # Convert the collation to a NumPy array and get its row and column labels first:
2668 matrix, witness_labels = self.to_distance_matrix(
2669 drop_constant=drop_constant, proportion=proportion, transform_matrix=transform_matrix, show_ext=show_ext
2670 )
2671 df = pd.DataFrame(matrix, index=witness_labels, columns=witness_labels)
2672 elif table_type == TableType.similarity:
2673 # Convert the collation to a NumPy array and get its row and column labels first:
2674 matrix, witness_labels = self.to_similarity_matrix(
2675 drop_constant=drop_constant, proportion=proportion, transform_matrix=transform_matrix, show_ext=show_ext
2676 )
2677 df = pd.DataFrame(matrix, index=witness_labels, columns=witness_labels)
2678 elif table_type == TableType.idf:
2679 # Convert the collation to a NumPy array and get its row and column labels first:
2680 matrix, witness_labels = self.to_idf_matrix(
2681 drop_constant=drop_constant,
2682 split_missing=split_missing,
2683 proportion=proportion,
2684 transform_matrix=transform_matrix,
2685 show_ext=show_ext,
2686 )
2687 df = pd.DataFrame(matrix, index=witness_labels, columns=witness_labels)
2688 elif table_type == TableType.mi:
2689 # Convert the collation to a NumPy array and get its row and column labels first:
2690 matrix, witness_labels = self.to_mi_matrix(
2691 drop_constant=drop_constant,
2692 split_missing=split_missing,
2693 proportion=proportion,
2694 transform_matrix=transform_matrix,
2695 show_ext=show_ext,
2696 )
2697 df = pd.DataFrame(matrix, index=witness_labels, columns=witness_labels)
2698 elif table_type == TableType.nexus:
2699 # Convert the collation to a NumPy array and get its row and column labels first:
2700 matrix, witness_labels, vu_labels = self.to_nexus_table(
2701 drop_constant=drop_constant, ambiguous_as_missing=ambiguous_as_missing
2702 )
2703 df = pd.DataFrame(matrix, index=witness_labels, columns=vu_labels)
2704 elif table_type == TableType.long:
2705 # Convert the collation to a long table and get its column labels first:
2706 long_table, column_labels = self.to_long_table(drop_constant=drop_constant)
2707 df = pd.DataFrame(long_table, columns=column_labels)
2708 return df
2710 def to_csv(
2711 self,
2712 file_addr: Union[Path, str],
2713 drop_constant: bool = False,
2714 ambiguous_as_missing: bool = False,
2715 proportion: bool = False,
2716 table_type: TableType = TableType.matrix,
2717 split_missing: SplitMissingType = None,
2718 transform_matrix: TransformMatrixType = None,
2719 show_ext: bool = False,
2720 **kwargs
2721 ):
2722 """Writes this Collation to a comma-separated value (CSV) file with the given address.
2724 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!
2726 Args:
2727 file_addr: A string representing the path to an output CSV file; the file type should be .csv.
2728 drop_constant: An optional flag indicating whether to ignore variation units with one substantive reading.
2729 Default value is False.
2730 ambiguous_as_missing: An optional flag indicating whether to treat all ambiguous states as missing data.
2731 Default value is False.
2732 proportion: An optional flag indicating whether or not to calculate distances as proportions over extant, unambiguous variation units.
2733 Default value is False.
2734 table_type: A TableType option indicating which type of tabular output to generate.
2735 Only applicable for tabular outputs.
2736 Default value is "matrix".
2737 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.
2738 If not specified, then missing data is ignored (i.e., all states are 0).
2739 If "uniform", then the contribution of 1 is divided evenly over all substantive readings.
2740 If "proportional", then the contribution of 1 is divided between the readings in proportion to their support among the witnesses that are not missing.
2741 Only applicable for table types "matrix" and "idf".
2742 transform_matrix: A TransformMatrixType option indicating how the columns of a witness-to-witness matrix output should be transformed.
2743 Only applicable for tabular outputs in which the rows and columns correspond to the witnesses in the collation.
2744 show_ext: An optional flag indicating whether each cell in the matrix
2745 should include the number of their extant, unambiguous variation units after the number of their disagreements/agreements.
2746 Only applicable for tabular output formats of type \"distance\" or \"similarity\".
2747 Default value is False.
2748 **kwargs: Keyword arguments for pandas.DataFrame.to_csv.
2749 """
2750 # Convert the collation to a Pandas DataFrame first:
2751 df = self.to_dataframe(
2752 drop_constant=drop_constant,
2753 ambiguous_as_missing=ambiguous_as_missing,
2754 proportion=proportion,
2755 table_type=table_type,
2756 split_missing=split_missing,
2757 show_ext=show_ext,
2758 transform_matrix=transform_matrix,
2759 )
2760 # Generate all parent folders for this file that don't already exist:
2761 Path(file_addr).parent.mkdir(parents=True, exist_ok=True)
2762 # Proceed based on the table type:
2763 if table_type == TableType.long:
2764 return df.to_csv(
2765 file_addr, encoding="utf-8-sig", index=False, **kwargs
2766 ) # add BOM to start of file so that Excel will know to read it as Unicode
2767 return df.to_csv(
2768 file_addr, encoding="utf-8-sig", **kwargs
2769 ) # add BOM to start of file so that Excel will know to read it as Unicode
2771 def to_excel(
2772 self,
2773 file_addr: Union[Path, str],
2774 drop_constant: bool = False,
2775 ambiguous_as_missing: bool = False,
2776 proportion: bool = False,
2777 table_type: TableType = TableType.matrix,
2778 split_missing: SplitMissingType = None,
2779 transform_matrix: TransformMatrixType = None,
2780 show_ext: bool = False,
2781 ):
2782 """Writes this Collation to an Excel (.xlsx) file with the given address.
2784 Since Pandas is deprecating its support for xlwt, specifying an output in old Excel (.xls) output is not recommended.
2786 Args:
2787 file_addr: A string representing the path to an output Excel file; the file type should be .xlsx.
2788 drop_constant: An optional flag indicating whether to ignore variation units with one substantive reading.
2789 Default value is False.
2790 ambiguous_as_missing: An optional flag indicating whether to treat all ambiguous states as missing data.
2791 Default value is False.
2792 proportion: An optional flag indicating whether or not to calculate distances as proportions over extant, unambiguous variation units.
2793 Default value is False.
2794 table_type: A TableType option indicating which type of tabular output to generate.
2795 Only applicable for tabular outputs.
2796 Default value is "matrix".
2797 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.
2798 If not specified, then missing data is ignored (i.e., all states are 0).
2799 If "uniform", then the contribution of 1 is divided evenly over all substantive readings.
2800 If "proportional", then the contribution of 1 is divided between the readings in proportion to their support among the witnesses that are not missing.
2801 Only applicable for table types "matrix" and "idf".
2802 transform_matrix: A TransformMatrixType option indicating how the columns of a witness-to-witness matrix output should be transformed.
2803 Only applicable for tabular outputs in which the rows and columns correspond to the witnesses in the collation.
2804 show_ext: An optional flag indicating whether each cell in the matrix
2805 should include the number of their extant, unambiguous variation units after the number of their disagreements/agreements.
2806 Only applicable for tabular output formats of type \"distance\" or \"similarity\".
2807 Default value is False.
2808 """
2809 # Convert the collation to a Pandas DataFrame first:
2810 df = self.to_dataframe(
2811 drop_constant=drop_constant,
2812 ambiguous_as_missing=ambiguous_as_missing,
2813 proportion=proportion,
2814 table_type=table_type,
2815 split_missing=split_missing,
2816 show_ext=show_ext,
2817 transform_matrix=transform_matrix,
2818 )
2819 # Generate all parent folders for this file that don't already exist:
2820 Path(file_addr).parent.mkdir(parents=True, exist_ok=True)
2821 # Proceed based on the table type:
2822 if table_type == TableType.long:
2823 return df.to_excel(file_addr, index=False)
2824 return df.to_excel(file_addr)
2826 def to_phylip_matrix(
2827 self,
2828 file_addr: Union[Path, str],
2829 drop_constant: bool = False,
2830 proportion: bool = False,
2831 table_type: TableType = TableType.distance,
2832 show_ext: bool = False,
2833 ):
2834 """Writes this Collation as a PHYLIP-formatted distance/similarity matrix to the file with the given address.
2836 Args:
2837 file_addr: A string representing the path to an output PHYLIP file; the file type should be .ph or .phy.
2838 drop_constant: An optional flag indicating whether to ignore variation units with one substantive reading.
2839 Default value is False.
2840 proportion: An optional flag indicating whether or not to calculate distances as proportions over extant, unambiguous variation units.
2841 Default value is False.
2842 table_type: A TableType option indicating which type of tabular output to generate.
2843 For PHYLIP-formatted outputs, distance and similarity matrices are the only supported table types.
2844 Default value is "distance".
2845 show_ext: An optional flag indicating whether each cell in the matrix
2846 should include the number of their extant, unambiguous variation units after the number of their disagreements/agreements.
2847 Only applicable for tabular output formats of type \"distance\" or \"similarity\".
2848 Default value is False.
2849 """
2850 # Convert the collation to a Pandas DataFrame first:
2851 matrix = None
2852 witness_labels = []
2853 # Proceed based on the table type:
2854 if table_type == TableType.distance:
2855 # Convert the collation to a NumPy array and get its row and column labels first:
2856 matrix, witness_labels = self.to_distance_matrix(
2857 drop_constant=drop_constant, proportion=proportion, show_ext=show_ext
2858 )
2859 elif table_type == TableType.similarity:
2860 # Convert the collation to a NumPy array and get its row and column labels first:
2861 matrix, witness_labels = self.to_similarity_matrix(
2862 drop_constant=drop_constant, proportion=proportion, show_ext=show_ext
2863 )
2864 # Generate all parent folders for this file that don't already exist:
2865 Path(file_addr).parent.mkdir(parents=True, exist_ok=True)
2866 with open(file_addr, "w", encoding="utf-8") as f:
2867 # The first line contains the number of taxa:
2868 f.write("%d\n" % len(witness_labels))
2869 # Every subsequent line contains a witness label, followed by the values in its row of the matrix:
2870 for i, wit_id in enumerate(witness_labels):
2871 wit_label = slugify(wit_id, lowercase=False, allow_unicode=True, separator='_')
2872 f.write("%s %s\n" % (wit_label, " ".join([str(v) for v in matrix[i]])))
2873 return
2875 def get_stemma_symbols(self):
2876 """Returns a list of one-character symbols needed to represent the states of all substantive readings in stemma format.
2878 The number of symbols equals the maximum number of substantive readings at any variation unit.
2880 Returns:
2881 A list of individual characters representing states in readings.
2882 """
2883 possible_symbols = (
2884 list(string.digits) + list(string.ascii_lowercase) + list(string.ascii_uppercase)
2885 ) # NOTE: the maximum number of symbols allowed in stemma format (other than "?" and "-") is 62
2886 # The number of symbols needed is equal to the length of the longest substantive reading vector:
2887 nsymbols = 0
2888 # If there are no witnesses, then no symbols are needed at all:
2889 if len(self.witnesses) == 0:
2890 return []
2891 wit_id = self.witnesses[0].id
2892 for rdg_support in self.readings_by_witness[wit_id]:
2893 nsymbols = max(nsymbols, len(rdg_support))
2894 stemma_symbols = possible_symbols[:nsymbols]
2895 return stemma_symbols
2897 def to_stemma(self, file_addr: Union[Path, str]):
2898 """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.
2900 Since this format does not support ambiguous states, all reading vectors with anything other than one nonzero entry will be interpreted as lacunose.
2901 If an interpGrp for weights is specified in the TEI XML collation, then the weights for the interp elements will be used as weights
2902 for the variation units that specify them in their ana attribute.
2904 Args:
2905 file_addr: A string representing the path to an output stemma prep file; the file should have no extension.
2906 The accompanying chron file will match this file name, except that it will have "_chron" appended to the end.
2907 drop_constant: An optional flag indicating whether to ignore variation units with one substantive reading.
2908 """
2909 # Populate a list of sites that will correspond to columns of the sequence alignment
2910 # (by default, constant sites are dropped):
2911 substantive_variation_unit_ids = [
2912 vu_id for vu_id in self.variation_unit_ids if len(self.substantive_readings_by_variation_unit_id[vu_id]) > 1
2913 ]
2914 substantive_variation_unit_ids_set = set(substantive_variation_unit_ids)
2915 substantive_variation_unit_reading_tuples_set = set(self.substantive_variation_unit_reading_tuples)
2916 # In a first pass, populate a dictionary mapping (variation unit index, reading index) tuples from the readings_by_witness dictionary
2917 # to the readings' texts:
2918 reading_texts_by_indices = {}
2919 for j, vu in enumerate(self.variation_units):
2920 if vu.id not in substantive_variation_unit_ids_set:
2921 continue
2922 k = 0
2923 for rdg in vu.readings:
2924 key = tuple([vu.id, rdg.id])
2925 if key not in substantive_variation_unit_reading_tuples_set:
2926 continue
2927 indices = tuple([j, k])
2928 reading_texts_by_indices[indices] = rdg.text
2929 k += 1
2930 # In a second pass, populate another dictionary mapping (variation unit index, reading index) tuples from the readings_by_witness dictionary
2931 # to the witnesses exclusively supporting those readings:
2932 reading_wits_by_indices = {}
2933 for indices in reading_texts_by_indices:
2934 reading_wits_by_indices[indices] = []
2935 for i, wit in enumerate(self.witnesses):
2936 for j, vu_id in enumerate(self.variation_unit_ids):
2937 if vu_id not in substantive_variation_unit_ids_set:
2938 continue
2939 rdg_support = self.readings_by_witness[wit.id][j]
2940 # If this witness does not exclusively support exactly one reading at this unit, then treat it as lacunose:
2941 if len([k for k, w in enumerate(rdg_support) if w > 0]) != 1:
2942 continue
2943 k = rdg_support.index(1)
2944 indices = tuple([j, k])
2945 reading_wits_by_indices[indices].append(wit.id)
2946 # In a third pass, write to the stemma file:
2947 symbols = self.get_stemma_symbols()
2948 Path(file_addr).parent.mkdir(
2949 parents=True, exist_ok=True
2950 ) # generate all parent folders for this file that don't already exist
2951 chron_file_addr = str(file_addr) + "_chron"
2952 with open(file_addr, "w", encoding="utf-8") as f:
2953 # Start with the witness list:
2954 f.write(
2955 "* %s ;\n\n"
2956 % " ".join(
2957 [slugify(wit.id, lowercase=False, allow_unicode=True, separator='_') for wit in self.witnesses]
2958 )
2959 )
2960 # f.write("^ %s\n\n" % chron_file_addr) #write the relative path to the chron file
2961 f.write(
2962 "^ %s\n\n" % ("." + os.sep + Path(chron_file_addr).name)
2963 ) # write the relative path to the chron file
2964 # Then add a line indicating that all witnesses are lacunose unless they are specified explicitly:
2965 f.write("= $? $* ;\n\n")
2966 with tqdm(total=len(self.variation_unit_ids)) as pbar:
2967 # Then proceed for each variation unit:
2968 for j, vu_id in enumerate(self.variation_unit_ids):
2969 if vu_id not in substantive_variation_unit_ids_set:
2970 pbar.update(1)
2971 continue
2972 # Print the variation unit ID first:
2973 f.write("@ %s\n" % vu_id)
2974 # In a first pass, print the texts of all readings enclosed in brackets:
2975 f.write("[ ")
2976 k = 0
2977 while True:
2978 indices = tuple([j, k])
2979 if indices not in reading_texts_by_indices:
2980 break
2981 text = slugify(
2982 reading_texts_by_indices[indices], lowercase=False, allow_unicode=True, separator='.'
2983 )
2984 # Denote omissions by en-dashes:
2985 if text == "":
2986 text = "\u2013"
2987 # The first reading should not be preceded by anything:
2988 if k == 0:
2989 f.write(text)
2990 f.write(" |")
2991 # Add the weight of this variation unit after the pipe by comparing its analysis categories to their weights:
2992 weight = 1
2993 vu = self.variation_units[j]
2994 if len(vu.analysis_categories) > 0:
2995 weight = int(
2996 sum(
2997 [
2998 self.weights_by_id[ana] if ana in self.weights_by_id else 1
2999 for ana in vu.analysis_categories
3000 ]
3001 )
3002 / len(vu.analysis_categories)
3003 )
3004 f.write("*%d" % weight)
3005 # Every subsequent reading should be preceded by a space:
3006 elif k > 0:
3007 f.write(" %s" % text)
3008 k += 1
3009 f.write(" ]\n")
3010 # In a second pass, print the indices and witnesses for all readings enclosed in angle brackets:
3011 k = 0
3012 f.write("\t< ")
3013 while True:
3014 indices = tuple([j, k])
3015 if indices not in reading_wits_by_indices:
3016 break
3017 rdg_symbol = symbols[k] # get the one-character alphanumeric code for this state
3018 wits = " ".join(reading_wits_by_indices[indices])
3019 # Open the variant reading support block with an angle bracket:
3020 if k == 0:
3021 f.write("%s %s" % (rdg_symbol, wits))
3022 # Open all subsequent variant reading support blocks with pipes on the next line:
3023 else:
3024 f.write("\n\t| %s %s" % (rdg_symbol, wits))
3025 k += 1
3026 f.write(" >\n")
3027 pbar.update(1)
3028 # In a fourth pass, write to the chron file:
3029 max_id_length = max(
3030 [len(slugify(wit.id, lowercase=False, allow_unicode=True, separator='_')) for wit in self.witnesses]
3031 )
3032 max_date_length = 0
3033 for wit in self.witnesses:
3034 if wit.date_range[0] is not None:
3035 max_date_length = max(max_date_length, len(str(wit.date_range[0])))
3036 if wit.date_range[1] is not None:
3037 max_date_length = max(max_date_length, len(str(wit.date_range[1])))
3038 # Attempt to get the minimum and maximum dates for witnesses; if we can't do this, then don't write a chron file:
3039 min_date = None
3040 max_date = None
3041 try:
3042 min_date = min([wit.date_range[0] for wit in self.witnesses if wit.date_range[0] is not None])
3043 max_date = max([wit.date_range[1] for wit in self.witnesses if wit.date_range[1] is not None])
3044 except Exception as e:
3045 print("WARNING: no witnesses have date ranges; no chron file will be written!")
3046 return
3047 with open(chron_file_addr, "w", encoding="utf-8") as f:
3048 for wit in self.witnesses:
3049 wit_label = slugify(wit.id, lowercase=False, allow_unicode=True, separator='_')
3050 f.write(wit_label)
3051 f.write(" " * (max_id_length - len(wit.id) + 1))
3052 # If either the lower bound on this witness's date is empty, then use the min and max dates over all witnesses as defaults:
3053 date_range = wit.date_range
3054 if date_range[0] is None:
3055 date_range = tuple([min_date, date_range[1]])
3056 # Then write the date range minimum, average, and maximum to the chron file:
3057 low_date = str(date_range[0])
3058 f.write(" " * (max_date_length - len(low_date) + 2))
3059 f.write(low_date)
3060 avg_date = str(int(((date_range[0] + date_range[1]) / 2)))
3061 f.write(" " * (max_date_length - len(str(avg_date)) + 2))
3062 f.write(avg_date)
3063 high_date = str(date_range[1])
3064 f.write(" " * (max_date_length - len(high_date) + 2))
3065 f.write(high_date)
3066 f.write("\n")
3067 return
3069 def to_file(
3070 self,
3071 file_addr: Union[Path, str],
3072 format: Format = None,
3073 drop_constant: bool = False,
3074 split_missing: SplitMissingType = None,
3075 char_state_labels: bool = True,
3076 frequency: bool = False,
3077 ambiguous_as_missing: bool = False,
3078 proportion: bool = False,
3079 calibrate_dates: bool = False,
3080 mrbayes: bool = False,
3081 clock_model: ClockModel = ClockModel.strict,
3082 ancestral_logger: AncestralLogger = AncestralLogger.state,
3083 table_type: TableType = TableType.matrix,
3084 transform_matrix: TransformMatrixType = None,
3085 show_ext: bool = False,
3086 seed: int = None,
3087 ):
3088 """Writes this Collation to the file with the given address.
3090 Args:
3091 file_addr (Union[Path, str]): The path to the output file.
3092 format (Format, optional): The desired output format.
3093 If None then it is infered from the file suffix.
3094 Defaults to None.
3095 drop_constant (bool, optional): An optional flag indicating whether to ignore variation units with one substantive reading.
3096 Default value is False.
3097 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.
3098 If not specified, then missing data is ignored (i.e., all states are 0).
3099 If "uniform", then the contribution of 1 is divided evenly over all substantive readings.
3100 If "proportional", then the contribution of 1 is divided between the readings in proportion to their support among the witnesses that are not missing.
3101 Only applicable for tabular outputs of type "matrix", "idf", "mean-idf", "mi", and "mean-mi".
3102 char_state_labels (bool, optional): An optional flag indicating whether to print
3103 the CharStateLabels block in NEXUS output.
3104 Default value is True.
3105 frequency (bool, optional): An optional flag indicating whether to use the StatesFormat=Frequency setting
3106 instead of the StatesFormat=StatesPresent setting
3107 (and thus represent all states with frequency vectors rather than symbols)
3108 in NEXUS output.
3109 Note that this setting is necessary to make use of certainty degrees assigned to multiple ambiguous states in the collation.
3110 Default value is False.
3111 ambiguous_as_missing (bool, optional): An optional flag indicating whether to treat all ambiguous states as missing data.
3112 If this flag is set, then only base symbols will be generated for the NEXUS file.
3113 It is only applied if the frequency option is False.
3114 Default value is False.
3115 proportion (bool, optional): An optional flag indicating whether to populate a distance matrix's cells
3116 with a proportion of disagreements to variation units where both witnesses are extant.
3117 It is only applied if the table_type option is "distance".
3118 Default value is False.
3119 calibrate_dates (bool, optional): An optional flag indicating whether to add an Assumptions block that specifies date distributions for witnesses
3120 in NEXUS output.
3121 This option is intended for inputs to BEAST 2.
3122 Default value is False.
3123 mrbayes (bool, optional): An optional flag indicating whether to add a MrBayes block that specifies model settings and age calibrations for witnesses
3124 in NEXUS output.
3125 This option is intended for inputs to MrBayes.
3126 Default value is False.
3127 clock_model (ClockModel, optional): A ClockModel option indicating which type of clock model to use.
3128 This option is intended for inputs to MrBayes and BEAST 2.
3129 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.
3130 Default value is "strict".
3131 ancestral_logger (AncestralLogger, optional): An AncestralLogger option indicating which class of logger (if any) to use for ancestral states.
3132 This option is intended for inputs to BEAST 2.
3133 table_type (TableType, optional): A TableType option indicating which type of tabular output to generate.
3134 Only applicable for tabular outputs and PHYLIP outputs.
3135 If the output is a PHYLIP file, then the type of tabular output must be "distance" or "similarity"; otherwise, it will be ignored.
3136 Default value is "matrix".
3137 transform_matrix (TransformMatrixType, optional): A TransformMatrixType option indicating how the columns of a witness-to-witness matrix output should be transformed.
3138 Only applicable for tabular outputs in which the rows and columns correspond to the witnesses in the collation.
3139 show_ext (bool, optional): An optional flag indicating whether each cell in the matrix
3140 should include the number of variation units where both witnesses are extant after the number of their disagreements/agreements.
3141 Only applicable for tabular output formats of type "distance" or "similarity".
3142 Default value is False.
3143 seed (optional, int): A seed for random number generation (for setting initial values of unspecified transcriptional rates in BEAST 2 XML output).
3144 """
3145 file_addr = Path(file_addr)
3146 format = format or Format.infer(
3147 file_addr.suffix
3148 ) # an exception will be raised here if the format or suffix is invalid
3150 if format == Format.NEXUS:
3151 return self.to_nexus(
3152 file_addr,
3153 drop_constant=drop_constant,
3154 char_state_labels=char_state_labels,
3155 frequency=frequency,
3156 ambiguous_as_missing=ambiguous_as_missing,
3157 calibrate_dates=calibrate_dates,
3158 mrbayes=mrbayes,
3159 clock_model=clock_model,
3160 )
3162 if format == format.HENNIG86:
3163 return self.to_hennig86(file_addr, drop_constant=drop_constant)
3165 if format == format.PHYLIP:
3166 if table_type in [TableType.distance, TableType.similarity]:
3167 return self.to_phylip_matrix(
3168 file_addr,
3169 drop_constant=drop_constant,
3170 proportion=proportion,
3171 table_type=table_type,
3172 show_ext=show_ext,
3173 )
3174 return self.to_phylip(file_addr, drop_constant=drop_constant)
3176 if format == format.FASTA:
3177 return self.to_fasta(file_addr, drop_constant=drop_constant)
3179 if format == format.BEAST:
3180 return self.to_beast(
3181 file_addr,
3182 drop_constant=drop_constant,
3183 clock_model=clock_model,
3184 ancestral_logger=ancestral_logger,
3185 seed=seed,
3186 )
3188 if format == Format.CSV:
3189 return self.to_csv(
3190 file_addr,
3191 drop_constant=drop_constant,
3192 ambiguous_as_missing=ambiguous_as_missing,
3193 proportion=proportion,
3194 table_type=table_type,
3195 split_missing=split_missing,
3196 transform_matrix=transform_matrix,
3197 show_ext=show_ext,
3198 )
3200 if format == Format.TSV:
3201 return self.to_csv(
3202 file_addr,
3203 drop_constant=drop_constant,
3204 ambiguous_as_missing=ambiguous_as_missing,
3205 proportion=proportion,
3206 table_type=table_type,
3207 split_missing=split_missing,
3208 transform_matrix=transform_matrix,
3209 show_ext=show_ext,
3210 sep="\t",
3211 )
3213 if format == Format.EXCEL:
3214 return self.to_excel(
3215 file_addr,
3216 drop_constant=drop_constant,
3217 ambiguous_as_missing=ambiguous_as_missing,
3218 proportion=proportion,
3219 table_type=table_type,
3220 split_missing=split_missing,
3221 transform_matrix=transform_matrix,
3222 show_ext=show_ext,
3223 )
3225 if format == Format.STEMMA:
3226 return self.to_stemma(file_addr)