#!/usr/bin/env python3 # -*- coding: utf-8 -*- # # Aaron LI # Created: 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 # # """ Calculate the central brightness excess value and ratio with respect to the fitted SBP model (i.e., single-beta model). NOTE: * excess value: brightness_observed - brightness_model_predicted * excess ratio: excess_value / brightness_observed """ import sys import argparse import json from collections import OrderedDict import numpy as np from configobj import ConfigObj from fit_sbp import make_model, make_sbpfit __version__ = "0.2.0" __date__ = "2016-05-06" 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 * fitted_model: fitted SBP model """ 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_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, } 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__)) 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) args = parser.parse_args() config = ConfigObj(args.config) modelname = config["model"] config_model = config[modelname] model = make_model(config, modelname=modelname) print("SBP model: %s" % model.long_name, file=sys.stderr) sbpfit_outfile = config.get("outfile") sbpfit_outfile = config_model.get("outfile", sbpfit_outfile) sbpfit_results = json.load(open(sbpfit_outfile), object_pairs_hook=OrderedDict) # Load fitted parameters for model 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, 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_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) open(args.outfile, "w").write(excess_data_json+"\n") if __name__ == "__main__": main() # vim: set ts=4 sw=4 tw=0 fenc=utf-8 ft=python: #