diff options
Diffstat (limited to 'astro')
| -rwxr-xr-x | astro/chandra/blanksky_add_time.py | 73 | ||||
| -rwxr-xr-x | astro/fitting/fit_betasbp_cut.py | 32 | ||||
| -rwxr-xr-x | astro/marx/marx_pntsrc.py | 95 |
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() + |