Completed
Push — master ( cc962b...240a5b )
by Andy
54s
created

papers.plot_cluster_comparison()   F

Complexity

Conditions 36

Size

Total Lines 198

Duplication

Lines 0
Ratio 0 %
Metric Value
cc 36
dl 0
loc 198
rs 2

How to fix   Long Method    Complexity   

Long Method

Small methods make your code easier to understand, in particular if combined with a good name. Besides, if your method is small, finding a good name is usually much easier.

For example, if you find yourself adding comments to a method's body, this is usually a good sign to extract the commented part to a new method, and use the comment as a starting point when coming up with a good name for this new method.

Commonly applied refactorings include:

Complexity

Complex classes like papers.plot_cluster_comparison() often do a lot of different things. To break such a class down, we need to identify a cohesive component within that class. A common approach to find such a component is to look for fields/methods that share the same prefixes, or suffixes.

Once you have determined the fields that belong together, you can apply the Extract Class refactoring. If the component makes sense as a sub-class, Extract Subclass is also a candidate, and is often faster.

1
"""
2
Plot some cluster abundance stuff compared to ASPCAP.
3
"""
4
5
6
import matplotlib
7
matplotlib.rcParams["text.usetex"] = True
8
9
10
import matplotlib.pyplot as plt
11
import numpy as np
12
13
from astropy.io import fits
14
from matplotlib.ticker import MaxNLocator
15
16
17
# x vs y
18
# membership criteria.
19
20
# C vs N
21
22
# O vs Na
23
# O vs Mg
24
# O vs Al
25
# O vs S
26
27
# C+N+O vs FE_H?
28
# Mg vs Al
29
30
31
32
33
def plot_cluster_comparison(data, cluster_name, membership, x_elements, 
34
    y_elements, used_cannon_for_target_selection=True, vel_lim=None,
35
    xlims=None, ylims=None):
36
    """
37
    membership should be same len as data
38
    """
39
40
    candidate_color, membership_color = ("#666666", "#3498DB")
41
    candidate_color, membership_color = ("#BBBBBB", "#3498DB")
42
    tc_suffix, aspcap_suffix = ("", "_ASPCAP")
43
    
44
    candidates = data["FIELD"] == cluster_name
45
46
    membership_kwds = {"s": 50, "lw": 1.5}
47
    candidate_kwds = {"s": 30, "marker": "+", "lw": 1.5}
48
49
    fig, axes = plt.subplots(6, 2, figsize=(5.1, 16))
50
    axes = np.array(axes).flatten()
51
52
    axes[0].set_visible(False)
53
    axes[1].set_visible(False)
54
55
    top_ax = plt.subplot(6, 1, 1)
56
57
58
    # Vhelio and FE_H_1 (our metallicity?)
59
    suffix = tc_suffix if used_cannon_for_target_selection else aspcap_suffix
60
    top_ax.scatter(
61
        data["VHELIO_AVG"][candidates], data["FE_H" + suffix][candidates], 
62
        facecolor=candidate_color, label=r"$\texttt{{FIELD = {0}}}$".format(cluster_name),
63
        **candidate_kwds)
64
    top_ax.scatter(
65
        data["VHELIO_AVG"][membership], data["FE_H" + suffix][membership], 
66
        facecolor=membership_color, **membership_kwds)
67
68
    N, M = len(data["VHELIO_AVG"][candidates]), len(data["VHELIO_AVG"][membership])
69
    top_ax.text(0.05, 0.95, r"${:,}$".format(N), color=candidate_color,
70
            verticalalignment="top", horizontalalignment="left",
71
            transform=top_ax.transAxes)  
72
    top_ax.text(0.05, 0.95 - 0.11, r"${:,}$".format(M), color=membership_color,
73
            verticalalignment="top", horizontalalignment="left",
74
            transform=top_ax.transAxes)  
75
76
    #top_ax.legend(frameon=True, fontsize=11, loc="upper left")
77
78
    top_ax.set_xlabel(r"$V_{\rm helio}$ $(\rm{km}$ $\rm{s}^{-1})$")
79
    if used_cannon_for_target_selection:
80
        top_ax.set_ylabel(r"$[\rm{Fe}/\rm{H}]$ $(\rm{The}$ $\rm{Cannon})$")
81
    else:
82
        top_ax.set_ylabel(r"$[\rm{Fe}/\rm{H}]$ $(\rm{ASPCAP})$")
83
84
    top_ax.set_title(r"$\rm{{{0}}}$ $\rm{{membership}}$ $\rm{{selection}}$".format(
85
        cluster_name))
86
87
    top_ax.xaxis.set_major_locator(MaxNLocator(4))
88
    top_ax.yaxis.set_major_locator(MaxNLocator(4))
89
90
91
    
92
    for j, (element_x, element_y) in enumerate(zip(x_elements, y_elements)):
93
94
        x_wrt_fe, y_wrt_fe = (True, True)
95
96
        if element_x.lower() == "fe":
97
            x_wrt_fe = False
98
99
        if element_y.lower() == "fe":
100
            y_wrt_fe = False
101
102
        # X/Y for The Cannon
103
        for i, (mask, color) \
104
        in enumerate(zip((candidates, membership), (candidate_color, membership_color))):
105
106
            if "," in element_x:
0 ignored issues
show
Duplication introduced by
This code seems to be duplicated in your project.

Duplicated code is one of the most pungent code smells. If you need to duplicate the same code in three or more different places, we strongly encourage you to look into extracting the code into a single class or operation.

You can also find more detailed suggestions in the “Code” section of your repository.

Loading history...
107
                x = 0
108
                for each in element_x.split(","):
109
                    x += data["{0}_H{1}".format(each.upper(), tc_suffix)]
110
                    if x_wrt_fe:
111
                        x = x - data["FE_H{}".format(tc_suffix)]
112
            else:
113
                x = data["{0}_H{1}".format(element_x.upper(), tc_suffix)]
114
                if x_wrt_fe:
115
                    x = x - data["FE_H{}".format(tc_suffix)]
116
117
            if "," in element_y:
0 ignored issues
show
Duplication introduced by
This code seems to be duplicated in your project.

Duplicated code is one of the most pungent code smells. If you need to duplicate the same code in three or more different places, we strongly encourage you to look into extracting the code into a single class or operation.

You can also find more detailed suggestions in the “Code” section of your repository.

Loading history...
118
                y = 0
119
                for each in element_y.split(","):
120
                    y += data["{0}_H{1}".format(each.upper(), tc_suffix)]
121
                    if y_wrt_fe:
122
                        y = y - data["FE_H{}".format(tc_suffix)]
123
            else:
124
                y = data["{0}_H{1}".format(element_y.upper(), tc_suffix)]
125
                if y_wrt_fe:
126
                    y = y - data["FE_H{}".format(tc_suffix)]
127
128
            kwds = candidate_kwds if i == 0 else membership_kwds
129
            axes[2*j + 2 + 1].scatter(x[mask], y[mask], facecolor=color, **kwds)
130
131
            # Quote the number of points.
132
            axes[2*j + 2 + 1].text(0.05, 0.95 - i * 0.10, r"${:,}$".format(len(x[mask])),
133
                color=color,
134
                verticalalignment="top", horizontalalignment="left",
135
                transform=axes[2*j + 2 + 1].transAxes)
136
137
138
        if xlims is None:
139
            tc_xlims = axes[2*j + 2 + 1].get_xlim()
140
            percent = 0.20 # 10%
141
            half_ptp = (np.ptp(tc_xlims) * (1 + percent))/2.
142
            tc_xlims = (np.mean(tc_xlims) - half_ptp, half_ptp + np.mean(tc_xlims))
143
144
        else:
145
            tc_xlims = xlims
146
147
        if ylims is None:
148
            tc_ylims = axes[2*j + 2 + 1].get_ylim()
149
            # Expand the scale just a little bit.
150
            percent = 0.20 # 10%
151
            half_ptp = (np.ptp(tc_ylims) * (1 + percent))/2.
152
            tc_ylims = (np.mean(tc_ylims) - half_ptp, half_ptp + np.mean(tc_ylims))
153
        else:
154
            tc_ylims = ylims
155
156
        # X/Y for ASPCAP.
157
        for i, (mask, color) \
158
        in enumerate(zip((candidates, membership), (candidate_color, membership_color))):
159
160
            if "," in element_x:
0 ignored issues
show
Duplication introduced by
This code seems to be duplicated in your project.

Duplicated code is one of the most pungent code smells. If you need to duplicate the same code in three or more different places, we strongly encourage you to look into extracting the code into a single class or operation.

You can also find more detailed suggestions in the “Code” section of your repository.

Loading history...
161
                x = 0
162
                for each in element_x.split(","):
163
                    x += data["{0}_H{1}".format(each.upper(), aspcap_suffix)]
164
                    if x_wrt_fe:
165
                        x = x - data["FE_H{}".format(aspcap_suffix)]
166
            else:
167
                x = data["{0}_H{1}".format(element_x.upper(), aspcap_suffix)]
168
                if x_wrt_fe:
169
                    x = x - data["FE_H{}".format(aspcap_suffix)]
170
171
172
            if "," in element_y:
0 ignored issues
show
Duplication introduced by
This code seems to be duplicated in your project.

Duplicated code is one of the most pungent code smells. If you need to duplicate the same code in three or more different places, we strongly encourage you to look into extracting the code into a single class or operation.

You can also find more detailed suggestions in the “Code” section of your repository.

Loading history...
173
                y = 0
174
                for each in element_y.split(","):
175
                    y += data["{0}_H{1}".format(each.upper(), aspcap_suffix)]
176
                    if y_wrt_fe:
177
                        y = y - data["FE_H{}".format(aspcap_suffix)]
178
            else:
179
                y = data["{0}_H{1}".format(element_y.upper(), aspcap_suffix)]
180
                if y_wrt_fe:
181
                    y = y - data["FE_H{}".format(aspcap_suffix)]
182
183
            kwds = candidate_kwds if i == 0 else membership_kwds
184
            axes[2*j + 2].scatter(x[mask], y[mask], facecolor=color, **kwds)
185
186
            N = sum((tc_xlims[1] > x[mask]) * (x[mask] > tc_xlims[0]) \
187
                  * (tc_ylims[1] > y[mask]) * (y[mask] > tc_ylims[0]))
188
            axes[2*j + 2].text(0.05, 0.95 - i * 0.10, r"${:,}$".format(N), color=color,
189
                verticalalignment="top", horizontalalignment="left",
190
                transform=axes[2*j + 2].transAxes)  
191
192
        
193
        if j == 0:
194
            axes[2*j + 2].set_title(r"${\rm ASPCAP}$", y=1.05)
195
            axes[2*j + 2 + 1].set_title(r"${\rm The}$ ${\rm Cannon}$", y=1.05)
196
    
197
198
        for ax in (axes[2*j + 2], axes[2*j + 2 + 1]):
199
            ax.set_xlim(tc_xlims)
200
            ax.set_ylim(tc_ylims)
201
202
            ax.xaxis.set_major_locator(MaxNLocator(4))
203
            ax.yaxis.set_major_locator(MaxNLocator(4))
204
205
            ax.set_xlabel(r"$[\rm{{{0}}}/\rm{{{1}}}]$".format(element_x.title(),
206
                "Fe" if x_wrt_fe else "H"))
207
        
208
        if "," in element_y:
209
            axes[2*j + 2].set_ylabel(r"$[(\rm{{{0}}})/{{{1}}}\rm{{{2}}}]$".format(
210
                element_y.replace(",", "+"), element_y.count(",") + 1,
211
                "Fe" if y_wrt_fe else "H"))
212
        else:
213
            axes[2*j + 2].set_ylabel(r"$[\rm{{{0}}}/\rm{{{1}}}]$".format(element_y.title(),
214
                    "Fe" if y_wrt_fe else "H"))
215
        axes[2*j + 2 + 1].yaxis.set_ticklabels([])
216
217
        
218
    for ax in axes[2:]:
219
        ax.set(adjustable='box-forced', aspect=np.ptp(ax.get_xlim())/np.ptp(ax.get_ylim()))
220
221
    fig.tight_layout()
222
223
    if vel_lim is not None:
224
        top_ax.set_xlim(vel_lim)
225
226
    fig.subplots_adjust(hspace=-0.0, bottom=0.03)
227
    pos = top_ax.get_position()
228
    top_ax.set_position([pos.x0, pos.y0 + 0.06, pos.width, pos.height - 0.06])
229
230
    return fig
231
232
233
if __name__ == "__main__":
234
235
    # To speed up development cycle..
236
    try:
237
        data
238
    except:
239
240
        data = fits.open(
241
            "/Users/arc/research/apogee/results-unregularized/results-unregularized-matched.fits")[1].data
242
243
        use = data["physical"]
244
        data = data[use]
245
246
247
    """
248
    M67 186
249
    N6791 173
250
    N2243 187
251
    N188 252
252
    N1333 380
253
    N6819 307
254
    N7789 327
255
256
    M35N2158 189
257
    M54SGRC1 386
258
    M5PAL5 317
259
    """
260
261
262
    # M13
263
    M13_members = (data["FIELD"] == "M13") \
264
                * (data["VHELIO_AVG"] > -265) * (data["VHELIO_AVG"] < -220) \
265
                * (data["FE_H"] < -1.2)
266
    M13_figure = plot_cluster_comparison(data, "M13", M13_members, 
267
        ["C", "O", "Mg", "Ca", "Fe"], ["N", "Na", "Al", "S", "C,N,O"],
268
        )
269
    M13_figure.savefig("M13_comparison.pdf", dpi=300)
270
271
272
    # M 92
273
    M92_members = (data["FIELD"] == "M92") \
274
                * (data["VHELIO_AVG"] > -140) * (data["VHELIO_AVG"] < -100) \
275
                * (data["FE_H"] < -1.7)
276
    M92_figure = plot_cluster_comparison(data, "M92", M92_members, 
277
        ["C", "O", "Mg", "Ca", "Fe"], ["N", "Na", "Al", "S", "C,N,O"],
278
        )
279
    M92_figure.savefig("M92_comparison.pdf", dpi=300)
280
281
282
    # M 53
283
    M53_members = (data["FIELD"] == "M53") \
284
                * (data["VHELIO_AVG"] > -80) * (data["VHELIO_AVG"] < -40) \
285
                * (data["FE_H"] < -1.5)
286
    M53_figure = plot_cluster_comparison(data, "M53", M53_members, 
287
        ["C", "O", "Mg", "Ca", "Fe"], ["N", "Na", "Al", "S", "C,N,O"],
288
        )
289
    M53_figure.savefig("M53_comparison.pdf", dpi=300)
290
291
292
    # M 15
293
    M15_members = (data["FIELD"] == "M15") \
294
                * (data["VHELIO_AVG"] > -130) * (data["VHELIO_AVG"] < 80) \
295
                * (data["FE_H"] < -1.7)
296
    
297
    M15_figure = plot_cluster_comparison(data, "M15", M15_members, 
298
        ["C", "O", "Mg", "Ca", "Fe"], ["N", "Na", "Al", "S", "C,N,O"],
299
        vel_lim=(-400, 150))
300
    M15_figure.savefig("M15_comparison.pdf", dpi=300)
301
        
302
303
    # M 3
304
    M3_members = (data["FIELD"] == "M3") \
305
               * (data["VHELIO_AVG"] > -165) * (data["VHELIO_AVG"] < -120) \
306
               * (data["FE_H"] > -1.65) * (data["FE_H"] < -1.2)
307
308
    M3_figure = plot_cluster_comparison(data, "M3", M3_members,
309
        ["C", "O", "Mg", "Ca", "Fe"], ["N", "Na", "Al", "S", "C,N,O"])
310
    M3_figure.savefig("M3_comparison.pdf", dpi=300)
311
312