diff options
Diffstat (limited to 'astro/fitting/fit_betasbp_cut.py')
| -rwxr-xr-x | astro/fitting/fit_betasbp_cut.py | 458 |
1 files changed, 458 insertions, 0 deletions
diff --git a/astro/fitting/fit_betasbp_cut.py b/astro/fitting/fit_betasbp_cut.py new file mode 100755 index 0000000..2d2931e --- /dev/null +++ b/astro/fitting/fit_betasbp_cut.py @@ -0,0 +1,458 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +# +# To fitting the given SBP data with the following beta model: +# s = s0 * pow((1.0+(r/rc)^2), 0.5-3*beta) + c +# And this tool supports the following two requirements for the fitting: +# (1) ignore the specified number of inner-most data points; +# (2) ignore the data points whose radius value less than the given value. +# +# Aaron LI +# 2015/05/29 +# +# Changelogs: +# v0.5.1, 2015/06/07, Aaron LI +# * fixed 'dof' calculation with 'n-p-1' +# * fixed 'par_eps' calculation/configuration +# v0.5.0, 2015/06/07, Aaron LI +# * added 'fit_model_bounds' using 'scipy.optimize.minimize' to +# perform function minimization with bounds +# * split the data cut section to function 'cut_data' +# * added argument 'options' to 'fit_model_bounds' +# v0.4.0, 2015/06/06, Aaron LI +# * replace getopt with 'argparse' +# * added 'get_parameter' to process model parameter initial value and bounds +# * support read parameter bounds from input file +# * added options '--s0', '--rc', '--beta', '--const' to get paramter +# initial values and bounds +# * renamed 'fit_beta_model' to 'fit_model', and added argument 'func' to +# support other models +# v0.3.0, 2015/06/02, Aaron LI +# * can output chi-squared and dof values +# * can output one standard deviation errors on the parameters +# v0.2.0, 2015/05/30: +# * Added option '-n' / '--no-radius' to ignore radius less than the +# given value. +# * Support read model initial parameter values from input file. +# +# TODO: +# * calculate fitting parameter's confidence interval / confidence bounds +# * to normalize fitting paramters to be the same order of magnitude +# for better minimization +# Ref: http://stackoverflow.com/questions/21369139/normalization-for-optimization-in-python +# + + +from __future__ import print_function, division + +__version__ = "0.5.1" +__date__ = "2015/06/07" + +import numpy as np +from scipy.optimize import curve_fit, minimize + +import os +import sys +import re +import argparse + + +# modes of to cut data +CUT_POINT = 'CUT_POINT' +CUT_RADIUS = 'CUT_RADIUS' +# default minimize method +MINIMIZE_METHOD = 'L-BFGS-B' + + +def beta_model(r, s0, rc, beta, c): + """ + SBP beta model, with a constant background. + """ + return s0 * np.power((1.0+(r/rc)**2), 0.5-3*beta) + c + + +def calc_chisq(func, xdata, ydata, yerrdata, *args): + """ + Calculate the chi-squared values for the given function according + to the provided data points. + + Arguments: + xdata: x values of data points + ydata: y values of data points + yerrdata: y standard deviation values + args: additional arguments for 'func' + + Return: + chi-squared value + """ + xdata = np.array(xdata) + ydata = np.array(ydata) + yerrdata = np.array(yerrdata) + return np.sum(((ydata - func(xdata, *args)) / yerrdata) ** 2) + + +def fit_model(func, xdata, ydata, yerrdata, p0): + """ + Fit the provided data with the beta model. + + Arguments: + p0: initial values for the parameters of beta model + + Return: + (popt, infodict) + popt: optimal values for the parameters + infodict: + * fvec: the function evaluated at the output parameters + * dof: degree of freedom + * chisq: chi-squared + * perr: one standard deviation errors on the parameters + """ + popt, pcov = curve_fit(func, xdata, ydata, p0=p0, sigma=yerrdata) + # the function evaluated at the output parameters + fvec = lambda x: func(x, *popt) + # degree of freedom + dof = len(xdata) - len(popt) - 1 + # chi squared + chisq = np.sum(((ydata - fvec(xdata)) / yerrdata) ** 2) + # one standard deviation errors on the parameters + perr = np.sqrt(np.diag(pcov)) + infodict = { + 'fvec': fvec, + 'dof': dof, + 'chisq': chisq, + 'perr': perr + } + return (popt, infodict) + + +def fit_model_bounds(func, xdata, ydata, yerrdata, p0=None, + bounds=None, options=None): + """ + Fit the provided data with the beta model, and apply parameter + bounds requirements. + Using 'scipy.optimize.minimize'. + + Arguments: + p0: initial values for the parameters of beta model + bounds: (min, max) pairs for each parameter bound + options: a dict of solver options (=> minimize: options) + + Return: + (popt, infodict) + popt: optimal values for the parameters + infodict: + * fvec: the function evaluated at the output parameters + * dof: degree of freedom + * chisq: chi-squared + * perr: one standard deviation errors on the parameters + """ + # objective function to be minimized, required format of 'f(x, *args)' + f = lambda p: calc_chisq(func, xdata, ydata, yerrdata, *p) + # minimize the given function using 'scipy.optimize.minimize' with bounds + res = minimize(f, p0, method=MINIMIZE_METHOD, bounds=bounds, + options=options) + popt = res.x + #print("DEBUG: minimization results:\n", res, file=sys.stderr) + + # check minimization results + if not res.success: + print("*** WARNING: minimization exited with error: ***\n" + \ + "*** %s ***" % res.message, file=sys.stderr) + + # the function evaluated at the output parameters + fvec = lambda x: func(x, *popt) + # degree of freedom + dof = len(xdata) - len(popt) - 1 + # chi squared + chisq = res.fun + # one standard deviation errors on the parameters + perr = popt * 0.0 # FIXME + infodict = { + 'fvec': fvec, + 'dof': dof, + 'chisq': chisq, + 'perr': perr + } + return (popt, infodict) + + +def cut_data(xdata, ydata, yerrdata, cutmode=CUT_POINT, cutvalue=0): + """ + Cut the given data with the provided cutmode and cutvalue, + return the cut data. + + Arguments: + xdata, ydata, yerrdata: input data (x, y, yerr) + cutmode: 'point' / 'radius'; ignore data by number of data points, + or by radius value less than the given value + cutvalue: the cut limit + + Return: + (xdata_cut, ydata_cut, yerrdata_cut) + """ + if cutmode == CUT_POINT: + xdata_cut = xdata[cutvalue:] + ydata_cut = ydata[cutvalue:] + yerrdata_cut = yerrdata[cutvalue:] + elif cutmode == CUT_RADIUS: + ii = xdata >= cutvalue + xdata_cut = xdata[ii] + ydata_cut = ydata[ii] + yerrdata_cut = yerrdata[ii] + else: + raise ValueError('Unknown cut mode: %s' % cutmode) + return (xdata_cut, ydata_cut, yerrdata_cut) + + +def get_parameter(pstring): + """ + Process the parameter string of the following format, and return + the initial value, lower limit, and upper limit. + "init_value" + "init_value lower upper" + "init_value,lower,upper" + If want to ignore the lower/upper limit, use 'None' (case-insensitive), + and the None is returned. + """ + parameters = pstring.replace(',', ' ').split() + if len(parameters) == 1: + init_value = float(parameters[0]) + return (init_value, None, None) + elif len(parameters) == 3: + init_value = float(parameters[0]) + if parameters[1].upper() == 'NONE': + lower_value = None + else: + lower_value = float(parameters[1]) + if parameters[2].upper() == 'NONE': + upper_value = None + else: + upper_value = float(parameters[2]) + return (init_value, lower_value, upper_value) + else: + raise ValueError('Invalid parameter format: %s' % pstring) + + +def main(): + # options + infile = None + outfilename= None + cutmode = CUT_POINT # ignore data by number of data points + cutvalue = 0 # do not ignore any data by default + # initial values for the four parameters of the beta model + s0_0 = 1.0e-7 + rc_0 = 10.0 + beta_0 = 0.6 + c_0 = 0.0 + # default bounds for the four parameters + s0_lower, s0_upper = None, None + rc_lower, rc_upper = None, None + beta_lower, beta_upper = None, None + c_lower, c_upper = None, None + + # parser for command line options and arguments + parser = argparse.ArgumentParser( + description="Fitting provided data with the beta model.", + epilog="Version: %s (%s)" % (__version__, __date__)) + parser.add_argument("-v", "--verbose", dest="verbose", + action="store_true", default=False, + help="show verbose/debug information (False)") + parser.add_argument("-V", "--version", action="version", + version="%(prog)s " + "%s (%s)" % (__version__, __date__)) + parser.add_argument("-i", "--infile", + dest="infile", required=True, + help="""input data file with the following 4 or 3 columns: + [radius radius_err brightness brightness_err], + [radius brightness brightness_err]. + Note: The initial values and lower/upper limits + for the beta models can also be provided with the + following syntax: + # s0_0 = init_value, lower_limit, upper_limit + # rc_0 = init_value, lower_limit, upper_limit + # beta_0 = init_value, lower_limit, upper_limit + # c_0 = init_value, lower_limit, upper_limit""") + parser.add_argument("-o", "--outfile", dest="outfilename", + help="output file to store the fitted data") + parser.add_argument("-b", "--bounds", dest="bounds", + action="store_true", default=False, + help="whether apply paramenter bounds (False)") + parser.add_argument("-c", "--cut-point", + dest="cut_point", metavar="N", type=int, default=0, + help="number of inner-most data points to be ignored") + parser.add_argument("-n", "--no-radius", + dest="cut_radius", metavar="RADIUS", type=float, default=0.0, + help="ignore data points with smaller radius") + parser.add_argument("--s0", dest="s0", + metavar="init_value,lower,upper", + help="initial value and lower/upper limits for parameter s0. " + \ + "Use 'none' (case-insensitive) to ignore the limit") + parser.add_argument("--rc", dest="rc", + metavar="init_value,lower,upper", + help="initial value and lower/upper limits for parameter rc. " + \ + "Use 'none' (case-insensitive) to ignore the limit") + parser.add_argument("--beta", dest="beta", + metavar="init_value,lower,upper", + help="initial value and lower/upper limits for parameter beta. " + \ + "Use 'none' (case-insensitive) to ignore the limit") + parser.add_argument("--const", dest="const", + metavar="init_value,lower,upper", + help="initial value and lower/upper limits for parameter const. " + \ + "Use 'none' (case-insensitive) to ignore the limit") + + args = parser.parse_args() + if args.outfilename: + outfile = open(args.outfilename, 'w') + else: + outfile = sys.stdout + # cut mode and value + if args.cut_point: + cutmode = CUT_POINT + cutvalue = args.cut_point + elif args.cut_radius: + cutmode = CUT_RADIUS + cutvalue = args.cut_radius + + if args.verbose: + print('DEBUG: apply parameter bounds: %s' % args.bounds, + file=sys.stderr) + print("DEBUG: cutmode: %s, cutvalue: %s" % (cutmode, cutvalue)) + + # input data list + r_data = [] + rerr_data = [] + s_data = [] + serr_data = [] + + # regex to match initial parameter names, blank line, and comment line + re_blank = re.compile(r'^\s*$') + re_comment = re.compile(r'^\s*#') + re_s0 = re.compile(r'^\s*#\s*s0_0\s*[:=]') + re_rc = re.compile(r'^\s*#\s*rc_0\s*[:=]') + re_beta = re.compile(r'^\s*#\s*beta_0\s*[:=]') + re_c = re.compile(r'^\s*#\s*c_0\s*[:=]') + for line in open(args.infile, 'r'): + if re_s0.match(line): + # read 's0_0': initial value for parameter 's0' + s0_pstring = re_s0.split(line)[1] + s0_0, s0_lower, s0_upper = get_parameter(s0_pstring) + elif re_rc.match(line): + # read 'rc_0': initial value for parameter 'rc' + rc_pstring = re_rc.split(line)[1] + rc_0, rc_lower, rc_upper = get_parameter(rc_pstring) + elif re_beta.match(line): + # read 'beta_0': initial value for parameter 'beta' + beta_pstring = re_beta.split(line)[1] + beta_0, beta_lower, beta_upper = get_parameter(beta_pstring) + elif re_c.match(line): + # read 'c_0': initial value for parameter 'c' + c_pstring = re_c.split(line)[1] + c_0, c_lower, c_upper = get_parameter(c_pstring) + elif re_blank.match(line): + # ignore blank line + continue + elif re_comment.match(line): + # ignore comment line + continue + else: + try: + r, rerr, s, serr = map(float, line.split()) + except ValueError: + try: + r, s, serr = map(float, line.split()) + rerr = 0.0 + except ValueError: + print('ERROR: unsupported input data format', + file=sys.stderr) + sys.exit(21) + r_data.append(r) + rerr_data.append(rerr) + s_data.append(s) + serr_data.append(serr) + + if args.verbose: + print('DEBUG: infile: s0_0 = %g (%s, %s)' % \ + (s0_0, s0_lower, s0_upper), file=sys.stderr) + print('DEBUG: infile: rc_0 = %g (%s, %s)' % \ + (rc_0, rc_lower, rc_upper), file=sys.stderr) + print('DEBUG: infile: beta_0 = %g (%s, %s)' % \ + (beta_0, beta_lower, beta_upper), file=sys.stderr) + print('DEBUG: infile: c_0 = %g (%s, %s)' % \ + (c_0, c_lower, c_upper), file=sys.stderr) + + # get parameter initial values and bounds from command line arguments + if args.s0: + s0_0, s0_lower, s0_upper = get_parameter(args.s0) + if args.rc: + rc_0, rc_lower, rc_upper = get_parameter(args.rc) + if args.beta: + beta_0, beta_lower, beta_upper = get_parameter(args.beta) + if args.const: + c_0, c_lower, c_upper = get_parameter(args.const) + + if args.verbose: + print('DEBUG: final: s0_0 = %g (%s, %s)' % \ + (s0_0, s0_lower, s0_upper), file=sys.stderr) + print('DEBUG: final: rc_0 = %g (%s, %s)' % \ + (rc_0, rc_lower, rc_upper), file=sys.stderr) + print('DEBUG: final: beta_0 = %g (%s, %s)' % \ + (beta_0, beta_lower, beta_upper), file=sys.stderr) + print('DEBUG: final: c_0 = %g (%s, %s)' % \ + (c_0, c_lower, c_upper), file=sys.stderr) + + # convert to numpy array + r_data = np.array(r_data) + rerr_data = np.array(rerr_data) + s_data = np.array(s_data) + serr_data = np.array(serr_data) + # cut data + r_data_cut, s_data_cut, serr_data_cut = cut_data(r_data, s_data, + serr_data, cutmode=cutmode, cutvalue=cutvalue) + + # model parameters + par_names = ["s0", "rc", "beta", "c"] + # initial values + par_0 = np.array([s0_0, rc_0, beta_0, c_0]) + # parameter bounds + par_bounds = np.array([(s0_lower, s0_upper), (rc_lower, rc_upper), + (beta_lower, beta_upper), (c_lower, c_upper)]) + # set eps for the parameters (required for the minimize method, + # otherwise error 'ABNORMAL_TERMINATION_IN_LNSRCH' occurs, which + # may due to the different order of magnitude of each parameters) + par_eps = np.absolute(par_0) * 1e-3 + par_eps[ par_eps<1e-15 ] = 1e-15 # set par_eps >= 1e-14 + if args.verbose: + print("DEBUG: parameters eps:\n", par_eps, file=sys.stderr) + + if args.bounds: + ## 'fit_model_bounds' to perform fitting with bounds + par_fit, infodict = fit_model_bounds(beta_model, r_data_cut, + s_data_cut, serr_data_cut, p0=par_0, bounds=par_bounds, + options={'eps': par_eps}) + else: + # 'fit_model' do not support parameter bounds + par_fit, infodict = fit_model(beta_model, r_data_cut, + s_data_cut, serr_data_cut, p0=par_0) + + fvec = infodict['fvec'] + dof = infodict['dof'] + chisq = infodict['chisq'] + perr = infodict['perr'] + + print("# beta-model fitting results:", file=outfile) + print("# s(r) = s0 * pow((1.0+(r/rc)^2), 0.5-3*beta) + c", file=outfile) + for i in range(len(par_names)): + print("# %s = %g +/- %g" % (par_names[i], par_fit[i], perr[i]), + file=outfile) + print("# chisq / dof = %g / %g = %g" % (chisq, dof, chisq/dof), + file=outfile) + print("# radius(input) brightness(fitted)", file=outfile) + for i in range(len(r_data)): + print("%g %g" % (r_data[i], fvec(r_data[i])), file=outfile) + + if args.outfilename: + outfile.close() + + +if __name__ == '__main__': + main() + |