From 9b647b3475a3021a1a6355ebbc2360003db5e8b3 Mon Sep 17 00:00:00 2001 From: "Sergey E. Koposov" Date: Fri, 15 Nov 2024 22:59:45 +0000 Subject: [PATCH] add NPIX_TOT add types to columndesc skip completely masked petals --- py/rvspecfit/desi/desi_fit.py | 67 ++++++++++++++++++++++------------- 1 file changed, 42 insertions(+), 25 deletions(-) diff --git a/py/rvspecfit/desi/desi_fit.py b/py/rvspecfit/desi/desi_fit.py index d107487..b685af5 100644 --- a/py/rvspecfit/desi/desi_fit.py +++ b/py/rvspecfit/desi/desi_fit.py @@ -21,6 +21,7 @@ import astropy.units as auni # noqa: E402 import numpy as np # noqa: E402 import scipy.stats # noqa: E402 +import scipy.sparse # noqa: E402 import rvspecfit # noqa: E402 from rvspecfit import fitter_ccf, vel_fit, spec_fit, utils, \ spec_inter # noqa: E402 @@ -308,15 +309,18 @@ def proc_onespec( outdict[name2] = fit_res['param'][name1] * unit outdict[name2 + '_ERR'] = fit_res['param_err'][name1] * unit + npixes = {} for i, curd in enumerate(specdata): if curd.name not in chisqs: chisqs[curd.name] = 0 chisqs_c[curd.name] = 0 chisqs[curd.name] += fit_res['chisq_array'][i] + chisqs[curd.name] += fit_res['npix_array'][i] chisqs_c[curd.name] += chisq_cont_array[i] outdict['CHISQ_TOT'] = sum(chisqs.values()) outdict['CHISQ_C_TOT'] = sum(chisqs_c.values()) + outdict['NPIX_TOT'] = sum(npixes.values()) for s in chisqs.keys(): outdict['CHISQ_%s' % s.replace('desi_', '').upper()] = chisqs[s] @@ -632,7 +636,7 @@ def resolution_mat_tocolumns(mat): def deconvolve_resolution_matrix(mat0, sig=0.451): width, npix = mat0.shape - pix_size = 0.8 # Angstrom, desi + pix_size = 0.8 # Angstrom DESI sig_pix = sig / pix_size xs = np.arange(width) gau_mat = np.array([ @@ -661,7 +665,6 @@ def deconvolve_resolution_matrix(mat0, sig=0.451): def construct_resolution_sparse_matrix(mat): width, npix = mat.shape w2 = width // 2 - from scipy.sparse import dia_matrix mat = mat.copy() deconvolve = True if deconvolve: @@ -678,7 +681,8 @@ def construct_resolution_sparse_matrix(mat): N2 = mat_rows[:w2 + 1 + i, j].sum() mat_rows[:, j] = mat_rows[:, j] / (N2 + (N2 == 0)) * mult mat = resolution_mat_tocolumns(mat_rows) - M = dia_matrix((mat, np.arange(w2, -w2 - 1, -1)), (npix, npix)) + M = scipy.sparse.dia_matrix((mat, np.arange(w2, -w2 - 1, -1)), + (npix, npix)) return M @@ -755,9 +759,13 @@ def get_specdata(waves, for s in setups: spec = fluxes[s][seqid] * 1 curivars = ivars[s][seqid] * 1 + badmask = (masks[s][seqid] > 0) medspec = np.nanmedian(spec) + if badmask.all(): + # skip bad data + continue if medspec == 0: - medspec = np.nanmedian(spec[spec > 0]) + medspec = np.nanmedian(spec[(spec > 0) & (~badmask)]) if not np.isfinite(medspec): medspec = np.nanmedian(np.abs(spec)) if not np.isfinite(medspec) or medspec == 0: @@ -769,7 +777,6 @@ def get_specdata(waves, dicroicmask = (waves[s] > 4300) & (waves[s] < 4450) else: dicroicmask = np.zeros(len(waves[s]), dtype=bool) - badmask = (masks[s][seqid] > 0) baderr = curivars <= 0 edge_mask = np.zeros(len(spec), dtype=bool) if use_resolution_matrix: @@ -815,7 +822,7 @@ def get_specdata(waves, def comment_filler(tab, desc): """ Fill comments in the FITS header """ for i, k in enumerate(tab.data.columns): - tab.header['TCOMM%d' % (i + 1)] = desc.get(k.name) or '' + tab.header['TCOMM%d' % (i + 1)] = desc.get(k.name)[1] or '' return tab @@ -831,31 +838,41 @@ def get_column_desc(setups): when we fitting a given set of configurations """ columnDesc = dict([ - ('VRAD', 'Radial velocity'), ('VRAD_ERR', 'Radial velocity error'), - ('VRAD_SKEW', 'Radial velocity posterior skewness'), - ('VRAD_KURT', 'Radial velocity posterior kurtosis'), - ('VSINI', 'Stellar rotation velocity'), - ('LOGG', 'Log of surface gravity'), ('TEFF', 'Effective temperature'), - ('FEH', '[Fe/H] from template fitting'), - ('ALPHAFE', '[alpha/Fe] from template fitting'), - ('LOGG_ERR', 'Log of surface gravity uncertainty'), - ('TEFF_ERR', 'Effective temperature uncertainty'), - ('FEH_ERR', '[Fe/H] uncertainty from template fitting'), - ('ALPHAFE_ERR', '[alpha/Fe] uncertainty from template fitting'), - ('CHISQ_TOT', 'Total chi-square for all arms'), - ('CHISQ_C_TOT', + ('VRAD', np.float32, 'Radial velocity'), + ('VRAD_ERR', np.float32, 'Radial velocity error'), + ('VRAD_SKEW', np.float32, 'Radial velocity posterior skewness'), + ('VRAD_KURT', np.float32, 'Radial velocity posterior kurtosis'), + ('VSINI', np.float32, 'Stellar rotation velocity'), + ('LOGG', np.float32, 'Log of surface gravity'), + ('TEFF', np.float32, 'Effective temperature'), + ('FEH', np.float32, '[Fe/H] from template fitting'), + ('ALPHAFE', np.float32, '[alpha/Fe] from template fitting'), + ('LOGG_ERR', np.float32, 'Log of surface gravity uncertainty'), + ('TEFF_ERR', np.float32, 'Effective temperature uncertainty'), + ('FEH_ERR', np.float32, '[Fe/H] uncertainty from template fitting'), + ('ALPHAFE_ERR', np.float32, + '[alpha/Fe] uncertainty from template fitting'), + ('CHISQ_TOT', np.float64, 'Total chi-square for all arms'), + ('NPIX_TOT', np.float64, 'Total number of unmasked pixels fitted'), + ('CHISQ_C_TOT', np.float64, 'Total chi-square for all arms for polynomial only fit'), - ('TARGETID', 'DESI targetid'), ('EXPID', 'DESI exposure id'), - ('SUCCESS', "Did we succeed or fail"), - ('RVS_WARN', "RVSpecFit warning flag"), ('RR_Z', 'Redrock redshift'), - ('RR_SPECTYPE', 'Redrock spectype') + ('VRAD_CCF', np.float32, 'Initial velocity from cross-correlation'), + ('TARGETID', np.int64, 'DESI targetid'), + ('EXPID', np.int64, 'DESI exposure id'), + ('SUCCESS', bool, "Did we succeed or fail"), + ('RVS_WARN', np.int64, "RVSpecFit warning flag"), + ('RR_Z', np.float64, 'Redrock redshift'), + ('RR_SPECTYPE', np.str, 'Redrock spectype') ]) for curs in setups: curs = curs.upper() - columnDesc['SN_%s' % curs] = ('Median S/N in the %s arm' % curs) - columnDesc['CHISQ_%s' % curs] = ('Chi-square in the %s arm' % curs) + columnDesc['SN_%s' % curs] = (np.float32, + 'Median S/N in the %s arm' % curs) + columnDesc['CHISQ_%s' % curs] = (np.float64, + 'Chi-square in the %s arm' % curs) columnDesc['CHISQ_C_%s' % curs] = ( + np.float64, 'Chi-square in the %s arm after fitting continuum only' % curs) return columnDesc