rename "sbp_fit.py" to "fit_sbp.py"

1 files changed, 0 insertions, 807 deletions

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     = <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: #