diff options
| -rwxr-xr-x[-rw-r--r--] | python/msvst_starlet.py | 0 | ||||
| -rwxr-xr-x[-rw-r--r--] | python/sbp_fit.py | 226 |
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__": |