#!/usr/bin/env python
'''
Calculate stand alone tsnr tables for a given desispec prod.

Standard invocation for ``daily`` specprod on ``perlmutter``::

    desi_tsnr_afterburner -o ${DESI_SPECTRO_REDUX}/${SPECPROD}/exposures-${SPECPROD}.fits --prod ${SPECPROD} \
                          --tile-completeness ${DESI_SPECTRO_REDUX}/${SPECPROD}/tiles-${SPECPROD}.csv \
                          --aux /global/cfs/cdirs/desi/survey/observations/SV1/sv1-tiles.fits \
                          --gfa-proc-dir /global/cfs/cdirs/desi/survey/GFA/ \
                          --compute-skymags --add-badexp --update --nproc 32 \
                          --nights $NIGHT &>> ${LOGLOCATION}/tsnr-${SPECPROD}-${NIGHT}.log

Use MPI::

    salloc -N 2 -C haswell -q interactive -t 00:30:00

    # 2 nodes, 2 mpi ranks reducing 2 nights (1 per rank) with 32 processes.
    srun -N 2 -n 2 -c 32 python desi_tsnr_afterburner
'''

import os,sys
import json
import glob
import itertools
import argparse
import astropy.io.fits as fits
import fitsio
import numpy as np
import multiprocessing
from astropy.table import Table,vstack

import yaml

from   desispec.io import read_sky
from   desispec.io import read_fiberflat
from   pathlib import Path
from   desispec.io.meta import findfile, specprod_root, faflavor2program
from   desispec.calibfinder import CalibFinder
from   desispec.io import read_frame
from   desispec.io import read_fibermap
from   desispec.io import read_table
from   desispec.io.fluxcalibration import read_flux_calibration
from   desispec.io.util import get_tempfilename, decode_camword, difference_camwords
from   desiutil.log import get_logger
from   desispec.tsnr import calc_tsnr2,tsnr2_to_efftime
# from   astropy.table import Table, vstack
from   desiutil.depend import getdep
from   desispec.tilecompleteness import read_gfa_data, compute_tile_completeness_table, merge_tile_completeness_table
from   desispec.skymag import compute_skymag
from   desispec.efftime import compute_efftime
from   desispec.parallel import  stdouterr_redirected, use_mpi
from   desispec.util import parse_int_args
from desispec.workflow.tableio import load_table


def parse(options=None):
    parser = argparse.ArgumentParser(
                description="Calculate template S/N ratio for exposures")
    parser.add_argument('-o','--outfile', type=str, default=None, required=False,
                        help = 'Output summary file')
    parser.add_argument('--update', action = 'store_true',
                        help = 'Update pre-existing output summary file (replace or append)')
    parser.add_argument('--details-dir', type=str, default = None, required=False,
                        help = 'Dir. to write per-exposure per-camera files with per-fiber tSNR details')
    parser.add_argument('--prod', type = str, default = None, required=False,
                        help = 'Path to input reduction, e.g. /global/cfs/cdirs/desi/spectro/redux/blanc/,  ' +
                               'or simply prod version, like blanc, but requires env. variable DESI_SPECTRO_REDUX. ' +
                               'Default is $DESI_SPECTRO_REDUX/$SPECPROD.')
    parser.add_argument('-c', '--cameras', type = str, default = None, required=False,
                        help = 'Cameras to reduce (comma separated)')
    parser.add_argument('-e', '--expids', type = str, default = None, required=False,
                        help = 'Comma separated list of exp ids to process')
    parser.add_argument('-n', '--nights', type = str, default = None, required=False,
                        help = 'Comma, or colon separated list of nights to process. ex: 20210501,20210502 or 20210501:20210531')
    parser.add_argument('--recompute', action='store_true',
                        help = 'Recompute TSNR values')
    parser.add_argument('--alpha_only', action='store_true',
                        help = 'Only compute the alpha for tsnr.')
    parser.add_argument('--nproc', type = int, default = 1,
                        help = 'Multiprocessing.')
    parser.add_argument('-t', '--tile-completeness', type = str, default = None, required=False,
                        help = 'Output tile completeness table')
    parser.add_argument('--aux', type = str, default = None, required=False, nargs="*",
                        help = 'Path to auxiliary tiles tables, like /global/cfs/cdirs/desi/survey/observations/SV1/sv1-tiles.fits ' +
                               '(optional, will not affect exposures after SV1)')
    parser.add_argument('--gfa-proc-dir', type = str, default = None, required=False,
                        help = 'Directory where the GFA processing from Aaron can be found, like /global/cfs/cdirs/desi/survey/GFA/')
    parser.add_argument('--skymags', type = str, default = None, required=False,
                        help = 'path to table with sky magnitudes. Expected keys=NIGHT,EXPID,SKY_MAG_G,SKY_MAG_R,SKY_MAG_Z')
    parser.add_argument('--compute-skymags', action='store_true',
                        help = 'recompute sky mags')
    parser.add_argument('--mpi', action='store_true',
                        help = 'Paralleize TSNR recomputation across multiple nodes using mpi4py. Requires args.update to be False.')
    parser.add_argument('--add-badexp', action='store_true',
                        help = "Include entries for known bad exposures that weren't processed")

    args = None
    if options is None:
        args = parser.parse_args()
    else:
        args = parser.parse_args(options)
    return args


def update_targ_info(entry, targ_in):
    """
    Updates and returns information from targeting dictionary derived from a header or table row.

    Args:
        entry, dict: the currently known information about an exposure
        targ_in, dict: dictionary with a single value per key (either built from header
                    or from a Table).

    Returns:
        entry, dict: the currently known information about an exposure updated with information about
                    SURVEY, GOALTYPE, FAPRGRM, FAFLAVOR, MINTFRAC, and GOALTIME if available in targ_in
    """
    targ_keys = list(targ_in.keys())
    targ_defs = {
        "SURVEY" : "unknown",
        "GOALTYPE" : "unknown",
        "FAPRGRM" : "unknown",
        "FAFLAVOR" : "unknown",
        "MINTFRAC" : 0.9,
        "GOALTIME" : 0.,
    }

    for key,default in targ_defs.items():
        ## If the key is defined to anything other than default, don't overwrite
        if key in entry.keys() and entry[key] != default:
            continue
        ## Otherwise if the key is in the target dict, assign that
        elif key in targ_keys:
            if isinstance(default,str):
                entry[key] = targ_in[key].strip().lower()
            else:
                entry[key] = targ_in[key]
        ## Otherwise if the key isn't in the entry, assign it the default
        elif key not in entry.keys():
            entry[key] = default

    if "FAFLAVOR" in targ_keys :
        faflavor=targ_in["FAFLAVOR"].strip().lower()

        if entry["FAPRGRM"] == "unknown"  :
            entry["FAPRGRM"] = faflavor.replace("sv1","").replace("sv2","").replace("cmx","")

        if entry["SURVEY"]=="unknown" :
            #- despite the name, "cmxlrgqso" and "cmxelg" are sv1 tiles
            if faflavor.find("sv1")>=0. or faflavor in ('cmxlrgqso', 'cmxelg') :
                entry["SURVEY"]="sv1"
            elif faflavor.find("sv2")>=0. :
                entry["SURVEY"]="sv2"
            elif faflavor.find("cmx")>=0. :
                entry["SURVEY"]="cmx"

    entry['GOALTYPE'] = faflavor2program(entry['FAFLAVOR'])
    if entry["GOALTYPE"]=="unknown" :
        if entry["FAPRGRM"].find("qso")>=0. or entry["FAPRGRM"].find("lrg")>=0. or \
           entry["FAPRGRM"].find("elg")>=0. or entry["FAPRGRM"].find("dark")>=0. :
            entry["GOALTYPE"]="dark"
        elif entry["FAPRGRM"].find("mws")>=0. or entry["FAPRGRM"].find("bgs")>=0. or entry["FAPRGRM"].find("bright")>=0. :
            entry["GOALTYPE"]="bright"

    if entry['FAPRGRM'].startswith('dith'):
        if entry['SURVEY'] == 'unknown':
            entry['SURVEY'] = 'cmx'

    if entry['GOALTYPE'] in ('unknown', 'other'):
        entry['GOALTYPE'] = faflavor2program(entry['FAFLAVOR'])

    if entry['GOALTYPE'] in ['dark1b', 'bright1b']:
        entry['GOALTYPE'] = entry['GOALTYPE'].replace('1b', '')

    return entry


def compute_tsnr_values(cframe_filename,cframe_hdulist,night,expid,camera,specprod_dir, alpha_only=False) :
    """
    Computes TSNR values
    Args:
       cframe_filename: str, cframe file path
       cframe_hdulist: astropy.fits.HDUlist object
       night: int
       expid: int
       camera: str
       specprod_dir: str, production directory
       alpha_only: bool, set to True to only recompute alpha

    Returns: astropy.table.Table obkect with TSNR values
    """
    log = get_logger()

    calib  = findfile('fluxcalib', night=night, expid=expid,
                      camera=camera, specprod_dir=specprod_dir,
                      readonly=True)

    hdr0 = cframe_hdulist[0].read_header()
    hdr1 = cframe_hdulist[1].read_header()

    if 'FIBERFLT' not in hdr0:
        log.critical('Failed to find FIBERFLT in cframe hdr of {:08d} on {}'.format(expid, night))
        log.info('Trying CalibFinder.')

        flat = CalibFinder([hdr0, hdr1]).findfile('FIBERFLAT')

    else:
        flat = hdr0['FIBERFLT']

    if 'SPECPROD' in flat:
        flat = flat.replace('SPECPROD', specprod_dir)
    if 'SPCALIB' in flat:
        hdr  = fitsio.read_header(cframe_filename)
        flat = flat.replace('SPCALIB', getdep(hdr, 'DESI_SPECTRO_CALIB'))

    if os.path.exists(flat):
        log.info('Retrieved {}'.format(flat))
    else:
        raise ValueError('Failed on flat retrieval for {} on {}.'.format(expid, night))

    iin = cframe_filename.replace('cframe', 'frame')
    sky = cframe_filename.replace('cframe', 'sky')
    psf = cframe_filename.replace('cframe', 'psf')

    cframe=read_frame(cframe_filename, skip_resolution=True)
    frame=read_frame(iin, skip_resolution=True)
    fiberflat=read_fiberflat(flat)
    fluxcalib=read_flux_calibration(calib)
    skymodel=read_sky(sky)

    results, alpha = calc_tsnr2(cframe, frame, fiberflat=fiberflat,
                                skymodel=skymodel, fluxcalib=fluxcalib, alpha_only=alpha_only)

    table=Table()
    for k in results:
        table[k] = results[k].astype(np.float32)
    table["TSNR2_ALPHA_"+camera[0].upper()] = np.repeat(alpha,len(frame.flux))

    return table


def update_table(table1,table2,keys) :
    """ Replace or append rows of table1 with content of table2 indexed by keys

    Args:
        table1: astropy.table.Table
        table2: astropy.table.Table
        keys: list of str

    Returns astropy.table.Table
    """

    log = get_logger()

    v1=table1[keys[0]]
    v2=table2[keys[0]]
    if len(keys)>1 : # I don't know how of a better generic way to create a joined index
        v1=v1.astype("str")
        v2=v2.astype("str")
        for k in keys[1:] :
            v1=[i+j for i,j in zip(v1,table1[k].astype(str))]
            v2=[i+j for i,j in zip(v2,table2[k].astype(str))]

    replace = np.isin(v1,v2)
    keep = ~replace

    log.info("keep {} entries, replace {}, add {}".format(np.sum(keep),np.sum(replace),len(table2)-np.sum(replace)))

    if np.sum(keep)>0 :
        for k in table2.dtype.names :
            if k not in table1.dtype.names :
                log.info("add new column {}".format(k))
                table1[k]=np.zeros(len(table1),dtype=table2[k].dtype)
        return vstack([table1[keep],table2])
    else :
        return table2

def compute_summary_tables(summary_rows, preexisting_tsnr2_expid_table,
                           preexisting_tsnr2_frame_table, specprod_dir,
                           add_badexp=True, nights=None):
    """Compute exposure and frame summary tables.

    Parameters
    ----------
    summary_rows : :class:`list`
        A list of :class:`dict` describing columns. In nightly operations, this
        is an empty list.
    preexisting_tsnr2_expid_table : :class:`~astropy.table.Table`
        Update this exposure table if the input is not ``None``. In nightly operations,
        this will always be set.
    preexisting_tsnr2_frame_table : :class:`~astropy.table.Table`
        Update this frame table if the input is not ``None``. In nightly operations,
        this will always be set.
    specprod_dir : :class:`str`
        Production directory.
    add_badexp : :class:`bool`, optional
        Add known bad exposures that weren't processed. Default is ``True``.
    nights : :class:`list`, optional
        List of nights for which to compute bad exposures summary.

    Returns
    -------
    :class:`tuple`
        A tuple containing the updated exposure and frames tables.

    Raises
    ------
    KeyError
        If the columns cannot be placed in the proper order.
    """
    log = get_logger()

    #- Create camera summary table; specify names to preserve column order
    ## Get copy of the columns, one to preserve the original cols and one to add any new cols encountered
    if len(summary_rows) == 0 and preexisting_tsnr2_expid_table is not None\
       and preexisting_tsnr2_frame_table is not None:
        cam_summary = Table(names=preexisting_tsnr2_frame_table.colnames,
                            dtype=preexisting_tsnr2_frame_table.dtype)
        cam_summary.meta['EXTNAME'] = 'FRAMES'
        exp_summary = Table(names=preexisting_tsnr2_expid_table.colnames,
                            dtype=preexisting_tsnr2_expid_table.dtype)
        exp_summary.meta['EXTNAME'] = 'EXPOSURES'
    else:
        colnames = list(summary_rows[0].keys())
        orig_colnames = list(summary_rows[0].keys())
        first_print = True
        for ii,row in enumerate(summary_rows):
            missing = [key for key in orig_colnames if key not in row.keys()]
            additional =  [key for key in row.keys() if key not in orig_colnames]
            if len(missing) > 0 or len(additional) > 0:
                ## Only print the column names once to reduce verbosity
                if first_print:
                    log.warning(f"First row of summary_rows had colnames={orig_colnames}")
                    first_print = False
                ## List missing columns compared to the first row (if any)
                if len(missing) > 0:
                    log.warning(f"Row {ii} differs from first, row is missing columns={missing}" )
                ## Could have both additional and missing columns, so this is independent check. List additional
                if len(additional) > 0:
                    log.warning(f"Row {ii} differs from first, row has additional columns={additional}" )

            ## Beyond the differences listed above, make sure to keep a complete list of column names
            ## for the table generation
            newcols =  [key for key in row.keys() if key not in colnames]
            if len(newcols) > 0:
                log.warning(f"Identified the following new columns: {newcols}. Adding to colnames.")
                colnames.extend(newcols)


        #- specify names to preserve column order
        masked_cam_summary = Table(rows=summary_rows, names=colnames)
        for col in masked_cam_summary.colnames:
            if isinstance(masked_cam_summary[col],Table.MaskedColumn):
                if col in ['NIGHT', 'EXPID', 'TILEID']:
                    log.warning(f'Column {col} identified as masked column, filling masked rows with value 0 .')
                    masked_cam_summary[col].set_fill_value(0)
                elif col in ['SURVEY', 'FAPRGRM', 'FAFLAVOR', 'GOALTYPE']:
                    log.warning(f'Column {col} identified as masked column, filling masked rows with value "unknown".')
                    masked_cam_summary[col].set_fill_value('unknown')
                else:
                    log.warning(f'Column {col} identified as masked column, filling masked rows with value 0.0.')
                    masked_cam_summary[col].set_fill_value(0.0)
        cam_summary = masked_cam_summary.filled()

        cam_summary.meta['EXTNAME'] = 'FRAMES'

        tsnr2_colnames = [x for x in cam_summary.colnames if x.startswith('TSNR2_')]

        #- which rows go with which petals:
        ispetal = dict()
        for petal in range(10):
            ii  = (cam_summary['CAMERA'] == 'b'+str(petal))
            ii |= (cam_summary['CAMERA'] == 'r'+str(petal))
            ii |= (cam_summary['CAMERA'] == 'z'+str(petal))
            ispetal[petal] = ii

        #- Distill into exposure summary per-petal
        rows = list()

        #- The unique key is the exposure Id
        for expid in np.unique(cam_summary['EXPID'] ) :
            ii = (cam_summary['EXPID'] == expid)

            row = dict()
            for k in ["NIGHT","EXPID","TILEID","TILERA","TILEDEC","MJD",
                      "AIRMASS","EBV","SEEING_ETC","EFFTIME_ETC","SURVEY","GOALTYPE",
                      "MINTFRAC","FAPRGRM","FAFLAVOR","GOALTIME"] :
                row[k]=cam_summary[k][ii][0]
            # mean
            row["EXPTIME"]=float(np.mean(cam_summary['EXPTIME'][ii]))

            #- Per EXPID TSNR2 is summed over cameras, averaged over petals
            for colname in tsnr2_colnames:
                tsnr2_petal = list()  #- TSNR^2 summed over cameras, per petal
                for petal in range(10):
                    jj = ii & ispetal[petal]
                    if np.any(jj):
                        tsnr2_petal.append(np.sum(cam_summary[colname][jj]))
                row[colname] = np.mean(tsnr2_petal) # mean of petals (each being median of valid fibers)

            rows.append(row)

        colnames = list(rows[0].keys())
        exp_summary = Table(rows=rows, names=colnames)
        exp_summary.meta['EXTNAME'] = 'EXPOSURES'

    if add_badexp:
        # check the production exposure table
        if nights is None:
            nights = np.unique(exp_summary["NIGHT"])
        for night in nights:
            pipe_exptab_filename = findfile("exposure_table", night=int(night), readonly=True)
            if not os.path.isfile(pipe_exptab_filename) :
                log.warning("cannot find prod. exposure table {}".format(pipe_exptab_filename))
                continue
            pipe_exptab = load_table(tablename=pipe_exptab_filename, tabletype="exptable", suppress_logging=True)
            missing = (~(np.isin(pipe_exptab["EXPID"],exp_summary["EXPID"])))&(pipe_exptab["TILEID"]>0)
            nmissing=np.sum(missing)
            if nmissing>0:
                for i in np.where(missing)[0] :
                    #- If the exposure was supposed to be processed but is unexpectedly
                    #- missing from prod, log error but don't add it to table because
                    #- we are in an ambiguous/error state of processing
                    if pipe_exptab["LASTSTEP"][i] == 'all':
                        tileid = pipe_exptab['TILEID'][i]
                        expid = pipe_exptab['EXPID'][i]
                        log.error(f'Tile {tileid} night {night} expid {expid} missing from prod; not adding to exposure table ')
                        continue

                    entry={}
                    entry["EXPID"]=pipe_exptab["EXPID"][i]
                    entry["NIGHT"]=night
                    entry["TILEID"]=pipe_exptab["TILEID"][i]
                    entry["MJD"] = pipe_exptab["MJD-OBS"][i]
                    entry["AIRMASS"]=0.
                    entry["EBV"]=0.
                    # EBVFAC from https://github.com/desihub/fiberassign/blob/466a99b926892be60c3c9679f521785e865ce315/py/fiberassign/fba_launch_io.py#L1658
                    if "EBVFAC" in pipe_exptab.dtype.names:
                        entry["EBV"] = 2.5 * np.log10(pipe_exptab["EBVFAC"][i]) / 2.165
                    entry["EXPTIME"]=pipe_exptab["EXPTIME"][i]
                    entry['SEEING_ETC']=0.
                    entry["EFFTIME_ETC"] = 0.
                    if "EFFTIME_ETC" in pipe_exptab.dtype.names:
                        entry["EFFTIME_ETC"] = pipe_exptab["EFFTIME_ETC"][i]

                    # SURVEY, GOALTYPE, FAPRGRM, FAFLAVOR, MINTFRAC, GOALTIME
                    targ_dict = {key : pipe_exptab[key][i] for key in pipe_exptab.dtype.names}
                    entry = update_targ_info(entry, targ_dict)
                    entry['PROGRAM'] = faflavor2program(entry['FAFLAVOR'])

                    # Get TILERA, TILEDEC from fiberassign file (*.fits or *.fits.gz)
                    # as of 20210624: SURVEY present in table_exposure since 20210611 only
                    dirname = os.path.dirname(findfile("raw", night, entry["EXPID"], readonly=True))
                    fns = glob.glob( os.path.join(dirname, "fiberassign-{:06d}.fits*".format(entry["TILEID"])) )
                    if len(fns) > 0:
                        fiberassign_header = fits.getheader(fns[0], 0)
                        targ_dict = {key : fiberassign_header[key] for key in fiberassign_header.keys()}
                        if 'SURVEY' not in targ_dict and 'FA_SURV' in targ_dict:
                            targ_dict['SURVEY'] = targ_dict['FA_SURV']

                        entry = update_targ_info(entry, targ_dict)
                        entry['TILERA'] = fiberassign_header['TILERA']
                        entry['TILEDEC'] = fiberassign_header['TILEDEC']
                    else:
                        log.error(f'No fiberassign*.fits* file in {dirname}')
                        entry["TILERA"] = -999.0
                        entry["TILEDEC"] = -999.0

                    for key in ("SKY_MAG_G_SPEC", "SKY_MAG_R_SPEC", "SKY_MAG_Z_SPEC"):
                        entry[key] = np.nan

                    vals=[]
                    for j,k in enumerate(exp_summary.dtype.names) :
                        if k in entry :
                            vals.append(entry[k])
                        else :
                            vals.append(0.)

                    log.warning("Adding missing exposure {}".format(vals))
                    exp_summary.add_row(vals)

                    # same for other table
                    for cam in decode_camword(difference_camwords(pipe_exptab['CAMWORD'][i],pipe_exptab["BADCAMWORD"][i])):
                        entry["CAMERA"]=cam
                        vals=[]
                        for j,k in enumerate(cam_summary.dtype.names) :
                            if k in entry :
                                vals.append(entry[k])
                            else :
                                vals.append(0.)
                        log.debug("Adding missing exposure.cam {}".format(vals))
                        cam_summary.add_row(vals)


    log.info("Add effective times from TSNR2 values")

    exp_summary['LRG_EFFTIME_DARK']   = tsnr2_to_efftime(exp_summary['TSNR2_LRG'],"LRG")
    exp_summary['ELG_EFFTIME_DARK']   = tsnr2_to_efftime(exp_summary['TSNR2_ELG'],"ELG")
    exp_summary['BGS_EFFTIME_BRIGHT'] = tsnr2_to_efftime(exp_summary['TSNR2_BGS'],"BGS")

    for program in ['DARK','BRIGHT','BACKUP']:
        if f'TSNR2_GPB{program}' in exp_summary.dtype.names:
            exp_summary[f'GPB_EFFTIME_{program}']   = tsnr2_to_efftime(exp_summary[f'TSNR2_GPB{program}'],f'GPB{program}')
        else:
            exp_summary[f'TSNR2_GPB{program}'] = 0.
            exp_summary[f'GPB_EFFTIME_{program}'] = 0.

    if 'TSNR2_LYA' in exp_summary.dtype.names:
        exp_summary['LYA_EFFTIME_DARK'] = tsnr2_to_efftime(exp_summary['TSNR2_LYA'],"LYA")

    # EFFTIME SPEC FOR GOALTYPE dark determined by TSNR2_ELG before Sept 1, 2021 and TSNR2_LRG after
    before_sept2021 = (exp_summary["NIGHT"]<20210901)
    after_sept2021 = ~before_sept2021
    exp_summary['EFFTIME_SPEC'] = np.zeros(len(exp_summary))
    exp_summary['EFFTIME_SPEC'][before_sept2021] = exp_summary['ELG_EFFTIME_DARK'][before_sept2021]
    exp_summary['EFFTIME_SPEC'][after_sept2021]  = exp_summary['LRG_EFFTIME_DARK'][after_sept2021]

    # EFFTIME SPEC FOR GOALTYPE bright determined by TSNR2_BGS.
    ii = (exp_summary['GOALTYPE']=='bright')
    exp_summary['EFFTIME_SPEC'][ii] = exp_summary['BGS_EFFTIME_BRIGHT'][ii]

    # EFFTIME SPEC FOR GOALTYPE backup determined by TSNR2_GPBBACKUP.
    # If no 'GPB_EFFTIME_BACKUP' then default to the BGS time
    ii = (exp_summary['GOALTYPE']=='backup')
    if 'GPB_EFFTIME_BACKUP' in exp_summary.dtype.names and np.any(ii) and np.any(exp_summary['EFFTIME_SPEC'][ii]>0.):
        exp_summary['EFFTIME_SPEC'][ii] = exp_summary['GPB_EFFTIME_BACKUP'][ii]
    else:
        exp_summary['EFFTIME_SPEC'][ii] = exp_summary['BGS_EFFTIME_BRIGHT'][ii]

    if preexisting_tsnr2_expid_table is not None :
        log.debug("Update to preexisting")

        # backward compatibility issue
        for k in ["SURVEY","GOALTYPE","FAPRGRM","FAFLAVOR"] :
            nentries=len(preexisting_tsnr2_expid_table)
            if k not in preexisting_tsnr2_expid_table.dtype.names :
                log.warning("adding column {}".format(k))
                preexisting_tsnr2_expid_table[k] = np.array(np.repeat("unknown",nentries),dtype='<U16')
        for k in ["GOALTIME","MINTFRAC","SEEING_ETC","EFFTIME_ETC","AIRMASS","EBV","LRG_EFFTIME_DARK"]:
            if k not in preexisting_tsnr2_expid_table.dtype.names :
                log.warning("adding column {}".format(k))
                preexisting_tsnr2_expid_table[k] = np.array(np.repeat(0.,nentries),dtype=float)

        exp_summary = update_table(preexisting_tsnr2_expid_table,exp_summary,["EXPID"])

    if preexisting_tsnr2_frame_table is not None :
        log.debug("Update to preexisting")
        cam_summary = update_table(preexisting_tsnr2_frame_table,cam_summary,["EXPID","CAMERA"])

    # sort table rows
    ii = np.argsort(exp_summary['EXPID'])
    exp_summary = exp_summary[ii]
    vals=["{}-{}".format(e,c) for e,c in zip(cam_summary['EXPID'],cam_summary['CAMERA'])]
    ii = np.argsort(vals)
    cam_summary = cam_summary[ii]

    if 'TSNR2_ALPHA' in exp_summary.dtype.names : exp_summary.remove_column('TSNR2_ALPHA')

    #- standardizing N>>1 sv1-era FAPRGRM/FAFLAVOR into dark/bright/backup/other
    exp_summary['PROGRAM'] = faflavor2program(exp_summary['FAFLAVOR'])

    neworder=['NIGHT','EXPID','TILEID','TILERA','TILEDEC','MJD',
              'SURVEY','PROGRAM','FAPRGRM','FAFLAVOR','EXPTIME','EFFTIME_SPEC','GOALTIME','GOALTYPE',
              'MINTFRAC','AIRMASS','EBV','SEEING_ETC','EFFTIME_ETC','TSNR2_ELG','TSNR2_QSO','TSNR2_LRG','TSNR2_LYA',
              'TSNR2_BGS','TSNR2_GPBDARK','TSNR2_GPBBRIGHT','TSNR2_GPBBACKUP','LRG_EFFTIME_DARK','ELG_EFFTIME_DARK',
              'BGS_EFFTIME_BRIGHT','LYA_EFFTIME_DARK','GPB_EFFTIME_DARK','GPB_EFFTIME_BRIGHT','GPB_EFFTIME_BACKUP',
              'TRANSPARENCY_GFA','SEEING_GFA','FIBER_FRACFLUX_GFA','FIBER_FRACFLUX_ELG_GFA','FIBER_FRACFLUX_BGS_GFA',
              'FIBERFAC_GFA','FIBERFAC_ELG_GFA','FIBERFAC_BGS_GFA','AIRMASS_GFA','SKY_MAG_AB_GFA','SKY_MAG_G_SPEC',
              'SKY_MAG_R_SPEC','SKY_MAG_Z_SPEC','EFFTIME_GFA','EFFTIME_DARK_GFA','EFFTIME_BRIGHT_GFA',
              'EFFTIME_BACKUP_GFA']

    if not np.all(np.isin(exp_summary.dtype.names, neworder)):
        missing = ~np.isin(exp_summary.dtype.names, neworder)
        msg = "Missing keys, '{}' in new order list!".format(np.array(exp_summary.dtype.names)[missing])
        log.critical(msg)
        log.critical("new order list:",sorted(neworder))
        log.critical("current table:",sorted(exp_summary.dtype.names))
        raise KeyError(msg)
    newtable=Table()
    newtable.meta=exp_summary.meta
    for k in neworder :
        if k in exp_summary.dtype.names :
            newtable[k]=exp_summary[k]
        else :
            newtable[k]=np.zeros(len(exp_summary))
    exp_summary=newtable

    return exp_summary, cam_summary


def write_summary_tables(tsnr2_expid_table,tsnr2_frame_table,output_fits_filename,output_csv_filename=None) :
    """Writes a summary fits file.

    Args:
        tsnr2_expid_table: astropy.table.Table
        tsnr2_frame_table: astropy.table.Table
        output_fits_filename: str, output fits filename
        output_csv_filename: str, output csv filename

    Returns nothing
    """
    log = get_logger()

    hdus = fits.HDUList()
    hdus.append(fits.convenience.table_to_hdu(tsnr2_expid_table))
    hdus.append(fits.convenience.table_to_hdu(tsnr2_frame_table))

    tmpfile = get_tempfilename(output_fits_filename)
    hdus.writeto(tmpfile, overwrite=True)
    os.rename(tmpfile,  output_fits_filename)

    log.info('Successfully wrote {}'.format(output_fits_filename))

    # also write first table as csv
    if output_csv_filename is not None :
        # make a copy before modifying precision of columns
        tsnr2_expid_table = tsnr2_expid_table.copy(copy_data=True)
        for k in tsnr2_expid_table.keys() :
            if k.find("TIME")>=0 or k.find("SNR2")>=0 :
                try :
                    tsnr2_expid_table[k]=np.around(tsnr2_expid_table[k].astype(float),1)
                except ValueError as e :
                    log.warning("cannot change number of digits for col {}: {}".format(k,e))
            if k.find("SEEING")>=0 or k.find("FWHM")>=0 or k.find("FIBER_FRACFLUX")>=0 or k.find("FIBERFAC")>=0 or \
               k.find("TRANS")>=0 or k.find("MAG")>=0 or k.find("AIRMASS")>=0 or k.find("EBV")>=0 :
                try :
                    tsnr2_expid_table[k]=np.around(tsnr2_expid_table[k].astype(float),3)
                except ValueError as e :
                    log.warning("cannot change number of digits for col {}: {}".format(k,e))
        tsnr2_expid_table.write(output_csv_filename,overwrite=True)
        log.info('Successfully wrote {}'.format(output_csv_filename))

def func(prod,night,expid,camera,recompute,alpha_only,details_dir) :
    """
    Compute TSNR for this night exposure camera (code in separate function for multiprocessing)

    Args:
        prod: str, production dir name
        night: int, night
        expid: int, expid
        camera: str, camera
        recompute: bool, recompute tsnr if true even if present in frame
        alpha_only: bool, only recompute alpha
        details_dir: str or None, save details per frame in this directory if not None

    Returns a dictionary with NIGHT,EXPID,TILEID,EXPTIME,CAMERA, and the TSNR2 values for this camera
    """
    log = get_logger()

    cframe_filename = findfile('cframe', night, expid, camera,
            specprod_dir=prod, readonly=True)
    if not os.path.isfile(cframe_filename) :
        return None

    cframe_hdulist = fitsio.FITS(cframe_filename)
    hdr    = cframe_hdulist[0].read_header()
    flavor = hdr['FLAVOR']
    if flavor != 'science':
        return None

    log.info("Processing {}".format(cframe_filename))

    prog   = hdr['PROGRAM']
    tileid = hdr['TILEID']

    table = None
    compute_tsnr = True

    #- check if already computed in cframe
    if "SCORES" in cframe_hdulist and not recompute:
        table = Table(cframe_hdulist["SCORES"].read())
        key = "TSNR2_ELG_"+camera[0].upper()
        if key in table.colnames :
            log.debug("Use TSNR values in cframe file")
            compute_tsnr = False

    #- check if already computed in details file
    table_output_filename = None
    if details_dir is not None and not recompute :
        table_output_filename= f'{details_dir}/{night}/{expid:08d}/tsnr-{camera}-{expid:08d}.fits'
        if os.path.isfile(table_output_filename):
            log.debug(f'Using previously generated {table_output_filename}')
            table = Table.read(table_output_filename)
            compute_tsnr = False

    #- compute tsnr
    if compute_tsnr :
        tsnr_table = compute_tsnr_values(cframe_filename,cframe_hdulist,night,expid,camera,specprod_dir=prod, alpha_only=alpha_only)

        if table is None :
            table = tsnr_table
        else :
            for k in tsnr_table.columns :
                table[k] = tsnr_table[k]

    fibermap = cframe_hdulist["FIBERMAP"].read()

    table['FIBER']       = fibermap['FIBER']
    table['TARGETID']    = fibermap['TARGETID']
    table.meta['NIGHT']  = night
    table.meta['EXPID']  = expid
    table.meta['TILEID'] = tileid
    table.meta['CAMERA'] = camera
    table.meta['EXTNAME'] = 'TSNR2'

    #- Write per-expid per-camera output if requested
    if compute_tsnr and table_output_filename is not None :
        Path(os.path.dirname(table_output_filename)).mkdir(parents=True,exist_ok=True)
        tmpfile = get_tempfilename(table_output_filename)
        table.write(tmpfile, format='fits', overwrite=True)
        os.rename(tmpfile, table_output_filename)
        log.debug('Successfully wrote {}.'.format(table_output_filename))

    #- Append to summary.
    entry = dict()
    entry['NIGHT'] = np.int32(night)
    entry['EXPID'] = np.int32(expid)
    entry['TILEID'] = np.int32(tileid)
    entry['TILERA'] = hdr['TILERA']
    entry['TILEDEC'] = hdr['TILEDEC']
    entry['MJD'] = hdr['MJD-OBS']
    entry['EXPTIME'] = np.float32(hdr['EXPTIME'])
    entry['AIRMASS'] = np.float32(hdr['AIRMASS'])
    entry['EBV'] = 0.
    etcfilename = findfile('etc', night=entry['NIGHT'], expid=entry['EXPID'],
                           readonly=True)
    if 'EBV' in fibermap.dtype.names:
        sel = (fibermap["EBV"] > 0)
        if sel.sum() > 0:
            entry['EBV'] = np.median(fibermap['EBV'][sel])
    if "ACQFWHM" in hdr :
        # FWHM of the guide star PSF, in arcseconds, measured in the initial acquisition image.
        entry['SEEING_ETC'] = np.float32(hdr['ACQFWHM'])

        log.debug('Retrieved ACQFWHM from frame hdr.')

    elif os.path.exists(etcfilename):
        with open(etcfilename) as etc_file:
            etcdata = json.load(etc_file)
        try:
            entry['SEEING_ETC'] = np.float32(etcdata['expinfo']['acq_fwhm'])
        except:
            entry['SEEING_ETC'] = np.float32(0.)

        log.debug('Retrieved ACQFWHM from etc json.')
    else :
        entry['SEEING_ETC'] = 0.

        log.debug('Failed to retrieved ACQFWHM from frame hdr.')

    ##  https://github.com/desihub/desietc/blob/main/header.md
    if "ETCTEFF" in hdr :
        entry['EFFTIME_ETC'] = np.float32(hdr['ETCTEFF'])

        log.debug('Retrieved ETCTEFF from frame hdr.')

    elif os.path.exists(etcfilename):
        with open(etcfilename) as etc_file:
            etcdata = json.load(etc_file)

        try:
            entry['EFFTIME_ETC'] = np.float32(etcdata['expinfo']['efftime'])
        except:
            entry['EFFTIME_ETC'] = np.float32(0.)

        log.debug('Retrieved ETCTEFF from etc json.')

    else:
        entry['EFFTIME_ETC'] = np.float32(0.)

        log.debug('Failed to retrieve ETCTEFF from etc json.')

    #- treat NaN the same as missing data, i.e. 0.0
    #- e.g. 20210516/00088872/etc-00088872.json etcdata['expinfo']['efftime']
    if np.isnan(entry['EFFTIME_ETC']):
        entry['EFFTIME_ETC'] = 0.0

    entry['CAMERA'] = camera
    for key in table.colnames:
        if key.startswith('TSNR2_'):
            #- TSNR2_{TRACER}_{BAND} -> TSNR2_{TRACER}
            short_key = '_'.join(key.split('_')[0:2])
            ok=(table[key]!=0)
            if np.sum(ok)>0 :
                entry[short_key]=np.median(table[key][ok]).astype(np.float32)
            else :
                entry[short_key]=0.
        if key.startswith('TSNR2_ALPHA'):
            entry['TSNR2_ALPHA'] = np.median(table[key]).astype(np.float32)

    #- Add information from targeting (SURVEY, GOALTYPE, FAPRGRM, FAFLAVOR, MINTFRAC, GOALTIME)
    fibermap_header = cframe_hdulist["FIBERMAP"].read_header()
    targ_dict = {key : fibermap_header[key] for key in fibermap_header.keys()}
    entry = update_targ_info(entry, targ_dict)

    cframe_hdulist.close()
    return entry

def _func(arg) :
    """
    Wrapper function to TSNR values for multiprocessing
    """
    return func(**arg)

def func2(night,expid,specprod_dir) :
    """
    Wrapper function to compute_skymag for multiprocessing
    """
    log = get_logger()
    mags = compute_skymag(night,expid,specprod_dir)
    entry = {'NIGHT':night,'EXPID':expid,'SKY_MAG_G':mags[0],'SKY_MAG_R':mags[1],'SKY_MAG_Z':mags[2]}
    log.info(entry)
    return(entry)

def _func2(arg) :
    """
    Wrapper function to compute_skymag for multiprocessing
    """
    return func2(**arg)

def add_skymags_columns(exposure_table,skymags_table) :
    """
    Add sky magnitudes to exposure table
    Args :
       exposure_table, astropy.table.Table with column EXPID
       skymags_table astropy.table.Table with columns EXPID,SKY_MAG_G,SKY_MAG_R,SKY_MAG_Z, AB mags/arcsec2 in decam bands
    returns nothing (updates input exposure table)
    """
    e2i = {e:i for i,e in enumerate(skymags_table["EXPID"])}
    jj=[]
    ii=[]
    for j,e in enumerate(exposure_table["EXPID"]) :
        if e in e2i :
            jj.append(j)
            ii.append(e2i[e])
        cols=["SKY_MAG_G","SKY_MAG_R","SKY_MAG_Z"]
        for col in cols :
            if col in skymags_table.dtype.names :
                col2=col+"_SPEC"
                if col2 not in exposure_table.dtype.names :
                    exposure_table[col2] = np.zeros(len(exposure_table),dtype=float)
                exposure_table[col2][jj] = skymags_table[col][ii]


def main():
    """Entry-point for command-line scripts.

    Returns
    -------
    :class:`int`
        An integer suitable for passing to :func:`sys.exit`.
    """
    log = get_logger()

    args=parse()
    if args.prod is None:
        args.prod = specprod_root()
    elif args.prod.find("/")<0 :
        args.prod = specprod_root(args.prod)

    if not args.outfile.endswith(".fits") :
        log.critical("Output filename '{}' is incorrect. It has to end with '.fits'.".format(args.outfile))
        return 1

    if args.mpi and use_mpi():
        from mpi4py import  MPI

        comm  = MPI.COMM_WORLD
        rank  = comm.Get_rank()
        size  = comm.Get_size()

        multinode = True

        if args.update:
            log.critical('--update is not supported for the --mpi option.  Remove --update.')
            raise RuntimeError()

        if not args.recompute:
            log.critical('need --recompute for the --mpi option.  Add --recompute.')
            raise RuntimeError()

        if (args.tile_completeness is not None):
            if os.path.isfile(args.tile_completeness):
                log.critical('Tile file exists.  Provide a new tile file path..')
                raise RuntimeError()

    else:
        comm  = None
        rank  = 0
        size  = 1

        multinode = False

    if rank == 0:
        log.info('outfile = {}'.format(args.outfile))
        if args.details_dir is not None :
            log.info('details-dir = {}'.format(args.details_dir))

        log.info('prod = {}'.format(args.prod))

    if args.cameras is None:
        petals  = np.arange(10).astype(str)
        cameras = [x[0] + x[1] for x in itertools.product(['b', 'r', 'z'], petals.astype(str))]
    else:
        cameras = args.cameras.split(',')

    if args.expids is not None:
        expids = [np.int(x) for x in args.expids.split(',')]
    else:
        expids = None

    if args.nights is None:
        dirnames = sorted(glob.glob('{}/exposures/*'.format(args.prod)))
        nights=[]
        for dirname in dirnames :
            try :
                night=int(os.path.basename(dirname))
                nights.append(night)
            except ValueError as e :
                log.warning("ignore {}".format(dirname))
    else :
        nights = parse_int_args(args.nights)

    if rank == 0:
        log.info('cameras = {}'.format(cameras))
        log.info("nights = {}".format(nights))
        if expids is not None :
            log.info('expids = {}'.format(expids))

    preexisting_tsnr2_expid_table = None
    preexisting_tsnr2_frame_table = None
    ## Check for existing files
    ## If opening returns a keyerror due to misnamed HDU, try the old naming convention
    if args.update and os.path.isfile(args.outfile) :
        log.info("Will append pre-existing table {}".format(args.outfile))
        try:
            preexisting_tsnr2_expid_table = read_table(args.outfile,"EXPOSURES")
        except KeyError:
            log.warning("Couldn't find the EXPOSURES HDU, looking for old format name TSNR2_EXPID")
            preexisting_tsnr2_expid_table = read_table(args.outfile,"TSNR2_EXPID")
        try:
            preexisting_tsnr2_frame_table = read_table(args.outfile,"FRAMES")
        except KeyError:
            log.warning("Couldn't find the FRAMES HDU, looking for old format name TSNR2_FRAME")
            preexisting_tsnr2_frame_table = read_table(args.outfile,"TSNR2_FRAME")

    # starting computing
    # one night at a time

    summary_rows  = list()

    def one_night(count, night, rows, tmpfilename):
        """Prepare input data for one night for passing to :func:`func`,
        which performs the actual TSNR2 calculations.

        Parameters
        ----------
        count : :class:`int`
            A MPI rank number. This may be unused inside this function.
        night : :class:`str`
            The night to analyze.
        rows : :class:`list`
            A list that will be passed to :func:`compute_summary_tables`.
        tmpfilename : :class:`str`
            A temporary filename.

        Returns
        -------
        :class:`int`
            An integer suitable for passing to :func:`sys.exit`.
        """
        dirnames = sorted(glob.glob('{}/exposures/{}/*'.format(args.prod,night)))
        night_expids=[]
        for dirname in dirnames :
            try :
                expid=int(os.path.basename(dirname))
                night_expids.append(expid)
            except ValueError as e :
                log.warning("ignore {}".format(dirname))
        if expids is not None :
            night_expids = np.intersect1d(expids,night_expids)
            if night_expids.size == 0 :
                return 0
        log.info("{} {}".format(night,night_expids))

        func_args = []
        for expid in night_expids :
            for camera in cameras:
                func_args.append({'prod':args.prod,'night':night,'expid':expid,'camera':camera,
                                  'recompute':args.recompute,'alpha_only':args.alpha_only,'details_dir':args.details_dir})

        if args.nproc == 1 :
            log.info(f"--nproc={args.nproc} so not multiprocessing")
            for func_arg in func_args :
                entry = func(**func_arg)
                if entry is not None :
                    rows.append(entry)
        else :
            log.info("Multiprocessing with {} procs".format(args.nproc))
            with multiprocessing.Pool(args.nproc) as pool:
                results  =  pool.map(_func, func_args)
                for entry in results :
                    if entry is not None :
                        rows.append(entry)

        # write result after every night.
        if len(rows) > 0 or args.add_badexp:
            try:
                tsnr2_expid_table, tsnr2_frame_table = compute_summary_tables(rows,
                                                                              preexisting_tsnr2_expid_table,
                                                                              preexisting_tsnr2_frame_table,
                                                                              args.prod, nights=[night],
                                                                              add_badexp=args.add_badexp)
            except KeyError:
                # Error messages have already been logged at this point.
                return 12
            write_summary_tables(tsnr2_expid_table, tsnr2_frame_table, output_fits_filename=tmpfilename)
            log.info("wrote {} entries in tmp file {}".format(len(rows), tmpfilename))

        return 0

    if comm is not None:
        comm.barrier()

        nights     = comm.bcast(nights, root=0)

        #- which nights should this rank process
        ranknights = nights[rank:len(nights):size]

        if len(ranknights) > 0:
            rank_summary_rows = list()

            log.info('Rank {:d} processes nights {} with nproc={}'.format(rank, ranknights, args.nproc))

            outdir = os.path.dirname(args.outfile)
            if outdir == "" :
                outdir = "."
            logfile = outdir+'/desi_tsnr_afterburner_{:d}_{:d}_{:d}.log'.format(ranknights[0], ranknights[-1], rank)

            if os.path.isfile(logfile):
                os.unlink(logfile)

            with stdouterr_redirected(to=logfile):
                for count, night in enumerate(ranknights):
                    tmpfilename = args.outfile.replace(".fits","_tmp_{:d}.fits".format(rank))

                    one_night_status = one_night(count, night, rank_summary_rows, tmpfilename)
                    #
                    # Not sure of the best way to exit in this case.
                    #
        else:
            log.warning('Rank {:d} has no nights to process'.format(rank))


        comm.barrier()

    else:
        for count, night in enumerate(nights):
            tmpfilename = args.outfile.replace(".fits", "_tmp.fits")

            one_night_status = one_night(count, night, summary_rows, tmpfilename)

            if one_night_status != 0:
                return one_night_status

    if (comm is None) | (rank == 0):
        log.info('Gathering TSNR tmp tables (over ranks: {}, {})'.format(multinode, size))

        tsnr2_expid = []
        tsnr2_frame = []

        for rank in range(size):
            if multinode:
                tmpfilename=args.outfile.replace('.fits', '_tmp_{:d}.fits'.format(rank))
            else:
                tmpfilename=args.outfile.replace('.fits', '_tmp.fits')

            try:
                tmp = fits.open(tmpfilename)

            except:
                log.warning('{} not found.'.format(tmpfilename))
                continue

            hdr = tmp[0].header

            tsnr2_expid.append(Table(tmp['EXPOSURES'].data))
            tsnr2_frame.append(Table(tmp['FRAMES'].data))

            # remove temporary file if successful
            if os.path.isfile(tmpfilename) :
                log.info("rm temporary file {}".format(tmpfilename))
                os.unlink(tmpfilename)

        tsnr2_expid_table = vstack(tsnr2_expid)
        tsnr2_frame_table = vstack(tsnr2_frame)

        tsnr2_expid_table.meta['EXTNAME'] = 'EXPOSURES'
        tsnr2_frame_table.meta['EXTNAME'] = 'FRAMES'

        if multinode:
            write_summary_tables(tsnr2_expid_table,tsnr2_frame_table,output_fits_filename=args.outfile.replace('.fits', '_tmp.fits'))

        if args.skymags is not None :
            skymags_table = Table.read(args.skymags)
            add_skymags_columns(tsnr2_expid_table,skymags_table)

        if args.compute_skymags :
            skymag_rows=[]
            for count,night in enumerate(nights) :
                dirnames = sorted(glob.glob('{}/exposures/{}/*'.format(args.prod,night)))
                night_expids=[]
                for dirname in dirnames :
                    try :
                        expid=int(os.path.basename(dirname))
                        night_expids.append(expid)
                    except ValueError as e :
                        log.warning("ignore {}".format(dirname))
                if expids is not None :
                    night_expids = np.intersect1d(expids,night_expids)
                    if night_expids.size == 0 :
                        continue
                log.info("{} {}".format(night,night_expids))

                func2_args = []
                for expid in night_expids :
                    func2_args.append({'night':night,'expid':expid,'specprod_dir':args.prod})

                if args.nproc == 1 :
                    log.info(f"--nproc={args.nproc} so not multiprocessing")
                    for func2_arg in func2_args :
                        entry = func2(**func2_arg)
                        if entry is not None :
                            skymag_rows.append(entry)
                else :
                    log.info("Multiprocessing with {} procs".format(args.nproc))
                    with multiprocessing.Pool(args.nproc) as pool:
                        results  =  pool.map(_func2, func2_args)
                        for entry in results :
                            if entry is not None :
                                skymag_rows.append(entry)

            if len(skymag_rows) > 0:
                log.info("Adding computed sky magnitudes to table.")
                colnames = list(skymag_rows[0].keys())
                skymags_table = Table(rows=skymag_rows, names=colnames)
                add_skymags_columns(tsnr2_expid_table,skymags_table)
            else:
                log.warning("No sky magnitudes to append to table.")

        gfa_table = None
        gfa_nights = list()

        if args.gfa_proc_dir is not None:
            try:
                gfa_table = read_gfa_data(args.gfa_proc_dir)
            except PermissionError as e:
                log.error(e)
                log.error("could not read some files in {}".format(args.gfa_proc_dir))
                args.gfa_proc_dir = None
            #
            # It is not unusual for GFA processing to be delayed or otherwise to be
            # out of sync with the current night. Therefore, when adding GFA
            # data to the exposures, make a note of any nights where the GFA
            # data changed, so this can be propagated to the tiles file.
            #
            if gfa_table is not None:
                #
                # Perform a join of gfa_table & tsnr2_expid_table.
                #
                e2i = {e: i for i, e in enumerate(gfa_table["EXPID"])}
                jj = []
                ii = []
                for j, e in enumerate(tsnr2_expid_table["EXPID"]):
                    if e in e2i:
                        jj.append(j)
                        ii.append(e2i[e])
                tsnr2_nights = tsnr2_expid_table["NIGHT"][jj]
                #
                # For daily reductions, all of the GFA columns should already exist.
                #
                cols = ("TRANSPARENCY", "FWHM_ASEC", "FIBER_FRACFLUX", "FIBER_FRACFLUX_ELG",
                        "FIBER_FRACFLUX_BGS", "FIBERFAC", "FIBERFAC_ELG", "FIBERFAC_BGS",
                        "SKY_MAG_AB", "AIRMASS")
                for col in cols:
                    if col in gfa_table.dtype.names:
                        if col == "FWHM_ASEC":
                            col2 = "SEEING_GFA"
                        else:
                            col2 = col + "_GFA"
                        if col2 not in tsnr2_expid_table.colnames:
                            tsnr2_expid_table[col2] = np.zeros(len(tsnr2_expid_table), dtype=float)
                        gfa_values = gfa_table[col][ii]
                        expid_values = tsnr2_expid_table[col2][jj]
                        changed_values = (expid_values != gfa_values) & np.isfinite(gfa_values)
                        gfa_nan = ~np.isfinite(gfa_values)
                        if changed_values.any():
                            expid_values[changed_values] = gfa_values[changed_values]
                            gfa_nights += tsnr2_nights[changed_values].tolist()
                            tsnr2_expid_table[col2][jj] = expid_values
                        if gfa_nan.any():
                            good_values = gfa_values.copy()
                            good_values[gfa_nan] = 0
                            #
                            # If we are working with exposure tables with historical data, chances
                            # are this correction has already been applied.
                            #
                            gfa_nan_changed_values = (tsnr2_expid_table[col2][jj] != good_values)
                            if gfa_nan_changed_values.any():
                                log.warning("gfa_table column '%s' contains NaN or other non-finite values. Invalid values will be replaced with zero (0).", col)
                                tsnr2_expid_table[col2][jj] = good_values
                                gfa_nights += tsnr2_nights[gfa_nan_changed_values].tolist()

            if args.skymags is not None or args.compute_skymags:
                #
                # Assuming NaN are replaced with zero above, compute_efftime() will return zero for the affected rows.
                # That is, there are no divide-by-zero or other errors induced by replacing NaN with zero.
                #
                efftime_dark, efftime_bright, efftime_backup = compute_efftime(tsnr2_expid_table[jj])
                efftime_columns = {"EFFTIME_DARK_GFA": efftime_dark,
                                   "EFFTIME_BRIGHT_GFA": efftime_bright,
                                   "EFFTIME_BACKUP_GFA": efftime_backup}
                goaltype = tsnr2_expid_table["GOALTYPE"][jj]
                for col in efftime_columns:
                    if col not in tsnr2_expid_table.colnames:
                        tsnr2_expid_table[col] = np.zeros(len(tsnr2_expid_table), dtype=float)
                    gfa_values = efftime_columns[col]
                    expid_values = tsnr2_expid_table[col][jj]
                    changed_values = (expid_values != gfa_values)
                    if changed_values.any():
                        expid_values[changed_values] = gfa_values[changed_values]
                        gfa_nights += tsnr2_nights[changed_values].tolist()
                        tsnr2_expid_table[col][jj] = expid_values
                        if col == "EFFTIME_DARK_GFA":
                            # tsnr2_expid_table["EFFTIME_GFA"] = tsnr2_expid_table["EFFTIME_DARK_GFA"]  # default
                            efftime_gfa = expid_values
                        elif col == "EFFTIME_BRIGHT_GFA":
                            # tsnr2_expid_table["EFFTIME_GFA"][goaltype == "bright"] = tsnr2_expid_table["EFFTIME_BRIGHT_GFA"][goaltype == "bright"]
                            efftime_gfa[goaltype == "bright"] = expid_values[goaltype == "bright"]
                        else:
                            # tsnr2_expid_table["EFFTIME_GFA"][goaltype == "backup"] = tsnr2_expid_table["EFFTIME_BACKUP_GFA"][goaltype == "backup"]
                            efftime_gfa[goaltype == "backup"] = expid_values[goaltype == "backup"]
                        tsnr2_expid_table["EFFTIME_GFA"][jj] = efftime_gfa

        gfa_nights = np.unique(np.array(gfa_nights))
        #
        # End of loop on nights.
        #
        if len(tsnr2_frame_table) == 0:
            log.error("no valid exposures added")
            return 1

        if args.tile_completeness is not None:
            # renaming for historical code reasons
            exposure_table = tsnr2_expid_table

            # get list of tiles to look at
            selection = np.repeat(True, len(exposure_table))
            if args.expids is not None or args.nights is not None:
                #
                # args.expids and args.nights should have already been parsed above.
                # Is there a specific reason to reparse them here?
                #
                if args.expids is not None:
                    expids = parse_int_args(args.expids)
                    selection &= np.isin(exposure_table["EXPID"], expids)
                if args.nights is not None:
                    nights = parse_int_args(args.nights)
                    selection &= np.isin(exposure_table["NIGHT"], nights)
                if len(gfa_nights) > 0:
                    selection &= np.isin(exposure_table["NIGHT"], gfa_nights)
                # consider only the tiles observed in this selection of exposures
                tiles=np.unique(exposure_table["TILEID"][selection])
                # keep all exposures of those tiles
                selection = np.isin(exposure_table["TILEID"],tiles)

            new_tile_table = compute_tile_completeness_table(exposure_table[selection],args.prod,auxiliary_table_filenames=args.aux)
            if os.path.isfile(args.tile_completeness) :
                previous_table = Table.read(args.tile_completeness)
                new_tile_table = merge_tile_completeness_table(previous_table, new_tile_table)

            head, foo = os.path.splitext(args.tile_completeness)
            write_formats = {'.fits': 'fits', '.csv': 'ascii.csv'}
            for ext in write_formats:
                new_tile_table.write(head + ext, format=write_formats[ext], overwrite=True)
                log.info("Wrote %s.", head + ext)

            # compute_tile_completeness_table may have found new GOALTIME info
            # from auxiliary_table_filename, so update exposures if needed
            for tileid, goaltime in new_tile_table['TILEID', 'GOALTIME']:
                ii = tsnr2_expid_table['TILEID'] == tileid
                if goaltime>0.0 and np.any(tsnr2_expid_table['GOALTIME'][ii] == 0.0):
                    tsnr2_expid_table['GOALTIME'][ii] = goaltime

                jj = tsnr2_frame_table['TILEID'] == tileid
                if goaltime>0.0 and np.any(tsnr2_frame_table['GOALTIME'][jj] == 0.0):
                    tsnr2_frame_table['GOALTIME'][jj] = goaltime

        output_csv_filename = os.path.splitext(args.outfile)[0] + '.csv'
        write_summary_tables(tsnr2_expid_table,tsnr2_frame_table,output_fits_filename=args.outfile,
                             output_csv_filename=output_csv_filename)

        # remove temporary file if successful
        if os.path.isfile(args.outfile) :
            tmpfilename=args.outfile.replace(".fits","_tmp.fits")
            if os.path.isfile(tmpfilename) :
                log.info("rm temporary file {}".format(tmpfilename))
                os.unlink(tmpfilename)
            tmpfilename=args.outfile.replace(".fits","_tmp.csv")
            if os.path.isfile(tmpfilename) :
                log.info("rm temporary file {}".format(tmpfilename))
                os.unlink(tmpfilename)
        else :
            log.error("no valid exposures added")

    return 0


if __name__ == '__main__':
    sys.exit(main())
