path: root/astro/fitting/fit_betasbp_cut.py

#!/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()