Renamed & updated fit_beta_model.py -> fit_betasbp_cut.py.

fit_betasbp_cut.py:
* added support ignore the data points whose radius value less than the given value.

2 files changed, 237 insertions, 133 deletions

diff --git a/astro/fitting/fit_beta_model.py b/astro/fitting/fit_beta_model.py
deleted file mode 100755
index aaf5781..0000000
--- a/astro/fitting/fit_beta_model.py
+++ /dev/null
@@ -1,133 +0,0 @@
-#!/usr/bin/env python3
-# -*- coding: utf-8 -*-
-#
-# Fitting data according to (SBP) beta model.
-#
-# Aaron LI
-# 2015/05/29
-
-
-import numpy as np
-from scipy.optimize import curve_fit
-
-import os
-import sys
-import getopt
-import re
-
-
-USAGE = """Usage:
-    %(prog)s [ -h -c -o outfile ] -i infile
-
-Required arguments:
-    -i, --infile
-        input data file with the following *four* columns:
-            r, rerr, s, serr
-    -o, --outfile
-        output file to store the fitted data
-        if not provided, then print results to screen
-
-Optional arguments:
-    -h, --help
-        print this usage
-    -c, --cut
-        accept an integer number
-        ignore the inner-most n data points
-""" % { 'prog': os.path.basename(sys.argv[0]) }
-
-
-def usage():
-    print(USAGE)
-
-
-def beta_model(r, s0, rc, beta, c):
-    """
-    SBP beta model.
-    """
-    return s0 * np.power((1.0+(r/rc)**2), 0.5-3*beta) + c
-
-
-def main():
-    infile = None
-    outfilename= None
-    cut = 0
-
-    try:
-        opts, args = getopt.getopt(sys.argv[1:], 'c:hi:o:',
-                ['cut=', 'help', 'infile=', 'outfile='])
-    except getopt.GetoptError as e:
-        print(e, file=sys.stderr)
-        usage()
-        sys.exit(2)
-
-    for opt, arg in opts:
-        if opt in ('-h', '--help'):
-            usage()
-            sys.exit(0)
-        elif opt in ('-i', '--infile'):
-            infile = arg
-        elif opt in ('-o', '--outfile'):
-            outfilename = arg
-        elif opt in ('-c', '--cut'):
-            cut = int(arg)
-        else:
-            assert False, 'unhandled option'
-
-    if not infile:
-        print('ERROR: --infile required', file=sys.stderr)
-        sys.exit(11)
-    if outfilename:
-        outfile = open(outfilename, 'w')
-    else:
-        outfile = sys.stdout
-
-    # input data list
-    r_data = []
-    rerr_data = []
-    s_data = []
-    serr_data = []
-
-    for line in open(infile, 'r'):
-        if re.match(r'^\s*#', line):
-            continue
-        if re.match(r'^\s*$', line):
-            continue
-        r, rerr, s, serr = map(float, line.split())
-        r_data.append(r)
-        rerr_data.append(rerr)
-        s_data.append(s)
-        serr_data.append(serr)
-
-    # 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)
-
-    # initial parameter values for beta model
-    s0_0 = 1.0e-7
-    rc_0 = 10.0
-    beta_0 = 0.6
-    c_0 = 0.0
-    par_0 = [s0_0, rc_0, beta_0, c_0]
-
-    par_fit, cov_fit = curve_fit(beta_model, r_data[cut:], s_data[cut:],
-            p0=par_0, sigma=serr_data[cut:])
-
-    s_fitted = beta_model(r_data, *par_fit)
-
-    print("# beta-model fitting results:", file=outfile)
-    print("# s0 * power((1.0+(r/rc)**2), 0.5-3*beta) + c", file=outfile)
-    print("# s0 = %g, rc = %g, beta = %g, c = %g" % tuple(par_fit),
-            file=outfile)
-    print("# radius(input)  brightness(fitted)", file=outfile)
-    for i in range(len(s_fitted)):
-        print("%g  %g" % (r_data[i], s_fitted[i]), file=outfile)
-
-    if not outfilename:
-        outfile.close()
-
-
-if __name__ == '__main__':
-    main()
-
diff --git a/astro/fitting/fit_betasbp_cut.py b/astro/fitting/fit_betasbp_cut.py
new file mode 100755
index 0000000..4b530f4
--- /dev/null
+++ b/astro/fitting/fit_betasbp_cut.py
@@ -0,0 +1,237 @@
+#!/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:
+# 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.
+#
+
+
+from __future__ import print_function, division
+
+import numpy as np
+from scipy.optimize import curve_fit
+
+import os
+import sys
+import getopt
+import re
+
+
+# modes of to cut data
+CUT_POINT = 'CUT_POINT'
+CUT_RADIUS = 'CUT_RADIUS'
+
+
+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 fit_beta_model(xdata, ydata, yerrdata, p0,
+        cutmode=CUT_POINT, cutvalue=0):
+    """
+    Fit the provided data with the beta model.
+
+    Arguments:
+        p0: initial values for the parameters of beta model
+        cutmode: 'point' / 'radius'; ignore data by number of data points,
+            or by radius value less than the given value
+        cutvalue: the cut limit
+
+    Return:
+        [fitted_par, fitted_cov, fitted_model_value]
+    """
+    if cutmode == CUT_POINT:
+        xdata2 = xdata[cutvalue:]
+        ydata2 = ydata[cutvalue:]
+        yerrdata2 = yerrdata[cutvalue:]
+    elif cutmode == CUT_RADIUS:
+        ii = xdata >= cutvalue
+        xdata2 = xdata[ii]
+        ydata2 = ydata[ii]
+        yerrdata2 = yerrdata[ii]
+    else:
+        raise ValueError('Unknown cut mode: %s' % cutmode)
+    par_fit, cov_fit = curve_fit(beta_model, xdata2, ydata2,
+            p0=p0, sigma=yerrdata2)
+    y_fit = beta_model(xdata, *par_fit)
+    return (par_fit, cov_fit, y_fit)
+
+
+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
+    # debug
+    verbose = False
+
+    try:
+        opts, args = getopt.getopt(sys.argv[1:], 'c:hi:n:o:v',
+                ['cut-point=', 'help', 'infile=', 'no-radius=',
+                 'outfile=', 'verbose'])
+    except getopt.GetoptError as e:
+        print(e, file=sys.stderr)
+        usage()
+        sys.exit(2)
+
+    for opt, arg in opts:
+        if opt in ('-h', '--help'):
+            usage()
+            sys.exit(0)
+        elif opt in ('-v', '--verbose'):
+            verbose = True
+        elif opt in ('-i', '--infile'):
+            infile = arg
+        elif opt in ('-o', '--outfile'):
+            outfilename = arg
+        elif opt in ('-c', '--cut-point'):
+            cutvalue = int(arg)
+            cutmode = CUT_POINT
+        elif opt in ('-n', '--no-radius'):
+            cutvalue = float(arg)
+            cutmode = CUT_RADIUS
+        else:
+            assert False, 'unhandled option'
+
+    if not infile:
+        print('ERROR: --infile required', file=sys.stderr)
+        sys.exit(11)
+    if outfilename:
+        outfile = open(outfilename, 'w')
+    else:
+        outfile = sys.stdout
+
+    # 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(infile, 'r'):
+        if re_s0.match(line):
+            # read 's0_0': initial value for parameter 's0'
+            s0_0 = float(re_s0.split(line)[1])
+        elif re_rc.match(line):
+            # read 'rc_0': initial value for parameter 'rc'
+            rc_0 = float(re_rc.split(line)[1])
+        elif re_beta.match(line):
+            # read 'beta_0': initial value for parameter 'beta'
+            beta_0 = float(re_beta.split(line)[1])
+        elif re_c.match(line):
+            # read 'c_0': initial value for parameter 'c'
+            c_0 = float(re_c.split(line)[1])
+        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)
+
+    # 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)
+
+    if verbose:
+        print('DEBUG: s0_0 = %g, rc_0 = %g, beta_0 = %g, c_0 = %g' % \
+                (s0_0, rc_0, beta_0, c_0), file=sys.stderr)
+
+    par_0 = [s0_0, rc_0, beta_0, c_0]
+
+    par_fit, cov_fit, s_fit = fit_beta_model(r_data, s_data, serr_data,
+            p0=par_0, cutmode=cutmode, cutvalue=cutvalue)
+
+    print("# beta-model fitting results:", file=outfile)
+    print("# s(r) = s0 * pow((1.0+(r/rc)^2), 0.5-3*beta) + c", file=outfile)
+    print("# s0 = %g\n# rc = %g\n# beta = %g\n# c = %g" % tuple(par_fit),
+            file=outfile)
+    print("# radius(input)  brightness(fitted)", file=outfile)
+    for i in range(len(s_fit)):
+        print("%g  %g" % (r_data[i], s_fit[i]), file=outfile)
+
+    if outfilename:
+        outfile.close()
+
+
+USAGE = """Usage:
+    %(prog)s [ -h -c -o outfile ] -i infile
+
+Required arguments:
+    -i, --infile
+        input data file with the following *four* or *three* columns:
+            r, rerr, s, serr
+            r, s, serr
+        Note: the initial values for beta model paramters can also be
+        provided with the following syntax:
+            # s0_0 = ??
+            # rc_0 = ??
+            # beta_0 = ??
+            # c_0 = ??
+    -o, --outfile
+        output file to store the fitted data
+        if not provided, then print results to screen
+
+Optional arguments:
+    -h, --help
+        print this usage
+    -c, --cut-point
+        accept an integer number (n)
+        ignore the inner-most n data points
+    -n, --no-radius
+        accept a float number (r)
+        ignore the data points whose radius is less than r
+""" % { 'prog': os.path.basename(sys.argv[0]) }
+
+
+def usage():
+    print(USAGE)
+
+
+if __name__ == '__main__':
+    main()
+
+