diff options
| -rwxr-xr-x | calc_pei.py | 126 |
1 files changed, 67 insertions, 59 deletions
diff --git a/calc_pei.py b/calc_pei.py index 37273d9..c8a79c1 100755 --- a/calc_pei.py +++ b/calc_pei.py @@ -3,7 +3,7 @@ # # Aaron LI # Created: 2016-04-29 -# Updated: 2016-05-01 +# Updated: 2016-05-02 # # TODO: # * to calculate the PEI error @@ -11,7 +11,7 @@ """ Calculate the power excess index (PEI), which is defined the area ratio of -the lower-left part with respect to the total rectangulr, which is further +the lower-left part with respect to the total rectangle, which is further defined by the radial power spectrum and the scale of 0.035R500 and 0.35R500, in the logarithmic space. @@ -19,11 +19,7 @@ Reference: Zhang, C., et al. 2016, ApJ """ -__version__ = "0.2.1" -__date__ = "2016-05-01" - -import sys import os import glob import argparse @@ -31,7 +27,6 @@ import json from collections import OrderedDict import numpy as np -import scipy as sp import scipy.interpolate import scipy.integrate @@ -43,46 +38,45 @@ import matplotlib.patches as patches from make_r500_regions import get_r500 +__version__ = "0.3.1" +__date__ = "2016-05-02" + plt.style.use("ggplot") def calc_pei(data, r500, interp_np=101): """ Calculate the power excess index (PEI), which is defined the area ratio - of the lower-left part with respect to the total rectangulr. + of the lower-left part with respect to the total rectangle. Arguments: * data: 3-column power spectrum data (frequency, power, power_err) * r500: R500 value in unit of the inverse of the above "frequency" * interp_np: number of data points interpolated to calculate PEI """ - freqs = data[:, 0] - psd1d = data[:, 1] - psd1d_err = data[:, 2] + freqs = data[:, 0] + psd1d = data[:, 1] # frequency values corresponding to 0.35R500 and 0.035R500 - f1 = 1.0 / (0.350 * r500) - f2 = 1.0 / (0.035 * r500) + pei_fmin = 1.0 / (0.350 * r500) + pei_fmax = 1.0 / (0.035 * r500) # switch to the logarithmic scale # XXX: how to deal with the errors mask = (freqs > 0.0) - x = np.log10(freqs[mask]) - y = np.log10(psd1d[mask]) - x1 = np.log10(f1) - x2 = np.log10(f2) + x = np.log10(freqs[mask]) + y = np.log10(psd1d[mask]) + pei_xmin = np.log10(pei_fmin) + pei_xmax = np.log10(pei_fmax) # interpolate the power spectrum - # FIXME: include a pair surrounding data points for interpolation - f_interp = sp.interpolate.interp1d(x, y, kind="cubic", assume_sorted=True) - y1 = f_interp(x1) - y2 = f_interp(x2) - if interp_np % 2 == 0: - # Simpson's rule requires an even number of intervals - interp_np += 1 - x_interp = np.linspace(x1, x2, num=interp_np) + f_interp = scipy.interpolate.interp1d(x, y, kind="linear", + assume_sorted=True) + x_interp = np.linspace(pei_xmin, pei_xmax, num=interp_np) y_interp = f_interp(x_interp) + pei_ymin = np.min(y_interp) + pei_ymax = np.max(y_interp) # calculate the PEI - area_total = abs(x1 - x2) * abs(y1 - y2) - area_below = sp.integrate.simps((y_interp-y2), x_interp) - pei_value = area_below / area_total + area_total = (pei_xmax - pei_xmin) * (pei_ymax - pei_ymin) + area_below = scipy.integrate.trapz((y_interp-pei_ymin), x_interp) + pei_value = area_below / area_total results = { "area_total": area_total, "area_below": area_below, @@ -93,15 +87,22 @@ def calc_pei(data, r500, interp_np=101): return (results, data_interp_log10) -def plot_pei(data, data_interp_log10, ax=None, fig=None): +def plot_pei(data, data_interp_log10, info={}, ax=None, fig=None): """ Make a plot to visualize the PEI rectangular. """ if ax is None: fig, ax = plt.subplots(1, 1) - freqs = data[:, 0] - psd1d = data[:, 1] + # prepare data + freqs = data[:, 0] + psd1d = data[:, 1] psd1d_err = data[:, 2] + x_interp = 10 ** data_interp_log10[:, 0] + y_interp = 10 ** data_interp_log10[:, 1] + pei_xmin = np.min(x_interp) + pei_xmax = np.max(x_interp) + pei_ymin = np.min(y_interp) + pei_ymax = np.max(y_interp) # mask = (freqs > 0.0) xmin = np.min(freqs[mask]) / 1.2 @@ -117,21 +118,25 @@ def plot_pei(data, data_interp_log10, ax=None, fig=None): ax.set_yscale("log") ax.set_xlim(xmin, xmax) ax.set_ylim(ymin, ymax) - ax.set_title("Radially Averaged Power Spectral Density") + ax.set_title("Power Spectral Density & PEI (%s; %d)" % + (info.get("name"), info.get("obsid"))) ax.set_xlabel(r"k (pixel$^{-1}$)") ax.set_ylabel("Power") + ax.text(x=xmax/1.1, y=ymax/1.2, + s="PEI = %.2f / %.2f = %.2f" % (info.get("area_below"), + info.get("area_total"), + info.get("pei")), + horizontalalignment="right", verticalalignment="top") # plot the interpolated data points and the PEI rectangle # credit: http://matplotlib.org/users/path_tutorial.html - x_interp = 10 ** data_interp_log10[:, 0] - y_interp = 10 ** data_interp_log10[:, 1] ax.plot(x_interp, y_interp, linestyle="--", marker="D", markersize=2, color="green", alpha=0.9) vertices = [ - (x_interp[0], y_interp[-1]), # left, bottom - (x_interp[0], y_interp[0]), # left, top - (x_interp[-1], y_interp[0]), # right, top - (x_interp[-1], y_interp[-1]), # right, bottom - (x_interp[0], y_interp[-1]), # ignored + (pei_xmin, pei_ymin), # left, bottom + (pei_xmin, pei_ymax), # left, top + (pei_xmax, pei_ymax), # right, top + (pei_xmax, pei_ymin), # right, bottom + (pei_xmin, pei_ymin), # ignored ] codes = [ Path.MOVETO, @@ -150,27 +155,30 @@ def plot_pei(data, data_interp_log10, ax=None, fig=None): def main(): # default arguments - default_infile = "psd.txt" - default_outfile = "pei.json" + default_infile = "psd.txt" + default_outfile = "pei.json" default_infojson = "../*_INFO.json" parser = argparse.ArgumentParser( description="Calculate the power excess index (PEI)", epilog="Version: %s (%s)" % (__version__, __date__)) parser.add_argument("-V", "--version", action="version", - version="%(prog)s " + "%s (%s)" % (__version__, __date__)) + version="%(prog)s " + "%s (%s)" % (__version__, + __date__)) parser.add_argument("-j", "--json", dest="json", required=False, - help="the *_INFO.json file (default: find %s)" % default_infojson) + help="the *_INFO.json file " + + "(default: find %s)" % default_infojson) parser.add_argument("-i", "--infile", dest="infile", - required=False, default=default_infile, - help="input data of the radial power spectrum " + \ - "(default: %s)" % default_infile) + required=False, default=default_infile, + help="input data of the radial power spectrum " + + "(default: %s)" % default_infile) parser.add_argument("-o", "--outfile", dest="outfile", - required=False, default=default_outfile, - help="output json file to save the results " + \ - "(default: %s)" % default_outfile) + required=False, default=default_outfile, + help="output json file to save the results " + + "(default: %s)" % default_outfile) parser.add_argument("-p", "--png", dest="png", default=None, - help="make PEI plot and save (default: same basename as outfile)") + help="make PEI plot and save " + + "(default: same basename as outfile)") args = parser.parse_args() if args.png is None: @@ -180,13 +188,13 @@ def main(): if args.json: info_json = args.json - json_str = open(info_json).read().rstrip().rstrip(",") - info = json.loads(json_str) - name = info["Source Name"] - obsid = int(info["Obs. ID"]) - r500 = get_r500(info) - r500_kpc = r500["r500_kpc"] - r500_pix = r500["r500_pix"] + json_str = open(info_json).read().rstrip().rstrip(",") + info = json.loads(json_str) + name = info["Source Name"] + obsid = int(info["Obs. ID"]) + r500 = get_r500(info) + r500_kpc = r500["r500_kpc"] + r500_pix = r500["r500_pix"] kpc_per_pix = r500["kpc_per_pix"] psd_data = np.loadtxt(args.infile) @@ -209,9 +217,9 @@ def main(): # Make and save a plot fig = Figure(figsize=(10, 8)) - canvas = FigureCanvas(fig) + FigureCanvas(fig) ax = fig.add_subplot(111) - plot_pei(psd_data, data_interp_log10, ax=ax, fig=fig) + plot_pei(psd_data, data_interp_log10, info=pei_data, ax=ax, fig=fig) fig.savefig(args.png, format="png", dpi=150) |