From bc7ed62dbb28499b565ecaa2464d37a278e8c2f4 Mon Sep 17 00:00:00 2001 From: Marcelo Alvarez Date: Mon, 5 Oct 2020 15:11:37 -0700 Subject: [PATCH 01/16] first commit --- py/desispec/scripts/specex.py | 11 ++++++++--- 1 file changed, 8 insertions(+), 3 deletions(-) diff --git a/py/desispec/scripts/specex.py b/py/desispec/scripts/specex.py index cfb23309f..8ca1e0156 100644 --- a/py/desispec/scripts/specex.py +++ b/py/desispec/scripts/specex.py @@ -17,6 +17,7 @@ from astropy.io import fits from desiutil.log import get_logger +import specex as spx modext = "so" if sys.platform == "darwin": @@ -117,7 +118,8 @@ def main(args, comm=None): specmin = int(args.specmin) nspec = int(args.nspec) bundlesize = int(args.bundlesize) - + print(bundlesize) + return specmax = specmin + nspec # Now we divide our spectra into bundles @@ -212,8 +214,10 @@ def main(args, comm=None): map(ct.addressof, arg_buffers) ] arg_pointers = (ct.POINTER(ct.c_char) * argc)(*addrlist) - retval = libspecex.cspecex_desi_psf_fit(argc, arg_pointers) - + retval=0 + #retval = libspecex.cspecex_desi_psf_fit(argc, arg_pointers) + retval = spx.specex_desi_psf_fit_main(com) + if retval != 0: comstr = " ".join(com) log.error("desi_psf_fit on process {} failed with return " @@ -233,6 +237,7 @@ def main(args, comm=None): inputs = [ "{}_{:02d}.fits".format(outroot, x) for x in bundles ] + args.disable_merge=True if args.disable_merge : log.info("don't merge") else : From addfeecd2a4a3aba5e537dd226d1989af32f6ebe Mon Sep 17 00:00:00 2001 From: Marcelo Alvarez Date: Thu, 8 Oct 2020 16:57:46 -0700 Subject: [PATCH 02/16] python-less spot-less write-less refactor --- py/desispec/scripts/specex.py | 23 ++++++++++++++++++----- 1 file changed, 18 insertions(+), 5 deletions(-) diff --git a/py/desispec/scripts/specex.py b/py/desispec/scripts/specex.py index 8ca1e0156..245c10433 100644 --- a/py/desispec/scripts/specex.py +++ b/py/desispec/scripts/specex.py @@ -118,8 +118,6 @@ def main(args, comm=None): specmin = int(args.specmin) nspec = int(args.nspec) bundlesize = int(args.bundlesize) - print(bundlesize) - return specmax = specmin + nspec # Now we divide our spectra into bundles @@ -214,10 +212,25 @@ def main(args, comm=None): map(ct.addressof, arg_buffers) ] arg_pointers = (ct.POINTER(ct.c_char) * argc)(*addrlist) - retval=0 - #retval = libspecex.cspecex_desi_psf_fit(argc, arg_pointers) - retval = spx.specex_desi_psf_fit_main(com) + # old way + # retval = libspecex.cspecex_desi_psf_fit(argc, arg_pointers) + + # new way + opts = spx.PyOptions() # input options + pyio = spx.PyIO() # IO options and methods + pypr = spx.PyPrior() # Gaussian priors + pymg = spx.PyImage() # image data + pyps = spx.PyPSF() # psf data + + opts.parse(com) # com is list of args + pyio.check_input_psf(opts) # set booleans for input psf + pypr.deal_with_priors(opts) # set Gaussian priors + pyio.read_img_data(opts,pymg) # read image data + pyio.read_psf_data(opts,pyps,pyio) # read psf data + retval = spx.specex_desi_psf_fit_main(opts,pyio,pypr,pymg,pyps) + + # check for success if retval != 0: comstr = " ".join(com) log.error("desi_psf_fit on process {} failed with return " From 25b202689bea72109ca355b38db9873a7c299e42 Mon Sep 17 00:00:00 2001 From: Marcelo Alvarez Date: Fri, 9 Oct 2020 18:04:54 -0700 Subject: [PATCH 03/16] python-less I/O refactor --- py/desispec/scripts/specex.py | 28 ++++++++++++++++------------ 1 file changed, 16 insertions(+), 12 deletions(-) diff --git a/py/desispec/scripts/specex.py b/py/desispec/scripts/specex.py index 245c10433..b928138c5 100644 --- a/py/desispec/scripts/specex.py +++ b/py/desispec/scripts/specex.py @@ -213,7 +213,7 @@ def main(args, comm=None): arg_pointers = (ct.POINTER(ct.c_char) * argc)(*addrlist) # old way - # retval = libspecex.cspecex_desi_psf_fit(argc, arg_pointers) + # rval = libspecex.cspecex_desi_psf_fit(argc, arg_pointers) # new way opts = spx.PyOptions() # input options @@ -221,22 +221,26 @@ def main(args, comm=None): pypr = spx.PyPrior() # Gaussian priors pymg = spx.PyImage() # image data pyps = spx.PyPSF() # psf data + pyft = spx.PyFitting() # psf fitting - opts.parse(com) # com is list of args - pyio.check_input_psf(opts) # set booleans for input psf - pypr.deal_with_priors(opts) # set Gaussian priors - pyio.read_img_data(opts,pymg) # read image data - pyio.read_psf_data(opts,pyps,pyio) # read psf data + rval = opts.parse(com) # com is list of args + rval = pyio.check_input_psf(opts) # set booleans for input psf + rval = pypr.deal_with_priors(opts) # set Gaussian priors + rval = pyio.read_img_data(opts,pymg) # read image + rval = pyio.read_psf_data(opts,pyps) # read psf + rval = pyft.fit_psf(opts,pyio,pypr,pymg,pyps) # fit psf - retval = spx.specex_desi_psf_fit_main(opts,pyio,pypr,pymg,pyps) - - # check for success - if retval != 0: + # check for fit_psf success + if rval != 0: comstr = " ".join(com) - log.error("desi_psf_fit on process {} failed with return " - "value {} running {}".format(rank, retval, comstr)) + log.error("fit_psf on process {} failed with return " + "value {} running {}".format(rank, rval, comstr)) failcount += 1 + rval = pyio.write_psf_data(opts,pyps) # write psf + rval = pyio.write_spots(opts,pyps) # write spots + + if comm is not None: from mpi4py import MPI failcount = comm.allreduce(failcount, op=MPI.SUM) From dee5c706811c80896fef38511c7720ae839443b7 Mon Sep 17 00:00:00 2001 From: Marcelo Alvarez Date: Wed, 14 Oct 2020 22:19:02 -0700 Subject: [PATCH 04/16] v0 of preproc image python read --- py/desispec/scripts/specex.py | 98 +++++++++++++++++++++++++++++++++-- 1 file changed, 94 insertions(+), 4 deletions(-) diff --git a/py/desispec/scripts/specex.py b/py/desispec/scripts/specex.py index b928138c5..9205b264b 100644 --- a/py/desispec/scripts/specex.py +++ b/py/desispec/scripts/specex.py @@ -63,6 +63,89 @@ ct.POINTER(ct.POINTER(ct.c_char)) ] +######################################### +# TEMPORARY SPECEX-DESISPEC I/O FUNCTIONS + +def meta2header(meta): + header = spx.MapStringString() + + for key in meta: + mkey = meta[key] + if type(mkey) == bool: + header[key]='F' + if mkey: header[key]='T' + elif type(mkey) == str: + mstr = mkey + if len(mstr) < 8: mstr = mstr.ljust(8, ' ') + header[key]="\'"+mstr+"\'" + else: + mstr = str(mkey) + if len(mstr) > 1: + if mstr[-2:]=='.0': mstr=mstr[:-1] + header[key]=mstr + + return header + +def read_desi_ppimage_spx(opts): + import desispec.io.image + + # read images from fits file + dsmg = desispec.io.image.read_image(opts.arc_image_filename) + + # set ivar=0 to pixels with mask!=0 + dsmg.ivar[dsmg.mask!=0] = 0.0 + + # convert from astropy.io.fits.header.Header dict-like object to + # std::map object via meta2header + # function above, changing formatting to match current use in + # specex + hdr = meta2header(dsmg.meta) + + # instantiate new specex::PyImage object with arrays and header + # from the preproc fits file. this object will be passed to the + # specex::PyFitting::fit_psf routine + pymg = spx.PyImage(dsmg.pix, dsmg.ivar, + dsmg.mask,dsmg.readnoise, + hdr) + + return pymg + +def compheaders(header1,header2): + if len(header1) != len(header2): + print('lengths not same') + return 1 + + for key in header1: + if header1[key] != header2[key]: + print('values for key ',key,' are ',header1[key],' ',header2[key]) + return 1 + + return 0 + +def comparrs(arr1,arr2,tag): + + if np.array_equal(arr1,arr2): return 0 + + if len(arr1) != len(arr2): + print(tag,' array length are different') + return 1 + + diff = np.abs(arr1-arr2) + print(tag,' arrays are different') + print(' min diff',diff.min()) + print(' max diff',diff.max()) + print(' avg diff',diff.mean()) + print(' arr1 avg',arr1.mean()) + print(' arr2 avg',arr2.mean()) + print(np.shape(arr1)) + print(len(diff[diff>0])) + print(arr1[diff>0],'\n') + print(arr2[diff>0],'\n') + + return 1 + +# END TEMPORARY I/O SUPPORT FUNCTIONS +######################################### def parse(options=None): parser = argparse.ArgumentParser(description="Estimate the PSF for " @@ -216,17 +299,25 @@ def main(args, comm=None): # rval = libspecex.cspecex_desi_psf_fit(argc, arg_pointers) # new way + specex_image_io = False + opts = spx.PyOptions() # input options pyio = spx.PyIO() # IO options and methods pypr = spx.PyPrior() # Gaussian priors - pymg = spx.PyImage() # image data pyps = spx.PyPSF() # psf data pyft = spx.PyFitting() # psf fitting rval = opts.parse(com) # com is list of args - rval = pyio.check_input_psf(opts) # set booleans for input psf + rval = pyio.check_input_psf(opts) # set input psf bools rval = pypr.deal_with_priors(opts) # set Gaussian priors - rval = pyio.read_img_data(opts,pymg) # read image + + if specex_image_io: + # read preproc images into pymg with specex (harp) + rval = pyio.read_img_data(opts,pymg) + else: + # read preproc images into pymg with desispec (astropy.io.fits) + pymg = read_desi_ppimage_spx(opts) + rval = pyio.read_psf_data(opts,pyps) # read psf rval = pyft.fit_psf(opts,pyio,pypr,pymg,pyps) # fit psf @@ -240,7 +331,6 @@ def main(args, comm=None): rval = pyio.write_psf_data(opts,pyps) # write psf rval = pyio.write_spots(opts,pyps) # write spots - if comm is not None: from mpi4py import MPI failcount = comm.allreduce(failcount, op=MPI.SUM) From 1d85c77ab6b082aaf7187ff09632f1a7a87bc091 Mon Sep 17 00:00:00 2001 From: Marcelo Alvarez Date: Thu, 15 Oct 2020 15:09:43 -0700 Subject: [PATCH 05/16] added info message --- py/desispec/scripts/specex.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/py/desispec/scripts/specex.py b/py/desispec/scripts/specex.py index 9205b264b..9a9448fad 100644 --- a/py/desispec/scripts/specex.py +++ b/py/desispec/scripts/specex.py @@ -243,6 +243,8 @@ def main(args, comm=None): if rank == 0: # Print parameters + log.info("specex: io_refactor") + time.sleep(5) log.info("specex: using {} processes".format(nproc)) log.info("specex: input image = {}".format(imgfile)) log.info("specex: input PSF = {}".format(inpsffile)) From 03cfc859912e5535ee2c0457fc78e5f927eaf4e1 Mon Sep 17 00:00:00 2001 From: Marcelo Alvarez Date: Mon, 9 Nov 2020 14:29:18 -0800 Subject: [PATCH 06/16] psf I/O testing --- py/desispec/scripts/specex.py | 55 +++++++++++++++++++++++++++++++++-- 1 file changed, 52 insertions(+), 3 deletions(-) diff --git a/py/desispec/scripts/specex.py b/py/desispec/scripts/specex.py index 9a9448fad..ff6ab1361 100644 --- a/py/desispec/scripts/specex.py +++ b/py/desispec/scripts/specex.py @@ -17,7 +17,9 @@ from astropy.io import fits from desiutil.log import get_logger + import specex as spx +import specter.psf modext = "so" if sys.platform == "darwin": @@ -66,6 +68,11 @@ ######################################### # TEMPORARY SPECEX-DESISPEC I/O FUNCTIONS +def compare_psfs(specter_psf, specexpy_psf, specex_psf): + print('specter nspec',specter_psf.nspec) + print('specex deg ',specex_psf.Degree()) + return + def meta2header(meta): header = spx.MapStringString() @@ -313,6 +320,9 @@ def main(args, comm=None): rval = pyio.check_input_psf(opts) # set input psf bools rval = pypr.deal_with_priors(opts) # set Gaussian priors + fiberkeys = spx.VectorInt() + wavekeys = spx.VectorDouble() + if specex_image_io: # read preproc images into pymg with specex (harp) rval = pyio.read_img_data(opts,pymg) @@ -320,9 +330,37 @@ def main(args, comm=None): # read preproc images into pymg with desispec (astropy.io.fits) pymg = read_desi_ppimage_spx(opts) - rval = pyio.read_psf_data(opts,pyps) # read psf - rval = pyft.fit_psf(opts,pyio,pypr,pymg,pyps) # fit psf + rval = pyio.read_psf_data(opts,pyps) # read psf (specex-py) + #################################### + # BEGIN TEST PSF I/O + + #spps = specter.psf.GaussHermitePSF( + # opts.input_psf_filename) # read psf (specter) + + #rval = pyio.write_psf_data(opts,pyps) # dummy write psf + #print(spx.get_pyps_image(pymg,pyps)) + #print(len(spx.get_pyps_image(pymg,pyps))) + + #return + + #print(spx.get_pyps_image(pyps)) + #return + # sxps = spx.GaussHermitePSF() # new PSF object (specex) + + #print(opts.output_fits_filename) + #opts.output_fits_filename=opts.output_fits_filename[:-5]+'_cp.fits' + #print(opts.output_fits_filename) + #print('writing psf') + #rval = pyio.write_psf_data(opts,pyps) # write psf + #compare_psfs(spps,pyps,sxps) + #return + + # END TEST PSF I/O + #################################### + + rval = pyft.fit_psf(opts,pyio,pypr,pymg,pyps) # fit psf + # check for fit_psf success if rval != 0: comstr = " ".join(com) @@ -330,8 +368,19 @@ def main(args, comm=None): "value {} running {}".format(rank, rval, comstr)) failcount += 1 - rval = pyio.write_psf_data(opts,pyps) # write psf + rval = pyio.write_psf_data(opts,pyps) # write psf rval = pyio.write_spots(opts,pyps) # write spots + + xtrace=spx.get_trace(pyps,'x') + ytrace=spx.get_trace(pyps,'y') + + fiberkeys = spx.VectorInt() + wavekeys = spx.VectorDouble() + spx.get_tracekeys(pyps,fiberkeys,wavekeys) + + print('first x entry',xtrace[0]) + print('first y entry',ytrace[0]) + print(fiberkeys,wavekeys) if comm is not None: from mpi4py import MPI From d6ade97ecb77b663934be1a05a673ebd345f8dc2 Mon Sep 17 00:00:00 2001 From: Marcelo Alvarez Date: Mon, 9 Nov 2020 15:54:17 -0800 Subject: [PATCH 07/16] psf table I/O refactor --- py/desispec/scripts/specex.py | 35 ++++------------------------------- 1 file changed, 4 insertions(+), 31 deletions(-) diff --git a/py/desispec/scripts/specex.py b/py/desispec/scripts/specex.py index ff6ab1361..911451aef 100644 --- a/py/desispec/scripts/specex.py +++ b/py/desispec/scripts/specex.py @@ -316,9 +316,9 @@ def main(args, comm=None): pyps = spx.PyPSF() # psf data pyft = spx.PyFitting() # psf fitting - rval = opts.parse(com) # com is list of args - rval = pyio.check_input_psf(opts) # set input psf bools - rval = pypr.deal_with_priors(opts) # set Gaussian priors + opts.parse(com) # com is list of args + pyio.check_input_psf(opts) # set input psf bools + pypr.deal_with_priors(opts) # set Gaussian priors fiberkeys = spx.VectorInt() wavekeys = spx.VectorDouble() @@ -330,34 +330,7 @@ def main(args, comm=None): # read preproc images into pymg with desispec (astropy.io.fits) pymg = read_desi_ppimage_spx(opts) - rval = pyio.read_psf_data(opts,pyps) # read psf (specex-py) - - #################################### - # BEGIN TEST PSF I/O - - #spps = specter.psf.GaussHermitePSF( - # opts.input_psf_filename) # read psf (specter) - - #rval = pyio.write_psf_data(opts,pyps) # dummy write psf - #print(spx.get_pyps_image(pymg,pyps)) - #print(len(spx.get_pyps_image(pymg,pyps))) - - #return - - #print(spx.get_pyps_image(pyps)) - #return - # sxps = spx.GaussHermitePSF() # new PSF object (specex) - - #print(opts.output_fits_filename) - #opts.output_fits_filename=opts.output_fits_filename[:-5]+'_cp.fits' - #print(opts.output_fits_filename) - #print('writing psf') - #rval = pyio.write_psf_data(opts,pyps) # write psf - #compare_psfs(spps,pyps,sxps) - #return - - # END TEST PSF I/O - #################################### + pyio.read_psf_data(opts,pyps) # read psf (specex-py) rval = pyft.fit_psf(opts,pyio,pypr,pymg,pyps) # fit psf From df6240d8a8815c6cd15e08d04da7e47abe23b9a9 Mon Sep 17 00:00:00 2001 From: Marcelo Alvarez Date: Mon, 16 Nov 2020 14:11:25 -0800 Subject: [PATCH 08/16] psfio.py added --- py/desispec/scripts/psfio.py | 98 +++++++++++++++++++++++++++++++++++ py/desispec/scripts/specex.py | 21 ++++---- 2 files changed, 107 insertions(+), 12 deletions(-) create mode 100644 py/desispec/scripts/psfio.py diff --git a/py/desispec/scripts/psfio.py b/py/desispec/scripts/psfio.py new file mode 100644 index 000000000..326a920a5 --- /dev/null +++ b/py/desispec/scripts/psfio.py @@ -0,0 +1,98 @@ +from astropy.io import fits + +import numpy as np +import specex as spx +import specter.psf + +def write_psf(pyps,opts): + + dotest = False + if not dotest: + + specter_psf = specter.psf.GaussHermitePSF(opts.output_fits_filename) + + xtrace=spx.get_trace(pyps,'x') + ytrace=spx.get_trace(pyps,'y') + + fiberkeys = spx.VectorInt() + wavekeys = spx.VectorDouble() + spx.get_tracekeys(pyps,fiberkeys,wavekeys) + + print(fiberkeys,wavekeys) + print(np.shape(xtrace)) + xtrace = np.reshape(xtrace,(500,7)) + ytrace = np.reshape(ytrace,(500,7)) + + else: + + spxpsf='/project/projectdirs/desi/users/malvarez/desi/psf-dev/fit-psf-r0-00055705_00.fits' + specter_psf = specter.psf.GaussHermitePSF(spxpsf) + + #x1 = np.linspace(0,1.,7) + #x2 = np.linspace(0,1.,500) + #xtrace, ytrace = np.meshgrid(x1,x2) + + fiberkeys = np.asarray([0,499]) + wavekeys = np.asarray([5543.,7826.]) + + hdul = fits.open(spxpsf) + xtrace = hdul[0].data + ytrace = hdul[1].data + print(np.shape(xtrace),np.shape(ytrace)) + + print('first x entry',xtrace[0][0]) + print('first y entry',ytrace[0][0]) + + fitsfile='foo.fits' + + # XTRACE AND HEADER + + hdu_xtrace = fits.PrimaryHDU(xtrace) + hdr_xtrace = hdu_xtrace.header + + hdr_xtrace['EXTNAME'] = 'XTRACE' + hdr_xtrace['FIBERMIN'] = fiberkeys[0] + hdr_xtrace['FIBERMAX'] = fiberkeys[1] + hdr_xtrace['WAVEMIN'] = wavekeys[0] + hdr_xtrace['WAVEMAX'] = wavekeys[1] + + hdr_xtrace.comments['EXTEND'] = 'FITS dataset may contain extensions' + hdr_xtrace.comments['SIMPLE'] = 'file does conform to FITS standard' + hdr_xtrace.comments['BITPIX'] = 'number of bits per data pixel' + hdr_xtrace.comments['NAXIS'] = 'number of data axes' + hdr_xtrace.comments['NAXIS1'] = 'length of data axis 1' + hdr_xtrace.comments['NAXIS2'] = 'length of data axis 2' + + hdr_xtrace.insert(6,('COMMENT', ' FITS (Flexible Image Transport System) format is defined in \'Astronomy and Astrophysics\', volume 376, page 359; bibcode: 2001A&A...376..359H')) + + hdr_xtrace['PSFTYPE'] = 'GAUSS-HERMITE' + hdr_xtrace['PSFVER'] = 3 + + hdr_xtrace.insert(16,('COMMENT', 'PSF generated by specex, https://github.com/desihub/specex')) + # YTRACE AND HEADER + + hdu_ytrace = fits.ImageHDU(ytrace) + hdr_ytrace = hdu_ytrace.header + + hdr_ytrace['PCOUNT'] = 0 + hdr_ytrace['GCOUNT'] = 1 + hdr_ytrace['EXTNAME'] = 'YTRACE' + hdr_ytrace['FIBERMIN'] = fiberkeys[0] + hdr_ytrace['FIBERMAX'] = fiberkeys[1] + hdr_ytrace['WAVEMIN'] = wavekeys[0] + hdr_ytrace['WAVEMAX'] = wavekeys[1] + + hdr_ytrace.comments['XTENSION'] = 'IMAGE extension' + hdr_ytrace.comments['PCOUNT'] = 'required keyword; must = 0' + hdr_ytrace.comments['GCOUNT'] = 'required keyword; must = 1' + hdr_ytrace.comments['BITPIX'] = 'number of bits per data pixel' + hdr_ytrace.comments['NAXIS'] = 'number of data axes' + hdr_ytrace.comments['NAXIS1'] = 'length of data axis 1' + hdr_ytrace.comments['NAXIS2'] = 'length of data axis 2' + + # WRITE HDUS + hdul = fits.HDUList([hdu_xtrace, hdu_ytrace]) + hdul.writeto(fitsfile,overwrite=True) + + return + diff --git a/py/desispec/scripts/specex.py b/py/desispec/scripts/specex.py index 911451aef..f9a8b3b19 100644 --- a/py/desispec/scripts/specex.py +++ b/py/desispec/scripts/specex.py @@ -20,6 +20,7 @@ import specex as spx import specter.psf +from . import psfio modext = "so" if sys.platform == "darwin": @@ -189,6 +190,10 @@ def parse(options=None): def main(args, comm=None): + #pyps = opts= 0 + #psfio.write_psf(pyps,opts) + #return + log = get_logger() imgfile = args.input_image @@ -315,7 +320,7 @@ def main(args, comm=None): pypr = spx.PyPrior() # Gaussian priors pyps = spx.PyPSF() # psf data pyft = spx.PyFitting() # psf fitting - + opts.parse(com) # com is list of args pyio.check_input_psf(opts) # set input psf bools pypr.deal_with_priors(opts) # set Gaussian priors @@ -344,17 +349,9 @@ def main(args, comm=None): rval = pyio.write_psf_data(opts,pyps) # write psf rval = pyio.write_spots(opts,pyps) # write spots - xtrace=spx.get_trace(pyps,'x') - ytrace=spx.get_trace(pyps,'y') - - fiberkeys = spx.VectorInt() - wavekeys = spx.VectorDouble() - spx.get_tracekeys(pyps,fiberkeys,wavekeys) - - print('first x entry',xtrace[0]) - print('first y entry',ytrace[0]) - print(fiberkeys,wavekeys) - + psfio.write_psf(pyps,opts) + return + if comm is not None: from mpi4py import MPI failcount = comm.allreduce(failcount, op=MPI.SUM) From a68a6cce2fd938cbcd5e79ff09e45aeba427cc3d Mon Sep 17 00:00:00 2001 From: Marcelo Alvarez Date: Mon, 23 Nov 2020 15:29:48 -0800 Subject: [PATCH 09/16] fitsio and pybind for table entries --- py/desispec/scripts/psfio.py | 218 +++++++++++++++++++++++++--------- py/desispec/scripts/specex.py | 11 +- 2 files changed, 167 insertions(+), 62 deletions(-) diff --git a/py/desispec/scripts/psfio.py b/py/desispec/scripts/psfio.py index 326a920a5..3068a3614 100644 --- a/py/desispec/scripts/psfio.py +++ b/py/desispec/scripts/psfio.py @@ -1,37 +1,99 @@ from astropy.io import fits +import fitsio +from fitsio import FITS,FITSHDR import numpy as np -import specex as spx import specter.psf -def write_psf(pyps,opts): +def write_psf(pyps,opts,dotest): - dotest = False if not dotest: + import specex as spx specter_psf = specter.psf.GaussHermitePSF(opts.output_fits_filename) + # SOME VARIABLES HARD CODED FOR DEVELOPMENT FOR LATER REPLACEMENT + # BY QUERIES TO specex.get_* METHODS + nrows=59 + ntab1=4 + ntab2=500 + + # GET TRACES xtrace=spx.get_trace(pyps,'x') ytrace=spx.get_trace(pyps,'y') + xtrace = np.reshape(xtrace,(500,7)) + ytrace = np.reshape(ytrace,(500,7)) + + # GET TRACE KEYS fiberkeys = spx.VectorInt() wavekeys = spx.VectorDouble() spx.get_tracekeys(pyps,fiberkeys,wavekeys) - print(fiberkeys,wavekeys) - print(np.shape(xtrace)) - xtrace = np.reshape(xtrace,(500,7)) - ytrace = np.reshape(ytrace,(500,7)) + mjd = spx.VectorLongLong() + plate_id = spx.VectorLongLong() + camera_id = spx.VectorString() + arc_exposure_id = spx.VectorLongLong() + NPIX_X = spx.VectorLongLong() + NPIX_Y = spx.VectorLongLong() + hSizeX = spx.VectorLongLong() + hSizeY = spx.VectorLongLong() + nparams_all = spx.VectorLongLong() + ncoeff = spx.VectorLongLong() + GHDEGX = spx.VectorLongLong() + GHDEGY = spx.VectorLongLong() + psf_error = spx.VectorDouble() + readout_noise = spx.VectorDouble() + gain = spx.VectorDouble() + + spx.get_tablekeys(pyps,mjd,plate_id,camera_id,arc_exposure_id, + NPIX_X,NPIX_Y,hSizeX,hSizeY, + nparams_all,ncoeff,GHDEGX,GHDEGY, + psf_error,readout_noise,gain) + + # GET TABLE + table_col0 = spx.VectorString() + table_col1 = spx.VectorDouble() + table_col2 = spx.VectorInt() + table_col3 = spx.VectorInt() + + spx.get_table(pyps,table_col0,table_col1,table_col2,table_col3) + + + # CONVERT COLUMNS 1-3 TO NUMPY ARRAYS OF THE CORRECT DIMENSIONS + col0 = table_col0 + col1 = np.zeros((nrows,ntab2,ntab1)) + col2 = np.zeros( nrows) + col3 = np.zeros( nrows) + + print(table_col0) + i = 0 + for r in np.arange(nrows): + for t2 in np.arange(ntab2): + for t1 in np.arange(ntab1): + col1[r,t2,t1] = table_col1[i] + i += 1 + col2[r] = table_col2[r] + col3[r] = table_col3[r] + + # LOAD TABLE INTO data ARRAY FOR WRITING + data = np.zeros(nrows, dtype=[('PARAM', 'U8'), + ('COEFF', 'f8', (ntab2,ntab1)), + ('LEGDEGX', 'i4'), + ('LEGDEGW', 'i4')]) + data['PARAM'] = col0 + data['COEFF'] = col1 + data['LEGDEGX'] = col2 + data['LEGDEGW'] = col3 + else: - spxpsf='/project/projectdirs/desi/users/malvarez/desi/psf-dev/fit-psf-r0-00055705_00.fits' + spxpsf=( + '/project/projectdirs/desi/users/malvarez/desi/'+ + 'psf-dev/fit-psf-r0-00055705_00.fits') specter_psf = specter.psf.GaussHermitePSF(spxpsf) - #x1 = np.linspace(0,1.,7) - #x2 = np.linspace(0,1.,500) - #xtrace, ytrace = np.meshgrid(x1,x2) - fiberkeys = np.asarray([0,499]) wavekeys = np.asarray([5543.,7826.]) @@ -40,59 +102,99 @@ def write_psf(pyps,opts): ytrace = hdul[1].data print(np.shape(xtrace),np.shape(ytrace)) - print('first x entry',xtrace[0][0]) - print('first y entry',ytrace[0][0]) + print('first x entry ',xtrace[0][0]) + print('first y entry ',ytrace[0][0]) + print('fiber,wave keys',fiberkeys,wavekeys) - fitsfile='foo.fits' + fitsfile = FITS('fio.fits','rw',clobber=True) - # XTRACE AND HEADER - - hdu_xtrace = fits.PrimaryHDU(xtrace) - hdr_xtrace = hdu_xtrace.header - - hdr_xtrace['EXTNAME'] = 'XTRACE' - hdr_xtrace['FIBERMIN'] = fiberkeys[0] - hdr_xtrace['FIBERMAX'] = fiberkeys[1] - hdr_xtrace['WAVEMIN'] = wavekeys[0] - hdr_xtrace['WAVEMAX'] = wavekeys[1] - - hdr_xtrace.comments['EXTEND'] = 'FITS dataset may contain extensions' - hdr_xtrace.comments['SIMPLE'] = 'file does conform to FITS standard' - hdr_xtrace.comments['BITPIX'] = 'number of bits per data pixel' - hdr_xtrace.comments['NAXIS'] = 'number of data axes' - hdr_xtrace.comments['NAXIS1'] = 'length of data axis 1' - hdr_xtrace.comments['NAXIS2'] = 'length of data axis 2' + # WRITE XTRACE + fitsfile.write(xtrace) - hdr_xtrace.insert(6,('COMMENT', ' FITS (Flexible Image Transport System) format is defined in \'Astronomy and Astrophysics\', volume 376, page 359; bibcode: 2001A&A...376..359H')) + fitsfile[0].write_key('EXTNAME','XTRACE','') + fitsfile[0].write_key('FIBERMIN',fiberkeys[0]) + fitsfile[0].write_key('FIBERMAX',fiberkeys[1]) + fitsfile[0].write_key('WAVEMIN',wavekeys[0]) + fitsfile[0].write_key('WAVEMAX',wavekeys[1]) + fitsfile[0].write_key('PSFTYPE','GAUSS-HERMITE') + fitsfile[0].write_key('PSFVER', 3) + fitsfile[0].write_comment('PSF generated by specex, https://github.com/desihub/specex') + + # WRITE YTRACE + fitsfile.write(ytrace) + + fitsfile[1].write_key('PCOUNT',0,'required keyword; must = 0') + fitsfile[1].write_key('GCOUNT',1,'required keyword; must = 1') + fitsfile[1].write_key('EXTNAME','YTRACE','') + fitsfile[1].write_key('FIBERMIN',fiberkeys[0],'') + fitsfile[1].write_key('FIBERMAX',fiberkeys[1],'') + fitsfile[1].write_key('WAVEMIN',wavekeys[0],'') + fitsfile[1].write_key('WAVEMAX',wavekeys[1],'') + + # WRITE TABLE + fitsfile.write(data) - hdr_xtrace['PSFTYPE'] = 'GAUSS-HERMITE' - hdr_xtrace['PSFVER'] = 3 + # WRITE REFORMATTED HEADER KEYS + fitsfile[2].write_key('TDIM3 ','( 1)','dimension') + fitsfile[2].write_key('TDIM4 ','( 1)','dimension') - hdr_xtrace.insert(16,('COMMENT', 'PSF generated by specex, https://github.com/desihub/specex')) - # YTRACE AND HEADER + # WRITE COMMENTS + fitsfile[2].write_comment('------------------------------------------------------------------------') + fitsfile[2].write_comment('PSF generated by specex, https://github.com/julienguy/spece') + fitsfile[2].write_comment('PSF fit date 2020-11-12') + fitsfile[2].write_comment('-') + fitsfile[2].write_comment('Each row of the table contains the data vector of one PSF parameter') + fitsfile[2].write_comment('The size of the vector is ((FIBERMAX-FIBERMIN+1)*(LEGDEG+1))') + fitsfile[2].write_comment('Description of the NPARAMS parameters :') + fitsfile[2].write_comment('GHSIGX : Sigma of first Gaussian along CCD columns for PSF core') + fitsfile[2].write_comment('GHSIGY : Sigma of first Gaussian along CCD rows for PSF core') + fitsfile[2].write_comment('GH-i-j : Hermite pol. coefficents, i along columns, j along rows,') + fitsfile[2].write_comment(' i is integer from 0 to GHDEGX, j is integer from 0 to GHDEGY,') + fitsfile[2].write_comment(' there are (GHDEGX+1)*(GHDEGY+1) such coefficents.') + fitsfile[2].write_comment('TAILAMP : Amplitude of PSF tail') + fitsfile[2].write_comment('TAILCORE : Size in pixels of PSF tail saturation in PSF core') + fitsfile[2].write_comment('TAILXSCA : Scaling apply to CCD coordinate along columns for PSF tail') + fitsfile[2].write_comment('TAILYSCA : Scaling apply to CCD coordinate along rows for PSF tail') + fitsfile[2].write_comment('TAILINDE : Asymptotic power law index of PSF tail') + fitsfile[2].write_comment('CONT : Continuum flux in arc image (not part of PSF)') + fitsfile[2].write_comment('- ') + fitsfile[2].write_comment('PSF_core(X,Y) = [ SUM_ij (GH-i-j)*HERM(i,X/GHSIGX)*HERM(j,Y/GHSIGX) ]') + fitsfile[2].write_comment(' *GAUS(X,GHSIGX)*GAUS(Y,GHSIGY)') + fitsfile[2].write_comment('- ') + fitsfile[2].write_comment('PSF_tail(X,Y) = TAILAMP*R^2/(TAILCORE^2+R^2)^(1+TAILINDE/2)') + fitsfile[2].write_comment(' with R^2=(X/TAILXSCA)^2+(Y/TAILYSCA)^2') + fitsfile[2].write_comment('- ') + fitsfile[2].write_comment('PSF_core is integrated in pixel') + fitsfile[2].write_comment('PSF_tail is not, it is evaluated at center of pixel') + fitsfile[2].write_comment('------------------------------------------------------------------------') - hdu_ytrace = fits.ImageHDU(ytrace) - hdr_ytrace = hdu_ytrace.header + # WRITE KEYS + fitsfile[2].write_key('EXTNAME','PSF',''); + fitsfile[2].write_key('PSFTYPE','GAUSS-HERMITE',''); + fitsfile[2].write_key('PSFVER','3',''); - hdr_ytrace['PCOUNT'] = 0 - hdr_ytrace['GCOUNT'] = 1 - hdr_ytrace['EXTNAME'] = 'YTRACE' - hdr_ytrace['FIBERMIN'] = fiberkeys[0] - hdr_ytrace['FIBERMAX'] = fiberkeys[1] - hdr_ytrace['WAVEMIN'] = wavekeys[0] - hdr_ytrace['WAVEMAX'] = wavekeys[1] - - hdr_ytrace.comments['XTENSION'] = 'IMAGE extension' - hdr_ytrace.comments['PCOUNT'] = 'required keyword; must = 0' - hdr_ytrace.comments['GCOUNT'] = 'required keyword; must = 1' - hdr_ytrace.comments['BITPIX'] = 'number of bits per data pixel' - hdr_ytrace.comments['NAXIS'] = 'number of data axes' - hdr_ytrace.comments['NAXIS1'] = 'length of data axis 1' - hdr_ytrace.comments['NAXIS2'] = 'length of data axis 2' + fitsfile[2].write_key('MJD',mjd[0],'MJD of arc lamp exposure'); + fitsfile[2].write_key('PLATEID',plate_id[0],'plate ID of arc lamp exposure'); + fitsfile[2].write_key('CAMERA',camera_id[0],'camera ID'); + fitsfile[2].write_key('ARCEXP',arc_exposure_id[0],'ID of arc lamp exposure used to fit PSF'); - # WRITE HDUS - hdul = fits.HDUList([hdu_xtrace, hdu_ytrace]) - hdul.writeto(fitsfile,overwrite=True) + fitsfile[2].write_key('NPIX_X',NPIX_X[0],'number of columns in input CCD image'); + fitsfile[2].write_key('NPIX_Y',NPIX_Y[0],'number of rows in input CCD image'); + fitsfile[2].write_key('HSIZEX',hSizeX[0],'Half size of PSF in fit, NX=2*HSIZEX+1'); + fitsfile[2].write_key('HSIZEY',hSizeY[0],'Half size of PSF in fit, NY=2*HSIZEY+1'); + fitsfile[2].write_key('FIBERMIN',fiberkeys[0],'first fiber (starting at 0)'); + fitsfile[2].write_key('FIBERMAX',fiberkeys[1],'last fiber (included)'); + fitsfile[2].write_key('NPARAMS',nparams_all[0],'number of PSF parameters'); + fitsfile[2].write_key('LEGDEG',(ncoeff[0]-1),'degree of Legendre pol.(wave) for parameters'); + fitsfile[2].write_key('GHDEGX',GHDEGX[0],'degree of Hermite polynomial along CCD columns'); + fitsfile[2].write_key('GHDEGY',GHDEGY[0],'degree of Hermite polynomial along CCD rows'); + fitsfile[2].write_key('WAVEMIN',wavekeys[0],'minimum wavelength (A), used for the Legendre polynomials'); + fitsfile[2].write_key('WAVEMAX',wavekeys[1],'maximum wavelength (A), used for the Legendre polynomials'); + fitsfile[2].write_key('PSFERROR',psf_error[0],'assumed PSF fractional error in chi2'); + fitsfile[2].write_key('READNOIS',readout_noise[0],'assumed read out noise in chi2'); + fitsfile[2].write_key('GAIN',gain[0],'assumed gain in chi2'); - return + # CLOSE FITS FILE + fitsfile.close() + return diff --git a/py/desispec/scripts/specex.py b/py/desispec/scripts/specex.py index f9a8b3b19..302311b66 100644 --- a/py/desispec/scripts/specex.py +++ b/py/desispec/scripts/specex.py @@ -320,8 +320,11 @@ def main(args, comm=None): pypr = spx.PyPrior() # Gaussian priors pyps = spx.PyPSF() # psf data pyft = spx.PyFitting() # psf fitting - - opts.parse(com) # com is list of args + + coms = spx.VectorString() + for strs in com: + coms.append(strs) + opts.parse(coms) # coms is list of args pyio.check_input_psf(opts) # set input psf bools pypr.deal_with_priors(opts) # set Gaussian priors @@ -346,10 +349,10 @@ def main(args, comm=None): "value {} running {}".format(rank, rval, comstr)) failcount += 1 - rval = pyio.write_psf_data(opts,pyps) # write psf + rval = pyio.write_psf_data(opts,pyps) # write psf (specex) + psfio.write_psf(pyps,opts,False) # write psf (fitsio) rval = pyio.write_spots(opts,pyps) # write spots - psfio.write_psf(pyps,opts) return if comm is not None: From 17c97cba8192733d8e14aa091c12eddbb1633729 Mon Sep 17 00:00:00 2001 From: Marcelo Alvarez Date: Mon, 23 Nov 2020 15:37:15 -0800 Subject: [PATCH 10/16] cleaning up --- py/desispec/scripts/psfio.py | 180 +++++++++++++++------------------- py/desispec/scripts/specex.py | 6 +- 2 files changed, 83 insertions(+), 103 deletions(-) diff --git a/py/desispec/scripts/psfio.py b/py/desispec/scripts/psfio.py index 3068a3614..05a8e62c6 100644 --- a/py/desispec/scripts/psfio.py +++ b/py/desispec/scripts/psfio.py @@ -5,107 +5,87 @@ import numpy as np import specter.psf -def write_psf(pyps,opts,dotest): - - if not dotest: - - import specex as spx - specter_psf = specter.psf.GaussHermitePSF(opts.output_fits_filename) - - # SOME VARIABLES HARD CODED FOR DEVELOPMENT FOR LATER REPLACEMENT - # BY QUERIES TO specex.get_* METHODS - nrows=59 - ntab1=4 - ntab2=500 - - # GET TRACES - xtrace=spx.get_trace(pyps,'x') - ytrace=spx.get_trace(pyps,'y') +def write_psf(pyps,opts): - xtrace = np.reshape(xtrace,(500,7)) - ytrace = np.reshape(ytrace,(500,7)) - - # GET TRACE KEYS - fiberkeys = spx.VectorInt() - wavekeys = spx.VectorDouble() - spx.get_tracekeys(pyps,fiberkeys,wavekeys) - - mjd = spx.VectorLongLong() - plate_id = spx.VectorLongLong() - camera_id = spx.VectorString() - arc_exposure_id = spx.VectorLongLong() - NPIX_X = spx.VectorLongLong() - NPIX_Y = spx.VectorLongLong() - hSizeX = spx.VectorLongLong() - hSizeY = spx.VectorLongLong() - nparams_all = spx.VectorLongLong() - ncoeff = spx.VectorLongLong() - GHDEGX = spx.VectorLongLong() - GHDEGY = spx.VectorLongLong() - psf_error = spx.VectorDouble() - readout_noise = spx.VectorDouble() - gain = spx.VectorDouble() - - spx.get_tablekeys(pyps,mjd,plate_id,camera_id,arc_exposure_id, - NPIX_X,NPIX_Y,hSizeX,hSizeY, - nparams_all,ncoeff,GHDEGX,GHDEGY, - psf_error,readout_noise,gain) - - # GET TABLE - table_col0 = spx.VectorString() - table_col1 = spx.VectorDouble() - table_col2 = spx.VectorInt() - table_col3 = spx.VectorInt() - - spx.get_table(pyps,table_col0,table_col1,table_col2,table_col3) - - - # CONVERT COLUMNS 1-3 TO NUMPY ARRAYS OF THE CORRECT DIMENSIONS - col0 = table_col0 - col1 = np.zeros((nrows,ntab2,ntab1)) - col2 = np.zeros( nrows) - col3 = np.zeros( nrows) - - print(table_col0) - i = 0 - for r in np.arange(nrows): - for t2 in np.arange(ntab2): - for t1 in np.arange(ntab1): - col1[r,t2,t1] = table_col1[i] - i += 1 - col2[r] = table_col2[r] - col3[r] = table_col3[r] - - # LOAD TABLE INTO data ARRAY FOR WRITING - data = np.zeros(nrows, dtype=[('PARAM', 'U8'), - ('COEFF', 'f8', (ntab2,ntab1)), - ('LEGDEGX', 'i4'), - ('LEGDEGW', 'i4')]) - - data['PARAM'] = col0 - data['COEFF'] = col1 - data['LEGDEGX'] = col2 - data['LEGDEGW'] = col3 - - else: - - spxpsf=( - '/project/projectdirs/desi/users/malvarez/desi/'+ - 'psf-dev/fit-psf-r0-00055705_00.fits') - specter_psf = specter.psf.GaussHermitePSF(spxpsf) - - fiberkeys = np.asarray([0,499]) - wavekeys = np.asarray([5543.,7826.]) - - hdul = fits.open(spxpsf) - xtrace = hdul[0].data - ytrace = hdul[1].data - print(np.shape(xtrace),np.shape(ytrace)) - - print('first x entry ',xtrace[0][0]) - print('first y entry ',ytrace[0][0]) - print('fiber,wave keys',fiberkeys,wavekeys) - + import specex as spx + specter_psf = specter.psf.GaussHermitePSF(opts.output_fits_filename) + + # SOME VARIABLES HARD CODED FOR DEVELOPMENT FOR LATER REPLACEMENT + # BY QUERIES TO specex.get_* METHODS + nrows=59 + ntab1=4 + ntab2=500 + + # GET TRACES + xtrace=spx.get_trace(pyps,'x') + ytrace=spx.get_trace(pyps,'y') + + xtrace = np.reshape(xtrace,(500,7)) + ytrace = np.reshape(ytrace,(500,7)) + + # GET TRACE KEYS + fiberkeys = spx.VectorInt() + wavekeys = spx.VectorDouble() + spx.get_tracekeys(pyps,fiberkeys,wavekeys) + + mjd = spx.VectorLongLong() + plate_id = spx.VectorLongLong() + camera_id = spx.VectorString() + arc_exposure_id = spx.VectorLongLong() + NPIX_X = spx.VectorLongLong() + NPIX_Y = spx.VectorLongLong() + hSizeX = spx.VectorLongLong() + hSizeY = spx.VectorLongLong() + nparams_all = spx.VectorLongLong() + ncoeff = spx.VectorLongLong() + GHDEGX = spx.VectorLongLong() + GHDEGY = spx.VectorLongLong() + psf_error = spx.VectorDouble() + readout_noise = spx.VectorDouble() + gain = spx.VectorDouble() + + spx.get_tablekeys(pyps,mjd,plate_id,camera_id,arc_exposure_id, + NPIX_X,NPIX_Y,hSizeX,hSizeY, + nparams_all,ncoeff,GHDEGX,GHDEGY, + psf_error,readout_noise,gain) + + # GET TABLE + table_col0 = spx.VectorString() + table_col1 = spx.VectorDouble() + table_col2 = spx.VectorInt() + table_col3 = spx.VectorInt() + + spx.get_table(pyps,table_col0,table_col1,table_col2,table_col3) + + + # CONVERT COLUMNS 1-3 TO NUMPY ARRAYS OF THE CORRECT DIMENSIONS + col0 = table_col0 + col1 = np.zeros((nrows,ntab2,ntab1)) + col2 = np.zeros( nrows) + col3 = np.zeros( nrows) + + print(table_col0) + i = 0 + for r in np.arange(nrows): + for t2 in np.arange(ntab2): + for t1 in np.arange(ntab1): + col1[r,t2,t1] = table_col1[i] + i += 1 + col2[r] = table_col2[r] + col3[r] = table_col3[r] + + # LOAD TABLE INTO data ARRAY FOR WRITING + data = np.zeros(nrows, dtype=[('PARAM', 'U8'), + ('COEFF', 'f8', (ntab2,ntab1)), + ('LEGDEGX', 'i4'), + ('LEGDEGW', 'i4')]) + + data['PARAM'] = col0 + data['COEFF'] = col1 + data['LEGDEGX'] = col2 + data['LEGDEGW'] = col3 + + # OPEN FITSFILE fitsfile = FITS('fio.fits','rw',clobber=True) # WRITE XTRACE diff --git a/py/desispec/scripts/specex.py b/py/desispec/scripts/specex.py index 302311b66..5945f6c2e 100644 --- a/py/desispec/scripts/specex.py +++ b/py/desispec/scripts/specex.py @@ -349,9 +349,9 @@ def main(args, comm=None): "value {} running {}".format(rank, rval, comstr)) failcount += 1 - rval = pyio.write_psf_data(opts,pyps) # write psf (specex) - psfio.write_psf(pyps,opts,False) # write psf (fitsio) - rval = pyio.write_spots(opts,pyps) # write spots + pyio.write_psf_data(opts,pyps) # write psf (specex) + psfio.write_psf(pyps,opts) # write psf (fitsio) + pyio.write_spots(opts,pyps) # write spots return From 05336dc84a3a6f613dceb629a8ddb873ed7ce081 Mon Sep 17 00:00:00 2001 From: Marcelo Alvarez Date: Mon, 23 Nov 2020 15:42:59 -0800 Subject: [PATCH 11/16] more cleanup --- py/desispec/scripts/specex.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/py/desispec/scripts/specex.py b/py/desispec/scripts/specex.py index 5945f6c2e..be9f1f59a 100644 --- a/py/desispec/scripts/specex.py +++ b/py/desispec/scripts/specex.py @@ -333,7 +333,7 @@ def main(args, comm=None): if specex_image_io: # read preproc images into pymg with specex (harp) - rval = pyio.read_img_data(opts,pymg) + pyio.read_img_data(opts,pymg) else: # read preproc images into pymg with desispec (astropy.io.fits) pymg = read_desi_ppimage_spx(opts) From 7962e1e19dc98266da893a2d928fc463445e3059 Mon Sep 17 00:00:00 2001 From: Marcelo Alvarez Date: Mon, 30 Nov 2020 13:04:25 -0800 Subject: [PATCH 12/16] 1st version pythono psf fitsio output --- py/desispec/scripts/psfio.py | 218 +++++++++++++++++------------------ 1 file changed, 109 insertions(+), 109 deletions(-) diff --git a/py/desispec/scripts/psfio.py b/py/desispec/scripts/psfio.py index 05a8e62c6..b0f29fdc2 100644 --- a/py/desispec/scripts/psfio.py +++ b/py/desispec/scripts/psfio.py @@ -3,125 +3,103 @@ from fitsio import FITS,FITSHDR import numpy as np -import specter.psf +import specex as spx -def write_psf(pyps,opts): +def write_psf(pyps,opts): - import specex as spx - specter_psf = specter.psf.GaussHermitePSF(opts.output_fits_filename) - - # SOME VARIABLES HARD CODED FOR DEVELOPMENT FOR LATER REPLACEMENT - # BY QUERIES TO specex.get_* METHODS - nrows=59 - ntab1=4 - ntab2=500 - - # GET TRACES + # initialize table writing + spx.tablewrite_init(pyps) + + # get traces xtrace=spx.get_trace(pyps,'x') ytrace=spx.get_trace(pyps,'y') + + xtrace = np.reshape(xtrace,(pyps.nfibers,pyps.trace_ncoeff)) + ytrace = np.reshape(ytrace,(pyps.nfibers,pyps.trace_ncoeff)) - xtrace = np.reshape(xtrace,(500,7)) - ytrace = np.reshape(ytrace,(500,7)) - - # GET TRACE KEYS - fiberkeys = spx.VectorInt() - wavekeys = spx.VectorDouble() - spx.get_tracekeys(pyps,fiberkeys,wavekeys) - - mjd = spx.VectorLongLong() - plate_id = spx.VectorLongLong() - camera_id = spx.VectorString() - arc_exposure_id = spx.VectorLongLong() - NPIX_X = spx.VectorLongLong() - NPIX_Y = spx.VectorLongLong() - hSizeX = spx.VectorLongLong() - hSizeY = spx.VectorLongLong() - nparams_all = spx.VectorLongLong() - ncoeff = spx.VectorLongLong() - GHDEGX = spx.VectorLongLong() - GHDEGY = spx.VectorLongLong() - psf_error = spx.VectorDouble() - readout_noise = spx.VectorDouble() - gain = spx.VectorDouble() - - spx.get_tablekeys(pyps,mjd,plate_id,camera_id,arc_exposure_id, - NPIX_X,NPIX_Y,hSizeX,hSizeY, - nparams_all,ncoeff,GHDEGX,GHDEGY, - psf_error,readout_noise,gain) - - # GET TABLE + # get table table_col0 = spx.VectorString() table_col1 = spx.VectorDouble() table_col2 = spx.VectorInt() table_col3 = spx.VectorInt() - spx.get_table(pyps,table_col0,table_col1,table_col2,table_col3) + table_bundle_id = spx.VectorInt() + table_bundle_ndata = spx.VectorInt() + table_bundle_nparams = spx.VectorInt() + table_bundle_chi2pdf = spx.VectorDouble() + spx.get_table(pyps,table_col0,table_col1,table_col2,table_col3, + table_bundle_id,table_bundle_ndata,table_bundle_nparams, + table_bundle_chi2pdf) - # CONVERT COLUMNS 1-3 TO NUMPY ARRAYS OF THE CORRECT DIMENSIONS + # copy columns to numpy arrays col0 = table_col0 - col1 = np.zeros((nrows,ntab2,ntab1)) - col2 = np.zeros( nrows) - col3 = np.zeros( nrows) - - print(table_col0) + col1 = np.zeros((pyps.table_nrows,pyps.nfibers,pyps.ncoeff)) + col2 = np.zeros( pyps.table_nrows) + col3 = np.zeros( pyps.table_nrows) + i = 0 - for r in np.arange(nrows): - for t2 in np.arange(ntab2): - for t1 in np.arange(ntab1): + for r in np.arange(pyps.table_nrows): + for t2 in np.arange(pyps.nfibers): + for t1 in np.arange(pyps.ncoeff): col1[r,t2,t1] = table_col1[i] i += 1 - col2[r] = table_col2[r] - col3[r] = table_col3[r] - - # LOAD TABLE INTO data ARRAY FOR WRITING - data = np.zeros(nrows, dtype=[('PARAM', 'U8'), - ('COEFF', 'f8', (ntab2,ntab1)), - ('LEGDEGX', 'i4'), - ('LEGDEGW', 'i4')]) + col2[r] = table_col2[r] + col3[r] = table_col3[r] + + # load table into data array for writing + data = np.zeros(pyps.table_nrows, + dtype=[('PARAM', 'U8'), + ('COEFF', 'f8', (pyps.nfibers,pyps.ncoeff)), + ('LEGDEGX', 'i4'), + ('LEGDEGW', 'i4')] + ) data['PARAM'] = col0 data['COEFF'] = col1 data['LEGDEGX'] = col2 data['LEGDEGW'] = col3 - - # OPEN FITSFILE + + # pad PARAM strings left justified with spaces to 8 characters + i=0 + for param in data['PARAM']: + data['PARAM'][i] = param.ljust(8,' ') + i += 1 + + # open fitsfile fitsfile = FITS('fio.fits','rw',clobber=True) - # WRITE XTRACE + # write xtrace fitsfile.write(xtrace) - fitsfile[0].write_key('EXTNAME','XTRACE','') - fitsfile[0].write_key('FIBERMIN',fiberkeys[0]) - fitsfile[0].write_key('FIBERMAX',fiberkeys[1]) - fitsfile[0].write_key('WAVEMIN',wavekeys[0]) - fitsfile[0].write_key('WAVEMAX',wavekeys[1]) - fitsfile[0].write_key('PSFTYPE','GAUSS-HERMITE') - fitsfile[0].write_key('PSFVER', 3) + fitsfile[0].write_key('EXTNAME', 'XTRACE','') + fitsfile[0].write_key('FIBERMIN', pyps.FIBERMIN) + fitsfile[0].write_key('FIBERMAX', pyps.FIBERMAX) + fitsfile[0].write_key('WAVEMIN', pyps.trace_WAVEMIN) + fitsfile[0].write_key('WAVEMAX', pyps.trace_WAVEMAX) + fitsfile[0].write_key('PSFTYPE', 'GAUSS-HERMITE') + fitsfile[0].write_key('PSFVER', 3) + fitsfile[0].write_comment('PSF generated by specex, https://github.com/desihub/specex') - # WRITE YTRACE + # write ytrace fitsfile.write(ytrace) - fitsfile[1].write_key('PCOUNT',0,'required keyword; must = 0') - fitsfile[1].write_key('GCOUNT',1,'required keyword; must = 1') - fitsfile[1].write_key('EXTNAME','YTRACE','') - fitsfile[1].write_key('FIBERMIN',fiberkeys[0],'') - fitsfile[1].write_key('FIBERMAX',fiberkeys[1],'') - fitsfile[1].write_key('WAVEMIN',wavekeys[0],'') - fitsfile[1].write_key('WAVEMAX',wavekeys[1],'') + fitsfile[1].write_key('PCOUNT', 0,'required keyword; must = 0') + fitsfile[1].write_key('GCOUNT', 1,'required keyword; must = 1') + fitsfile[1].write_key('EXTNAME', 'YTRACE') + fitsfile[1].write_key('FIBERMIN', pyps.FIBERMIN) + fitsfile[1].write_key('FIBERMAX', pyps.FIBERMAX) + fitsfile[1].write_key('WAVEMIN', pyps.trace_WAVEMIN) + fitsfile[1].write_key('WAVEMAX', pyps.trace_WAVEMAX) - # WRITE TABLE + # write table fitsfile.write(data) - # WRITE REFORMATTED HEADER KEYS - fitsfile[2].write_key('TDIM3 ','( 1)','dimension') - fitsfile[2].write_key('TDIM4 ','( 1)','dimension') - - # WRITE COMMENTS + # write comments fitsfile[2].write_comment('------------------------------------------------------------------------') - fitsfile[2].write_comment('PSF generated by specex, https://github.com/julienguy/spece') - fitsfile[2].write_comment('PSF fit date 2020-11-12') + fitsfile[2].write_comment('PSF generated by specex, https://github.com/julienguy/specex') + fitsfile[2].write_comment('PSF fit date 2020-11-12') # HARDCODED FOR DEV fitsfile[2].write_comment('-') fitsfile[2].write_comment('Each row of the table contains the data vector of one PSF parameter') fitsfile[2].write_comment('The size of the vector is ((FIBERMAX-FIBERMIN+1)*(LEGDEG+1))') @@ -148,33 +126,55 @@ def write_psf(pyps,opts): fitsfile[2].write_comment('PSF_tail is not, it is evaluated at center of pixel') fitsfile[2].write_comment('------------------------------------------------------------------------') - # WRITE KEYS + # write keys fitsfile[2].write_key('EXTNAME','PSF',''); fitsfile[2].write_key('PSFTYPE','GAUSS-HERMITE',''); fitsfile[2].write_key('PSFVER','3',''); - fitsfile[2].write_key('MJD',mjd[0],'MJD of arc lamp exposure'); - fitsfile[2].write_key('PLATEID',plate_id[0],'plate ID of arc lamp exposure'); - fitsfile[2].write_key('CAMERA',camera_id[0],'camera ID'); - fitsfile[2].write_key('ARCEXP',arc_exposure_id[0],'ID of arc lamp exposure used to fit PSF'); - - fitsfile[2].write_key('NPIX_X',NPIX_X[0],'number of columns in input CCD image'); - fitsfile[2].write_key('NPIX_Y',NPIX_Y[0],'number of rows in input CCD image'); - fitsfile[2].write_key('HSIZEX',hSizeX[0],'Half size of PSF in fit, NX=2*HSIZEX+1'); - fitsfile[2].write_key('HSIZEY',hSizeY[0],'Half size of PSF in fit, NY=2*HSIZEY+1'); - fitsfile[2].write_key('FIBERMIN',fiberkeys[0],'first fiber (starting at 0)'); - fitsfile[2].write_key('FIBERMAX',fiberkeys[1],'last fiber (included)'); - fitsfile[2].write_key('NPARAMS',nparams_all[0],'number of PSF parameters'); - fitsfile[2].write_key('LEGDEG',(ncoeff[0]-1),'degree of Legendre pol.(wave) for parameters'); - fitsfile[2].write_key('GHDEGX',GHDEGX[0],'degree of Hermite polynomial along CCD columns'); - fitsfile[2].write_key('GHDEGY',GHDEGY[0],'degree of Hermite polynomial along CCD rows'); - fitsfile[2].write_key('WAVEMIN',wavekeys[0],'minimum wavelength (A), used for the Legendre polynomials'); - fitsfile[2].write_key('WAVEMAX',wavekeys[1],'maximum wavelength (A), used for the Legendre polynomials'); - fitsfile[2].write_key('PSFERROR',psf_error[0],'assumed PSF fractional error in chi2'); - fitsfile[2].write_key('READNOIS',readout_noise[0],'assumed read out noise in chi2'); - fitsfile[2].write_key('GAIN',gain[0],'assumed gain in chi2'); - - # CLOSE FITS FILE + fitsfile[2].write_key('MJD',pyps.mjd,'MJD of arc lamp exposure'); + fitsfile[2].write_key('PLATEID',pyps.plate_id,'plate ID of arc lamp exposure'); + fitsfile[2].write_key('CAMERA',pyps.camera_id,'camera ID'); + fitsfile[2].write_key('ARCEXP',pyps.arc_exposure_id,'ID of arc lamp exposure used to fit PSF'); + + fitsfile[2].write_key('NPIX_X',pyps.NPIX_X,'number of columns in input CCD image'); + fitsfile[2].write_key('NPIX_Y',pyps.NPIX_Y,'number of rows in input CCD image'); + fitsfile[2].write_key('HSIZEX',pyps.hSizeX,'Half size of PSF in fit, NX=2*HSIZEX+1'); + fitsfile[2].write_key('HSIZEY',pyps.hSizeY,'Half size of PSF in fit, NY=2*HSIZEY+1'); + fitsfile[2].write_key('FIBERMIN',pyps.FIBERMIN,'first fiber (starting at 0)'); + fitsfile[2].write_key('FIBERMAX',pyps.FIBERMAX,'last fiber (included)'); + fitsfile[2].write_key('NPARAMS',pyps.nparams_all,'number of PSF parameters'); + fitsfile[2].write_key('LEGDEG',(pyps.ncoeff-1),'degree of Legendre pol.(wave) for parameters'); + fitsfile[2].write_key('GHDEGX',pyps.GHDEGX,'degree of Hermite polynomial along CCD columns'); + fitsfile[2].write_key('GHDEGY',pyps.GHDEGY,'degree of Hermite polynomial along CCD rows'); + fitsfile[2].write_key('WAVEMIN',pyps.table_WAVEMIN,'minimum wavelength (A), used for the Legendre polynomials'); + fitsfile[2].write_key('WAVEMAX',pyps.table_WAVEMAX,'maximum wavelength (A), used for the Legendre polynomials'); + fitsfile[2].write_key('PSFERROR',pyps.psf_error,'assumed PSF fractional error in chi2'); + fitsfile[2].write_key('READNOIS',pyps.readout_noise,'assumed read out noise in chi2'); + fitsfile[2].write_key('GAIN',pyps.gain,'assumed gain in chi2'); + + i=0 + for bid in table_bundle_id: + ndata = table_bundle_ndata[i] + nparams = table_bundle_nparams[i] + chi2pdf = table_bundle_chi2pdf[i] + i += 1 + + keybase = 'B'+str(bid).rjust(2,'0') + print(keybase) + # chi2 + fitsfile[2].write_key( + keybase+'RCHI2',chi2pdf,'best fit chi2/ndf for fiber bundle ' + +str(bid)) + # ndata + fitsfile[2].write_key( + keybase+'NDATA',ndata,'number of pixels in fit for fiber bundle ' + +str(bid)) + # chi2 + fitsfile[2].write_key( + keybase+'NPAR ',nparams,'number of parameters in fit for fiber bundle ' + +str(bid)) + + # close fits file fitsfile.close() return From 8c4ea405efc556b4a68d26f8af97e7e75b543a80 Mon Sep 17 00:00:00 2001 From: Marcelo Alvarez Date: Mon, 30 Nov 2020 13:06:25 -0800 Subject: [PATCH 13/16] 1st version python psf fitsio output --- py/desispec/scripts/psfio.py | 1 + 1 file changed, 1 insertion(+) diff --git a/py/desispec/scripts/psfio.py b/py/desispec/scripts/psfio.py index b0f29fdc2..b13b92497 100644 --- a/py/desispec/scripts/psfio.py +++ b/py/desispec/scripts/psfio.py @@ -161,6 +161,7 @@ def write_psf(pyps,opts): keybase = 'B'+str(bid).rjust(2,'0') print(keybase) + # chi2 fitsfile[2].write_key( keybase+'RCHI2',chi2pdf,'best fit chi2/ndf for fiber bundle ' From 864b856d538efe1d8b1e9269a9f0c6ef72a2e091 Mon Sep 17 00:00:00 2001 From: Marcelo Alvarez Date: Mon, 30 Nov 2020 14:46:29 -0800 Subject: [PATCH 14/16] match harp header format --- py/desispec/scripts/psfio.py | 29 ++++++++++++++++++++++++++--- 1 file changed, 26 insertions(+), 3 deletions(-) diff --git a/py/desispec/scripts/psfio.py b/py/desispec/scripts/psfio.py index b13b92497..34774cfb8 100644 --- a/py/desispec/scripts/psfio.py +++ b/py/desispec/scripts/psfio.py @@ -1,6 +1,7 @@ from astropy.io import fits import fitsio from fitsio import FITS,FITSHDR +from datetime import datetime import numpy as np import specex as spx @@ -67,7 +68,8 @@ def write_psf(pyps,opts): i += 1 # open fitsfile - fitsfile = FITS('fio.fits','rw',clobber=True) + fitsfilename=opts.output_fits_filename + fitsfile = FITS(fitsfilename,'rw',clobber=True) # write xtrace fitsfile.write(xtrace) @@ -96,10 +98,14 @@ def write_psf(pyps,opts): # write table fitsfile.write(data) + # write extra keys to match harp behaviour + fitsfile[2].write_key('TDIM3', '( 1) ', 'dimension') + fitsfile[2].write_key('TDIM4', '( 1) ', 'dimension') + # write comments fitsfile[2].write_comment('------------------------------------------------------------------------') fitsfile[2].write_comment('PSF generated by specex, https://github.com/julienguy/specex') - fitsfile[2].write_comment('PSF fit date 2020-11-12') # HARDCODED FOR DEV + fitsfile[2].write_comment('PSF fit date '+datetime.today().strftime('%Y-%m-%d')) fitsfile[2].write_comment('-') fitsfile[2].write_comment('Each row of the table contains the data vector of one PSF parameter') fitsfile[2].write_comment('The size of the vector is ((FIBERMAX-FIBERMIN+1)*(LEGDEG+1))') @@ -160,7 +166,6 @@ def write_psf(pyps,opts): i += 1 keybase = 'B'+str(bid).rjust(2,'0') - print(keybase) # chi2 fitsfile[2].write_key( @@ -178,4 +183,22 @@ def write_psf(pyps,opts): # close fits file fitsfile.close() + # final attempts to make header identical to specex version using harp + fin = open(fitsfilename,"rb") + data = fin.read() + namein = str(pyps.ncoeff)+','+str(pyps.nfibers)+') \' ' + nameout = ' '+str(pyps.ncoeff)+', '+str(pyps.nfibers)+')\' / dimension' + data = data.replace( + bytes(namein,encoding='utf8'), + bytes(nameout,encoding='utf8') + ) + # b'4,500) \' ', + # b' 4, 500)\' / dimension' + #) + fin.close() + + fin = open(fitsfilename,"wb") + fin.write(data) + fin.close() + return From 4a70f1c189b9a4c2529bfd76c9cf62119f8a713a Mon Sep 17 00:00:00 2001 From: Marcelo Alvarez Date: Tue, 1 Dec 2020 13:33:19 -0800 Subject: [PATCH 15/16] rename write_psf_data to prepare_psf --- py/desispec/scripts/specex.py | 56 +++++++++++++---------------------- 1 file changed, 20 insertions(+), 36 deletions(-) diff --git a/py/desispec/scripts/specex.py b/py/desispec/scripts/specex.py index be9f1f59a..3040dcf9c 100644 --- a/py/desispec/scripts/specex.py +++ b/py/desispec/scripts/specex.py @@ -302,59 +302,43 @@ def main(args, comm=None): log.debug("proc {} calling {}".format(rank, " ".join(com))) - argc = len(com) - arg_buffers = [ct.create_string_buffer(com[i].encode('ascii')) \ - for i in range(argc)] - addrlist = [ ct.cast(x, ct.POINTER(ct.c_char)) for x in \ - map(ct.addressof, arg_buffers) ] - arg_pointers = (ct.POINTER(ct.c_char) * argc)(*addrlist) - # old way + + # argc = len(com) + # arg_buffers = [ct.create_string_buffer(com[i].encode('ascii')) \ + # for i in range(argc)] + # addrlist = [ ct.cast(x, ct.POINTER(ct.c_char)) for x in \ + # map(ct.addressof, arg_buffers) ] + # arg_pointers = (ct.POINTER(ct.c_char) * argc)(*addrlist) # rval = libspecex.cspecex_desi_psf_fit(argc, arg_pointers) # new way - specex_image_io = False + # instantiate specex c++ objects exposed to python opts = spx.PyOptions() # input options pyio = spx.PyIO() # IO options and methods pypr = spx.PyPrior() # Gaussian priors pyps = spx.PyPSF() # psf data pyft = spx.PyFitting() # psf fitting - coms = spx.VectorString() + # copy com to opaque pybind VectorString object args + spxargs = spx.VectorString() for strs in com: - coms.append(strs) - opts.parse(coms) # coms is list of args + spxargs.append(strs) + + opts.parse(spxargs) # parse args pyio.check_input_psf(opts) # set input psf bools pypr.deal_with_priors(opts) # set Gaussian priors - - fiberkeys = spx.VectorInt() - wavekeys = spx.VectorDouble() - if specex_image_io: - # read preproc images into pymg with specex (harp) - pyio.read_img_data(opts,pymg) - else: - # read preproc images into pymg with desispec (astropy.io.fits) - pymg = read_desi_ppimage_spx(opts) - - pyio.read_psf_data(opts,pyps) # read psf (specex-py) + pymg = read_desi_ppimage_spx(opts) # read preproc images (desispec) + pyio.read_psf_data(opts,pyps) # read psf (specex) - rval = pyft.fit_psf(opts,pyio,pypr,pymg,pyps) # fit psf + pyft.fit_psf(opts,pyio,pypr,pymg,pyps) # fit psf (specex) - # check for fit_psf success - if rval != 0: - comstr = " ".join(com) - log.error("fit_psf on process {} failed with return " - "value {} running {}".format(rank, rval, comstr)) - failcount += 1 + pyio.prepare_psf(opts,pyps) # prepare psf (specex) + psfio.write_psf(pyps,opts) # write psf (fitsio) + pyio.write_spots(opts,pyps) # write spots - pyio.write_psf_data(opts,pyps) # write psf (specex) - psfio.write_psf(pyps,opts) # write psf (fitsio) - pyio.write_spots(opts,pyps) # write spots - - return - if comm is not None: from mpi4py import MPI failcount = comm.allreduce(failcount, op=MPI.SUM) @@ -368,7 +352,7 @@ def main(args, comm=None): inputs = [ "{}_{:02d}.fits".format(outroot, x) for x in bundles ] - args.disable_merge=True + args.disable_merge=False if args.disable_merge : log.info("don't merge") else : From 19fda1e14d98161d7598a735dd89b06d1c6a5dab Mon Sep 17 00:00:00 2001 From: Marcelo Alvarez Date: Mon, 14 Dec 2020 22:29:15 -0800 Subject: [PATCH 16/16] specex.py --- py/desispec/scripts/psfio.py | 4 +--- py/desispec/scripts/specex.py | 5 ++--- 2 files changed, 3 insertions(+), 6 deletions(-) diff --git a/py/desispec/scripts/psfio.py b/py/desispec/scripts/psfio.py index 34774cfb8..3f241c326 100644 --- a/py/desispec/scripts/psfio.py +++ b/py/desispec/scripts/psfio.py @@ -192,9 +192,7 @@ def write_psf(pyps,opts): bytes(namein,encoding='utf8'), bytes(nameout,encoding='utf8') ) - # b'4,500) \' ', - # b' 4, 500)\' / dimension' - #) + fin.close() fin = open(fitsfilename,"wb") diff --git a/py/desispec/scripts/specex.py b/py/desispec/scripts/specex.py index 3040dcf9c..fd25b582f 100644 --- a/py/desispec/scripts/specex.py +++ b/py/desispec/scripts/specex.py @@ -352,16 +352,15 @@ def main(args, comm=None): inputs = [ "{}_{:02d}.fits".format(outroot, x) for x in bundles ] - args.disable_merge=False if args.disable_merge : log.info("don't merge") else : #- Empirically it appears that files written by one rank sometimes #- aren't fully buffer-flushed and closed before getting here, #- despite the MPI allreduce barrier. Pause to let I/O catch up. - log.info('HACK: taking a 20 sec pause before merging') + log.info('HACK: taking a 2 sec pause before merging') sys.stdout.flush() - time.sleep(20.) + time.sleep(2.) merge_psf(inputs,outfits)