aboutsummaryrefslogtreecommitdiffstats
path: root/python/fit_sbp.py
diff options
context:
space:
mode:
Diffstat (limited to 'python/fit_sbp.py')
-rwxr-xr-xpython/fit_sbp.py807
1 files changed, 807 insertions, 0 deletions
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 = <NAME>
+obsid = <OBSID>
+r500_pix = <R500_PIX>
+r500_kpc = <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: #