diff options
Diffstat (limited to 'astro/fitting')
| -rwxr-xr-x | astro/fitting/fit_betasbp_cut.py | 458 |
1 files changed, 0 insertions, 458 deletions
diff --git a/astro/fitting/fit_betasbp_cut.py b/astro/fitting/fit_betasbp_cut.py deleted file mode 100755 index 2d2931e..0000000 --- a/astro/fitting/fit_betasbp_cut.py +++ /dev/null @@ -1,458 +0,0 @@ -#!/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() - |