diff options
| -rwxr-xr-x | calc_sbp_excess.py | 156 |
1 files changed, 123 insertions, 33 deletions
diff --git a/calc_sbp_excess.py b/calc_sbp_excess.py index 2b472ca..6bf91b2 100755 --- a/calc_sbp_excess.py +++ b/calc_sbp_excess.py @@ -3,11 +3,14 @@ # # Aaron LI # Created: 2016-04-26 -# Updated: 2016-04-26 +# Updated: 2016-05-06 +# +# Changelog: +# 2016-05-06: +# * Add function `estimate_excess_error()` to estimate uncertainty +# * update according to `sbp_fit` renamed to `fit_sbp` +# * PEP8 fixes # -# TODO: -# * to estimate errors for the excess value and ratio (Monte Carlo: -# repeatedly fit the randomized SBP data and calculate excess) # """ @@ -19,9 +22,6 @@ NOTE: * excess ratio: excess_value / brightness_observed """ -__version__ = "0.1.0" -__date__ = "2016-04-26" - import sys import argparse @@ -31,51 +31,125 @@ from collections import OrderedDict import numpy as np from configobj import ConfigObj -from sbp_fit import make_model +from fit_sbp import make_model, make_sbpfit +__version__ = "0.2.0" +__date__ = "2016-05-06" -def calc_excess(data, model): + +def calc_excess(data, fitted_model): """ Calculate the central brightness excess value and ratio with respect to the fitted SBP model (single-beta). + TODO/XXX: + * whether to interpolate the SBP? + Arguments: * data: raw 4-column SBP data - * model: fitted SBP model + * fitted_model: fitted SBP model """ - radii = data[:, 0] - radii_err = data[:, 1] - brightness = data[:, 2] - brightness_model = model.f(radii) - rin = radii - radii_err - rout = radii + radii_err + radii = data[:, 0] + radii_err = data[:, 1] + brightness = data[:, 2] + brightness_model = fitted_model.f(radii) + rin = radii - radii_err + rout = radii + radii_err areas = np.pi * (rout**2 - rin**2) - bsum_obs = np.sum(brightness * areas) + bsum_obs = np.sum(brightness * areas) bsum_model = np.sum(brightness_model * areas) excess_value = bsum_obs - bsum_model excess_ratio = excess_value / bsum_obs excess = { - "brightness_obs": bsum_obs, - "brightness_model": bsum_model, - "excess_value": excess_value, - "excess_ratio": excess_ratio, + "brightness_obs": bsum_obs, + "brightness_model": bsum_model, + "excess_value": excess_value, + "excess_ratio": excess_ratio, } return excess +def estimate_excess_error(data, sbpfit, mctimes, verbose=False): + """ + Estimate the uncertainty of central excess value by Monte Carlo method. + + XXX/TODO: + * whether also consider the uncertainty of R500? + + Arguments: + * data: 4-column SBP data (radius, r_err, brightness, brightness_err) + * sbpfit: `SbpFit` object used to perform SBP fitting + * mctimes: number of Monte Carlo iterations + """ + brightness = data[:, 2] + brightness_err = data[:, 3] + params = sbpfit.dump_params() + ev_results = [] + er_results = [] + if verbose: + print("Estimating excess uncertainty by Monte Carlo " + + "(%d times) ..." % mctimes, end="", flush=True) + for i in range(mctimes): + if verbose and i % 100 == 0: + print("%d..." % i, end="", flush=True) + # randomize SBP data + brightness_rand = np.random.normal(brightness, scale=brightness_err) + sbpdata_rand = data.copy() + sbpdata_rand[:, 2] = brightness_rand + # load randomized data and perform SBP fit + sbpfit.reset(keep_source=True) + sbpfit.load_data(sbpdata_rand, keep_mask=True) + sbpfit.load_params(params) + sbpfit.fit() + model_rand = sbpfit.get_model() + excess = calc_excess(data=sbpdata_rand, fitted_model=model_rand) + ev_results.append(excess["excess_value"]) + er_results.append(excess["excess_ratio"]) + if verbose: + print("DONE!", flush=True) + ev_mean = np.mean(ev_results) + ev_std = np.std(ev_results) + ev_q25, ev_median, ev_q75 = np.percentile(ev_results, q=(25, 50, 75)) + er_mean = np.mean(er_results) + er_std = np.std(er_results) + er_q25, er_median, er_q75 = np.percentile(er_results, q=(25, 50, 75)) + results = { + "excess_value_mean": ev_mean, + "excess_value_median": ev_median, + "excess_value_q25": ev_q25, + "excess_value_q75": ev_q75, + "excess_value_std": ev_std, + "excess_ratio_mean": er_mean, + "excess_ratio_median": er_median, + "excess_ratio_q25": er_q25, + "excess_ratio_q75": er_q75, + "excess_ratio_std": er_std, + } + return results + + def main(): # default arguments default_outfile = "excess.json" + default_mctimes = 1000 parser = argparse.ArgumentParser( description="Calculate the central brightness excess value", 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("-v", "--verbose", dest="verbose", + required=False, action="store_true", + help="show verbose information") + parser.add_argument("-m", "--mctimes", dest="mctimes", required=False, + type=int, default=default_mctimes, + help="number of MC times to estimate excess error " + + "(default: %d)" % default_mctimes) parser.add_argument("config", help="Config file for SBP fitting") parser.add_argument("outfile", nargs="?", default=default_outfile, - help="Output json file to save the results " + \ - "(default: %s)" % default_outfile) + help="Output json file to save the results " + + "(default: %s)" % default_outfile) args = parser.parse_args() config = ConfigObj(args.config) @@ -83,7 +157,7 @@ def main(): config_model = config[modelname] model = make_model(config, modelname=modelname) - print("SBP model: %s" % model.name, file=sys.stderr) + print("SBP model: %s" % model.long_name, file=sys.stderr) sbpfit_outfile = config.get("outfile") sbpfit_outfile = config_model.get("outfile", sbpfit_outfile) @@ -94,17 +168,33 @@ def main(): for pname, pvalue in sbpfit_results["params"].items(): model.set_param(name=pname, value=pvalue[0]) + sbpfit = make_sbpfit(config, model=model) + sbpdata = np.loadtxt(config["sbpfile"]) - excess = calc_excess(data=sbpdata, model=model) + excess = calc_excess(data=sbpdata, fitted_model=model) + excess_err = estimate_excess_error(data=sbpdata, sbpfit=sbpfit, + mctimes=args.mctimes, + verbose=args.verbose) excess_data = OrderedDict([ - ("name", config["name"]), - ("obsid", int(config["obsid"])), - ("model", modelname), - ("brightness_obs", excess["brightness_obs"]), - ("brightness_model", excess["brightness_model"]), - ("excess_value", excess["excess_value"]), - ("excess_ratio", excess["excess_ratio"]), + ("name", config["name"]), + ("obsid", int(config["obsid"])), + ("model", modelname), + ("brightness_obs", excess["brightness_obs"]), + ("brightness_model", excess["brightness_model"]), + ("excess_value", excess["excess_value"]), + ("excess_value_mean", excess_err["excess_value_mean"]), + ("excess_value_median", excess_err["excess_value_median"]), + ("excess_value_q25", excess_err["excess_value_q25"]), + ("excess_value_q75", excess_err["excess_value_q75"]), + ("excess_value_std", excess_err["excess_value_std"]), + ("excess_ratio", excess["excess_ratio"]), + ("excess_ratio_mean", excess_err["excess_ratio_mean"]), + ("excess_ratio_median", excess_err["excess_ratio_median"]), + ("excess_ratio_q25", excess_err["excess_ratio_q25"]), + ("excess_ratio_q75", excess_err["excess_ratio_q75"]), + ("excess_ratio_std", excess_err["excess_ratio_std"]), + ("mc_times", args.mctimes), ]) excess_data_json = json.dumps(excess_data, indent=2) print(excess_data_json) |