Updated astro tools.

* fit_betasbp_cut.py
  - fixed 'dof' calculation with 'n-p-1'
  - fixed 'par_eps' calculation/configuration
  - 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'
* marx_pntsrc.py
  Run MARX simulation on a given list of point sources, merge the
  output simulation results, and finally convert into FITS image.
* blanksky_add_time.py
  Add a time column for the chandra blanksky event file.
  The time data are generated with a uniform distribution
  between TSTART and TSTOP.

3 files changed, 188 insertions, 12 deletions

diff --git a/astro/chandra/blanksky_add_time.py b/astro/chandra/blanksky_add_time.py
new file mode 100755
index 0000000..a4d0cb6
--- /dev/null
+++ b/astro/chandra/blanksky_add_time.py
@@ -0,0 +1,73 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+#
+# Aaron LI
+# 2015/06/16
+
+"""
+Add a time column for the chandra blanksky event file.
+The time data are generated with a uniform distribution
+between TSTART and TSTOP.
+"""
+
+__version__ = "0.1.0"
+__date__ = "2015/06/16"
+
+import sys
+import argparse
+
+import numpy as np
+try:
+    from astropy.io import fits
+except ImportError:
+    try:
+        import pyfits as fits
+    except ImportError:
+        raise ImportError("cannot import 'astropy.io.fits' or 'pyfits'")
+
+
+def add_time_column(fitsfile, blockname="EVENTS"):
+    """
+    Add a time column to the specified block of the input fits file.
+    The time data are generated with a uniform distribution
+    between TSTART and TSTOP.
+
+    Return:
+        A fits object with the new time column.
+    """
+    if isinstance(fitsfile, str):
+        fitsfile = fits.open(fitsfile)
+    table = fitsfile[blockname]
+    tstart = table.header["TSTART"]
+    tstop = table.header["TSTOP"]
+    counts = len(table.data)
+    time_data = np.random.uniform(tstart, tstop, counts)
+    time_col = fits.Column(name="time", format="1D", unit="s", array=time_data)
+    newtable = fits.BinTableHDU.from_columns(
+            fits.ColDefs([time_col]) + table.columns)
+    fitsfile[blockname].data = newtable.data
+    return fitsfile
+
+
+def main():
+    parser = argparse.ArgumentParser(
+            description="Add a time column for Chandra blanksky event file.")
+    parser.add_argument("-V", "--version", action="version",
+            version="%(prog)s " + "%s (%s)" % (__version__, __date__))
+    parser.add_argument("infile", help="input chandra blanksky file")
+    parser.add_argument("outfile", nargs="?", default=None,
+            help="modified blanksky file. IN-PLACE modification if omitted.")
+    parser.add_argument("-C", "--clobber", dest="clobber",
+            action="store_true", help="overwrite output file if exists")
+    args = parser.parse_args()
+
+    newfits = add_time_column(args.infile)
+    if args.outfile:
+        newfits.writeto(args.outfile, clobber=args.clobber)
+    else:
+        newfits.writeto(args.infile, clobber=True)
+
+
+if __name__ == "__main__":
+    main()
+
diff --git a/astro/fitting/fit_betasbp_cut.py b/astro/fitting/fit_betasbp_cut.py
index 9dcf730..2d2931e 100755
--- a/astro/fitting/fit_betasbp_cut.py
+++ b/astro/fitting/fit_betasbp_cut.py
@@ -11,6 +11,9 @@
 # 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
@@ -33,7 +36,7 @@
 #   * Support read model initial parameter values from input file.
 #
 # TODO:
-# * calculate fitting parameter's standard deviation
+# * 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
@@ -42,7 +45,7 @@
 
 from __future__ import print_function, division
 
-__version__ = "0.5.0"
+__version__ = "0.5.1"
 __date__ = "2015/06/07"
 
 import numpy as np
@@ -108,7 +111,7 @@ def fit_model(func, xdata, ydata, yerrdata, p0):
     # the function evaluated at the output parameters
     fvec = lambda x: func(x, *popt)
     # degree of freedom
-    dof = len(xdata) - len(popt)
+    dof = len(xdata) - len(popt) - 1
     # chi squared
     chisq = np.sum(((ydata - fvec(xdata)) / yerrdata) ** 2)
     # one standard deviation errors on the parameters
@@ -122,10 +125,12 @@ def fit_model(func, xdata, ydata, yerrdata, p0):
     return (popt, infodict)
 
 
-def fit_model_bounds(func, xdata, ydata, yerrdata, p0=None, bounds=None,
-        options=None):
+def fit_model_bounds(func, xdata, ydata, yerrdata, p0=None,
+        bounds=None, options=None):
     """
-    Fit the provided data with the beta model.
+    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
@@ -142,12 +147,12 @@ def fit_model_bounds(func, xdata, ydata, yerrdata, p0=None, bounds=None,
             * perr: one standard deviation errors on the parameters
     """
     # objective function to be minimized, required format of 'f(x, *args)'
-    f = lambda x: calc_chisq(func, xdata, ydata, yerrdata, *x)
+    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)
+    #print("DEBUG: minimization results:\n", res, file=sys.stderr)
 
     # check minimization results
     if not res.success:
@@ -157,7 +162,7 @@ def fit_model_bounds(func, xdata, ydata, yerrdata, p0=None, bounds=None,
     # the function evaluated at the output parameters
     fvec = lambda x: func(x, *popt)
     # degree of freedom
-    dof = len(xdata) - len(popt)
+    dof = len(xdata) - len(popt) - 1
     # chi squared
     chisq = res.fun
     # one standard deviation errors on the parameters
@@ -408,12 +413,15 @@ def main():
     # initial values
     par_0 = np.array([s0_0, rc_0, beta_0, c_0])
     # parameter bounds
-    par_bounds = [(s0_lower, s0_upper), (rc_lower, rc_upper),
-                  (beta_lower, beta_upper), (c_lower, c_upper)]
+    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 = par_0 * 1e-4
+    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
diff --git a/astro/marx/marx_pntsrc.py b/astro/marx/marx_pntsrc.py
new file mode 100755
index 0000000..31ffe49
--- /dev/null
+++ b/astro/marx/marx_pntsrc.py
@@ -0,0 +1,95 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+#
+# Aaron LI
+# 2015/06/16
+
+"""
+Run MARX simulation on a given list of point sources, merge the
+output simulation results, and finally convert into FITS image.
+"""
+
+__version__ = "0.1.0"
+__date__ = "2015/06/16"
+
+import sys
+import argparse
+import subprocess
+import re
+
+
+def marx_pntsrc(pfile, ra, dec, flux, outdir):
+    """
+    Run MARX simulation for the provided point source.
+    """
+    cmd = "marx @@%(pfile)s SourceRA=%(ra)s " % {"pfile": pfile, "ra": ra} + \
+            "SourceDEC=%(dec)s SourceFlux=%(flux)s OutputDir=%(outdir)s" % \
+            {"dec": dec, "flux": flux, "outdir": outdir}
+    print("CMD: %s" % cmd, file=sys.stderr)
+    subprocess.call(cmd, shell=True)
+
+
+def marxcat(indirs, outdir):
+    """
+    Concatenate a list of MARX simulation results.
+    """
+    if isinstance(indirs, list):
+        marxdirs = " ".join(indirs)
+    elif isinstance(indirs, str):
+        marxdirs = indirs
+    else:
+        raise ValueError("invalid indirs type: %s" % indirs)
+    cmd = "marxcat %(marxdirs)s %(outdir)s" % \
+            {"marxdirs": marxdirs, "outdir": outdir}
+    print("CMD: %s" % cmd, file=sys.stderr)
+    subprocess.call(cmd, shell=True)
+
+
+def marx2fits(indir, outfile, params=""):
+    """
+    Convert the results of MARX simulation into FITS image.
+    """
+    cmd = "marx2fits %(params)s %(indir)s %(outfile)s" % \
+            {"params": params, "indir": indir, "outfile": outfile}
+    print("CMD: %s" % cmd, file=sys.stderr)
+    subprocess.call(cmd, shell=True)
+
+
+def main():
+    parser = argparse.ArgumentParser(
+            description="Run MARX on a given list of point sources")
+    parser.add_argument("-V", "--version", action="version",
+            version="%(prog)s " + "%s (%s)" % (__version__, __date__))
+    parser.add_argument("pfile", help="marx paramter file")
+    parser.add_argument("srclist", help="point source list file")
+    parser.add_argument("outprefix", help="prefix of output directories")
+    args = parser.parse_args()
+
+    outdirs = []
+    i = 0
+    for line in open(args.srclist, "r"):
+        if re.match(r"^\s*$", line):
+            # skip blank line
+            continue
+        elif re.match(r"^\s*#", line):
+            # skip comment line
+            continue
+        i += 1
+        ra, dec, flux = map(float, line.split())
+        print("INFO: ra = %g, dec = %g, flux = %g" % (ra, dec, flux),
+                file=sys.stderr)
+        outdir = "%sp%02d" % (args.outprefix, i)
+        print("INFO: outdir = %s" % outdir, file=sys.stderr)
+        outdirs.append(outdir)
+        marx_pntsrc(args.pfile, ra, dec, flux, outdir)
+    # merge results
+    merged = args.outprefix + "merged"
+    marxcat(outdirs, merged)
+    # convert to FITS image
+    merged_fits = merged + ".fits"
+    marx2fits(merged, merged_fits)
+
+
+if __name__ == "__main__":
+    main()
+