aboutsummaryrefslogtreecommitdiffstats
diff options
context:
space:
mode:
authorAaron LI <aaronly.me@outlook.com>2016-04-21 09:48:48 +0800
committerAaron LI <aaronly.me@outlook.com>2016-04-21 09:48:48 +0800
commitc11cdf28d7b188a68969a760280df31a0566dd01 (patch)
tree4dd4b6059e373261ca5f18e8a58429ce929f94b7
parent35e480b73ea90cfd6ce115b590ea72db5175a0aa (diff)
downloadatoolbox-c11cdf28d7b188a68969a760280df31a0566dd01.tar.bz2
sbp_fit.py: support pix/kpc units; update config syntax
-rwxr-xr-x[-rw-r--r--]python/msvst_starlet.py0
-rwxr-xr-x[-rw-r--r--]python/sbp_fit.py226
2 files changed, 166 insertions, 60 deletions
diff --git a/python/msvst_starlet.py b/python/msvst_starlet.py
index b90adf9..b90adf9 100644..100755
--- a/python/msvst_starlet.py
+++ b/python/msvst_starlet.py
diff --git a/python/sbp_fit.py b/python/sbp_fit.py
index c9a67a5..db7719d 100644..100755
--- a/python/sbp_fit.py
+++ b/python/sbp_fit.py
@@ -3,9 +3,14 @@
#
# Aaron LI
# Created: 2016-03-13
-# Updated: 2016-04-05
+# Updated: 2016-04-20
#
# Changelogs:
+# 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:
@@ -25,6 +30,7 @@
#
# TODO:
# * to allow fit the outer beta component, then fix it, and fit the inner one
+# * to also plot another X axis with unit (R500) (also mark R500 value)
#
"""
@@ -37,42 +43,63 @@ or the double-beta model:
Sample config file:
-------------------------------------------------
-name = <NAME>
-obsid = <OBSID>
+name = <NAME>
+obsid = <OBSID>
+r500_pix = <R500_pixel>
+r500_kpc = <R500_kpc>
-sbpfile = sbprofile.txt
-
-model = sbeta
-outfile = sbpfit_sbeta.txt
-imgfile = sbpfit_sbeta.png
+sbpfile = sbprofile.txt
+# unit of radius: pix (default) or kpc
+unit = pixel
+# sbp model: "sbeta" or "dbeta"
+model = sbeta
#model = dbeta
-#outfile = sbpfit_dbeta.txt
-#imgfile = sbpfit_dbeta.png
-# data range to be ignored during fitting
+# output file to store the fitting results
+outfile = sbpfit.txt
+# 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]
-# name = initial, lower, upper
-s0 = 1.0e-8, 0.0, 1.0e-6
-rc = 30.0, 1.0, 1.0e4
-beta = 0.7, 0.3, 1.1
-bkg = 1.0e-9, 0.0, 1.0e-7
+# model-related options (OVERRIDE the upper level options)
+outfile = sbpfit_sbeta.txt
+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, 1.0, 1.0e4
+ #rc = 30.0, 1.0, 1.0e4, FIXED
+ beta = 0.7, 0.3, 1.1
+ bkg = 1.0e-9, 0.0, 1.0e-7
+
[dbeta]
-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, 1.0, 5.0e2
-beta2 = 0.7, 0.3, 1.1
-bkg = 1.0e-9, 0.0, 1.0e-7
+outfile = sbpfit_dbeta.txt
+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, 1.0, 5.0e2
+ beta2 = 0.7, 0.3, 1.1
+ bkg = 1.0e-9, 0.0, 1.0e-7
-------------------------------------------------
"""
-__version__ = "0.5.1"
-__date__ = "2016-04-05"
+__version__ = "0.6.0"
+__date__ = "2016-04-20"
import numpy as np
@@ -397,62 +424,104 @@ 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,
- name=None, obsid=None):
+ 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.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
- self.name = name
- self.obsid = obsid
+ 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):
+ def set_source(self, name, obsid=None, r500_pix=None, r500_kpc=None):
self.name = name
self.obsid = obsid
+ self.r500_pix = r500_pix
+ self.r500_kpc = r500_kpc
+ try:
+ self.kpc_per_pix = r500_kpc / r500_pix
+ except (TypeError, ZeroDivisionError):
+ self.kpc_per_pix = -1
- def load_data(self, xdata, ydata, xerr, yerr):
+ def load_data(self, xdata, ydata, xerr, yerr, xunit="pix"):
self.xdata = xdata
self.ydata = ydata
self.xerr = xerr
self.yerr = yerr
- self.mask = np.ones(xdata.shape, dtype=np.bool)
+ 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):
+ 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 xmin is None:
- xmin = min(self.xdata)
- if xmax is None:
- xmax = max(self.xdata)
+ 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):
+ 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 xmin is None:
- xmin = min(self.xdata)
- if xmax is None:
- xmax = max(self.xdata)
+ 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 set_residual(self):
def f_residual(params):
if self.yerr is None:
@@ -571,11 +640,10 @@ class SbpFit:
return (fig, ax)
-def make_model(config):
+def make_model(config, modelname):
"""
Make the model with parameters set according to the config.
"""
- modelname = config["model"]
if modelname == "sbeta":
# single-beta model
model = FitModelSBeta()
@@ -585,7 +653,7 @@ def make_model(config):
else:
raise ValueError("Invalid model")
# set initial values and bounds for the model parameters
- params = config.get(modelname)
+ params = config[modelname]["params"]
for p, value in params.items():
variable = True
if len(value) == 4 and value[3].upper() in ["FIXED", "FALSE"]:
@@ -604,38 +672,76 @@ def main():
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)
- # fit model
- model = make_model(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])
- sbpfit.set_source(config["name"], obsid=config.get("obsid"))
+ 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(config["outfile"], "w") as outfile:
- sbpfit.report(outfile=outfile)
+ with open(outfile, "w") as ofile:
+ sbpfit.report(outfile=ofile)
# make and save a plot
fig = Figure()
canvas = FigureCanvas(fig)
ax = fig.add_subplot(111)
sbpfit.plot(ax=ax, fig=fig)
- fig.savefig(config["imgfile"])
+ fig.savefig(imgfile)
if __name__ == "__main__":