source_modelling.rupture_propagation

  1"""
  2Rupture Propagation Module
  3
  4This module provides functions for computing likely rupture paths from
  5information about the distances between faults.
  6
  7Functions:
  8    - shaw_dieterich_distance_model:
  9      Compute fault jump probabilities using the Shaw-Dieterich distance model.
 10    - prune_distance_graph: Prune the distance graph based on a cutoff value.
 11    - probability_graph:
 12      Convert a graph of distances between faults into a graph of jump
 13      probabilities using the Shaw-Dieterich model.
 14    - probabilistic_minimum_spanning_tree: Generate a probabilistic minimum spanning tree.
 15
 16Typing Aliases:
 17    - DistanceGraph: A graph representing distances between faults.
 18    - ProbabilityGraph: A graph representing jump probabilities between faults.
 19    - RuptureCausalityTree: A tree representing the causality of ruptures between faults.
 20"""
 21
 22from collections import defaultdict, namedtuple
 23from typing import Generator, Optional
 24
 25import numpy as np
 26from qcore import coordinates
 27
 28from source_modelling import sources
 29
 30DistanceGraph = dict[str, dict[str, int]]
 31ProbabilityGraph = defaultdict[str, dict[str, float]]
 32RuptureCausalityTree = dict[str, Optional[str]]
 33
 34JumpPair = namedtuple("JumpPair", ["from_point", "to_point"])
 35
 36
 37def shaw_dieterich_distance_model(distance: float, d0: float, delta: float) -> float:
 38    """
 39    Compute fault jump probabilities using the Shaw-Dieterich distance model[0].
 40
 41    Parameters
 42    ----------
 43    distance : float
 44        The distance between two faults.
 45    d0 : float
 46        The characteristic distance parameter.
 47    delta : float
 48        The characteristic slip distance parameter.
 49
 50    Returns
 51    -------
 52    float
 53        The calculated probability.
 54
 55    References
 56    ----------
 57    [0]: Shaw, B. E., & Dieterich, J. H. (2007). Probabilities for jumping fault
 58         segment stepovers. Geophysical Research Letters, 34(1).
 59    """
 60    return min(1, np.exp(-(distance - delta) / d0))
 61
 62
 63def prune_distance_graph(distances: DistanceGraph, cutoff: int) -> DistanceGraph:
 64    """
 65    Prune the distance graph based on a cutoff value.
 66
 67    Parameters
 68    ----------
 69    distances : DistanceGraph
 70        The graph of distances between faults.
 71    cutoff : int
 72        The cutoff distance in metres.
 73
 74    Returns
 75    -------
 76    DistanceGraph
 77        A copy of the input distance graph, keeping only edges that are less
 78        than the cutoff.
 79    """
 80    return {
 81        fault_u: {
 82            fault_v: distance
 83            for fault_v, distance in neighbours_fault_u.items()
 84            if distance < cutoff
 85        }
 86        for fault_u, neighbours_fault_u in distances.items()
 87    }
 88
 89
 90def probability_graph(
 91    distances: DistanceGraph, d0: float = 3, delta: float = 1
 92) -> ProbabilityGraph:
 93    """
 94    Convert a graph of distances between faults into a graph of jump
 95    probablities using the Shaw-Dieterich model.
 96
 97    Parameters
 98    ----------
 99    distances : DistanceGraph
100        The distance graph between faults.
101    d0 : float, optional
102        The d0 parameter for the Shaw_Dieterich model. See `shaw_dieterich_distance_model`.
103    delta : float, optional
104        The delta parameter for the Shaw_Dieterich model. See `shaw_dieterich_distance_model`.
105
106    Returns
107    -------
108    ProbabilityGraph
109        The graph with faults as vertices. Each edge (fault_u, fault_v)
110        has a log-probability -p as a weight. The log-probability -p here
111        is the negative of the log-probability a rupture propogates from
112        fault_u to fault_v, relative to the probability it propogates to
113        any of the other neighbours of fault_u.
114    """
115    probabilities_raw = {
116        fault_u: {
117            fault_v: shaw_dieterich_distance_model(distance / 1000, d0, delta)
118            for fault_v, distance in neighbours_fault_u.items()
119        }
120        for fault_u, neighbours_fault_u in distances.items()
121    }
122    probabilities_log = defaultdict(dict)
123    for fault_u, neighbours_fault_u in probabilities_raw.items():
124        normalising_constant = sum(neighbours_fault_u.values())
125        if normalising_constant == 0:
126            for fault_v, _ in neighbours_fault_u.items():
127                probabilities_log[fault_u][fault_v] = -np.log(
128                    1 / len(neighbours_fault_u)
129                )
130        for fault_v, prob in neighbours_fault_u.items():
131            probabilities_log[fault_u][fault_v] = -np.log(prob / normalising_constant)
132    return probabilities_log
133
134
135def probabilistic_minimum_spanning_tree(
136    probability_graph: ProbabilityGraph, initial_fault: str
137) -> RuptureCausalityTree:
138    """
139    Generate a probabilistic minimum spanning tree.
140
141    The minimum spanning tree of the probability graph represents rupture
142    causality tree with the highest likelihood of occuring assuming that
143    rupture jumps are independent.
144
145    NOTE: While the overall probability of the minimum spanning tree is high,
146    the paths from the initial fault may not be the most likely.
147
148    Parameters
149    ----------
150    probability_graph : ProbabilityGraph
151        The probability graph.
152    initial_fault : str
153        The initial fault.
154
155    Returns
156    -------
157    RuptureCausalityTree
158        The probabilistic minimum spanning tree.
159    """
160    rupture_causality_tree = {initial_fault: None}
161    path_probabilities = defaultdict(lambda: np.inf)
162    path_probabilities[initial_fault] = 0
163    processed_faults = set()
164    for _ in range(len(probability_graph)):
165        current_fault = min(
166            (fault for fault in probability_graph if fault not in processed_faults),
167            key=lambda fault: path_probabilities[fault],
168        )
169        processed_faults.add(current_fault)
170        for fault_neighbour, probability in probability_graph[current_fault].items():
171            if (
172                fault_neighbour not in processed_faults
173                and probability < path_probabilities[fault_neighbour]
174            ):
175                path_probabilities[fault_neighbour] = probability
176                rupture_causality_tree[fault_neighbour] = current_fault
177    return rupture_causality_tree
178
179
180def distance_between(
181    source_a: sources.IsSource,
182    source_b: sources.IsSource,
183    source_a_point: np.ndarray,
184    source_b_point: np.ndarray,
185) -> float:
186    global_point_a = source_a.fault_coordinates_to_wgs_depth_coordinates(source_a_point)
187    global_point_b = source_b.fault_coordinates_to_wgs_depth_coordinates(source_b_point)
188    return coordinates.distance_between_wgs_depth_coordinates(
189        global_point_a, global_point_b
190    )
191
192
193def estimate_most_likely_rupture_propagation(
194    sources_map: dict[str, sources.IsSource],
195    initial_source: str,
196    jump_impossibility_limit_distance: int = 15000,
197) -> RuptureCausalityTree:
198    distance_graph = {
199        source_a_name: {
200            source_b_name: distance_between(
201                source_a,
202                source_b,
203                *sources.closest_point_between_sources(source_a, source_b),
204            )
205            for source_b_name, source_b in sources_map.items()
206            if source_a_name != source_b_name
207        }
208        for source_a_name, source_a in sources_map.items()
209    }
210    distance_graph = prune_distance_graph(
211        distance_graph, jump_impossibility_limit_distance
212    )
213    jump_probability_graph = probability_graph(distance_graph)
214    return probabilistic_minimum_spanning_tree(
215        jump_probability_graph, initial_fault=initial_source
216    )
217
218
219def jump_points_from_rupture_tree(
220    source_map: dict[str, sources.IsSource],
221    rupture_causality_tree: RuptureCausalityTree,
222) -> dict[str, JumpPair]:
223    jump_points = {}
224    for source, parent in rupture_causality_tree.items():
225        if parent is None:
226            continue
227        source_point, parent_point = sources.closest_point_between_sources(
228            source_map[source], source_map[parent]
229        )
230        jump_points[source] = JumpPair(parent_point, source_point)
231    return jump_points
232
233
234def tree_nodes_in_order(
235    tree: dict[str, str],
236) -> Generator[str, None, None]:
237    """Generate faults in topologically sorted order.
238
239    Parameters
240    ----------
241    faults : list[RealisationFault]
242        List of RealisationFault objects.
243
244    Yields
245    ------
246    RealisationFault
247        The next fault in the topologically sorted order.
248    """
249    tree_child_map = defaultdict(list)
250    for cur, parent in tree.items():
251        if parent:
252            tree_child_map[parent].append(cur)
253
254    def in_order_traversal(
255        node: str,
256    ) -> Generator[str, None, None]:
257        yield node
258        for child in tree_child_map[node]:
259            yield from in_order_traversal(child)
260
261    initial_fault = next(cur for cur, parent in tree.items() if not parent)
262    yield from in_order_traversal(initial_fault)