From 5268684f61e712f852219d6dd327f907bc236bfa Mon Sep 17 00:00:00 2001 From: Aaron LI Date: Wed, 4 May 2016 17:37:20 +0800 Subject: rename "sbp_fit.py" to "fit_sbp.py" --- python/fit_sbp.py | 807 ++++++++++++++++++++++++++++++++++++++++++++++++++++++ python/sbp_fit.py | 807 ------------------------------------------------------ 2 files changed, 807 insertions(+), 807 deletions(-) create mode 100755 python/fit_sbp.py delete mode 100755 python/sbp_fit.py diff --git a/python/fit_sbp.py b/python/fit_sbp.py new file mode 100755 index 0000000..c22e0c8 --- /dev/null +++ b/python/fit_sbp.py @@ -0,0 +1,807 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +# +# Aaron LI +# Created: 2016-03-13 +# Updated: 2016-04-26 +# +# Changelogs: +# 2016-04-26: +# * Reorder some methods of classes 'FitModelSBeta' and 'FitModelDBeta' +# * Change the output file extension from ".txt" to ".json" +# 2016-04-21: +# * Plot another X axis with unit "r500", with R500 values marked +# * Adjust output image size/resolution +# 2016-04-20: +# * Support "pix" and "kpc" units +# * Allow ignore data w.r.t R500 value +# * Major changes to the config syntax +# * Add commandline argument to select the sbp model +# 2016-04-05: +# * Allow fix parameters +# 2016-03-31: +# * Remove `ci_report()' +# * Add `make_results()' to orgnize all results as s Python dictionary +# * Report results as json string +# 2016-03-28: +# * Add `main()', `make_model()' +# * Use `configobj' to handle configurations +# * Save fit results and plot +# * Add `ci_report()' +# 2016-03-14: +# * Refactor classes `FitModelSBeta' and `FitModelDBeta' +# * Add matplotlib plot support +# * Add `ignore_data()' and `notice_data()' support +# * Add classes `FitModelSBetaNorm' and `FitModelDBetaNorm' +# +# TODO: +# * to allow fit the outer beta component, then fix it, and fit the inner one +# * to integrate basic information of config file to the output json +# * to output the ignored radius range in the same unit as input sbp data +# + +""" +Fit the surface brightness profile (SBP) with the single-beta model: + s(r) = s0 * [1.0 + (r/rc)^2] ^ (0.5-3*beta) + bkg +or the double-beta model: + s(r) = s01 * [1.0 + (r/rc1)^2] ^ (0.5-3*beta1) + + s02 * [1.0 + (r/rc2)^2] ^ (0.5-3*beta2) + bkg + + +Sample config file: +------------------------------------------------- +name = +obsid = +r500_pix = +r500_kpc = + +sbpfile = sbprofile.txt +# unit of radius: pix (default) or kpc +unit = pixel + +# sbp model: "sbeta" or "dbeta" +model = sbeta +#model = dbeta + +# output file to store the fitting results +outfile = sbpfit.json +# output file to save the fitting plot +imgfile = sbpfit.png + +# data range to be ignored during fitting (same unit as the above "unit") +#ignore = 0.0-20.0, +# specify the ignore range w.r.t R500 ("r500_pix" or "r500_kpc" required) +#ignore_r500 = 0.0-0.15, + +[sbeta] +# model-related options (OVERRIDE the upper level options) +outfile = sbpfit_sbeta.json +imgfile = sbpfit_sbeta.png +#ignore = 0.0-20.0, +#ignore_r500 = 0.0-0.15, + [[params]] + # model parameters + # name = initial, lower, upper, variable (FIXED/False to fix the parameter) + s0 = 1.0e-8, 0.0, 1.0e-6 + rc = 30.0, 5.0, 1.0e4 + #rc = 30.0, 5.0, 1.0e4, FIXED + beta = 0.7, 0.3, 1.1 + bkg = 1.0e-10, 0.0, 1.0e-8 + + +[dbeta] +outfile = sbpfit_dbeta.json +imgfile = sbpfit_dbeta.png +#ignore = 0.0-20.0, +#ignore_r500 = 0.0-0.15, + [[params]] + s01 = 1.0e-8, 0.0, 1.0e-6 + rc1 = 50.0, 10.0, 1.0e4 + beta1 = 0.7, 0.3, 1.1 + s02 = 1.0e-8, 0.0, 1.0e-6 + rc2 = 30.0, 2.0, 5.0e2 + beta2 = 0.7, 0.3, 1.1 + bkg = 1.0e-10, 0.0, 1.0e-8 +------------------------------------------------- +""" + +__version__ = "0.6.2" +__date__ = "2016-04-26" + + +import os +import sys +import re +import argparse +import json +from collections import OrderedDict + +import numpy as np +import lmfit +import matplotlib.pyplot as plt +from matplotlib.backends.backend_agg import FigureCanvasAgg as FigureCanvas +from matplotlib.figure import Figure +from configobj import ConfigObj + + +plt.style.use("ggplot") + + +class FitModel: + """ + Meta-class of the fitting model. + + The supplied `func' should have the following syntax: + y = f(x, params) + where the `params' is `lmfit.Parameters' instance which contains all + the model parameters to be fitted, and should be provided as well. + """ + def __init__(self, name=None, func=None, params=lmfit.Parameters()): + self.name = name + self.func = func + self.params = params + + def f(self, x): + return self.func(x, self.params) + + def get_param(self, name=None): + """ + Return the requested `Parameter' object or the whole + `Parameters' object of no name supplied. + """ + try: + return self.params[name] + except KeyError: + return self.params + + def set_param(self, name, *args, **kwargs): + """ + Set the properties of the specified parameter. + """ + param = self.params[name] + param.set(*args, **kwargs) + + def plot(self, params, xdata, ax): + """ + Plot the fitted model. + """ + f_fitted = lambda x: self.func(x, params) + ydata = f_fitted(xdata) + ax.plot(xdata, ydata, 'k-') + +class FitModelSBeta(FitModel): + """ + The single-beta model to be fitted. + Single-beta model, with a constant background. + """ + params = lmfit.Parameters() + params.add_many( # (name, value, vary, min, max, expr) + ("s0", 1.0e-8, True, 0.0, 1.0e-6, None), + ("rc", 30.0, True, 1.0, 1.0e4, None), + ("beta", 0.7, True, 0.3, 1.1, None), + ("bkg", 1.0e-9, True, 0.0, 1.0e-7, None)) + + def __init__(self): + super(self.__class__, self).__init__(name="Single-beta", + func=self.sbeta, params=self.params) + + @staticmethod + def sbeta(r, params): + parvals = params.valuesdict() + s0 = parvals["s0"] + rc = parvals["rc"] + beta = parvals["beta"] + bkg = parvals["bkg"] + return s0 * np.power((1 + (r/rc)**2), (0.5 - 3*beta)) + bkg + + def plot(self, params, xdata, ax): + """ + Plot the fitted model, as well as the fitted parameters. + """ + super(self.__class__, self).plot(params, xdata, ax) + ydata = self.sbeta(xdata, params) + # fitted paramters + ax.vlines(x=params["rc"].value, ymin=min(ydata), ymax=max(ydata), + linestyles="dashed") + ax.hlines(y=params["bkg"].value, xmin=min(xdata), xmax=max(xdata), + linestyles="dashed") + ax.text(x=params["rc"].value, y=min(ydata), + s="beta: %.2f\nrc: %.2f" % (params["beta"].value, + params["rc"].value)) + ax.text(x=min(xdata), y=min(ydata), + s="bkg: %.3e" % params["bkg"].value, + verticalalignment="top") + + +class FitModelDBeta(FitModel): + """ + The double-beta model to be fitted. + Double-beta model, with a constant background. + + NOTE: + the first beta component (s01, rc1, beta1) describes the main and + outer SBP; while the second beta component (s02, rc2, beta2) accounts + for the central brightness excess. + """ + params = lmfit.Parameters() + params.add("s01", value=1.0e-8, min=0.0, max=1.0e-6) + params.add("rc1", value=50.0, min=10.0, max=1.0e4) + params.add("beta1", value=0.7, min=0.3, max=1.1) + #params.add("df_s0", value=1.0e-8, min=0.0, max=1.0e-6) + #params.add("s02", expr="s01 + df_s0") + params.add("s02", value=1.0e-8, min=0.0, max=1.0e-6) + #params.add("df_rc", value=30.0, min=0.0, max=1.0e4) + #params.add("rc2", expr="rc1 - df_rc") + params.add("rc2", value=20.0, min=1.0, max=5.0e2) + params.add("beta2", value=0.7, min=0.3, max=1.1) + params.add("bkg", value=1.0e-9, min=0.0, max=1.0e-7) + + def __init__(self): + super(self.__class__, self).__init__(name="Double-beta", + func=self.dbeta, params=self.params) + + @classmethod + def dbeta(self, r, params): + return self.beta1(r, params) + self.beta2(r, params) + + @staticmethod + def beta1(r, params): + """ + This beta component describes the main/outer part of the SBP. + """ + parvals = params.valuesdict() + s01 = parvals["s01"] + rc1 = parvals["rc1"] + beta1 = parvals["beta1"] + bkg = parvals["bkg"] + return s01 * np.power((1 + (r/rc1)**2), (0.5 - 3*beta1)) + bkg + + @staticmethod + def beta2(r, params): + """ + This beta component describes the central/excess part of the SBP. + """ + parvals = params.valuesdict() + s02 = parvals["s02"] + rc2 = parvals["rc2"] + beta2 = parvals["beta2"] + return s02 * np.power((1 + (r/rc2)**2), (0.5 - 3*beta2)) + + def plot(self, params, xdata, ax): + """ + Plot the fitted model, and each beta component, + as well as the fitted parameters. + """ + super(self.__class__, self).plot(params, xdata, ax) + beta1_ydata = self.beta1(xdata, params) + beta2_ydata = self.beta2(xdata, params) + ax.plot(xdata, beta1_ydata, 'b-.') + ax.plot(xdata, beta2_ydata, 'b-.') + # fitted paramters + ydata = beta1_ydata + beta2_ydata + ax.vlines(x=params["rc1"].value, ymin=min(ydata), ymax=max(ydata), + linestyles="dashed") + ax.vlines(x=params["rc2"].value, ymin=min(ydata), ymax=max(ydata), + linestyles="dashed") + ax.hlines(y=params["bkg"].value, xmin=min(xdata), xmax=max(xdata), + linestyles="dashed") + ax.text(x=params["rc1"].value, y=min(ydata), + s="beta1: %.2f\nrc1: %.2f" % (params["beta1"].value, + params["rc1"].value)) + ax.text(x=params["rc2"].value, y=min(ydata), + s="beta2: %.2f\nrc2: %.2f" % (params["beta2"].value, + params["rc2"].value)) + ax.text(x=min(xdata), y=min(ydata), + s="bkg: %.3e" % params["bkg"].value, + verticalalignment="top") + + +class FitModelSBetaNorm(FitModel): + """ + The single-beta model to be fitted. + Single-beta model, with a constant background. + Normalized the `s0' and `bkg' parameters by take the logarithm. + """ + params = lmfit.Parameters() + params.add_many( # (name, value, vary, min, max, expr) + ("log10_s0", -8.0, True, -12.0, -6.0, None), + ("rc", 30.0, True, 1.0, 1.0e4, None), + ("beta", 0.7, True, 0.3, 1.1, None), + ("log10_bkg", -9.0, True, -12.0, -7.0, None)) + + @staticmethod + def sbeta(r, params): + parvals = params.valuesdict() + s0 = 10 ** parvals["log10_s0"] + rc = parvals["rc"] + beta = parvals["beta"] + bkg = 10 ** parvals["log10_bkg"] + return s0 * np.power((1 + (r/rc)**2), (0.5 - 3*beta)) + bkg + + def __init__(self): + super(self.__class__, self).__init__(name="Single-beta", + func=self.sbeta, params=self.params) + + def plot(self, params, xdata, ax): + """ + Plot the fitted model, as well as the fitted parameters. + """ + super(self.__class__, self).plot(params, xdata, ax) + ydata = self.sbeta(xdata, params) + # fitted paramters + ax.vlines(x=params["rc"].value, ymin=min(ydata), ymax=max(ydata), + linestyles="dashed") + ax.hlines(y=(10 ** params["bkg"].value), xmin=min(xdata), + xmax=max(xdata), linestyles="dashed") + ax.text(x=params["rc"].value, y=min(ydata), + s="beta: %.2f\nrc: %.2f" % (params["beta"].value, + params["rc"].value)) + ax.text(x=min(xdata), y=min(ydata), + s="bkg: %.3e" % (10 ** params["bkg"].value), + verticalalignment="top") + + +class FitModelDBetaNorm(FitModel): + """ + The double-beta model to be fitted. + Double-beta model, with a constant background. + Normalized the `s01', `s02' and `bkg' parameters by take the logarithm. + + NOTE: + the first beta component (s01, rc1, beta1) describes the main and + outer SBP; while the second beta component (s02, rc2, beta2) accounts + for the central brightness excess. + """ + params = lmfit.Parameters() + params.add("log10_s01", value=-8.0, min=-12.0, max=-6.0) + params.add("rc1", value=50.0, min=10.0, max=1.0e4) + params.add("beta1", value=0.7, min=0.3, max=1.1) + #params.add("df_s0", value=1.0e-8, min=0.0, max=1.0e-6) + #params.add("s02", expr="s01 + df_s0") + params.add("log10_s02", value=-8.0, min=-12.0, max=-6.0) + #params.add("df_rc", value=30.0, min=0.0, max=1.0e4) + #params.add("rc2", expr="rc1 - df_rc") + params.add("rc2", value=20.0, min=1.0, max=5.0e2) + params.add("beta2", value=0.7, min=0.3, max=1.1) + params.add("log10_bkg", value=-9.0, min=-12.0, max=-7.0) + + @staticmethod + def beta1(r, params): + """ + This beta component describes the main/outer part of the SBP. + """ + parvals = params.valuesdict() + s01 = 10 ** parvals["log10_s01"] + rc1 = parvals["rc1"] + beta1 = parvals["beta1"] + bkg = 10 ** parvals["log10_bkg"] + return s01 * np.power((1 + (r/rc1)**2), (0.5 - 3*beta1)) + bkg + + @staticmethod + def beta2(r, params): + """ + This beta component describes the central/excess part of the SBP. + """ + parvals = params.valuesdict() + s02 = 10 ** parvals["log10_s02"] + rc2 = parvals["rc2"] + beta2 = parvals["beta2"] + return s02 * np.power((1 + (r/rc2)**2), (0.5 - 3*beta2)) + + @classmethod + def dbeta(self, r, params): + return self.beta1(r, params) + self.beta2(r, params) + + def __init__(self): + super(self.__class__, self).__init__(name="Double-beta", + func=self.dbeta, params=self.params) + + def plot(self, params, xdata, ax): + """ + Plot the fitted model, and each beta component, + as well as the fitted parameters. + """ + super(self.__class__, self).plot(params, xdata, ax) + beta1_ydata = self.beta1(xdata, params) + beta2_ydata = self.beta2(xdata, params) + ax.plot(xdata, beta1_ydata, 'b-.') + ax.plot(xdata, beta2_ydata, 'b-.') + # fitted paramters + ydata = beta1_ydata + beta2_ydata + ax.vlines(x=params["log10_rc1"].value, ymin=min(ydata), ymax=max(ydata), + linestyles="dashed") + ax.vlines(x=params["rc2"].value, ymin=min(ydata), ymax=max(ydata), + linestyles="dashed") + ax.hlines(y=(10 ** params["bkg"].value), xmin=min(xdata), + xmax=max(xdata), linestyles="dashed") + ax.text(x=params["rc1"].value, y=min(ydata), + s="beta1: %.2f\nrc1: %.2f" % (params["beta1"].value, + params["rc1"].value)) + ax.text(x=params["rc2"].value, y=min(ydata), + s="beta2: %.2f\nrc2: %.2f" % (params["beta2"].value, + params["rc2"].value)) + ax.text(x=min(xdata), y=min(ydata), + s="bkg: %.3e" % (10 ** params["bkg"].value), + verticalalignment="top") + + +class SbpFit: + """ + Class to handle the SBP fitting with single-/double-beta model. + """ + def __init__(self, model, method="lbfgsb", + xdata=None, ydata=None, xerr=None, yerr=None, xunit="pix", + name=None, obsid=None, r500_pix=None, r500_kpc=None): + self.method = method + self.model = model + self.load_data(xdata=xdata, ydata=ydata, xerr=xerr, yerr=yerr, + xunit=xunit) + self.set_source(name=name, obsid=obsid, r500_pix=r500_pix, + r500_kpc=r500_kpc) + + def set_source(self, name, obsid=None, r500_pix=None, r500_kpc=None): + self.name = name + try: + self.obsid = int(obsid) + except TypeError: + self.obsid = None + try: + self.r500_pix = float(r500_pix) + except TypeError: + self.r500_pix = None + try: + self.r500_kpc = float(r500_kpc) + except TypeError: + self.r500_kpc = None + try: + self.kpc_per_pix = self.r500_kpc / self.r500_pix + except (TypeError, ZeroDivisionError): + self.kpc_per_pix = -1 + + def load_data(self, xdata, ydata, xerr, yerr, xunit="pix"): + self.xdata = xdata + self.ydata = ydata + self.xerr = xerr + self.yerr = yerr + if xdata is not None: + self.mask = np.ones(xdata.shape, dtype=np.bool) + else: + self.mask = None + if xunit.lower() in ["pix", "pixel"]: + self.xunit = "pix" + elif xunit.lower() == "kpc": + self.xunit = "kpc" + else: + raise ValueError("invalid xunit: %s" % xunit) + + def ignore_data(self, xmin=None, xmax=None, unit=None): + """ + Ignore the data points within range [xmin, xmax]. + If xmin is None, then xmin=min(xdata); + if xmax is None, then xmax=max(xdata). + + if unit is None, then assume the same unit as `self.xunit'. + """ + if unit is None: + unit = self.xunit + if xmin is not None: + xmin = self.convert_unit(xmin, unit=unit) + else: + xmin = np.min(self.xdata) + if xmax is not None: + xmax = self.convert_unit(xmax, unit=unit) + else: + xmax = np.max(self.xdata) + ignore_idx = np.logical_and(self.xdata >= xmin, self.xdata <= xmax) + self.mask[ignore_idx] = False + # reset `f_residual' + self.f_residual = None + + def notice_data(self, xmin=None, xmax=None, unit=None): + """ + Notice the data points within range [xmin, xmax]. + If xmin is None, then xmin=min(xdata); + if xmax is None, then xmax=max(xdata). + + if unit is None, then assume the same unit as `self.xunit'. + """ + if unit is None: + unit = self.xunit + if xmin is not None: + xmin = self.convert_unit(xmin, unit=unit) + else: + xmin = np.min(self.xdata) + if xmax is not None: + xmax = self.convert_unit(xmax, unit=unit) + else: + xmax = np.max(self.xdata) + notice_idx = np.logical_and(self.xdata >= xmin, self.xdata <= xmax) + self.mask[notice_idx] = True + # reset `f_residual' + self.f_residual = None + + def convert_unit(self, x, unit): + """ + Convert the value x in given unit to be the unit `self.xunit' + """ + if unit == self.xunit: + return x + elif (unit == "pix") and (self.xunit == "kpc"): + return (x / self.r500_pix * self.r500_kpc) + elif (unit == "kpc") and (self.xunit == "pix"): + return (x / self.r500_kpc * self.r500_pix) + elif (unit == "r500") and (self.xunit == "pix"): + return (x * self.r500_pix) + elif (unit == "r500") and (self.xunit == "kpc"): + return (x * self.r500_kpc) + else: + raise ValueError("invalid units: %s vs. %s" % (unit, self.xunit)) + + def convert_to_r500(self, x, unit=None): + """ + Convert the value x in given unit to be in unit "r500". + """ + if unit is None: + unit = self.xunit + if unit == "r500": + return x + elif unit == "pix": + return (x / self.r500_pix) + elif unit == "kpc": + return (x / self.r500_kpc) + else: + raise ValueError("invalid unit: %s" % unit) + + def set_residual(self): + def f_residual(params): + if self.yerr is None: + return self.model.func(self.xdata[self.mask], params) - \ + self.ydata + else: + return (self.model.func(self.xdata[self.mask], params) - \ + self.ydata[self.mask]) / self.yerr[self.mask] + self.f_residual = f_residual + + def fit(self, method=None): + if method is None: + method = self.method + if not hasattr(self, "f_residual") or self.f_residual is None: + self.set_residual() + self.fitter = lmfit.Minimizer(self.f_residual, self.model.params) + self.fitted = self.fitter.minimize(method=method) + self.fitted_model = lambda x: self.model.func(x, self.fitted.params) + + def calc_ci(self, sigmas=[0.68, 0.90]): + # `conf_interval' requires the fitted results have valid `stderr', + # so we need to re-fit the model with method `leastsq'. + fitted = self.fitter.minimize(method="leastsq", + params=self.fitted.params) + self.ci, self.trace = lmfit.conf_interval(self.fitter, fitted, + sigmas=sigmas, trace=True) + + def make_results(self): + """ + Make the `self.results' dictionary which contains all the fitting + results as well as the confidence intervals. + """ + fitted = self.fitted + self.results = OrderedDict() + ## fitting results + self.results.update( + nfev = fitted.nfev, + ndata = fitted.ndata, + nvarys = fitted.nvarys, # number of varible paramters + nfree = fitted.nfree, # degree of freem + chisqr = fitted.chisqr, + redchi = fitted.redchi, + aic = fitted.aic, + bic = fitted.bic) + params = fitted.params + pnames = list(params.keys()) + pvalues = OrderedDict() + for pn in pnames: + par = params.get(pn) + pvalues[pn] = [par.value, par.min, par.max, par.vary] + self.results["params"] = pvalues + ## confidence intervals + if hasattr(self, "ci") and self.ci is not None: + ci = self.ci + ci_values = OrderedDict() + ci_sigmas = [ "ci%02d" % (v[0]*100) for v in ci.get(pnames[0]) ] + ci_names = sorted(list(set(ci_sigmas))) + ci_idx = { k: [] for k in ci_names } + for cn, idx in zip(ci_sigmas, range(len(ci_sigmas))): + ci_idx[cn].append(idx) + # parameters ci + for pn in pnames: + ci_pv = OrderedDict() + pv = [ v[1] for v in ci.get(pn) ] + # best + pv_best = pv[ ci_idx["ci00"][0] ] + ci_pv["best"] = pv_best + # ci of each sigma + pv2 = [ v-pv_best for v in pv ] + for cn in ci_names[1:]: + ci_pv[cn] = [ pv2[idx] for idx in ci_idx[cn] ] + ci_values[pn] = ci_pv + self.results["ci"] = ci_values + + def report(self, outfile=sys.stdout): + if not hasattr(self, "results") or self.results is None: + self.make_results() + jd = json.dumps(self.results, indent=2) + print(jd, file=outfile) + + def plot(self, ax=None, fig=None, r500_axis=True): + """ + Arguments: + * r500_axis: whether to add a second X axis in unit "r500" + """ + if ax is None: + fig, ax = plt.subplots(1, 1) + # noticed data points + eb = ax.errorbar(self.xdata[self.mask], self.ydata[self.mask], + xerr=self.xerr[self.mask], yerr=self.yerr[self.mask], + fmt="none") + # ignored data points + ignore_mask = np.logical_not(self.mask) + if np.sum(ignore_mask) > 0: + eb = ax.errorbar(self.xdata[ignore_mask], self.ydata[ignore_mask], + xerr=self.xerr[ignore_mask], yerr=self.yerr[ignore_mask], + fmt="none") + eb[-1][0].set_linestyle("-.") + # fitted model + xmax = self.xdata[-1] + self.xerr[-1] + xpred = np.power(10, np.linspace(0, np.log10(xmax), 2*len(self.xdata))) + ypred = self.fitted_model(xpred) + ymin = min(min(self.ydata), min(ypred)) + ymax = max(max(self.ydata), max(ypred)) + self.model.plot(params=self.fitted.params, xdata=xpred, ax=ax) + ax.set_xscale("log") + ax.set_yscale("log") + ax.set_xlim(1.0, xmax) + ax.set_ylim(ymin/1.2, ymax*1.2) + name = self.name + if self.obsid is not None: + name += "; %s" % self.obsid + ax.set_title("Fitted Surface Brightness Profile (%s)" % name) + ax.set_xlabel("Radius (%s)" % self.xunit) + ax.set_ylabel(r"Surface Brightness (photons/cm$^2$/pixel$^2$/s)") + ax.text(x=xmax, y=ymax, + s="redchi: %.2f / %.2f = %.2f" % (self.fitted.chisqr, + self.fitted.nfree, self.fitted.chisqr/self.fitted.nfree), + horizontalalignment="right", verticalalignment="top") + plot_ret = [fig, ax] + if r500_axis: + # Add a second X-axis with labels in unit "r500" + # Credit: https://stackoverflow.com/a/28192477/4856091 + try: + ax.title.set_position([0.5, 1.1]) # raise title position + ax2 = ax.twiny() + # NOTE: the ORDER of the following lines MATTERS + ax2.set_xscale(ax.get_xscale()) + ax2_ticks = ax.get_xticks() + ax2.set_xticks(ax2_ticks) + ax2.set_xbound(ax.get_xbound()) + ax2.set_xticklabels([ "%.2g" % self.convert_to_r500(x) + for x in ax2_ticks ]) + ax2.set_xlabel("Radius (r500; r500 = %s pix = %s kpc)" % (\ + self.r500_pix, self.r500_kpc)) + ax2.grid(False) + plot_ret.append(ax2) + except ValueError: + # cannot convert X values to unit "r500" + pass + # automatically adjust layout + fig.tight_layout() + return plot_ret + + +def make_model(config, modelname): + """ + Make the model with parameters set according to the config. + """ + if modelname == "sbeta": + # single-beta model + model = FitModelSBeta() + elif modelname == "dbeta": + # double-beta model + model = FitModelDBeta() + else: + raise ValueError("Invalid model: %s" % modelname) + # set initial values and bounds for the model parameters + params = config[modelname]["params"] + for p, value in params.items(): + variable = True + if len(value) == 4 and value[3].upper() in ["FIXED", "FALSE"]: + variable = False + model.set_param(name=p, value=float(value[0]), + min=float(value[1]), max=float(value[2]), vary=variable) + return model + + +def main(): + # parser for command line options and arguments + parser = argparse.ArgumentParser( + description="Fit surface brightness profile with " + \ + "single-/double-beta model", + epilog="Version: %s (%s)" % (__version__, __date__)) + parser.add_argument("-V", "--version", action="version", + version="%(prog)s " + "%s (%s)" % (__version__, __date__)) + parser.add_argument("config", help="Config file for SBP fitting") + # exclusive argument group for model selection + grp_model = parser.add_mutually_exclusive_group(required=False) + grp_model.add_argument("-s", "--sbeta", dest="sbeta", + action="store_true", help="single-beta model for SBP") + grp_model.add_argument("-d", "--dbeta", dest="dbeta", + action="store_true", help="double-beta model for SBP") + # + args = parser.parse_args() + + config = ConfigObj(args.config) + + # determine the model name + if args.sbeta: + modelname = "sbeta" + elif args.dbeta: + modelname = "dbeta" + else: + modelname = config["model"] + + config_model = config[modelname] + # determine the "outfile" and "imgfile" + outfile = config.get("outfile") + outfile = config_model.get("outfile", outfile) + imgfile = config.get("imgfile") + imgfile = config_model.get("imgfile", imgfile) + + # SBP fitting model + model = make_model(config, modelname=modelname) + + # sbp data and fit object + sbpdata = np.loadtxt(config["sbpfile"]) + sbpfit = SbpFit(model=model, xdata=sbpdata[:, 0], xerr=sbpdata[:, 1], + ydata=sbpdata[:, 2], yerr=sbpdata[:, 3], + xunit=config.get("unit", "pix")) + sbpfit.set_source(name=config["name"], obsid=config.get("obsid"), + r500_pix=config.get("r500_pix"), r500_kpc=config.get("r500_kpc")) + + # apply data range ignorance + if "ignore" in config.keys(): + for ig in config.as_list("ignore"): + xmin, xmax = map(float, ig.split("-")) + sbpfit.ignore_data(xmin=xmin, xmax=xmax) + if "ignore_r500" in config.keys(): + for ig in config.as_list("ignore_r500"): + xmin, xmax = map(float, ig.split("-")) + sbpfit.ignore_data(xmin=xmin, xmax=xmax, unit="r500") + + # apply additional data range ignorance specified within model section + if "ignore" in config_model.keys(): + for ig in config_model.as_list("ignore"): + xmin, xmax = map(float, ig.split("-")) + sbpfit.ignore_data(xmin=xmin, xmax=xmax) + if "ignore_r500" in config_model.keys(): + for ig in config_model.as_list("ignore_r500"): + xmin, xmax = map(float, ig.split("-")) + sbpfit.ignore_data(xmin=xmin, xmax=xmax, unit="r500") + + # fit and calculate confidence intervals + sbpfit.fit() + sbpfit.calc_ci() + sbpfit.report() + with open(outfile, "w") as ofile: + sbpfit.report(outfile=ofile) + + # make and save a plot + fig = Figure(figsize=(10, 8)) + canvas = FigureCanvas(fig) + ax = fig.add_subplot(111) + sbpfit.plot(ax=ax, fig=fig, r500_axis=True) + fig.savefig(imgfile, dpi=150) + + +if __name__ == "__main__": + main() + +# vim: set ts=4 sw=4 tw=0 fenc=utf-8 ft=python: # diff --git a/python/sbp_fit.py b/python/sbp_fit.py deleted file mode 100755 index c22e0c8..0000000 --- a/python/sbp_fit.py +++ /dev/null @@ -1,807 +0,0 @@ -#!/usr/bin/env python3 -# -*- coding: utf-8 -*- -# -# Aaron LI -# Created: 2016-03-13 -# Updated: 2016-04-26 -# -# Changelogs: -# 2016-04-26: -# * Reorder some methods of classes 'FitModelSBeta' and 'FitModelDBeta' -# * Change the output file extension from ".txt" to ".json" -# 2016-04-21: -# * Plot another X axis with unit "r500", with R500 values marked -# * Adjust output image size/resolution -# 2016-04-20: -# * Support "pix" and "kpc" units -# * Allow ignore data w.r.t R500 value -# * Major changes to the config syntax -# * Add commandline argument to select the sbp model -# 2016-04-05: -# * Allow fix parameters -# 2016-03-31: -# * Remove `ci_report()' -# * Add `make_results()' to orgnize all results as s Python dictionary -# * Report results as json string -# 2016-03-28: -# * Add `main()', `make_model()' -# * Use `configobj' to handle configurations -# * Save fit results and plot -# * Add `ci_report()' -# 2016-03-14: -# * Refactor classes `FitModelSBeta' and `FitModelDBeta' -# * Add matplotlib plot support -# * Add `ignore_data()' and `notice_data()' support -# * Add classes `FitModelSBetaNorm' and `FitModelDBetaNorm' -# -# TODO: -# * to allow fit the outer beta component, then fix it, and fit the inner one -# * to integrate basic information of config file to the output json -# * to output the ignored radius range in the same unit as input sbp data -# - -""" -Fit the surface brightness profile (SBP) with the single-beta model: - s(r) = s0 * [1.0 + (r/rc)^2] ^ (0.5-3*beta) + bkg -or the double-beta model: - s(r) = s01 * [1.0 + (r/rc1)^2] ^ (0.5-3*beta1) + - s02 * [1.0 + (r/rc2)^2] ^ (0.5-3*beta2) + bkg - - -Sample config file: -------------------------------------------------- -name = -obsid = -r500_pix = -r500_kpc = - -sbpfile = sbprofile.txt -# unit of radius: pix (default) or kpc -unit = pixel - -# sbp model: "sbeta" or "dbeta" -model = sbeta -#model = dbeta - -# output file to store the fitting results -outfile = sbpfit.json -# output file to save the fitting plot -imgfile = sbpfit.png - -# data range to be ignored during fitting (same unit as the above "unit") -#ignore = 0.0-20.0, -# specify the ignore range w.r.t R500 ("r500_pix" or "r500_kpc" required) -#ignore_r500 = 0.0-0.15, - -[sbeta] -# model-related options (OVERRIDE the upper level options) -outfile = sbpfit_sbeta.json -imgfile = sbpfit_sbeta.png -#ignore = 0.0-20.0, -#ignore_r500 = 0.0-0.15, - [[params]] - # model parameters - # name = initial, lower, upper, variable (FIXED/False to fix the parameter) - s0 = 1.0e-8, 0.0, 1.0e-6 - rc = 30.0, 5.0, 1.0e4 - #rc = 30.0, 5.0, 1.0e4, FIXED - beta = 0.7, 0.3, 1.1 - bkg = 1.0e-10, 0.0, 1.0e-8 - - -[dbeta] -outfile = sbpfit_dbeta.json -imgfile = sbpfit_dbeta.png -#ignore = 0.0-20.0, -#ignore_r500 = 0.0-0.15, - [[params]] - s01 = 1.0e-8, 0.0, 1.0e-6 - rc1 = 50.0, 10.0, 1.0e4 - beta1 = 0.7, 0.3, 1.1 - s02 = 1.0e-8, 0.0, 1.0e-6 - rc2 = 30.0, 2.0, 5.0e2 - beta2 = 0.7, 0.3, 1.1 - bkg = 1.0e-10, 0.0, 1.0e-8 -------------------------------------------------- -""" - -__version__ = "0.6.2" -__date__ = "2016-04-26" - - -import os -import sys -import re -import argparse -import json -from collections import OrderedDict - -import numpy as np -import lmfit -import matplotlib.pyplot as plt -from matplotlib.backends.backend_agg import FigureCanvasAgg as FigureCanvas -from matplotlib.figure import Figure -from configobj import ConfigObj - - -plt.style.use("ggplot") - - -class FitModel: - """ - Meta-class of the fitting model. - - The supplied `func' should have the following syntax: - y = f(x, params) - where the `params' is `lmfit.Parameters' instance which contains all - the model parameters to be fitted, and should be provided as well. - """ - def __init__(self, name=None, func=None, params=lmfit.Parameters()): - self.name = name - self.func = func - self.params = params - - def f(self, x): - return self.func(x, self.params) - - def get_param(self, name=None): - """ - Return the requested `Parameter' object or the whole - `Parameters' object of no name supplied. - """ - try: - return self.params[name] - except KeyError: - return self.params - - def set_param(self, name, *args, **kwargs): - """ - Set the properties of the specified parameter. - """ - param = self.params[name] - param.set(*args, **kwargs) - - def plot(self, params, xdata, ax): - """ - Plot the fitted model. - """ - f_fitted = lambda x: self.func(x, params) - ydata = f_fitted(xdata) - ax.plot(xdata, ydata, 'k-') - -class FitModelSBeta(FitModel): - """ - The single-beta model to be fitted. - Single-beta model, with a constant background. - """ - params = lmfit.Parameters() - params.add_many( # (name, value, vary, min, max, expr) - ("s0", 1.0e-8, True, 0.0, 1.0e-6, None), - ("rc", 30.0, True, 1.0, 1.0e4, None), - ("beta", 0.7, True, 0.3, 1.1, None), - ("bkg", 1.0e-9, True, 0.0, 1.0e-7, None)) - - def __init__(self): - super(self.__class__, self).__init__(name="Single-beta", - func=self.sbeta, params=self.params) - - @staticmethod - def sbeta(r, params): - parvals = params.valuesdict() - s0 = parvals["s0"] - rc = parvals["rc"] - beta = parvals["beta"] - bkg = parvals["bkg"] - return s0 * np.power((1 + (r/rc)**2), (0.5 - 3*beta)) + bkg - - def plot(self, params, xdata, ax): - """ - Plot the fitted model, as well as the fitted parameters. - """ - super(self.__class__, self).plot(params, xdata, ax) - ydata = self.sbeta(xdata, params) - # fitted paramters - ax.vlines(x=params["rc"].value, ymin=min(ydata), ymax=max(ydata), - linestyles="dashed") - ax.hlines(y=params["bkg"].value, xmin=min(xdata), xmax=max(xdata), - linestyles="dashed") - ax.text(x=params["rc"].value, y=min(ydata), - s="beta: %.2f\nrc: %.2f" % (params["beta"].value, - params["rc"].value)) - ax.text(x=min(xdata), y=min(ydata), - s="bkg: %.3e" % params["bkg"].value, - verticalalignment="top") - - -class FitModelDBeta(FitModel): - """ - The double-beta model to be fitted. - Double-beta model, with a constant background. - - NOTE: - the first beta component (s01, rc1, beta1) describes the main and - outer SBP; while the second beta component (s02, rc2, beta2) accounts - for the central brightness excess. - """ - params = lmfit.Parameters() - params.add("s01", value=1.0e-8, min=0.0, max=1.0e-6) - params.add("rc1", value=50.0, min=10.0, max=1.0e4) - params.add("beta1", value=0.7, min=0.3, max=1.1) - #params.add("df_s0", value=1.0e-8, min=0.0, max=1.0e-6) - #params.add("s02", expr="s01 + df_s0") - params.add("s02", value=1.0e-8, min=0.0, max=1.0e-6) - #params.add("df_rc", value=30.0, min=0.0, max=1.0e4) - #params.add("rc2", expr="rc1 - df_rc") - params.add("rc2", value=20.0, min=1.0, max=5.0e2) - params.add("beta2", value=0.7, min=0.3, max=1.1) - params.add("bkg", value=1.0e-9, min=0.0, max=1.0e-7) - - def __init__(self): - super(self.__class__, self).__init__(name="Double-beta", - func=self.dbeta, params=self.params) - - @classmethod - def dbeta(self, r, params): - return self.beta1(r, params) + self.beta2(r, params) - - @staticmethod - def beta1(r, params): - """ - This beta component describes the main/outer part of the SBP. - """ - parvals = params.valuesdict() - s01 = parvals["s01"] - rc1 = parvals["rc1"] - beta1 = parvals["beta1"] - bkg = parvals["bkg"] - return s01 * np.power((1 + (r/rc1)**2), (0.5 - 3*beta1)) + bkg - - @staticmethod - def beta2(r, params): - """ - This beta component describes the central/excess part of the SBP. - """ - parvals = params.valuesdict() - s02 = parvals["s02"] - rc2 = parvals["rc2"] - beta2 = parvals["beta2"] - return s02 * np.power((1 + (r/rc2)**2), (0.5 - 3*beta2)) - - def plot(self, params, xdata, ax): - """ - Plot the fitted model, and each beta component, - as well as the fitted parameters. - """ - super(self.__class__, self).plot(params, xdata, ax) - beta1_ydata = self.beta1(xdata, params) - beta2_ydata = self.beta2(xdata, params) - ax.plot(xdata, beta1_ydata, 'b-.') - ax.plot(xdata, beta2_ydata, 'b-.') - # fitted paramters - ydata = beta1_ydata + beta2_ydata - ax.vlines(x=params["rc1"].value, ymin=min(ydata), ymax=max(ydata), - linestyles="dashed") - ax.vlines(x=params["rc2"].value, ymin=min(ydata), ymax=max(ydata), - linestyles="dashed") - ax.hlines(y=params["bkg"].value, xmin=min(xdata), xmax=max(xdata), - linestyles="dashed") - ax.text(x=params["rc1"].value, y=min(ydata), - s="beta1: %.2f\nrc1: %.2f" % (params["beta1"].value, - params["rc1"].value)) - ax.text(x=params["rc2"].value, y=min(ydata), - s="beta2: %.2f\nrc2: %.2f" % (params["beta2"].value, - params["rc2"].value)) - ax.text(x=min(xdata), y=min(ydata), - s="bkg: %.3e" % params["bkg"].value, - verticalalignment="top") - - -class FitModelSBetaNorm(FitModel): - """ - The single-beta model to be fitted. - Single-beta model, with a constant background. - Normalized the `s0' and `bkg' parameters by take the logarithm. - """ - params = lmfit.Parameters() - params.add_many( # (name, value, vary, min, max, expr) - ("log10_s0", -8.0, True, -12.0, -6.0, None), - ("rc", 30.0, True, 1.0, 1.0e4, None), - ("beta", 0.7, True, 0.3, 1.1, None), - ("log10_bkg", -9.0, True, -12.0, -7.0, None)) - - @staticmethod - def sbeta(r, params): - parvals = params.valuesdict() - s0 = 10 ** parvals["log10_s0"] - rc = parvals["rc"] - beta = parvals["beta"] - bkg = 10 ** parvals["log10_bkg"] - return s0 * np.power((1 + (r/rc)**2), (0.5 - 3*beta)) + bkg - - def __init__(self): - super(self.__class__, self).__init__(name="Single-beta", - func=self.sbeta, params=self.params) - - def plot(self, params, xdata, ax): - """ - Plot the fitted model, as well as the fitted parameters. - """ - super(self.__class__, self).plot(params, xdata, ax) - ydata = self.sbeta(xdata, params) - # fitted paramters - ax.vlines(x=params["rc"].value, ymin=min(ydata), ymax=max(ydata), - linestyles="dashed") - ax.hlines(y=(10 ** params["bkg"].value), xmin=min(xdata), - xmax=max(xdata), linestyles="dashed") - ax.text(x=params["rc"].value, y=min(ydata), - s="beta: %.2f\nrc: %.2f" % (params["beta"].value, - params["rc"].value)) - ax.text(x=min(xdata), y=min(ydata), - s="bkg: %.3e" % (10 ** params["bkg"].value), - verticalalignment="top") - - -class FitModelDBetaNorm(FitModel): - """ - The double-beta model to be fitted. - Double-beta model, with a constant background. - Normalized the `s01', `s02' and `bkg' parameters by take the logarithm. - - NOTE: - the first beta component (s01, rc1, beta1) describes the main and - outer SBP; while the second beta component (s02, rc2, beta2) accounts - for the central brightness excess. - """ - params = lmfit.Parameters() - params.add("log10_s01", value=-8.0, min=-12.0, max=-6.0) - params.add("rc1", value=50.0, min=10.0, max=1.0e4) - params.add("beta1", value=0.7, min=0.3, max=1.1) - #params.add("df_s0", value=1.0e-8, min=0.0, max=1.0e-6) - #params.add("s02", expr="s01 + df_s0") - params.add("log10_s02", value=-8.0, min=-12.0, max=-6.0) - #params.add("df_rc", value=30.0, min=0.0, max=1.0e4) - #params.add("rc2", expr="rc1 - df_rc") - params.add("rc2", value=20.0, min=1.0, max=5.0e2) - params.add("beta2", value=0.7, min=0.3, max=1.1) - params.add("log10_bkg", value=-9.0, min=-12.0, max=-7.0) - - @staticmethod - def beta1(r, params): - """ - This beta component describes the main/outer part of the SBP. - """ - parvals = params.valuesdict() - s01 = 10 ** parvals["log10_s01"] - rc1 = parvals["rc1"] - beta1 = parvals["beta1"] - bkg = 10 ** parvals["log10_bkg"] - return s01 * np.power((1 + (r/rc1)**2), (0.5 - 3*beta1)) + bkg - - @staticmethod - def beta2(r, params): - """ - This beta component describes the central/excess part of the SBP. - """ - parvals = params.valuesdict() - s02 = 10 ** parvals["log10_s02"] - rc2 = parvals["rc2"] - beta2 = parvals["beta2"] - return s02 * np.power((1 + (r/rc2)**2), (0.5 - 3*beta2)) - - @classmethod - def dbeta(self, r, params): - return self.beta1(r, params) + self.beta2(r, params) - - def __init__(self): - super(self.__class__, self).__init__(name="Double-beta", - func=self.dbeta, params=self.params) - - def plot(self, params, xdata, ax): - """ - Plot the fitted model, and each beta component, - as well as the fitted parameters. - """ - super(self.__class__, self).plot(params, xdata, ax) - beta1_ydata = self.beta1(xdata, params) - beta2_ydata = self.beta2(xdata, params) - ax.plot(xdata, beta1_ydata, 'b-.') - ax.plot(xdata, beta2_ydata, 'b-.') - # fitted paramters - ydata = beta1_ydata + beta2_ydata - ax.vlines(x=params["log10_rc1"].value, ymin=min(ydata), ymax=max(ydata), - linestyles="dashed") - ax.vlines(x=params["rc2"].value, ymin=min(ydata), ymax=max(ydata), - linestyles="dashed") - ax.hlines(y=(10 ** params["bkg"].value), xmin=min(xdata), - xmax=max(xdata), linestyles="dashed") - ax.text(x=params["rc1"].value, y=min(ydata), - s="beta1: %.2f\nrc1: %.2f" % (params["beta1"].value, - params["rc1"].value)) - ax.text(x=params["rc2"].value, y=min(ydata), - s="beta2: %.2f\nrc2: %.2f" % (params["beta2"].value, - params["rc2"].value)) - ax.text(x=min(xdata), y=min(ydata), - s="bkg: %.3e" % (10 ** params["bkg"].value), - verticalalignment="top") - - -class SbpFit: - """ - Class to handle the SBP fitting with single-/double-beta model. - """ - def __init__(self, model, method="lbfgsb", - xdata=None, ydata=None, xerr=None, yerr=None, xunit="pix", - name=None, obsid=None, r500_pix=None, r500_kpc=None): - self.method = method - self.model = model - self.load_data(xdata=xdata, ydata=ydata, xerr=xerr, yerr=yerr, - xunit=xunit) - self.set_source(name=name, obsid=obsid, r500_pix=r500_pix, - r500_kpc=r500_kpc) - - def set_source(self, name, obsid=None, r500_pix=None, r500_kpc=None): - self.name = name - try: - self.obsid = int(obsid) - except TypeError: - self.obsid = None - try: - self.r500_pix = float(r500_pix) - except TypeError: - self.r500_pix = None - try: - self.r500_kpc = float(r500_kpc) - except TypeError: - self.r500_kpc = None - try: - self.kpc_per_pix = self.r500_kpc / self.r500_pix - except (TypeError, ZeroDivisionError): - self.kpc_per_pix = -1 - - def load_data(self, xdata, ydata, xerr, yerr, xunit="pix"): - self.xdata = xdata - self.ydata = ydata - self.xerr = xerr - self.yerr = yerr - if xdata is not None: - self.mask = np.ones(xdata.shape, dtype=np.bool) - else: - self.mask = None - if xunit.lower() in ["pix", "pixel"]: - self.xunit = "pix" - elif xunit.lower() == "kpc": - self.xunit = "kpc" - else: - raise ValueError("invalid xunit: %s" % xunit) - - def ignore_data(self, xmin=None, xmax=None, unit=None): - """ - Ignore the data points within range [xmin, xmax]. - If xmin is None, then xmin=min(xdata); - if xmax is None, then xmax=max(xdata). - - if unit is None, then assume the same unit as `self.xunit'. - """ - if unit is None: - unit = self.xunit - if xmin is not None: - xmin = self.convert_unit(xmin, unit=unit) - else: - xmin = np.min(self.xdata) - if xmax is not None: - xmax = self.convert_unit(xmax, unit=unit) - else: - xmax = np.max(self.xdata) - ignore_idx = np.logical_and(self.xdata >= xmin, self.xdata <= xmax) - self.mask[ignore_idx] = False - # reset `f_residual' - self.f_residual = None - - def notice_data(self, xmin=None, xmax=None, unit=None): - """ - Notice the data points within range [xmin, xmax]. - If xmin is None, then xmin=min(xdata); - if xmax is None, then xmax=max(xdata). - - if unit is None, then assume the same unit as `self.xunit'. - """ - if unit is None: - unit = self.xunit - if xmin is not None: - xmin = self.convert_unit(xmin, unit=unit) - else: - xmin = np.min(self.xdata) - if xmax is not None: - xmax = self.convert_unit(xmax, unit=unit) - else: - xmax = np.max(self.xdata) - notice_idx = np.logical_and(self.xdata >= xmin, self.xdata <= xmax) - self.mask[notice_idx] = True - # reset `f_residual' - self.f_residual = None - - def convert_unit(self, x, unit): - """ - Convert the value x in given unit to be the unit `self.xunit' - """ - if unit == self.xunit: - return x - elif (unit == "pix") and (self.xunit == "kpc"): - return (x / self.r500_pix * self.r500_kpc) - elif (unit == "kpc") and (self.xunit == "pix"): - return (x / self.r500_kpc * self.r500_pix) - elif (unit == "r500") and (self.xunit == "pix"): - return (x * self.r500_pix) - elif (unit == "r500") and (self.xunit == "kpc"): - return (x * self.r500_kpc) - else: - raise ValueError("invalid units: %s vs. %s" % (unit, self.xunit)) - - def convert_to_r500(self, x, unit=None): - """ - Convert the value x in given unit to be in unit "r500". - """ - if unit is None: - unit = self.xunit - if unit == "r500": - return x - elif unit == "pix": - return (x / self.r500_pix) - elif unit == "kpc": - return (x / self.r500_kpc) - else: - raise ValueError("invalid unit: %s" % unit) - - def set_residual(self): - def f_residual(params): - if self.yerr is None: - return self.model.func(self.xdata[self.mask], params) - \ - self.ydata - else: - return (self.model.func(self.xdata[self.mask], params) - \ - self.ydata[self.mask]) / self.yerr[self.mask] - self.f_residual = f_residual - - def fit(self, method=None): - if method is None: - method = self.method - if not hasattr(self, "f_residual") or self.f_residual is None: - self.set_residual() - self.fitter = lmfit.Minimizer(self.f_residual, self.model.params) - self.fitted = self.fitter.minimize(method=method) - self.fitted_model = lambda x: self.model.func(x, self.fitted.params) - - def calc_ci(self, sigmas=[0.68, 0.90]): - # `conf_interval' requires the fitted results have valid `stderr', - # so we need to re-fit the model with method `leastsq'. - fitted = self.fitter.minimize(method="leastsq", - params=self.fitted.params) - self.ci, self.trace = lmfit.conf_interval(self.fitter, fitted, - sigmas=sigmas, trace=True) - - def make_results(self): - """ - Make the `self.results' dictionary which contains all the fitting - results as well as the confidence intervals. - """ - fitted = self.fitted - self.results = OrderedDict() - ## fitting results - self.results.update( - nfev = fitted.nfev, - ndata = fitted.ndata, - nvarys = fitted.nvarys, # number of varible paramters - nfree = fitted.nfree, # degree of freem - chisqr = fitted.chisqr, - redchi = fitted.redchi, - aic = fitted.aic, - bic = fitted.bic) - params = fitted.params - pnames = list(params.keys()) - pvalues = OrderedDict() - for pn in pnames: - par = params.get(pn) - pvalues[pn] = [par.value, par.min, par.max, par.vary] - self.results["params"] = pvalues - ## confidence intervals - if hasattr(self, "ci") and self.ci is not None: - ci = self.ci - ci_values = OrderedDict() - ci_sigmas = [ "ci%02d" % (v[0]*100) for v in ci.get(pnames[0]) ] - ci_names = sorted(list(set(ci_sigmas))) - ci_idx = { k: [] for k in ci_names } - for cn, idx in zip(ci_sigmas, range(len(ci_sigmas))): - ci_idx[cn].append(idx) - # parameters ci - for pn in pnames: - ci_pv = OrderedDict() - pv = [ v[1] for v in ci.get(pn) ] - # best - pv_best = pv[ ci_idx["ci00"][0] ] - ci_pv["best"] = pv_best - # ci of each sigma - pv2 = [ v-pv_best for v in pv ] - for cn in ci_names[1:]: - ci_pv[cn] = [ pv2[idx] for idx in ci_idx[cn] ] - ci_values[pn] = ci_pv - self.results["ci"] = ci_values - - def report(self, outfile=sys.stdout): - if not hasattr(self, "results") or self.results is None: - self.make_results() - jd = json.dumps(self.results, indent=2) - print(jd, file=outfile) - - def plot(self, ax=None, fig=None, r500_axis=True): - """ - Arguments: - * r500_axis: whether to add a second X axis in unit "r500" - """ - if ax is None: - fig, ax = plt.subplots(1, 1) - # noticed data points - eb = ax.errorbar(self.xdata[self.mask], self.ydata[self.mask], - xerr=self.xerr[self.mask], yerr=self.yerr[self.mask], - fmt="none") - # ignored data points - ignore_mask = np.logical_not(self.mask) - if np.sum(ignore_mask) > 0: - eb = ax.errorbar(self.xdata[ignore_mask], self.ydata[ignore_mask], - xerr=self.xerr[ignore_mask], yerr=self.yerr[ignore_mask], - fmt="none") - eb[-1][0].set_linestyle("-.") - # fitted model - xmax = self.xdata[-1] + self.xerr[-1] - xpred = np.power(10, np.linspace(0, np.log10(xmax), 2*len(self.xdata))) - ypred = self.fitted_model(xpred) - ymin = min(min(self.ydata), min(ypred)) - ymax = max(max(self.ydata), max(ypred)) - self.model.plot(params=self.fitted.params, xdata=xpred, ax=ax) - ax.set_xscale("log") - ax.set_yscale("log") - ax.set_xlim(1.0, xmax) - ax.set_ylim(ymin/1.2, ymax*1.2) - name = self.name - if self.obsid is not None: - name += "; %s" % self.obsid - ax.set_title("Fitted Surface Brightness Profile (%s)" % name) - ax.set_xlabel("Radius (%s)" % self.xunit) - ax.set_ylabel(r"Surface Brightness (photons/cm$^2$/pixel$^2$/s)") - ax.text(x=xmax, y=ymax, - s="redchi: %.2f / %.2f = %.2f" % (self.fitted.chisqr, - self.fitted.nfree, self.fitted.chisqr/self.fitted.nfree), - horizontalalignment="right", verticalalignment="top") - plot_ret = [fig, ax] - if r500_axis: - # Add a second X-axis with labels in unit "r500" - # Credit: https://stackoverflow.com/a/28192477/4856091 - try: - ax.title.set_position([0.5, 1.1]) # raise title position - ax2 = ax.twiny() - # NOTE: the ORDER of the following lines MATTERS - ax2.set_xscale(ax.get_xscale()) - ax2_ticks = ax.get_xticks() - ax2.set_xticks(ax2_ticks) - ax2.set_xbound(ax.get_xbound()) - ax2.set_xticklabels([ "%.2g" % self.convert_to_r500(x) - for x in ax2_ticks ]) - ax2.set_xlabel("Radius (r500; r500 = %s pix = %s kpc)" % (\ - self.r500_pix, self.r500_kpc)) - ax2.grid(False) - plot_ret.append(ax2) - except ValueError: - # cannot convert X values to unit "r500" - pass - # automatically adjust layout - fig.tight_layout() - return plot_ret - - -def make_model(config, modelname): - """ - Make the model with parameters set according to the config. - """ - if modelname == "sbeta": - # single-beta model - model = FitModelSBeta() - elif modelname == "dbeta": - # double-beta model - model = FitModelDBeta() - else: - raise ValueError("Invalid model: %s" % modelname) - # set initial values and bounds for the model parameters - params = config[modelname]["params"] - for p, value in params.items(): - variable = True - if len(value) == 4 and value[3].upper() in ["FIXED", "FALSE"]: - variable = False - model.set_param(name=p, value=float(value[0]), - min=float(value[1]), max=float(value[2]), vary=variable) - return model - - -def main(): - # parser for command line options and arguments - parser = argparse.ArgumentParser( - description="Fit surface brightness profile with " + \ - "single-/double-beta model", - epilog="Version: %s (%s)" % (__version__, __date__)) - parser.add_argument("-V", "--version", action="version", - version="%(prog)s " + "%s (%s)" % (__version__, __date__)) - parser.add_argument("config", help="Config file for SBP fitting") - # exclusive argument group for model selection - grp_model = parser.add_mutually_exclusive_group(required=False) - grp_model.add_argument("-s", "--sbeta", dest="sbeta", - action="store_true", help="single-beta model for SBP") - grp_model.add_argument("-d", "--dbeta", dest="dbeta", - action="store_true", help="double-beta model for SBP") - # - args = parser.parse_args() - - config = ConfigObj(args.config) - - # determine the model name - if args.sbeta: - modelname = "sbeta" - elif args.dbeta: - modelname = "dbeta" - else: - modelname = config["model"] - - config_model = config[modelname] - # determine the "outfile" and "imgfile" - outfile = config.get("outfile") - outfile = config_model.get("outfile", outfile) - imgfile = config.get("imgfile") - imgfile = config_model.get("imgfile", imgfile) - - # SBP fitting model - model = make_model(config, modelname=modelname) - - # sbp data and fit object - sbpdata = np.loadtxt(config["sbpfile"]) - sbpfit = SbpFit(model=model, xdata=sbpdata[:, 0], xerr=sbpdata[:, 1], - ydata=sbpdata[:, 2], yerr=sbpdata[:, 3], - xunit=config.get("unit", "pix")) - sbpfit.set_source(name=config["name"], obsid=config.get("obsid"), - r500_pix=config.get("r500_pix"), r500_kpc=config.get("r500_kpc")) - - # apply data range ignorance - if "ignore" in config.keys(): - for ig in config.as_list("ignore"): - xmin, xmax = map(float, ig.split("-")) - sbpfit.ignore_data(xmin=xmin, xmax=xmax) - if "ignore_r500" in config.keys(): - for ig in config.as_list("ignore_r500"): - xmin, xmax = map(float, ig.split("-")) - sbpfit.ignore_data(xmin=xmin, xmax=xmax, unit="r500") - - # apply additional data range ignorance specified within model section - if "ignore" in config_model.keys(): - for ig in config_model.as_list("ignore"): - xmin, xmax = map(float, ig.split("-")) - sbpfit.ignore_data(xmin=xmin, xmax=xmax) - if "ignore_r500" in config_model.keys(): - for ig in config_model.as_list("ignore_r500"): - xmin, xmax = map(float, ig.split("-")) - sbpfit.ignore_data(xmin=xmin, xmax=xmax, unit="r500") - - # fit and calculate confidence intervals - sbpfit.fit() - sbpfit.calc_ci() - sbpfit.report() - with open(outfile, "w") as ofile: - sbpfit.report(outfile=ofile) - - # make and save a plot - fig = Figure(figsize=(10, 8)) - canvas = FigureCanvas(fig) - ax = fig.add_subplot(111) - sbpfit.plot(ax=ax, fig=fig, r500_axis=True) - fig.savefig(imgfile, dpi=150) - - -if __name__ == "__main__": - main() - -# vim: set ts=4 sw=4 tw=0 fenc=utf-8 ft=python: # -- cgit v1.2.2