diff options
| -rwxr-xr-x | deproject_sbp.py | 204 | ||||
| -rw-r--r-- | fitting_models.py | 206 |
2 files changed, 211 insertions, 199 deletions
diff --git a/deproject_sbp.py b/deproject_sbp.py index 313e586..3e60556 100755 --- a/deproject_sbp.py +++ b/deproject_sbp.py @@ -2,9 +2,12 @@ # # Weitian LI # Created: 2016-06-10 -# Updated: 2016-06-25 +# Updated: 2016-06-26 # # Change logs: +# 2016-06-26: +# * Split classes 'FitModel', 'ABModel' and 'PLCModel' into separate +# module 'fitting_models.py' # 2016-06-25: # * Use 'InterpolatedUnivariateSpline' instead of 'interp1d' # 2016-06-24: @@ -155,7 +158,6 @@ import numpy as np import pandas as pd import scipy.optimize as optimize import scipy.interpolate as interpolate -import lmfit import matplotlib.pyplot as plt from matplotlib.backends.backend_agg import FigureCanvasAgg as FigureCanvas from matplotlib.figure import Figure @@ -163,207 +165,11 @@ from configobj import ConfigObj from astro_params import AstroParams, ChandraPixel from projection import Projection +from fitting_models import ABModel, PLCModel plt.style.use("ggplot") -class FitModel: - """ - Base/Meta class for model fitting, with data and parameters scaling. - """ - name = "" - params = lmfit.Parameters() - # optimization method - fit_method = "lbfgsb" - # whether the 'ydata' and 'yerr' to be scaled in order to reduce - # the dynamical range for a more stable fitting - scale = False - scale_factor = 1.0 - - def __init__(self, fit_method="lbfgsb", params=None, scale=True): - self.fit_method = fit_method - if params is not None: - self.load_params(params) - self.scale = scale - - @staticmethod - def model(x, params): - pass - - def f(self, x): - return self.model(x, self.params) * self.scale_factor - - def load_data(self, xdata, ydata=None, xerr=None, yerr=None, - update_params=False): - if xdata.ndim == 2 and xdata.shape[1] == 4: - # 4-column data - self.xdata = xdata[:, 0].copy() - self.xerr = xdata[:, 1].copy() - self.ydata = xdata[:, 2].copy() - self.yerr = xdata[:, 3].copy() - else: - self.xdata = np.array(xdata) - self.ydata = np.array(ydata) - self.xerr = np.array(xerr) - self.yerr = np.array(yerr) - self.scale_data(update_params=update_params) - - def scale_data(self, update_params=False): - """ - Scale the ydata and yerr to reduce their dynamical ranges, - for a more stable model fitting. - """ - if self.scale: - y_min = np.min(self.ydata) - y_max = np.max(self.ydata) - self.scale_factor = np.exp(np.log(y_min*y_max) / 2) - self.ydata /= self.scale_factor - self.yerr /= self.scale_factor - if update_params: - self.scale_params() - - def scale_params(self): - """ - Scale the paramters' min/max values accordingly. - """ - pass - - def f_residual(self, params): - if self.yerr is None: - return self.model(self.xdata, params) - self.ydata - else: - return (self.model(self.xdata, params) - self.ydata) / self.yerr - - def fit(self, method=None): - if method is None: - method = self.fit_method - self.fitter = lmfit.Minimizer(self.f_residual, self.params) - self.fitted = self.fitter.minimize(method=method) - self.load_params(self.fitted.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 dump_params(self, serialize=True): - """ - Dump the current values/settings for all model parameters, - and these dumped results can be later loaded by 'load_params()'. - """ - if serialize: - return self.params.dumps() - else: - return self.params.copy() - - def load_params(self, params): - """ - Load the provided parameters values/settings. - """ - if isinstance(params, lmfit.parameter.Parameters): - self.params = params.copy() - else: - p = lmfit.parameter.Parameters() - p.loads(params) - self.params = p - - -class ABModel(FitModel): - """ - AB model is a modified beta model, which can roughly fit both - centrally peaked and cored models, e.g., central excess emission. - - This model is used here to constrain the deprojected 3D gas density - profile, in order to require it is smooth enough. - - References: - [1] Pratt & Arnaud, 2002, A&A, 394, 375; eq.(2) - [2] ref.[1], eq.(10) - """ - name = "AB model" - # model parameters - params = lmfit.Parameters() - params.add_many( # (name, value, vary, min, max, expr) - ("A", 1.0e-9, True, 0.0, 1.0e-5, None), - ("alpha", 0.7, True, 0.1, 1.1, None), - ("rc", 30.0, True, 1.0, 1.0e4, None), - ("beta", 0.7, True, 0.3, 1.1, None)) - - def scale_params(self): - A_min = 1.0 - A_max = np.max(self.ydata) - self.set_param("A", value=(A_min+A_max)*0.5, - min=A_min, max=A_max) - - @staticmethod - def model(x, params): - parvals = params.valuesdict() - A = parvals["A"] - alpha = parvals["alpha"] - rc = parvals["rc"] - beta = parvals["beta"] - return (A * np.power(x/rc, -alpha) * - np.power((1 + (x/rc)**2), -1.5*beta + 0.5*alpha)) - - -class PLCModel(FitModel): - """ - PLC model consists of a powerlaw and an constant, that is used - to fit/approximate the outer SBP. - Therefore, the fitted constant is used to subtract the uniform - background from the SBP, and the fitted powerlaw index is used - to extrapolate the SBP in order to mitigate the deprojection - errors due to FoV limit. - - NOTE: - I think the uniform background (i.e., by fitting the whole or - core-excluded SBP) should be subtracted from the SBP first, then - adopt this PLCModel to fit the outer part of SBP, with the 'bkg' - parameter fixed at zero. - """ - name = "PLC model" - # model parameters - params = lmfit.Parameters() - params.add_many( # (name, value, vary, min, max, expr) - ("A", 1.0e-9, True, 0.0, 1.0e-5, None), - ("rmin", 30.0, False, 1.0, 1.0e4, None), - ("alpha", 1.6, True, 0.4, 2.8, None), - ("bkg", 0.0, False, 0.0, 1.0e-5, None)) - - def load_data(self, xdata, ydata=None, xerr=None, yerr=None, - update_params=False): - super().load_data(xdata=xdata, ydata=ydata, xerr=xerr, yerr=yerr, - update_params=update_params) - self.set_param("rmin", value=np.min(xdata), vary=False) - - def scale_params(self): - ymin = np.min(self.ydata) - ymax = np.max(self.ydata) - self.set_param("A", value=ymax, min=ymax/10.0, max=ymax*10.0) - self.set_param("bkg", value=ymin, min=0.0, max=ymin) - - @staticmethod - def model(x, params): - parvals = params.valuesdict() - A = parvals["A"] - rmin = parvals["rmin"] - alpha = parvals["alpha"] - bkg = parvals["bkg"] - return A * np.power(x/rmin, -alpha) + bkg - - class DeprojectSBP: """ Deproject the observed SBP to derive the 3D emission measure (EM) diff --git a/fitting_models.py b/fitting_models.py new file mode 100644 index 0000000..a864a92 --- /dev/null +++ b/fitting_models.py @@ -0,0 +1,206 @@ +# -*- mode: python; coding: utf-8 -*- +# +# Weitan LI +# Created: 2016-06-26 +# Updated: 2016-06-26 +# + +import numpy as np +import lmfit + + +class FittingModel: + """ + Base/Meta class for model fitting, with data and parameters scaling. + """ + name = "" + params = lmfit.Parameters() + # optimization method + fit_method = "lbfgsb" + # whether the 'ydata' and 'yerr' to be scaled in order to reduce + # the dynamical range for a more stable fitting + scale = False + scale_factor = 1.0 + + def __init__(self, fit_method="lbfgsb", params=None, scale=True): + self.fit_method = fit_method + if params is not None: + self.load_params(params) + self.scale = scale + + @staticmethod + def model(x, params): + pass + + def f(self, x): + return self.model(x, self.params) * self.scale_factor + + def load_data(self, xdata, ydata=None, xerr=None, yerr=None, + update_params=False): + if xdata.ndim == 2 and xdata.shape[1] == 4: + # 4-column data + self.xdata = xdata[:, 0].copy() + self.xerr = xdata[:, 1].copy() + self.ydata = xdata[:, 2].copy() + self.yerr = xdata[:, 3].copy() + else: + self.xdata = np.array(xdata) + self.ydata = np.array(ydata) + self.xerr = np.array(xerr) + self.yerr = np.array(yerr) + self.scale_data(update_params=update_params) + + def scale_data(self, update_params=False): + """ + Scale the ydata and yerr to reduce their dynamical ranges, + for a more stable model fitting. + """ + if self.scale: + y_min = np.min(self.ydata) + y_max = np.max(self.ydata) + self.scale_factor = np.exp(np.log(y_min*y_max) / 2) + self.ydata /= self.scale_factor + self.yerr /= self.scale_factor + if update_params: + self.scale_params() + + def scale_params(self): + """ + Scale the paramters' min/max values accordingly. + """ + pass + + def f_residual(self, params): + if self.yerr is None: + return self.model(self.xdata, params) - self.ydata + else: + return (self.model(self.xdata, params) - self.ydata) / self.yerr + + def fit(self, method=None): + if method is None: + method = self.fit_method + self.fitter = lmfit.Minimizer(self.f_residual, self.params) + self.fitted = self.fitter.minimize(method=method) + self.load_params(self.fitted.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 dump_params(self, serialize=True): + """ + Dump the current values/settings for all model parameters, + and these dumped results can be later loaded by 'load_params()'. + """ + if serialize: + return self.params.dumps() + else: + return self.params.copy() + + def load_params(self, params): + """ + Load the provided parameters values/settings. + """ + if isinstance(params, lmfit.parameter.Parameters): + self.params = params.copy() + else: + p = lmfit.parameter.Parameters() + p.loads(params) + self.params = p + + +class ABModel(FittingModel): + """ + AB model is a modified beta model, which can roughly fit both + centrally peaked and cored models, e.g., central excess emission. + + This model is used here to constrain the deprojected 3D gas density + profile, in order to require it is smooth enough. + + References: + [1] Pratt & Arnaud, 2002, A&A, 394, 375; eq.(2) + [2] Croston et al. 2006, A&A, 459, 1007-1019; eq.(10) + """ + name = "AB model" + # model parameters + params = lmfit.Parameters() + params.add_many( # (name, value, vary, min, max, expr) + ("A", 1.0e-9, True, 0.0, 1.0e-5, None), + ("alpha", 0.7, True, 0.1, 1.1, None), + ("rc", 30.0, True, 1.0, 1.0e4, None), + ("beta", 0.7, True, 0.3, 1.1, None)) + + def scale_params(self): + A_min = 1.0 + A_max = np.max(self.ydata) + self.set_param("A", value=(A_min+A_max)*0.5, + min=A_min, max=A_max) + + @staticmethod + def model(x, params): + parvals = params.valuesdict() + A = parvals["A"] + alpha = parvals["alpha"] + rc = parvals["rc"] + beta = parvals["beta"] + return (A * np.power(x/rc, -alpha) * + np.power((1 + (x/rc)**2), -1.5*beta + 0.5*alpha)) + + +class PLCModel(FittingModel): + """ + PLC model consists of a powerlaw and an constant, that is used + to fit/approximate the outer SBP. + Therefore, the fitted constant is used to subtract the uniform + background from the SBP, and the fitted powerlaw index is used + to extrapolate the SBP in order to mitigate the deprojection + errors due to FoV limit. + + NOTE: + I think the uniform background (i.e., by fitting the whole or + core-excluded SBP) should be subtracted from the SBP first, then + adopt this PLCModel to fit the outer part of SBP, with the 'bkg' + parameter fixed at zero. + """ + name = "PLC model" + # model parameters + params = lmfit.Parameters() + params.add_many( # (name, value, vary, min, max, expr) + ("A", 1.0e-9, True, 0.0, 1.0e-5, None), + ("rmin", 30.0, False, 1.0, 1.0e4, None), + ("alpha", 1.6, True, 0.4, 2.8, None), + ("bkg", 0.0, False, 0.0, 1.0e-5, None)) + + def load_data(self, xdata, ydata=None, xerr=None, yerr=None, + update_params=False): + super().load_data(xdata=xdata, ydata=ydata, xerr=xerr, yerr=yerr, + update_params=update_params) + self.set_param("rmin", value=np.min(xdata), vary=False) + + def scale_params(self): + ymin = np.min(self.ydata) + ymax = np.max(self.ydata) + self.set_param("A", value=ymax, min=ymax/10.0, max=ymax*10.0) + self.set_param("bkg", value=ymin, min=0.0, max=ymin) + + @staticmethod + def model(x, params): + parvals = params.valuesdict() + A = parvals["A"] + rmin = parvals["rmin"] + alpha = parvals["alpha"] + bkg = parvals["bkg"] + return A * np.power(x/rmin, -alpha) + bkg |