From e3741533770c98bec2936edc44b892e8c4d63ed9 Mon Sep 17 00:00:00 2001 From: binho Date: Tue, 20 Feb 2024 12:29:19 -0600 Subject: [PATCH] Start event series code/example/test --- examples/test_event_series.pj | 2 + .../trees_maps_files/geosse_dummy_tree2.tre | 1 + .../geosse_dummy_tree2_maps.tsv | 18 ++ .../geosse_dummy_tree2_tip_states.tsv | 5 + src/phylojunction/functionality/stoch_map.py | 157 ++++++++++++++++-- 5 files changed, 167 insertions(+), 16 deletions(-) create mode 100644 examples/test_event_series.pj create mode 100644 examples/trees_maps_files/geosse_dummy_tree2.tre create mode 100644 examples/trees_maps_files/geosse_dummy_tree2_maps.tsv create mode 100644 examples/trees_maps_files/geosse_dummy_tree2_tip_states.tsv diff --git a/examples/test_event_series.pj b/examples/test_event_series.pj new file mode 100644 index 0000000..07a6a32 --- /dev/null +++ b/examples/test_event_series.pj @@ -0,0 +1,2 @@ +trs <- read_tree(n=1, file_path="examples/trees_maps_files/geosse_dummy_tree2.tre", node_name_attr="index") +mapped_trs := map_attr(tree=trs, fun="smap", maps_file_path="examples/trees_maps_files/geosse_dummy_tree2_maps.tsv", tip_attr_file_path="examples/trees_maps_files/geosse_dummy_tree2_tip_states.tsv", attr_name="state", n_regions=2, geo="true") \ No newline at end of file diff --git a/examples/trees_maps_files/geosse_dummy_tree2.tre b/examples/trees_maps_files/geosse_dummy_tree2.tre new file mode 100644 index 0000000..1416043 --- /dev/null +++ b/examples/trees_maps_files/geosse_dummy_tree2.tre @@ -0,0 +1 @@ +((A:[&index=1]4.0,B:[&index=2]4.0):[&index=6]1.0,((C:[&index=3]1.0,D:[&index=4]1.0):[&index=7]2.0,E:[&index=5]3.0):[&index=8]2.0)[&index=9]; \ No newline at end of file diff --git a/examples/trees_maps_files/geosse_dummy_tree2_maps.tsv b/examples/trees_maps_files/geosse_dummy_tree2_maps.tsv new file mode 100644 index 0000000..775e256 --- /dev/null +++ b/examples/trees_maps_files/geosse_dummy_tree2_maps.tsv @@ -0,0 +1,18 @@ +iteration node_index branch_start_time branch_end_time start_state end_state transition_time transition_type parent_index child1_index child2_index +1 1 4 0 0 0 NA no_change 6 NA NA +1 2 4 0 1 1 NA no_change 6 NA NA +1 3 1 0 0 0 NA no_change 7 NA NA +1 4 1 0 1 1 NA no_change 7 NA NA +1 5 3 0 0 0 NA no_change 8 NA NA +1 6 5 4 0 2 4.75 anagenetic 9 1 2 +1 6 4 4 2 2 NA no_change 9 1 2 +1 6 4 4 2 0 4 cladogenetic 9 1 2 +1 6 4 4 2 1 4 cladogenetic 9 1 2 +1 7 3 1 1 2 2.75 anagenetic 8 3 4 +1 7 1 1 2 2 NA no_change 8 3 4 +1 7 1 1 2 0 1 cladogenetic 8 3 4 +1 7 1 1 2 1 1 cladogenetic 8 3 4 +1 8 5 3 0 2 3.25 anagenetic 9 5 7 +1 8 3 3 2 2 NA no_change 9 5 7 +1 8 3 3 2 0 4 cladogenetic 9 5 7 +1 8 3 3 2 1 4 cladogenetic 9 5 7 diff --git a/examples/trees_maps_files/geosse_dummy_tree2_tip_states.tsv b/examples/trees_maps_files/geosse_dummy_tree2_tip_states.tsv new file mode 100644 index 0000000..3539a62 --- /dev/null +++ b/examples/trees_maps_files/geosse_dummy_tree2_tip_states.tsv @@ -0,0 +1,5 @@ +nd1 0 +nd2 1 +nd3 0 +nd4 1 +nd5 0 \ No newline at end of file diff --git a/src/phylojunction/functionality/stoch_map.py b/src/phylojunction/functionality/stoch_map.py index ab36cc3..4de3703 100644 --- a/src/phylojunction/functionality/stoch_map.py +++ b/src/phylojunction/functionality/stoch_map.py @@ -1,12 +1,9 @@ -import sys import os import copy import enum -import math -import numpy as np import typing as ty -import matplotlib import matplotlib.pyplot as plt # type: ignore +from abc import ABC, abstractmethod # pj imports import phylojunction.data.tree as pjt @@ -20,13 +17,138 @@ __email__ = "f.mendes@wustl.edu" -class StochMap(): +class Hypothesis(enum.Enum): + VICARIANCE = 0 + FOUNDER_EVENT = 1 + + +class EvolRelevantEvent(ABC): + """Evolution-relevant event. + + Examples of evolution-relevant events are those represented by + stochastic maps (e.g., range contraction, range expansion, + extinction, dispersal) or the (paleogeographic) appearance and + disappearance of barriers. + + Parameters: + n_chars (int): Number of characters involved in event. + age (float): Age of event (age of present moment is 0.0). + time (float, optional): Time of event (0.0 at origin or + root). Defaults to None. + """ + + _n_chars: int + _age: float + _time: float + + @abstractmethod + def __init__(self, + n_chars: int, + age: ty.Optional[float], + time: ty.Optional[float] = None) -> None: + + self._n_chars = n_chars + self._age = age + self._time = time + + + @property + def age(self) -> float: + return self._age + + @property + def time(self) -> float: + return self._time + + +class EvolRelevantEventSeries: + """Series of evolution-relevant events. + + This class has member methods that interrogate a series of event, + truncating it depending on what the user specifies. + + """ + + _event_list: ty.List[EvolRelevantEvent] + _n_events: int + _series_type: str + _trunc_event_list: ty.List[EvolRelevantEvent] + _trunc_n_events: int + _supported_hyp: Hypothesis + + def __init__(self, + event_list: ty.List[EvolRelevantEvent], + series_type: str) -> None: + + self._series_type = series_type + + self.health_check() + + self._event_list = event_list + self._n_events = len(event_list) + + def health_check(self) -> None: + """Confirm validity of event list given series type.""" + + if self._series_type == "speciation": + last_event = self._event_list[-1] + + if not isinstance(last_event, RangeSplitOrBirth): + raise ec.SequenceInvalidLastError(type(last_event)) + + def truncate_upstream(self, truncator_fn: ty.Callable) -> None: + """Deep-copy series and truncate copy upstream. + + Truncation of an event series is done according to some + user-specified function. Truncation may be carried out, for + example, at the first event that destabilizes a range (i.e., + the first "destabilizer" event). + + Side-effects: + (i) Populate self._trunc_event_list + (ii) Populate self._trunc_n_events + + Parameters: + truncator_fn (function): Method of static Truncate class + that executes truncation. + """ + + self._trunc_event_list = truncator_fn( + copy.deepcopy(self._event_list)) + + self._trunc_n_events = len(self._trunc_event_list) + + + @property + def event_list(self) -> ty.List[EvolRelevantEvent]: + return self._event_list + + @property + def n_events(self) -> int: + return self._n_events + + @property + def series_type(self) -> str: + return self._series_type + + @property + def supported_hyp(self) -> Hypothesis: + return self._supported_hyp + + @supported_hyp.setter + def supported_hyp(self, a_hyp: Hypothesis) -> None: + if (a_hyp == Hypothesis.VICARIANCE or \ + a_hyp == Hypothesis.FOUNDER_EVENT) and \ + self._series_type != "speciation": + raise ec.SequenceHypothesisTypeError(a_hyp, self._series_type) + + self._supported_hyp = a_hyp + + +class StochMap(EvolRelevantEvent): """Parent class for a single stochastic map. Parameters: - n_chars (int): Number of characters (e.g., atomic regions). - age (float): Age of stochastic map (0.0 at present). - time (float): Time of stochastic map (0.0 at origin or root). from_state (int): Integer representation of source state (i.e., source biogeographic range). from_state_bit_patt (str): Bit pattern representation of @@ -51,9 +173,6 @@ class StochMap(): focal node. """ - n_chars: int - age: float - time: float from_state: int from_state_bit_patt: str to_state: int @@ -79,9 +198,10 @@ def __init__(self, to_state2: ty.Optional[int] = None, to_state2_bit_patt: ty.Optional[str] = None, child2_node_name: ty.Optional[str] = None) -> None: - - self.age = age - self.time = time + + # initializing EvolRelevantEvent members + super().__init__(n_chars, age, time=time) + self.from_state = from_state self.from_state_bit_patt = from_state_bit_patt self.to_state = to_state @@ -89,7 +209,6 @@ def __init__(self, self.focal_node_name = focal_node_name self.parent_node_name = parent_node_name self.child_node_name = child_node_name - self.n_chars = n_chars # if cladogenetic range split self.child2_node_name = child2_node_name @@ -1274,7 +1393,13 @@ def iter_mcmc_generations(self) -> ty.Iterator[StochMapsOnTree]: yield self.stoch_maps_tree_dict[idx] +################################## +# Event series parsing functions # +################################## +def truncate_at_split_relevant_event(event_series: EvolRelevantEventSeries): + def _is_split_relevant_event(): + pass if __name__ == "__main__": @@ -1325,7 +1450,7 @@ def iter_mcmc_generations(self) -> ty.Iterator[StochMapsOnTree]: node_states_file_path="examples/trees_maps_files/geosse_dummy_tree1_tip_states.tsv", stoch_map_attr_name="state") - fig = matplotlib.pyplot.figure() + fig = plt.figure() ax = fig.add_axes([0.25, 0.2, 0.5, 0.6]) ax.patch.set_alpha(0.0)