aboutsummaryrefslogtreecommitdiffstats
path: root/fg21sim/extragalactic/clusters/mergertree.py
blob: 1df7b2fb57b6227eaa3a66e0b5a0fa0b4c1c84e6 (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
# Copyright (c) 2017-2019 Weitian LI <wt@liwt.net>
# MIT License

"""
Merger tree that represents the merging history of a cluster using
the binary tree data structure.
"""

import operator as op
import logging

from ...utils.io import pickle_dump, pickle_load

logger = logging.getLogger(__name__)


class MergerTree:
    """
    A binary tree that represents the cluster merging history.

    Description
    -----------
                   Merged (M0, z0, age0)
                   ~~~~~~~~~~~~~~~~~~~~~
                   |                   |
        Main (M1, z1, age1)    Sub (M2, z2, age2)
        ~~~~~~~~~~~~~~~~~~~    ~~~~~~~~~~~~~~~~~~

    * "Merged" is the merged cluster from "Main" and "Sub" at redshift
      z1 (=z2) or cosmic time age1 (=age2).
      M0 = M1 + M2, M1 > M2
    * If "Sub" is missing, then this is an accretion event (not a merger).

    Parameters
    ----------
    data : any (e.g., a dictionary)
        Data (e.g., mass, redshift, age) associated with this tree node.
    main, sub : `~MergerTree`
        Links to the main and sub (optional) clusters between which the
        merger happens.
        The ``sub`` cluster may be missing, which is regarded as an
        accretion event rather than a merger.
    """
    def __init__(self, data, main=None, sub=None):
        self.data = data
        self.main = main
        self.sub = sub

    def itermain(self):
        """
        Iterate by tracing the main cluster.
        """
        maintree = self
        subtree = None
        while maintree:
            main = maintree.data
            if subtree is not None:
                sub = subtree.data
            else:
                sub = {}
            yield (main, sub)
            subtree = maintree.sub
            maintree = maintree.main

    @property
    def lmain(self):
        """
        Return the length (i.e., number of events) along the main cluster.
        """
        return len(list(self.itermain()))

    def imain(self, idx):
        """
        Trace the main cluster to locate the idx-th (0-based) event, and
        return its data (both the main and sub clusters).

        Parameters
        ----------
        idx : int
            The event index (0-based) along the main cluster to be retrieved.

        Returns
        -------
        (main, sub) : dict
            The data dictionaries of the main and sub clusters.  The ``sub``
            dictionary may be ``None`` if it is an accretion event.
        """
        for main, sub in self.itermain():
            if idx == 0:
                return (main, sub)
            idx -= 1
        raise IndexError("index out of range: %d" % idx)


def save_mtree(mtree, outfile, clobber=False):
    """
    Pickle the merger tree data and save to file.
    """
    pickle_dump(mtree, outfile=outfile, clobber=clobber)
    logger.info("Saved merger tree to file: {0}".format(outfile))


def read_mtree(infile):
    mtree = pickle_load(infile)
    logger.info("Loaded merger tree from file: {0}".format(infile))
    return mtree


def get_history(mtree, merger_only=False):
    """
    Extract all the formation history (merger and accretion events).

    Parameters
    ----------
    mtree : `~MergerTree`
        The simulated merger tree from which to extract the history.
        Default: ``self.mtree``
    merger_only : bool, optional
        If ``True``, only extract the merger events.

    Returns
    -------
    evlist : list[event]
        List of events with each element being a dictionary of the
        event properties.
    """
    evlist = []
    for main, sub in mtree.itermain():
        z, age, M_main = op.itemgetter("z", "age", "mass")(main)
        if sub:
            # merger
            M_sub = sub["mass"]
            R_mass = M_main / M_sub
        else:
            # accretion
            if merger_only:
                continue
            M_sub, R_mass = None, None

        evlist.append({
            "z": z,
            "age": age,
            "M_main": M_main,
            "M_sub": M_sub,
            "R_mass": R_mass,
        })

    return evlist


def show_mtree(mtree):
    """
    Trace the main cluster and show its formation history.
    """
    parent = None
    for i, (main, sub) in enumerate(mtree.itermain()):
        info = "%2d: " % i
        z = main["z"]
        age = main["age"]
        info += "<z=%.3f; t=%5.2f> " % (z, age)
        mass = main["mass"]
        if sub:
            # merger event
            Msub = sub["mass"]
            Rmass = mass / Msub
            info += "[%.3e @@@ %.3e] (Rm=%5.1f)" % (mass, Msub, Rmass)
        elif parent:
            # accretion event
            dM = parent["mass"] - mass
            info += " %.3e  +  %.3e            " % (mass, dM)
        else:
            # root cluster
            info += " %.3e" % mass
        if parent:
            info += " <dz=%.3f/dt=%.2f>" % (z-parent["z"], parent["age"]-age)
        parent = main
        print(info)


def plot_mtree(mtree, outfile, figsize=(12, 8)):
    """
    Plot the cluster merger tree.

    XXX: Need to speed up this function.

    Parameters
    ----------
    mtree : `~MergerTree`
        The merger tree to be plotted
    outfile : str
        Output filename to save the plotted figure
    figsize : tuple
        The (width, height) of the plotting figure
    """
    def _plot(tree, ax):
        if tree is None:
            return
        if tree.main is None:
            # Only plot a point for current tree node
            x = [tree.data["age"]]
            y = [tree.data["mass"]]
            ax.plot(x, y, marker="o", markersize=1.5, color="black",
                    linestyle=None)
            return

        # Plot a point for current tree node
        x = [tree.data["age"]]
        y = [tree.data["mass"]]
        ax.plot(x, y, marker="o", markersize=1.5, color="black",
                linestyle=None)
        # Plot a line from current tree node to its main node
        x = [tree.data["age"], tree.main.data["age"]]
        y = [tree.data["mass"], tree.main.data["mass"]]
        ax.plot(x, y, color="blue")
        if tree.sub:
            # Plot a line between main and sub nodes
            x = [tree.main.data["age"], tree.sub.data["age"]]
            y = [tree.main.data["mass"], tree.sub.data["mass"]]
            ax.plot(x, y, color="green", linewidth=1, alpha=0.8)

        # Recursively plot the descendant nodes
        _plot(tree.main, ax)
        _plot(tree.sub, ax)

    from matplotlib.figure import Figure
    from matplotlib.backends.backend_agg import FigureCanvasAgg as FigureCanvas

    fig = Figure(figsize=figsize)
    canvas = FigureCanvas(fig)
    ax = fig.add_subplot(1, 1, 1)
    print("Plotting the merger tree, may take a while ...")
    _plot(mtree, ax=ax)
    ax.set_xlabel("Cosmic time [Gyr]")
    ax.set_ylabel("Mass [Msun]")
    ax.set_xlim((0, mtree.data["age"]))
    ax.set_ylim((0, mtree.data["mass"]))
    fig.tight_layout()
    canvas.print_figure(outfile)
    print("Saved plot to file: %s" % outfile)