From 62a928b67b06a8eea296f2a84e9f5e24e488b8f7 Mon Sep 17 00:00:00 2001 From: Kevin Ludwick Date: Wed, 20 Nov 2024 23:26:42 -0600 Subject: [PATCH 01/26] photon counting, fixed an err estimation from converting with k gain, fixed a minor error in combine.py almost finished with unit tests for photon counting --- corgidrp/combine.py | 2 +- corgidrp/darks.py | 2 +- corgidrp/l2a_to_l2b.py | 10 +- corgidrp/mocks.py | 81 ++++ corgidrp/photon_counting.py | 368 +++++++++++++++++++ corgidrp/recipe_templates/l1_to_l2b_pc.json | 55 +++ corgidrp/recipe_templates/l2a_to_l2b_pc.json | 29 ++ corgidrp/walker.py | 4 +- tests/test_kgain.py | 3 + tests/test_photon_counting.py | 66 ++++ 10 files changed, 614 insertions(+), 6 deletions(-) create mode 100644 corgidrp/photon_counting.py create mode 100644 corgidrp/recipe_templates/l1_to_l2b_pc.json create mode 100644 corgidrp/recipe_templates/l2a_to_l2b_pc.json create mode 100644 tests/test_photon_counting.py diff --git a/corgidrp/combine.py b/corgidrp/combine.py index 6bc38b9e..7a7c1abb 100644 --- a/corgidrp/combine.py +++ b/corgidrp/combine.py @@ -91,7 +91,7 @@ def combine_subexposures(input_dataset, num_frames_per_group=None, collapse="mea pri_hdr = input_dataset[num_frames_per_group*i].pri_hdr.copy() ext_hdr = input_dataset[num_frames_per_group*i].ext_hdr.copy() err_hdr = input_dataset[num_frames_per_group*i].err_hdr.copy() - dq_hdr = input_dataset[num_frames_per_group*i].err_hdr.copy() + dq_hdr = input_dataset[num_frames_per_group*i].dq_hdr.copy() hdulist = input_dataset[num_frames_per_group*i].hdu_list.copy() new_image = data.Image(data_collapse, pri_hdr=pri_hdr, ext_hdr=ext_hdr, err=err_collapse, dq=dq_collapse, err_hdr=err_hdr, dq_hdr=dq_hdr, input_hdulist=hdulist) diff --git a/corgidrp/darks.py b/corgidrp/darks.py index 93c53f82..bb490d60 100644 --- a/corgidrp/darks.py +++ b/corgidrp/darks.py @@ -111,7 +111,7 @@ def mean_combine(image_list, bpmap_list, err=False): # Divide sum_im by map_im only where map_im is not equal to 0 (i.e., # not masked). # Where map_im is equal to 0, set combined_im to zero - comb_image = np.divide(sum_im, map_im, out=np.zeros_like(sum_im), + comb_image = np.divide(sum_im, np.sqrt(map_im), out=np.zeros_like(sum_im), where=map_im != 0) # Mask any value that was never mapped (aka masked in every frame) diff --git a/corgidrp/l2a_to_l2b.py b/corgidrp/l2a_to_l2b.py index 78bfd980..cefb4e81 100644 --- a/corgidrp/l2a_to_l2b.py +++ b/corgidrp/l2a_to_l2b.py @@ -251,13 +251,17 @@ def convert_to_electrons(input_dataset, k_gain): kgain_error = kgain_dataset.all_err kgain = k_gain.value #extract from caldb + kgainerr = k_gain.error[0] kgain_cube *= kgain - kgain_error *= kgain - + kgain_err_up = kgain_error*(kgain+kgainerr) + kgain_err_low = kgain_error*(kgain-kgainerr) + kgain_error = np.max([kgain_err_up - kgain, kgain - kgain_err_low], axis=0) + history_msg = "data converted to detected EM electrons by kgain {0}".format(str(kgain)) # update the output dataset with this converted data and update the history - kgain_dataset.update_after_processing_step(history_msg, new_all_data=kgain_cube, new_all_err=kgain_error, header_entries = {"BUNIT":"detected EM electrons", "KGAIN":kgain}) + kgain_dataset.update_after_processing_step(history_msg, new_all_data=kgain_cube, new_all_err=kgain_error, header_entries = {"BUNIT":"detected EM electrons", "KGAIN":kgain, + "KGAIN_ER": k_gain.error[0], "RN":k_gain.ext_hdr['RN'], "RN_ERR":k_gain.ext_hdr["RN_ERR"]}) return kgain_dataset def em_gain_division(input_dataset): diff --git a/corgidrp/mocks.py b/corgidrp/mocks.py index cf278aae..c17cbd18 100644 --- a/corgidrp/mocks.py +++ b/corgidrp/mocks.py @@ -1,6 +1,7 @@ import os from pathlib import Path import numpy as np +import warnings import scipy.ndimage import pandas as pd import astropy.io.fits as fits @@ -1857,3 +1858,83 @@ def add_defect(sch_imgs, prob, ori, temp): hdul.writeto(str(filename)[:-4]+'_'+str(mult_counter)+'.fits', overwrite = True) else: hdul.writeto(filename, overwrite = True) + +def create_photon_countable_frames(Nbrights=30, Ndarks=40, EMgain=5000, kgain=7, exptime=0.1, cosmic_rate=0): + '''This creates mock Dataset containing frames with large gain and short exposure time, illuminated and dark frames. + Used for unit tests for photon counting. + + Args: + Nbrights (int): number of illuminated frames to simulate + Ndarks (int): number of dark frames to simulate + EMgain (float): EM gain + kgain (float): k gain (e-/DN) + exptime (float): exposure time (in s) + cosmic_rate (float) : simulated cosmic rays incidence, hits/cm^2/s + + Returns: + dataset (corgidrp.data.Dataset): Dataset containing both the illuminated and dark frames + ill_mean (float): mean electron count value simulated in the illuminated frames + dark_mean (float): mean electron count value simulated in the dark frames + ''' + pix_row = 1024 #number of rows and number of columns + fluxmap = np.ones((pix_row,pix_row)) #photon flux map, photons/s + + emccd = EMCCDDetect( + em_gain=EMgain, + full_well_image=60000., # e- + full_well_serial=100000., # e- + dark_current=8.33e-4, # e-/pix/s + cic=0.01, # e-/pix/frame + read_noise=100., # e-/pix/frame + bias=0, # e-; 0 for simplicity + qe=0.9, # quantum efficiency, e-/photon + cr_rate=cosmic_rate, # cosmic rays incidence, hits/cm^2/s + pixel_pitch=13e-6, # m + eperdn=kgain, + nbits=64, # number of ADU bits + numel_gain_register=604 #number of gain register elements + ) + + thresh = emccd.em_gain/10 # threshold + + if np.average(exptime*fluxmap) > 0.1: + warnings.warn('average # of photons/pixel is > 0.1. Decrease frame ' + 'time to get lower average # of photons/pixel.') + + if emccd.read_noise <=0: + warnings.warn('read noise should be greater than 0 for effective ' + 'photon counting') + if thresh < 4*emccd.read_noise: + warnings.warn('thresh should be at least 4 or 5 times read_noise for ' + 'accurate photon counting') + + avg_ph_flux = np.mean(fluxmap) + # theoretical electron flux for brights + ill_mean = avg_ph_flux*emccd.qe*exptime + emccd.dark_current*exptime + emccd.cic + # theoretical electron flux for darks + dark_mean = emccd.dark_current*exptime + emccd.cic + + frame_e_list = [] + frame_e_dark_list = [] + prihdr, exthdr = create_default_headers() + for i in range(Nbrights): + # Simulate bright + frame_dn = emccd.sim_full_frame(fluxmap, exptime) + frame = data.Image(frame_dn, pri_hdr=prihdr, ext_hdr=exthdr) + frame.ext_hdr['CMDGAIN'] = EMgain + frame.ext_hdr['EXPTIME'] = exptime + frame.pri_hdr["VISTYPE"] = "TDEMO" + frame_e_list.append(frame) + + for i in range(Ndarks): + # Simulate dark + frame_dn_dark = emccd.sim_full_frame(np.zeros_like(fluxmap), exptime) + frame_dark = data.Image(frame_dn_dark, pri_hdr=prihdr.copy(), ext_hdr=exthdr.copy()) + frame_dark.ext_hdr['CMDGAIN'] = EMgain + frame_dark.ext_hdr['EXPTIME'] = exptime + frame_dark.pri_hdr["VISTYPE"] = "DARK" + frame_e_dark_list.append(frame_dark) + + dataset = data.Dataset(frame_e_list+frame_e_dark_list) + + return dataset, ill_mean, dark_mean \ No newline at end of file diff --git a/corgidrp/photon_counting.py b/corgidrp/photon_counting.py new file mode 100644 index 00000000..6629f211 --- /dev/null +++ b/corgidrp/photon_counting.py @@ -0,0 +1,368 @@ +"""Photon count a stack of analog images and return a mean expected electron +count array with photometric corrections. + +""" +import warnings +import numpy as np +import corgidrp.data as data + +def varL23(g, L, T, N): + '''Expected variance after photometric correction. + See https://doi.org/10.1117/1.JATIS.9.4.048006 for details. + + Parameters + ---------- + g (scalar): EM gain + L (2-D array): mean expected number of electrons + T (scalar): threshold + N (2-D array): number of frames + + Returns + ------- + variance from photon counting and the photometric correction + + ''' + Const = 6/(6 + L*(6 + L*(3 + L))) + eThresh = (Const*(np.e**(-T/g)*L*(2*g**2*(6 + L*(3 + L)) + + 2*g*L*(3 + L)*T + L**2*T**2))/(12*g**2)) + std_dev = np.sqrt(N * eThresh * (1-eThresh)) + + return (std_dev)**2*(((np.e**((T/g)))/N) + + 2*((np.e**((2*T)/g)*(g - T))/(2*g*N**2))*(N*eThresh) + + 3*(((np.e**(((3*T)/g)))*(4*g**2 - 8*g*T + 5*T**2))/( + 12*g**2*N**3))*(N*eThresh)**2)**2 + +class PhotonCountException(Exception): + """Exception class for photon_count module.""" + + +def photon_count(e_image, thresh): + """Convert an analog image into a photon-counted image. + + Parameters + ---------- + e_image : array_like, float + Analog image (e-). + thresh : float + Photon counting threshold (e-). Values > thresh will be assigned a 1, + values <= thresh will be assigned a 0. + + Returns + ------- + pc_image : array_like, float + Output digital image in units of photons. + + B Nemati and S Miller - UAH - 06-Aug-2018 + + """ + # Check if input is an array/castable to one + e_image = np.array(e_image).astype(float) + if len(e_image.shape) == 0: + raise PhotonCountException('e_image must have length > 0') + + pc_image = np.zeros(e_image.shape, dtype=int) + pc_image[e_image > thresh] = 1 + + return pc_image + + +def get_pc_mean(input_dataset, detector_params, niter=2): + """Take a stack of images, illuminated and dark frames of the same exposure + time, k gain, read noise, and EM gain, and return the mean expected value per + pixel. The frames are taken in photon-counting mode, which means high EM + gain and very low exposure time. The number of darks can be different from + the number of illuminated frames. + + The frames in each stack should be SCI image (1024x1024) frames that have + had some of the L2b steps applied: + + - have had their bias subtracted + - have had masks made for cosmic rays + - have been corrected for nonlinearity (These first 3 steps make the frames L2a.) + - have been frame-selected (to weed out bad frames) + - have been converted from DN to e- + - have been desmeared if desmearing is appropriate. Under normal + circumstances, these frames should not be desmeared. The only time desmearing + would be useful is in the unexpected case that, for example, + dark current is so high that it stands far above other noise that is + not smeared upon readout, such as clock-induced charge, + fixed-pattern noise, and read noise. For such low-exposure frames, though, + this really should not be a concern. + + This algorithm will photon count each frame in each of the 2 sets (illuminated and dark) individually, + then co-add the photon counted frames. The co-added frame is then averaged + and corrected for thresholding and coincidence loss, returning the mean + expected array in units of photoelectrons. The threshold is + determined by "T_factor" stored in DetectorParams. The dark mean expected array + is then subtracted from the illuminated mean expected array. + + Parameters + ---------- + input_dataset (corgidrp.data.Dataset): + This is an instance of corgidrp.data.Dataset containing the + photon-counted illuminated frames as well as the + photon-counted dark frames, and all the frames must have the same + exposure time, EM gain, k gain, and read noise. + detector_params: (corgidrp.data.DetectorParams) + A calibration file storing detector calibration values. + niter : int, optional + Number of Newton's method iterations (used for photometric + correction). Defaults to 2. + + Returns + ------- + mean_expected : array_like + Mean expected photoelectron array (dark-subtracted) + + References + ---------- + [1] https://www.spiedigitallibrary.org/conference-proceedings-of-spie/11443/114435F/Photon-counting-and-precision-photometry-for-the-Roman-Space-Telescope/10.1117/12.2575983.full + [2] https://doi.org/10.1117/1.JATIS.9.4.048006 + + B Nemati, S Miller - UAH - 13-Dec-2020 + Kevin Ludwick - UAH - 2023 + + """ + test_dataset, _ = input_dataset.copy().split_dataset(exthdr_keywords=['EXPTIME', 'CMDGAIN', 'KGAIN', 'RN']) + if len(test_dataset) > 1: + raise PhotonCountException('All frames must have the same exposure time, ' + 'commanded EM gain, and k gain.') + datasets, vals = test_dataset[0].split_dataset(prihdr_keywords=['VISTYPE']) + if len(vals) != 2: + raise PhotonCountException('There should only be 2 datasets, and one should have \'VISTYPE\' = \'DARK\'.') + dark_dataset = datasets[vals.index('DARK')] + ill_dataset = datasets[1-vals.index('DARK')] + + pc_means = [] + errs = [] + dqs = [] + for dataset in [ill_dataset, dark_dataset]: + # getting number of read noise standard deviations at which to set the + # photon-counting threshold + T_factor = detector_params.params['T_factor'] + # now get threshold to use for photon-counting + thresh = T_factor*dataset.frames[0].ext_hdr['RN'] + if thresh < 0: + raise PhotonCountException('thresh must be nonnegative') + if not isinstance(niter, (int, np.integer)) or niter < 1: + raise PhotonCountException('niter must be an integer greater than ' + '0') + try: # if EM gain measured directly from frame TODO change hdr name if necessary + em_gain = dataset.frames[0].ext_hdr['EMGAIN_M'] + except: + try: # use applied EM gain if available + em_gain = dataset.frames[0].ext_hdr['EMGAIN_A'] + except: # use commanded gain otherwise + em_gain = dataset.frames[0].ext_hdr['CMDGAIN'] + if thresh >= em_gain: + warnings.warn('thresh should be less than em_gain for effective ' + 'photon counting') + + # Photon count stack of frames + frames_pc = photon_count(dataset.all_data, thresh) + # frames_pc = np.array([photon_count(frame, thresh) for frame in dataset.all_data]) XXX + bool_map = dataset.all_dq.astype(bool).astype(float) + bool_map[bool_map > 0] = np.nan + bool_map[bool_map == 0] = 1 + nframes = np.nansum(bool_map, axis=0) + frames_pc_masked = frames_pc * bool_map + # Co-add frames + frame_pc_coadded = np.nansum(frames_pc_masked, axis=0) + + # Correct for thresholding and coincidence loss; any pixel masked all the + # way through the stack will give NaN + mean_expected = corr_photon_count(frame_pc_coadded, nframes, thresh, + em_gain, niter) + + ##### error calculation + # combine all err frames to get standard error, which is due to inherent error in + # in the data before any photon counting is done + err_sum = np.sum(dataset.all_err**2, axis=0) + bool_map_sum = np.nansum(bool_map, axis=0) + comb_err = np.divide(np.sqrt(err_sum), np.sqrt(bool_map_sum), out=np.zeros_like(err_sum), + where=bool_map_sum != 0) + # expected error from photon counting + pc_variance = varL23(em_gain,mean_expected,thresh,nframes) + # now combine this error with the statistical error coming from the process of photon counting + total_err = np.sqrt(comb_err**2 + pc_variance) + total_dq = (bool_map_sum == 0).astype(int) # 1 where flagged, 0 otherwise + pc_means.append(mean_expected) + errs.append(total_err) + dqs.append(total_dq) + + # now subtract the dark PC mean + combined_pc_mean = pc_means[0] - pc_means[1] + combined_err = np.sqrt(errs[0]**2 + errs[1]**2) + combined_dq = np.bitwise_or(dqs[0], dqs[1]) + pri_hdr = ill_dataset[0].pri_hdr.copy() + ext_hdr = ill_dataset[0].ext_hdr.copy() + err_hdr = ill_dataset[0].err_hdr.copy() + dq_hdr = ill_dataset[0].dq_hdr.copy() + hdulist = ill_dataset[0].hdu_list.copy() + new_image = data.Image(combined_pc_mean, pri_hdr=pri_hdr, ext_hdr=ext_hdr, err=combined_err, dq=combined_dq, err_hdr=err_hdr, + dq_hdr=dq_hdr, input_hdulist=hdulist) + + new_image.filename = ill_dataset[0].filename.split('.')[0]+'_pc.fits' + new_image._record_parent_filenames(input_dataset) + pc_dataset = data.Dataset([new_image]) + pc_dataset.update_after_processing_step("Photon-counted {0} frames using T_factor = {1} and niter={2}".format(len(input_dataset), T_factor, niter)) + + return pc_dataset + +def corr_photon_count(nobs, nfr, t, g, niter=2): + """Correct photon counted images. + + Parameters + ---------- + nobs : array_like + Number of observations (Co-added photon counted frame). + nfr : int + Number of coadded frames. + t : float + Photon counting threshold. + g : float + EM gain. + niter : int, optional + Number of Newton's method iterations. Defaults to 2. + + Returns + ------- + lam : array_like + Mean expeted detected electrons per pixel (lambda). + + """ + # Get an approximate value of lambda for the first guess of the Newton fit + lam0 = calc_lam_approx(nobs, nfr, t, g) + + # Use Newton's method to converge at a value for lambda + lam = lam_newton_fit(nobs, nfr, t, g, lam0, niter) + + return lam + + +def calc_lam_approx(nobs, nfr, t, g): + """Approximate lambda calculation. + + This will calculate the first order approximation of lambda, and for values + that are out of bounds (e.g. from statistical fluctuations) it will revert + to the zeroth order. + + Parameters + ---------- + nobs : array_like + Number of observations (Co-added photon counted frame). + nfr : int + Number of coadded frames. + t : float + Photon counting threshold. + g : float + EM gain used when taking images. + + Returns + ------- + array_like + Mean expeted (lambda). + + """ + # First step of equation (before taking log) + init = 1 - (nobs/nfr) * np.exp(t/g) + # Mask out all values less than or equal to 0 + lam_m = np.zeros_like(init).astype(bool) + lam_m[init > 0] = True + + # Use the first order approximation on all values greater than zero + lam1 = np.zeros_like(init) + lam1[lam_m] = -np.log(init[lam_m]) + + # For any values less than zero, revert to the zeroth order approximation + lam0 = nobs / nfr + lam1[~lam_m] = lam0[~lam_m] + + return lam1 + + +def lam_newton_fit(nobs, nfr, t, g, lam0, niter): + """Newton fit for finding lambda. + + Parameters + ---------- + nobs : array_like + Number of observations (Co-added photon counted frame). + nfr : int + Number of coadded frames. + t : float + Photon counting threshold. + g : float + EM gain used when taking images. + lam0 : array_like + Initial guess for lambda. + niter : int + Number of Newton's fit iterations to take. + + Returns + ------- + lam : array_like + Mean expeted (lambda). + + """ + # Mask out zero values to avoid divide by zero + lam_est_m = np.ma.masked_where(lam0 == 0, lam0) + nobs_m = np.ma.masked_where(nobs == 0, nobs) + + # Iterate Newton's method + for i in range(niter): + func = _calc_func(nobs_m, nfr, t, g, lam_est_m) + dfunc = _calc_dfunc(nfr, t, g, lam_est_m) + lam_est_m -= func / dfunc + + if lam_est_m.min() < 0: + raise PhotonCountException('negative number of photon counts; ' + 'try decreasing the frametime') + + # Fill zero values back in + lam = lam_est_m.filled(0) + + return lam + + +def _calc_func(nobs, nfr, t, g, lam): + """Objective function for lambda.""" + e_thresh = ( + np.exp(-t/g) + * ( + t**2 * lam**2 + + 2*g * t * lam * (3 + lam) + + 2*g**2 * (6 + 3*lam + lam**2) + ) + / (2*g**2 * (6 + 3*lam + lam**2)) + ) + + e_coinloss = (1 - np.exp(-lam)) / lam + + #if (lam * nfr * e_thresh * e_coinloss).any() > nobs.any(): + # warnings.warn('Input photon flux is too high; decrease frametime') + # This warning isn't necessary; could have a negative func but still + # close enough to 0 for Newton's method + func = lam * nfr * e_thresh * e_coinloss - nobs + + return func + + +def _calc_dfunc(nfr, t, g, lam): + """Derivative wrt lambda of objective function.""" + dfunc = ( + (np.exp(-t/g - lam) * nfr) + / (2 * g**2 * (6 + 3*lam + lam**2)**2) + * ( + 2*g**2 * (6 + 3*lam + lam**2)**2 + + t**2 * lam * ( + -12 + 3*lam + 3*lam**2 + lam**3 + 3*np.exp(lam)*(4 + lam) + ) + + 2*g*t * ( + -18 + 6*lam + 15*lam**2 + 6*lam**3 + lam**4 + + 6*np.exp(lam)*(3 + 2*lam) + ) + ) + ) + + return dfunc \ No newline at end of file diff --git a/corgidrp/recipe_templates/l1_to_l2b_pc.json b/corgidrp/recipe_templates/l1_to_l2b_pc.json new file mode 100644 index 00000000..b35144ba --- /dev/null +++ b/corgidrp/recipe_templates/l1_to_l2b_pc.json @@ -0,0 +1,55 @@ +{ + "name" : "l1_to_l2b_pc", + "template" : true, + "drpconfig" : { + "track_individual_errors" : false + }, + "inputs" : [], + "outputdir" : "", + "steps" : [ + { + "name" : "prescan_biassub", + "calibs" : { + "DetectorNoiseMaps" : "AUTOMATIC" + }, + "keywords" : { + "return_full_frame" : false + } + }, + { + "name" : "detect_cosmic_rays", + "calibs" : { + "DetectorParams" : "AUTOMATIC", + "KGain" : "AUTOMATIC" + } + }, + { + "name" : "correct_nonlinearity", + "calibs" : { + "NonLinearityCalibration" : "AUTOMATIC" + } + }, + { + "name" : "update_to_l2a" + }, + + { + "name" : "frame_select" + }, + { + "name" : "convert_to_electrons", + "calibs" : { + "KGain" : "AUTOMATIC" + } + }, + { + "name": "get_pc_mean", + "calibs" : { + "DetectorParams" : "AUTOMATIC" + } + }, + { + "name" : "save" + } + ] +} \ No newline at end of file diff --git a/corgidrp/recipe_templates/l2a_to_l2b_pc.json b/corgidrp/recipe_templates/l2a_to_l2b_pc.json new file mode 100644 index 00000000..2f8b0236 --- /dev/null +++ b/corgidrp/recipe_templates/l2a_to_l2b_pc.json @@ -0,0 +1,29 @@ +{ + "name" : "l1_to_l2b_pc", + "template" : true, + "drpconfig" : { + "track_individual_errors" : false + }, + "inputs" : [], + "outputdir" : "", + "steps" : [ + { + "name" : "frame_select" + }, + { + "name" : "convert_to_electrons", + "calibs" : { + "KGain" : "AUTOMATIC" + } + }, + { + "name": "get_pc_mean", + "calibs" : { + "DetectorParams" : "AUTOMATIC" + } + }, + { + "name" : "save" + } + ] +} \ No newline at end of file diff --git a/corgidrp/walker.py b/corgidrp/walker.py index 25f9b441..d7a1774d 100644 --- a/corgidrp/walker.py +++ b/corgidrp/walker.py @@ -11,6 +11,7 @@ import corgidrp.caldb as caldb import corgidrp.l1_to_l2a import corgidrp.l2a_to_l2b +import corgidrp.photon_counting import corgidrp.pump_trap_calibration import corgidrp.calibrate_nonlin import corgidrp.detector @@ -41,7 +42,8 @@ "create_onsky_flatfield" : corgidrp.detector.create_onsky_flatfield, "combine_subexposures" : corgidrp.combine.combine_subexposures, "build_trad_dark" : corgidrp.darks.build_trad_dark, - "sort_pupilimg_frames" : corgidrp.sorting.sort_pupilimg_frames + "sort_pupilimg_frames" : corgidrp.sorting.sort_pupilimg_frames, + "get_pc_mean" : corgidrp.photon_counting.get_pc_mean } recipe_dir = os.path.join(os.path.dirname(__file__), "recipe_templates") diff --git a/tests/test_kgain.py b/tests/test_kgain.py index 990ada5a..6313a58b 100644 --- a/tests/test_kgain.py +++ b/tests/test_kgain.py @@ -84,6 +84,9 @@ def test_kgain(): assert gain_dataset[0].ext_hdr["BUNIT"] == "detected EM electrons" assert gain_dataset[0].err_hdr["BUNIT"] == "detected EM electrons" assert gain_dataset[0].ext_hdr["KGAIN"] == k_gain + assert gain_dataset[0].ext_hdr["KGAIN_ER"] == kgain.error + assert gain_dataset[0].ext_hdr["RN"] > 0 + assert gain_dataset[0].ext_hdr["RN_ERR"] > 0 assert("converted" in str(gain_dataset[0].ext_hdr["HISTORY"])) if __name__ == "__main__": diff --git a/tests/test_photon_counting.py b/tests/test_photon_counting.py new file mode 100644 index 00000000..10e974f7 --- /dev/null +++ b/tests/test_photon_counting.py @@ -0,0 +1,66 @@ +import pytest +import os +import astropy.time as time +import corgidrp.mocks as mocks +import corgidrp +from astropy.io import fits +import numpy as np +import corgidrp.walker as walker +import corgidrp.caldb as caldb +import corgidrp.data as data + +dataset, ill_mean, dark_mean = mocks.create_photon_countable_frames(Nbrights=20, Ndarks=30, cosmic_rate=0) +thisfile_dir = os.path.dirname(__file__) # this file's folder +outputdir = os.path.join(thisfile_dir, 'simdata', 'pc_test_data') +if not os.path.exists(outputdir): + os.mkdir(outputdir) +# empty out directory of any previous files +for f in os.listdir(outputdir): + os.remove(os.path.join(outputdir,f)) +dataset.save(outputdir, ['pc_frame_{0}.fits'.format(i) for i in range(len(dataset))]) +l1_data_filelist = [] +for f in os.listdir(outputdir): + l1_data_filelist.append(os.path.join(outputdir, f)) + +def test_expected_results(): + '''Results are as expected theoretically. Also runs raw frames through pre-processing pipeline.''' + + this_caldb = caldb.CalDB() # connection to cal DB + # remove other KGain calibrations that may exist in case they don't have the added header keywords + for i in range(len(this_caldb._db['Type'])): + if this_caldb._db['Type'][i] == 'KGain': + this_caldb._db = this_caldb._db.drop(i) + this_caldb.save() + # KGain + kgain_val = 7 # default value used in mocks.create_photon_countable_frames() + pri_hdr, ext_hdr = mocks.create_default_headers() + ext_hdr["DRPCTIME"] = time.Time.now().isot + ext_hdr['DRPVERSN'] = corgidrp.__version__ + mock_input_dataset = data.Dataset(l1_data_filelist) + kgain = data.KGain(np.array([[kgain_val]]), pri_hdr=pri_hdr, ext_hdr=ext_hdr, + input_dataset=mock_input_dataset) + # add in keywords that didn't make it into mock_kgain.fits, using values used in mocks.create_photon_countable_frames() + kgain.ext_hdr['RN'] = 100 + kgain.ext_hdr['RN_ERR'] = 0 + kgain.save(filedir=outputdir, filename="mock_kgain.fits") + this_caldb.create_entry(kgain) + walker.walk_corgidrp(l1_data_filelist, '', outputdir, template="l1_to_l2b_pc.json") + # get photon-counted frame + for f in os.listdir(outputdir): + if f.endswith('_pc.fits'): + pc_filename = f + pc_file = os.path.join(outputdir, pc_filename) + pc_frame = fits.getdata(pc_file) + pc_ext_hdr = fits.getheader(pc_file, 1) + assert np.isclose(pc_frame.mean(), ill_mean - dark_mean) + assert 'niter=2' in pc_ext_hdr["HISTORY"][-1] + assert 'T_factor=5' in pc_ext_hdr["HISTORY"][-1] + + # empty out directory now + for f in os.listdir(outputdir): + os.remove(os.path.join(outputdir,f)) + + #XXX negative value test, others from PhotonCount? + #XXX change pc error: should come from doing pc on upper and lower bound of pixel values; same goes for trad and synthesized dark calibration +if __name__ == '__main__': + test_expected_results() \ No newline at end of file From f181ccea55f81cb83f0d20e2d3ff4f42bf5ef9e6 Mon Sep 17 00:00:00 2001 From: Kevin Ludwick Date: Fri, 22 Nov 2024 20:59:26 -0600 Subject: [PATCH 02/26] unit tests done, a few other small improvements --- corgidrp/darks.py | 16 +-- corgidrp/data.py | 2 + corgidrp/l2a_to_l2b.py | 7 +- corgidrp/mocks.py | 19 ++- corgidrp/photon_counting.py | 116 ++++++++++--------- corgidrp/recipe_templates/l2a_to_l2b_pc.json | 2 +- tests/test_calibrate_darks_lsq.py | 2 +- tests/test_kgain.py | 4 +- tests/test_mean_combine.py | 6 +- tests/test_photon_counting.py | 77 ++++++++---- 10 files changed, 152 insertions(+), 99 deletions(-) diff --git a/corgidrp/darks.py b/corgidrp/darks.py index bb490d60..dfb33ff4 100644 --- a/corgidrp/darks.py +++ b/corgidrp/darks.py @@ -105,14 +105,14 @@ def mean_combine(image_list, bpmap_list, err=False): sum_im += masked map_im += (im_m.mask == False).astype(int) - if err: # sqrt of sum of sigma**2 terms - sum_im = np.sqrt(sum_im) - # Divide sum_im by map_im only where map_im is not equal to 0 (i.e., # not masked). # Where map_im is equal to 0, set combined_im to zero - comb_image = np.divide(sum_im, np.sqrt(map_im), out=np.zeros_like(sum_im), + comb_image = np.divide(sum_im, map_im, out=np.zeros_like(sum_im), where=map_im != 0) + + if err: # (sqrt of sum of sigma**2 terms)/sqrt(n) + comb_image = np.sqrt(comb_image) # Mask any value that was never mapped (aka masked in every frame) comb_bpmap = (map_im == 0).astype(int) @@ -329,7 +329,7 @@ def calibrate_darks_lsq(dataset, detector_params, detector_regions=None): output Dark's dq after assigning these pixels a flag value of 256. They should have large err values. The pixels that are masked for EVERY frame in all sub-stacks - but 4 (or less) are assigned a flag value of + but 3 (or less) are assigned a flag value of 1, which falls under the category of "Bad pixel - unspecified reason". These pixels would have no reliability for dark subtraction. @@ -415,7 +415,7 @@ def calibrate_darks_lsq(dataset, detector_params, detector_regions=None): output Dark's dq after assigning these pixels a flag value of 256. They should have large err values. The pixels that are masked for EVERY frame in all sub-stacks - but 4 (or less) are assigned a flag value of + but 3 (or less) are assigned a flag value of 1, which falls under the category of "Bad pixel - unspecified reason". These pixels would have no reliability for dark subtraction. FPN_std_map : array-like (full frame) @@ -522,7 +522,7 @@ def calibrate_darks_lsq(dataset, detector_params, detector_regions=None): # flag value of 256; unreliable pixels, large err output_dq = (unreliable_pix_map >= len(datasets)-3).astype(int)*256 # flag value of 1 for those that are masked all the way through for all - # but 4 (or less) stacks; this overwrites the flag value of 256 that was assigned to + # but 3 (or less) stacks; this overwrites the flag value of 256 that was assigned to # these pixels in previous line unfittable_ind = np.where(unfittable_pix_map >= len(datasets)-3) output_dq[unfittable_ind] = 1 @@ -577,7 +577,7 @@ def calibrate_darks_lsq(dataset, detector_params, detector_regions=None): # input data error comes from .err arrays; could use this for error bars # in input data for weighted least squares, but we'll just simply get the # std error and add it in quadrature to least squares fit standard dev - stacks_err = np.sqrt(np.sum(mean_err_stack**2, axis=0))/len(mean_err_stack) + stacks_err = np.sqrt(np.sum(mean_err_stack**2, axis=0)/len(mean_err_stack)) # matrix to be used for least squares and covariance matrix X = np.array([np.ones([len(EMgain_arr)]), EMgain_arr, EMgain_arr*exptime_arr]).T diff --git a/corgidrp/data.py b/corgidrp/data.py index 32eff9ef..d18aa99e 100644 --- a/corgidrp/data.py +++ b/corgidrp/data.py @@ -104,6 +104,8 @@ def update_after_processing_step(self, history_entry, new_all_data=None, new_all if new_all_err.shape[-2:] != self.all_err.shape[-2:] or new_all_err.shape[0] != self.all_err.shape[0]: raise ValueError("The shape of new_all_err is {0}, whereas we are expecting {1}".format(new_all_err.shape, self.all_err.shape)) self.all_err = new_all_err + for i in range(len(self.frames)): + self.frames[i].err = self.all_err[i] if new_all_dq is not None: if new_all_dq.shape != self.all_dq.shape: raise ValueError("The shape of new_all_dq is {0}, whereas we are expecting {1}".format(new_all_dq.shape, self.all_dq.shape)) diff --git a/corgidrp/l2a_to_l2b.py b/corgidrp/l2a_to_l2b.py index cefb4e81..a3db3bef 100644 --- a/corgidrp/l2a_to_l2b.py +++ b/corgidrp/l2a_to_l2b.py @@ -248,14 +248,13 @@ def convert_to_electrons(input_dataset, k_gain): # you should make a copy the dataset to start kgain_dataset = input_dataset.copy() kgain_cube = kgain_dataset.all_data - kgain_error = kgain_dataset.all_err kgain = k_gain.value #extract from caldb kgainerr = k_gain.error[0] + # x = c*kgain, where c (counts beforehand) and kgain both have error, so do propogation of error due to the product of 2 independent sources + #[:,0:1,:,:] to get same dimensions as input_dataset.all_err + kgain_error = (np.abs(kgain_cube*kgain)*np.sqrt((kgain_dataset.all_err/kgain_cube)**2 + (kgainerr/kgain)**2))[:,0:1,:,:] kgain_cube *= kgain - kgain_err_up = kgain_error*(kgain+kgainerr) - kgain_err_low = kgain_error*(kgain-kgainerr) - kgain_error = np.max([kgain_err_up - kgain, kgain - kgain_err_low], axis=0) history_msg = "data converted to detected EM electrons by kgain {0}".format(str(kgain)) diff --git a/corgidrp/mocks.py b/corgidrp/mocks.py index c17cbd18..e485183c 100644 --- a/corgidrp/mocks.py +++ b/corgidrp/mocks.py @@ -167,7 +167,7 @@ def create_synthesized_master_dark_calib(detector_areas): # image area, including "shielded" rows and cols: imrows, imcols, imr0c0 = imaging_area_geom('SCI', detector_areas) prerows, precols, prer0c0 = unpack_geom('SCI', 'prescan', detector_areas) - + np.random.seed(4567) frame_list = [] for i in range(len(EMgain_arr)): for l in range(N): #number of frames to produce @@ -1859,7 +1859,7 @@ def add_defect(sch_imgs, prob, ori, temp): else: hdul.writeto(filename, overwrite = True) -def create_photon_countable_frames(Nbrights=30, Ndarks=40, EMgain=5000, kgain=7, exptime=0.1, cosmic_rate=0): +def create_photon_countable_frames(Nbrights=30, Ndarks=40, EMgain=5000, kgain=7, exptime=0.1, cosmic_rate=0, full_frame=True): '''This creates mock Dataset containing frames with large gain and short exposure time, illuminated and dark frames. Used for unit tests for photon counting. @@ -1870,12 +1870,13 @@ def create_photon_countable_frames(Nbrights=30, Ndarks=40, EMgain=5000, kgain=7, kgain (float): k gain (e-/DN) exptime (float): exposure time (in s) cosmic_rate (float) : simulated cosmic rays incidence, hits/cm^2/s - + full_frame (bool) : If True, simulated frames are SCI full frames. If False, 5x5 images are simulated. Defaults to True. Returns: dataset (corgidrp.data.Dataset): Dataset containing both the illuminated and dark frames ill_mean (float): mean electron count value simulated in the illuminated frames dark_mean (float): mean electron count value simulated in the dark frames ''' + np.random.seed(1234) pix_row = 1024 #number of rows and number of columns fluxmap = np.ones((pix_row,pix_row)) #photon flux map, photons/s @@ -1886,7 +1887,7 @@ def create_photon_countable_frames(Nbrights=30, Ndarks=40, EMgain=5000, kgain=7, dark_current=8.33e-4, # e-/pix/s cic=0.01, # e-/pix/frame read_noise=100., # e-/pix/frame - bias=0, # e-; 0 for simplicity + bias=2000, # e- qe=0.9, # quantum efficiency, e-/photon cr_rate=cosmic_rate, # cosmic rays incidence, hits/cm^2/s pixel_pitch=13e-6, # m @@ -1919,7 +1920,10 @@ def create_photon_countable_frames(Nbrights=30, Ndarks=40, EMgain=5000, kgain=7, prihdr, exthdr = create_default_headers() for i in range(Nbrights): # Simulate bright - frame_dn = emccd.sim_full_frame(fluxmap, exptime) + if full_frame: + frame_dn = emccd.sim_full_frame(fluxmap, exptime) + else: + frame_dn = emccd.sim_sub_frame(fluxmap[:5,:5], exptime) frame = data.Image(frame_dn, pri_hdr=prihdr, ext_hdr=exthdr) frame.ext_hdr['CMDGAIN'] = EMgain frame.ext_hdr['EXPTIME'] = exptime @@ -1928,7 +1932,10 @@ def create_photon_countable_frames(Nbrights=30, Ndarks=40, EMgain=5000, kgain=7, for i in range(Ndarks): # Simulate dark - frame_dn_dark = emccd.sim_full_frame(np.zeros_like(fluxmap), exptime) + if full_frame: + frame_dn_dark = emccd.sim_full_frame(np.zeros_like(fluxmap), exptime) + else: + frame_dn_dark = emccd.sim_sub_frame(np.zeros_like(fluxmap[:5,:5]), exptime) frame_dark = data.Image(frame_dn_dark, pri_hdr=prihdr.copy(), ext_hdr=exthdr.copy()) frame_dark.ext_hdr['CMDGAIN'] = EMgain frame_dark.ext_hdr['EXPTIME'] = exptime diff --git a/corgidrp/photon_counting.py b/corgidrp/photon_counting.py index 6629f211..42e0e63d 100644 --- a/corgidrp/photon_counting.py +++ b/corgidrp/photon_counting.py @@ -134,14 +134,20 @@ def get_pc_mean(input_dataset, detector_params, niter=2): ill_dataset = datasets[1-vals.index('DARK')] pc_means = [] + ups = [] + lows = [] + t = [] errs = [] dqs = [] + pc_means_up = [] + pc_means_low = [] for dataset in [ill_dataset, dark_dataset]: # getting number of read noise standard deviations at which to set the # photon-counting threshold T_factor = detector_params.params['T_factor'] # now get threshold to use for photon-counting - thresh = T_factor*dataset.frames[0].ext_hdr['RN'] + read_noise = dataset.frames[0].ext_hdr['RN'] + thresh = T_factor*read_noise if thresh < 0: raise PhotonCountException('thresh must be nonnegative') if not isinstance(niter, (int, np.integer)) or niter < 1: @@ -157,41 +163,58 @@ def get_pc_mean(input_dataset, detector_params, niter=2): if thresh >= em_gain: warnings.warn('thresh should be less than em_gain for effective ' 'photon counting') + if np.average(dataset.all_data) > 0.1: + warnings.warn('average # of photons/pixel is > 0.1. Decrease frame ' + 'time to get lower average # of photons/pixel.') + if read_noise <=0: + warnings.warn('read noise should be greater than 0 for effective ' + 'photon counting') + if thresh < 4*read_noise: + warnings.warn('thresh should be at least 4 or 5 times read_noise for ' + 'accurate photon counting') # Photon count stack of frames frames_pc = photon_count(dataset.all_data, thresh) - # frames_pc = np.array([photon_count(frame, thresh) for frame in dataset.all_data]) XXX bool_map = dataset.all_dq.astype(bool).astype(float) bool_map[bool_map > 0] = np.nan bool_map[bool_map == 0] = 1 nframes = np.nansum(bool_map, axis=0) + # upper and lower bounds for PC (for getting accurate err) + frames_pc_up = photon_count(dataset.all_data+dataset.all_err[:,0], thresh) + frames_pc_low = photon_count(dataset.all_data-dataset.all_err[:,0], thresh) frames_pc_masked = frames_pc * bool_map + frames_pc_masked_up = frames_pc_up * bool_map + frames_pc_masked_low = frames_pc_low * bool_map # Co-add frames frame_pc_coadded = np.nansum(frames_pc_masked, axis=0) + frame_pc_coadded_up = np.nansum(frames_pc_masked_up, axis=0) + frame_pc_coadded_low = np.nansum(frames_pc_masked_low, axis=0) # Correct for thresholding and coincidence loss; any pixel masked all the - # way through the stack will give NaN + # way through the stack may give NaN, but it should be masked in lam_newton_fit(); + # and it doesn't matter anyways since its DQ value will be 1 (it will be masked when the + # bad pixel correction is run, which comes after this photon-counting step) mean_expected = corr_photon_count(frame_pc_coadded, nframes, thresh, em_gain, niter) - - ##### error calculation - # combine all err frames to get standard error, which is due to inherent error in - # in the data before any photon counting is done - err_sum = np.sum(dataset.all_err**2, axis=0) - bool_map_sum = np.nansum(bool_map, axis=0) - comb_err = np.divide(np.sqrt(err_sum), np.sqrt(bool_map_sum), out=np.zeros_like(err_sum), - where=bool_map_sum != 0) - # expected error from photon counting + mean_expected_up = corr_photon_count(frame_pc_coadded_up, nframes, thresh, + em_gain, niter) + mean_expected_low = corr_photon_count(frame_pc_coadded_low, nframes, thresh, + em_gain, niter) + ##### error calculation: accounts for err coming from input dataset and + # statistical error from the photon-counting and photometric correction process. + # expected error from photon counting (biggest source from the actual values, not + # mean_expected_up or mean_expected_low): pc_variance = varL23(em_gain,mean_expected,thresh,nframes) - # now combine this error with the statistical error coming from the process of photon counting - total_err = np.sqrt(comb_err**2 + pc_variance) - total_dq = (bool_map_sum == 0).astype(int) # 1 where flagged, 0 otherwise + up = mean_expected_up + pc_variance + low = mean_expected_low - pc_variance + errs.append(np.max([up - mean_expected, mean_expected - low], axis=0)) + dq = (nframes == 0).astype(int) pc_means.append(mean_expected) - errs.append(total_err) - dqs.append(total_dq) + dqs.append(dq) # now subtract the dark PC mean combined_pc_mean = pc_means[0] - pc_means[1] + combined_pc_mean[combined_pc_mean<0] = 0 combined_err = np.sqrt(errs[0]**2 + errs[1]**2) combined_dq = np.bitwise_or(dqs[0], dqs[1]) pri_hdr = ill_dataset[0].pri_hdr.copy() @@ -205,7 +228,7 @@ def get_pc_mean(input_dataset, detector_params, niter=2): new_image.filename = ill_dataset[0].filename.split('.')[0]+'_pc.fits' new_image._record_parent_filenames(input_dataset) pc_dataset = data.Dataset([new_image]) - pc_dataset.update_after_processing_step("Photon-counted {0} frames using T_factor = {1} and niter={2}".format(len(input_dataset), T_factor, niter)) + pc_dataset.update_after_processing_step("Photon-counted {0} frames using T_factor={1} and niter={2}".format(len(input_dataset), T_factor, niter)) return pc_dataset @@ -215,11 +238,11 @@ def corr_photon_count(nobs, nfr, t, g, niter=2): Parameters ---------- nobs : array_like - Number of observations (Co-added photon counted frame). + Number of observations (Co-added photon-counted frame). nfr : int - Number of coadded frames. + Number of coadded frames, accounting for masked pixels in the frames. t : float - Photon counting threshold. + Photon-counting threshold. g : float EM gain. niter : int, optional @@ -228,7 +251,7 @@ def corr_photon_count(nobs, nfr, t, g, niter=2): Returns ------- lam : array_like - Mean expeted detected electrons per pixel (lambda). + Mean expeted electrons per pixel (lambda). """ # Get an approximate value of lambda for the first guess of the Newton fit @@ -261,7 +284,7 @@ def calc_lam_approx(nobs, nfr, t, g): Returns ------- array_like - Mean expeted (lambda). + Mean expected (lambda). """ # First step of equation (before taking log) @@ -302,12 +325,12 @@ def lam_newton_fit(nobs, nfr, t, g, lam0, niter): Returns ------- lam : array_like - Mean expeted (lambda). + Mean expected (lambda). """ # Mask out zero values to avoid divide by zero - lam_est_m = np.ma.masked_where(lam0 == 0, lam0) - nobs_m = np.ma.masked_where(nobs == 0, nobs) + lam_est_m = np.ma.masked_array(lam0, mask=(lam0==0)) + nobs_m = np.ma.masked_array(nobs, mask=(nobs==0)) # Iterate Newton's method for i in range(niter): @@ -315,7 +338,7 @@ def lam_newton_fit(nobs, nfr, t, g, lam0, niter): dfunc = _calc_dfunc(nfr, t, g, lam_est_m) lam_est_m -= func / dfunc - if lam_est_m.min() < 0: + if np.nanmin(lam_est_m.data) < 0: raise PhotonCountException('negative number of photon counts; ' 'try decreasing the frametime') @@ -327,42 +350,25 @@ def lam_newton_fit(nobs, nfr, t, g, lam0, niter): def _calc_func(nobs, nfr, t, g, lam): """Objective function for lambda.""" - e_thresh = ( - np.exp(-t/g) - * ( - t**2 * lam**2 - + 2*g * t * lam * (3 + lam) - + 2*g**2 * (6 + 3*lam + lam**2) - ) - / (2*g**2 * (6 + 3*lam + lam**2)) - ) - - e_coinloss = (1 - np.exp(-lam)) / lam - - #if (lam * nfr * e_thresh * e_coinloss).any() > nobs.any(): + epsilon_prime = (lam*(2*g**2*(6 + lam*(3 + lam)) + 2*g*lam*(3 + lam)*t + + lam**2*t**2))/(2.*np.e**(t/g)*g**2*(6 + lam*(6 + lam*(3 + lam)))) + + #if (nfr * epsilon_prime).any() > nobs.any(): # warnings.warn('Input photon flux is too high; decrease frametime') # This warning isn't necessary; could have a negative func but still # close enough to 0 for Newton's method - func = lam * nfr * e_thresh * e_coinloss - nobs + func = nfr * epsilon_prime - nobs return func def _calc_dfunc(nfr, t, g, lam): """Derivative wrt lambda of objective function.""" - dfunc = ( - (np.exp(-t/g - lam) * nfr) - / (2 * g**2 * (6 + 3*lam + lam**2)**2) - * ( - 2*g**2 * (6 + 3*lam + lam**2)**2 - + t**2 * lam * ( - -12 + 3*lam + 3*lam**2 + lam**3 + 3*np.exp(lam)*(4 + lam) - ) - + 2*g*t * ( - -18 + 6*lam + 15*lam**2 + 6*lam**3 + lam**4 - + 6*np.exp(lam)*(3 + 2*lam) - ) - ) - ) + dfunc = (lam*nfr*(2*g**2*(3 + 2*lam) + 2*g*lam*t + 2*g*(3 + lam)*t + + 2*lam*t**2))/(2.*np.e**(t/g)*g**2*(6 + lam*(6 + lam*(3 + lam)))) - (lam*(6 + + lam*(3 + lam) + lam*(3 + 2*lam))*nfr* + (2*g**2*(6 + lam*(3 + lam)) + 2*g*lam*(3 + lam)*t + lam**2*t**2))/(2.*np.e**(t/g)*g**2*(6 + + lam*(6 + lam*(3 + lam)))**2) + (nfr*(2*g**2*(6 + lam*(3 + lam)) + 2*g*lam*(3 + lam)*t + + lam**2*t**2))/(2.*np.e**(t/g)*g**2*(6 + lam*(6 + lam*(3 + lam)))) return dfunc \ No newline at end of file diff --git a/corgidrp/recipe_templates/l2a_to_l2b_pc.json b/corgidrp/recipe_templates/l2a_to_l2b_pc.json index 2f8b0236..82550c2e 100644 --- a/corgidrp/recipe_templates/l2a_to_l2b_pc.json +++ b/corgidrp/recipe_templates/l2a_to_l2b_pc.json @@ -1,5 +1,5 @@ { - "name" : "l1_to_l2b_pc", + "name" : "l2a_to_l2b_pc", "template" : true, "drpconfig" : { "track_individual_errors" : false diff --git a/tests/test_calibrate_darks_lsq.py b/tests/test_calibrate_darks_lsq.py index 1c22310a..2e85a144 100644 --- a/tests/test_calibrate_darks_lsq.py +++ b/tests/test_calibrate_darks_lsq.py @@ -207,6 +207,7 @@ def test_mean_num(): if __name__ == '__main__': + test_mean_num() test_expected_results_sub() test_sub_stack_len() test_g_arr_unique() @@ -214,7 +215,6 @@ def test_mean_num(): test_t_arr_unique() test_t_gtr_0() test_k_gtr_0() - test_mean_num() diff --git a/tests/test_kgain.py b/tests/test_kgain.py index 6313a58b..9da40a3d 100644 --- a/tests/test_kgain.py +++ b/tests/test_kgain.py @@ -17,6 +17,8 @@ def test_kgain(): dat = np.ones([1024,1024]) * 2 err = np.ones([1,1024,1024]) * 0.5 ptc = np.ones([2,1024]) + exthd["RN"] = 100 + exthd["RN_ERR"] = 2 ptc_hdr = fits.Header() image1 = data.Image(dat,pri_hdr = prhd, ext_hdr = exthd, err = err) image2 = image1.copy() @@ -84,7 +86,7 @@ def test_kgain(): assert gain_dataset[0].ext_hdr["BUNIT"] == "detected EM electrons" assert gain_dataset[0].err_hdr["BUNIT"] == "detected EM electrons" assert gain_dataset[0].ext_hdr["KGAIN"] == k_gain - assert gain_dataset[0].ext_hdr["KGAIN_ER"] == kgain.error + assert gain_dataset[0].ext_hdr["KGAIN_ER"] == kgain.error[0] assert gain_dataset[0].ext_hdr["RN"] > 0 assert gain_dataset[0].ext_hdr["RN_ERR"] > 0 assert("converted" in str(gain_dataset[0].ext_hdr["HISTORY"])) diff --git a/tests/test_mean_combine.py b/tests/test_mean_combine.py index 84dec301..7b8b0cbd 100644 --- a/tests/test_mean_combine.py +++ b/tests/test_mean_combine.py @@ -29,12 +29,12 @@ def test_mean_im(): - """Verify method calculates mean image and erorr term.""" + """Verify method calculates mean image and error term.""" tol = 1e-13 check_combined_im = np.mean(check_ims, axis=0) check_combined_im_err = np.sqrt(np.sum(np.array(check_ims)**2, - axis=0))/len(check_ims) + axis=0)/len(check_ims)) # For the pixel that is masked throughout check_combined_im[all_masks_i] = 0 check_combined_im_err[all_masks_i] = 0 @@ -42,7 +42,7 @@ def test_mean_im(): # For pixel that is only masked once check_combined_im[single_mask_i] = np.mean(unmasked_vals) check_combined_im_err[single_mask_i] = \ - np.sqrt(np.sum(unmasked_vals**2))/unmasked_vals.size + np.sqrt(np.sum(unmasked_vals**2)/unmasked_vals.size) combined_im, _, _, _ = mean_combine(check_ims, check_masks) combined_im_err, _, _, _ = mean_combine(check_ims, diff --git a/tests/test_photon_counting.py b/tests/test_photon_counting.py index 10e974f7..c0c20149 100644 --- a/tests/test_photon_counting.py +++ b/tests/test_photon_counting.py @@ -8,23 +8,25 @@ import corgidrp.walker as walker import corgidrp.caldb as caldb import corgidrp.data as data +from corgidrp.photon_counting import get_pc_mean, photon_count, PhotonCountException -dataset, ill_mean, dark_mean = mocks.create_photon_countable_frames(Nbrights=20, Ndarks=30, cosmic_rate=0) -thisfile_dir = os.path.dirname(__file__) # this file's folder -outputdir = os.path.join(thisfile_dir, 'simdata', 'pc_test_data') -if not os.path.exists(outputdir): - os.mkdir(outputdir) -# empty out directory of any previous files -for f in os.listdir(outputdir): - os.remove(os.path.join(outputdir,f)) -dataset.save(outputdir, ['pc_frame_{0}.fits'.format(i) for i in range(len(dataset))]) -l1_data_filelist = [] -for f in os.listdir(outputdir): - l1_data_filelist.append(os.path.join(outputdir, f)) +detector_params = data.DetectorParams({}) def test_expected_results(): '''Results are as expected theoretically. Also runs raw frames through pre-processing pipeline.''' - + dataset, ill_mean, dark_mean = mocks.create_photon_countable_frames(Nbrights=20, Ndarks=30, cosmic_rate=0) + thisfile_dir = os.path.dirname(__file__) # this file's folder + outputdir = os.path.join(thisfile_dir, 'simdata', 'pc_test_data') + if not os.path.exists(outputdir): + os.mkdir(outputdir) + # empty out directory of any previous files + for f in os.listdir(outputdir): + os.remove(os.path.join(outputdir,f)) + dataset.save(outputdir, ['pc_frame_{0}.fits'.format(i) for i in range(len(dataset))]) + l1_data_filelist = [] + for f in os.listdir(outputdir): + l1_data_filelist.append(os.path.join(outputdir, f)) + this_caldb = caldb.CalDB() # connection to cal DB # remove other KGain calibrations that may exist in case they don't have the added header keywords for i in range(len(this_caldb._db['Type'])): @@ -51,16 +53,51 @@ def test_expected_results(): pc_filename = f pc_file = os.path.join(outputdir, pc_filename) pc_frame = fits.getdata(pc_file) + pc_frame_err = fits.getdata(pc_file, 'ERR') pc_ext_hdr = fits.getheader(pc_file, 1) - assert np.isclose(pc_frame.mean(), ill_mean - dark_mean) + # more frames (which would take longer to run) would give an even better agreement than the 5% agreement below + assert np.isclose(pc_frame.mean(), ill_mean - dark_mean, rtol=0.05) assert 'niter=2' in pc_ext_hdr["HISTORY"][-1] assert 'T_factor=5' in pc_ext_hdr["HISTORY"][-1] + assert pc_frame_err.min() >= 0 + assert pc_frame_err[0][3,3] >= 0 + + +def test_negative(): + """Values at or below the + threshold are photon-counted as 0, including negative values.""" + fr = np.array([-3,2,0]) + pc_fr = photon_count(fr, thresh=2) + assert (pc_fr == np.array([0,0,0])).all() - # empty out directory now - for f in os.listdir(outputdir): - os.remove(os.path.join(outputdir,f)) - #XXX negative value test, others from PhotonCount? - #XXX change pc error: should come from doing pc on upper and lower bound of pixel values; same goes for trad and synthesized dark calibration +def test_masking_increases_err(): + '''Test that a pixel that is masked heavily (reducing the number of usable frames for that pixel) has a bigger err than the average pixel. + And if masking is all the way through for a certain pixel, that pixel will + be flagged in the DQ. Also, make sure there is failure when niter<1.''' + # exposure time too long to get reasonable photon-counted result (after photometric correction) + dataset_err, _, _ = mocks.create_photon_countable_frames(Nbrights=60, Ndarks=60, cosmic_rate=0, full_frame=False) + # instead of running through walker, just do the pre-processing steps simply + # using kgain=7 and bias=2000 and read noise = 100, from mocks.create_photon_countable_frames() + dataset_err.all_data = dataset_err.all_data.astype(float)*7 - 2000. + # masked for half of the bright and half of the dark frames + dataset_err.all_dq[30:90,3,3] = 1 + dataset_err.all_dq[:, 2, 2] = 1 #masked all the way through for (2,2) + for f in dataset_err.frames: + f.ext_hdr['RN'] = 100 + f.ext_hdr['KGAIN'] = 7 + pc_dataset_err = get_pc_mean(dataset_err, detector_params) + # the DQ for that pixel should be 1 + assert pc_dataset_err[0].dq[2,2] == 1 + # err for (3,3) is above the 95th percentile of error: + assert pc_dataset_err[0].err[0][3,3]>np.nanpercentile(pc_dataset_err[0].err,95) + # also when niter<1, exception + with pytest.raises(PhotonCountException): + get_pc_mean(dataset_err, detector_params, niter=0) + if __name__ == '__main__': - test_expected_results() \ No newline at end of file + test_masking_increases_err() + test_negative() + test_expected_results() + + \ No newline at end of file From 7ef22a959001f23c3d67dfe2b79c622d10e6644c Mon Sep 17 00:00:00 2001 From: Kevin Ludwick Date: Fri, 22 Nov 2024 21:21:03 -0600 Subject: [PATCH 03/26] linting --- corgidrp/photon_counting.py | 170 ++++++++++++++++-------------------- 1 file changed, 75 insertions(+), 95 deletions(-) diff --git a/corgidrp/photon_counting.py b/corgidrp/photon_counting.py index 42e0e63d..aca89df4 100644 --- a/corgidrp/photon_counting.py +++ b/corgidrp/photon_counting.py @@ -10,16 +10,14 @@ def varL23(g, L, T, N): '''Expected variance after photometric correction. See https://doi.org/10.1117/1.JATIS.9.4.048006 for details. - Parameters - ---------- - g (scalar): EM gain - L (2-D array): mean expected number of electrons - T (scalar): threshold - N (2-D array): number of frames + Args: + g (scalar): EM gain + L (2-D array): mean expected number of electrons + T (scalar): threshold + N (2-D array): number of frames - Returns - ------- - variance from photon counting and the photometric correction + Returns: + (float): variance from photon counting and the photometric correction ''' Const = 6/(6 + L*(6 + L*(3 + L))) @@ -39,18 +37,13 @@ class PhotonCountException(Exception): def photon_count(e_image, thresh): """Convert an analog image into a photon-counted image. - Parameters - ---------- - e_image : array_like, float - Analog image (e-). - thresh : float - Photon counting threshold (e-). Values > thresh will be assigned a 1, - values <= thresh will be assigned a 0. + Args: + e_image (array_like, float): Analog image (e-). + thresh (float): Photon counting threshold (e-). Values > thresh will be assigned a 1, + values <= thresh will be assigned a 0. - Returns - ------- - pc_image : array_like, float - Output digital image in units of photons. + Returns: + pc_image (array_like, float): Output digital image in units of photons. B Nemati and S Miller - UAH - 06-Aug-2018 @@ -96,23 +89,17 @@ def get_pc_mean(input_dataset, detector_params, niter=2): determined by "T_factor" stored in DetectorParams. The dark mean expected array is then subtracted from the illuminated mean expected array. - Parameters - ---------- - input_dataset (corgidrp.data.Dataset): - This is an instance of corgidrp.data.Dataset containing the - photon-counted illuminated frames as well as the - photon-counted dark frames, and all the frames must have the same - exposure time, EM gain, k gain, and read noise. - detector_params: (corgidrp.data.DetectorParams) - A calibration file storing detector calibration values. - niter : int, optional - Number of Newton's method iterations (used for photometric - correction). Defaults to 2. - - Returns - ------- - mean_expected : array_like - Mean expected photoelectron array (dark-subtracted) + Args: + input_dataset (corgidrp.data.Dataset): This is an instance of corgidrp.data.Dataset containing the + photon-counted illuminated frames as well as the + photon-counted dark frames, and all the frames must have the same + exposure time, EM gain, k gain, and read noise. + detector_params (corgidrp.data.DetectorParams): A calibration file storing detector calibration values. + niter (int, optional): Number of Newton's method iterations (used for photometric + correction). Defaults to 2. + + Returns: + pc_dataset (corgidrp.data.Dataset): Contains mean expected photoelectron array (dark-subtracted) References ---------- @@ -134,13 +121,8 @@ def get_pc_mean(input_dataset, detector_params, niter=2): ill_dataset = datasets[1-vals.index('DARK')] pc_means = [] - ups = [] - lows = [] - t = [] errs = [] dqs = [] - pc_means_up = [] - pc_means_low = [] for dataset in [ill_dataset, dark_dataset]: # getting number of read noise standard deviations at which to set the # photon-counting threshold @@ -235,23 +217,15 @@ def get_pc_mean(input_dataset, detector_params, niter=2): def corr_photon_count(nobs, nfr, t, g, niter=2): """Correct photon counted images. - Parameters - ---------- - nobs : array_like - Number of observations (Co-added photon-counted frame). - nfr : int - Number of coadded frames, accounting for masked pixels in the frames. - t : float - Photon-counting threshold. - g : float - EM gain. - niter : int, optional - Number of Newton's method iterations. Defaults to 2. - - Returns - ------- - lam : array_like - Mean expeted electrons per pixel (lambda). + Args: + nobs (array_like): Number of observations (Co-added photon-counted frame). + nfr (int): Number of coadded frames, accounting for masked pixels in the frames. + t (float): Photon-counting threshold. + g (float): EM gain. + niter (int, optional): Number of Newton's method iterations. Defaults to 2. + + Returns: + lam (array_like): Mean expeted electrons per pixel (lambda). """ # Get an approximate value of lambda for the first guess of the Newton fit @@ -270,21 +244,14 @@ def calc_lam_approx(nobs, nfr, t, g): that are out of bounds (e.g. from statistical fluctuations) it will revert to the zeroth order. - Parameters - ---------- - nobs : array_like - Number of observations (Co-added photon counted frame). - nfr : int - Number of coadded frames. - t : float - Photon counting threshold. - g : float - EM gain used when taking images. - - Returns - ------- - array_like - Mean expected (lambda). + Args: + nobs (array_like): Number of observations (Co-added photon counted frame). + nfr (int): Number of coadded frames. + t (float): Photon counting threshold. + g (float): EM gain used when taking images. + + Returns: + lam1 (array_like): Mean expected (lambda). """ # First step of equation (before taking log) @@ -307,25 +274,16 @@ def calc_lam_approx(nobs, nfr, t, g): def lam_newton_fit(nobs, nfr, t, g, lam0, niter): """Newton fit for finding lambda. - Parameters - ---------- - nobs : array_like - Number of observations (Co-added photon counted frame). - nfr : int - Number of coadded frames. - t : float - Photon counting threshold. - g : float - EM gain used when taking images. - lam0 : array_like - Initial guess for lambda. - niter : int - Number of Newton's fit iterations to take. - - Returns - ------- - lam : array_like - Mean expected (lambda). + Args: + nobs (array_like): Number of observations (Co-added photon counted frame). + nfr (int): Number of coadded frames. + t (float): Photon counting threshold. + g (float): EM gain used when taking images. + lam0 (array_like): Initial guess for lambda. + niter (int): Number of Newton's fit iterations to take. + + Returns: + lam (array_like): Mean expected (lambda). """ # Mask out zero values to avoid divide by zero @@ -349,7 +307,19 @@ def lam_newton_fit(nobs, nfr, t, g, lam0, niter): def _calc_func(nobs, nfr, t, g, lam): - """Objective function for lambda.""" + """Objective function for lambda for Newton's method for all applying photometric correction. + + Args: + nobs (array-like): number of frames per pixel that passed the threshold + nfr (array-like): number of unmasked frames per pixel total + t (float): threshold for photon counting + g (float): EM gain + lam (array-like): estimated mean expected electron count + + Returns: + func (array-like): objective function + + """ epsilon_prime = (lam*(2*g**2*(6 + lam*(3 + lam)) + 2*g*lam*(3 + lam)*t + lam**2*t**2))/(2.*np.e**(t/g)*g**2*(6 + lam*(6 + lam*(3 + lam)))) @@ -363,7 +333,17 @@ def _calc_func(nobs, nfr, t, g, lam): def _calc_dfunc(nfr, t, g, lam): - """Derivative wrt lambda of objective function.""" + """Derivative with respect to lambda of objective function. + + Args: + nfr (array-like): number of unmasked frames per pixel total + t (float): threshold for photon counting + g (float): EM gain + lam (array-like): estimated mean expected electron count + + Returns: + dfunc (array-like): derivative with respect to lambda of objective function + """ dfunc = (lam*nfr*(2*g**2*(3 + 2*lam) + 2*g*lam*t + 2*g*(3 + lam)*t + 2*lam*t**2))/(2.*np.e**(t/g)*g**2*(6 + lam*(6 + lam*(3 + lam)))) - (lam*(6 + lam*(3 + lam) + lam*(3 + 2*lam))*nfr* From 01482b6bdda026718c9873cd665c6fc024f16b63 Mon Sep 17 00:00:00 2001 From: Kevin Ludwick Date: Fri, 22 Nov 2024 21:27:48 -0600 Subject: [PATCH 04/26] lint --- corgidrp/mocks.py | 1 + 1 file changed, 1 insertion(+) diff --git a/corgidrp/mocks.py b/corgidrp/mocks.py index e485183c..b15a573e 100644 --- a/corgidrp/mocks.py +++ b/corgidrp/mocks.py @@ -1871,6 +1871,7 @@ def create_photon_countable_frames(Nbrights=30, Ndarks=40, EMgain=5000, kgain=7, exptime (float): exposure time (in s) cosmic_rate (float) : simulated cosmic rays incidence, hits/cm^2/s full_frame (bool) : If True, simulated frames are SCI full frames. If False, 5x5 images are simulated. Defaults to True. + Returns: dataset (corgidrp.data.Dataset): Dataset containing both the illuminated and dark frames ill_mean (float): mean electron count value simulated in the illuminated frames From 0d9f6f8803efbee327cfbe1e4977571deedc301a Mon Sep 17 00:00:00 2001 From: Kevin Ludwick Date: Fri, 22 Nov 2024 21:37:11 -0600 Subject: [PATCH 05/26] lint --- corgidrp/mocks.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/corgidrp/mocks.py b/corgidrp/mocks.py index b15a573e..b6857a19 100644 --- a/corgidrp/mocks.py +++ b/corgidrp/mocks.py @@ -1869,8 +1869,8 @@ def create_photon_countable_frames(Nbrights=30, Ndarks=40, EMgain=5000, kgain=7, EMgain (float): EM gain kgain (float): k gain (e-/DN) exptime (float): exposure time (in s) - cosmic_rate (float) : simulated cosmic rays incidence, hits/cm^2/s - full_frame (bool) : If True, simulated frames are SCI full frames. If False, 5x5 images are simulated. Defaults to True. + cosmic_rate: (float) simulated cosmic rays incidence, hits/cm^2/s + full_frame: (bool) If True, simulated frames are SCI full frames. If False, 5x5 images are simulated. Defaults to True. Returns: dataset (corgidrp.data.Dataset): Dataset containing both the illuminated and dark frames From 948ccd243a62b1422f4781bbe2ee10f6d2f6b62e Mon Sep 17 00:00:00 2001 From: Kevin Ludwick Date: Fri, 22 Nov 2024 23:06:17 -0600 Subject: [PATCH 06/26] eliminated potential division by 0 --- corgidrp/l2a_to_l2b.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/corgidrp/l2a_to_l2b.py b/corgidrp/l2a_to_l2b.py index a3db3bef..13d0ef5a 100644 --- a/corgidrp/l2a_to_l2b.py +++ b/corgidrp/l2a_to_l2b.py @@ -253,7 +253,7 @@ def convert_to_electrons(input_dataset, k_gain): kgainerr = k_gain.error[0] # x = c*kgain, where c (counts beforehand) and kgain both have error, so do propogation of error due to the product of 2 independent sources #[:,0:1,:,:] to get same dimensions as input_dataset.all_err - kgain_error = (np.abs(kgain_cube*kgain)*np.sqrt((kgain_dataset.all_err/kgain_cube)**2 + (kgainerr/kgain)**2))[:,0:1,:,:] + kgain_error = (np.sqrt((kgain*kgain_dataset.all_err)**2 + (kgain_cube*kgainerr)**2))[:,0:1,:,:] kgain_cube *= kgain history_msg = "data converted to detected EM electrons by kgain {0}".format(str(kgain)) From 8cabe716b06445278a029fb4677f97c0e4436f70 Mon Sep 17 00:00:00 2001 From: Kevin Ludwick Date: Fri, 22 Nov 2024 23:21:54 -0600 Subject: [PATCH 07/26] make a detectornoisemaps cal file for CI's sake so that I can leave the "OPTIONAL" out of the recipe in prescan subtraction since we do want the calibrated bias offset ideally --- tests/test_photon_counting.py | 16 ++++++++++++++++ 1 file changed, 16 insertions(+) diff --git a/tests/test_photon_counting.py b/tests/test_photon_counting.py index c0c20149..25df5f3d 100644 --- a/tests/test_photon_counting.py +++ b/tests/test_photon_counting.py @@ -8,6 +8,7 @@ import corgidrp.walker as walker import corgidrp.caldb as caldb import corgidrp.data as data +import corgidrp.detector as detector from corgidrp.photon_counting import get_pc_mean, photon_count, PhotonCountException detector_params = data.DetectorParams({}) @@ -46,6 +47,21 @@ def test_expected_results(): kgain.ext_hdr['RN_ERR'] = 0 kgain.save(filedir=outputdir, filename="mock_kgain.fits") this_caldb.create_entry(kgain) + + # NoiseMap + noise_map_dat = np.zeros((3, detector.detector_areas['SCI']['frame_rows'], detector.detector_areas['SCI']['frame_cols'])) + noise_map_noise = np.zeros([1,] + list(noise_map_dat.shape)) + noise_map_dq = np.zeros(noise_map_dat.shape, dtype=int) + err_hdr = fits.Header() + err_hdr['BUNIT'] = 'detected electrons' + ext_hdr['B_O'] = 0 + ext_hdr['B_O_ERR'] = 0 + noise_map = data.DetectorNoiseMaps(noise_map_dat, pri_hdr=pri_hdr, ext_hdr=ext_hdr, + input_dataset=mock_input_dataset, err=noise_map_noise, + dq = noise_map_dq, err_hdr=err_hdr) + noise_map.save(filedir=outputdir, filename="mock_detnoisemaps.fits") + this_caldb.create_entry(noise_map) + walker.walk_corgidrp(l1_data_filelist, '', outputdir, template="l1_to_l2b_pc.json") # get photon-counted frame for f in os.listdir(outputdir): From 87f6d95d148e67f141753302ade5546db54154ca Mon Sep 17 00:00:00 2001 From: Kevin Ludwick Date: Sat, 23 Nov 2024 10:33:44 -0600 Subject: [PATCH 08/26] had to add nonlin calib as well --- tests/test_photon_counting.py | 20 ++++++++++++++++++++ 1 file changed, 20 insertions(+) diff --git a/tests/test_photon_counting.py b/tests/test_photon_counting.py index 25df5f3d..8bbc32c2 100644 --- a/tests/test_photon_counting.py +++ b/tests/test_photon_counting.py @@ -34,6 +34,7 @@ def test_expected_results(): if this_caldb._db['Type'][i] == 'KGain': this_caldb._db = this_caldb._db.drop(i) this_caldb.save() + # KGain kgain_val = 7 # default value used in mocks.create_photon_countable_frames() pri_hdr, ext_hdr = mocks.create_default_headers() @@ -62,6 +63,22 @@ def test_expected_results(): noise_map.save(filedir=outputdir, filename="mock_detnoisemaps.fits") this_caldb.create_entry(noise_map) + # Make a non-linearity correction calibration file + input_non_linearity_filename = "nonlin_table_TVAC.txt" + input_non_linearity_path = os.path.join(os.path.dirname(__file__), "test_data", input_non_linearity_filename) + test_non_linearity_filename = input_non_linearity_filename.split(".")[0] + ".fits" + nonlin_fits_filepath = os.path.join(os.path.dirname(__file__), "test_data", test_non_linearity_filename) + tvac_nonlin_data = np.genfromtxt(input_non_linearity_path, delimiter=",") + + pri_hdr, ext_hdr = mocks.create_default_headers() + new_nonlinearity = data.NonLinearityCalibration(tvac_nonlin_data,pri_hdr=pri_hdr,ext_hdr=ext_hdr,input_dataset = dummy_dataset) + new_nonlinearity.filename = nonlin_fits_filepath + # fake the headers because this frame doesn't have the proper headers + new_nonlinearity.pri_hdr = pri_hdr + new_nonlinearity.ext_hdr = ext_hdr + this_caldb.create_entry(new_nonlinearity) + + walker.walk_corgidrp(l1_data_filelist, '', outputdir, template="l1_to_l2b_pc.json") # get photon-counted frame for f in os.listdir(outputdir): @@ -78,6 +95,9 @@ def test_expected_results(): assert pc_frame_err.min() >= 0 assert pc_frame_err[0][3,3] >= 0 + this_caldb.remove_entry(kgain) + this_caldb.remove_entry(noise_map) + this_caldb.remove_entry(new_nonlinearity) def test_negative(): """Values at or below the From 569ae0cd67160dbb143ae12cc5af766e6f251128 Mon Sep 17 00:00:00 2001 From: Kevin Ludwick Date: Sat, 23 Nov 2024 11:00:26 -0600 Subject: [PATCH 09/26] oops, forgot the dummy dataset --- tests/test_photon_counting.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/tests/test_photon_counting.py b/tests/test_photon_counting.py index 8bbc32c2..638df5d0 100644 --- a/tests/test_photon_counting.py +++ b/tests/test_photon_counting.py @@ -63,6 +63,10 @@ def test_expected_results(): noise_map.save(filedir=outputdir, filename="mock_detnoisemaps.fits") this_caldb.create_entry(noise_map) + # Nonlinearity calibration + ### Create a dummy non-linearity file #### + #Create a mock dataset because it is a required input when creating a NonLinearityCalibration + dummy_dataset = mocks.create_prescan_files() # Make a non-linearity correction calibration file input_non_linearity_filename = "nonlin_table_TVAC.txt" input_non_linearity_path = os.path.join(os.path.dirname(__file__), "test_data", input_non_linearity_filename) From e995bf372ac38d53afc7b574bde3b699610aaf6a Mon Sep 17 00:00:00 2001 From: Kevin Ludwick Date: Sat, 23 Nov 2024 17:20:23 -0600 Subject: [PATCH 10/26] fewer full frames simulated in hopes that the CI will not cancel the tests when it runs --- corgidrp/mocks.py | 6 +++--- corgidrp/photon_counting.py | 1 - tests/test_photon_counting.py | 14 ++++++++------ 3 files changed, 11 insertions(+), 10 deletions(-) diff --git a/corgidrp/mocks.py b/corgidrp/mocks.py index b6857a19..6ed993d5 100644 --- a/corgidrp/mocks.py +++ b/corgidrp/mocks.py @@ -1870,7 +1870,7 @@ def create_photon_countable_frames(Nbrights=30, Ndarks=40, EMgain=5000, kgain=7, kgain (float): k gain (e-/DN) exptime (float): exposure time (in s) cosmic_rate: (float) simulated cosmic rays incidence, hits/cm^2/s - full_frame: (bool) If True, simulated frames are SCI full frames. If False, 5x5 images are simulated. Defaults to True. + full_frame: (bool) If True, simulated frames are SCI full frames. If False, 50x50 images are simulated. Defaults to True. Returns: dataset (corgidrp.data.Dataset): Dataset containing both the illuminated and dark frames @@ -1924,7 +1924,7 @@ def create_photon_countable_frames(Nbrights=30, Ndarks=40, EMgain=5000, kgain=7, if full_frame: frame_dn = emccd.sim_full_frame(fluxmap, exptime) else: - frame_dn = emccd.sim_sub_frame(fluxmap[:5,:5], exptime) + frame_dn = emccd.sim_sub_frame(fluxmap[:50,:50], exptime) frame = data.Image(frame_dn, pri_hdr=prihdr, ext_hdr=exthdr) frame.ext_hdr['CMDGAIN'] = EMgain frame.ext_hdr['EXPTIME'] = exptime @@ -1936,7 +1936,7 @@ def create_photon_countable_frames(Nbrights=30, Ndarks=40, EMgain=5000, kgain=7, if full_frame: frame_dn_dark = emccd.sim_full_frame(np.zeros_like(fluxmap), exptime) else: - frame_dn_dark = emccd.sim_sub_frame(np.zeros_like(fluxmap[:5,:5]), exptime) + frame_dn_dark = emccd.sim_sub_frame(np.zeros_like(fluxmap[:50,:50]), exptime) frame_dark = data.Image(frame_dn_dark, pri_hdr=prihdr.copy(), ext_hdr=exthdr.copy()) frame_dark.ext_hdr['CMDGAIN'] = EMgain frame_dark.ext_hdr['EXPTIME'] = exptime diff --git a/corgidrp/photon_counting.py b/corgidrp/photon_counting.py index aca89df4..b6e2b705 100644 --- a/corgidrp/photon_counting.py +++ b/corgidrp/photon_counting.py @@ -305,7 +305,6 @@ def lam_newton_fit(nobs, nfr, t, g, lam0, niter): return lam - def _calc_func(nobs, nfr, t, g, lam): """Objective function for lambda for Newton's method for all applying photometric correction. diff --git a/tests/test_photon_counting.py b/tests/test_photon_counting.py index 638df5d0..aa497d9a 100644 --- a/tests/test_photon_counting.py +++ b/tests/test_photon_counting.py @@ -15,7 +15,7 @@ def test_expected_results(): '''Results are as expected theoretically. Also runs raw frames through pre-processing pipeline.''' - dataset, ill_mean, dark_mean = mocks.create_photon_countable_frames(Nbrights=20, Ndarks=30, cosmic_rate=0) + dataset, ill_mean, dark_mean = mocks.create_photon_countable_frames(Nbrights=5, Ndarks=6, cosmic_rate=0) thisfile_dir = os.path.dirname(__file__) # this file's folder outputdir = os.path.join(thisfile_dir, 'simdata', 'pc_test_data') if not os.path.exists(outputdir): @@ -82,7 +82,6 @@ def test_expected_results(): new_nonlinearity.ext_hdr = ext_hdr this_caldb.create_entry(new_nonlinearity) - walker.walk_corgidrp(l1_data_filelist, '', outputdir, template="l1_to_l2b_pc.json") # get photon-counted frame for f in os.listdir(outputdir): @@ -114,12 +113,14 @@ def test_negative(): def test_masking_increases_err(): '''Test that a pixel that is masked heavily (reducing the number of usable frames for that pixel) has a bigger err than the average pixel. And if masking is all the way through for a certain pixel, that pixel will - be flagged in the DQ. Also, make sure there is failure when niter<1.''' + be flagged in the DQ. Also, make sure there is failure when niter<1. Also, very good agreement with theoretically expected value since we use more frames.''' # exposure time too long to get reasonable photon-counted result (after photometric correction) - dataset_err, _, _ = mocks.create_photon_countable_frames(Nbrights=60, Ndarks=60, cosmic_rate=0, full_frame=False) + dataset_err, ill_mean, dark_mean = mocks.create_photon_countable_frames(Nbrights=60, Ndarks=60, cosmic_rate=0, full_frame=False) # instead of running through walker, just do the pre-processing steps simply - # using kgain=7 and bias=2000 and read noise = 100, from mocks.create_photon_countable_frames() - dataset_err.all_data = dataset_err.all_data.astype(float)*7 - 2000. + # using kgain=7 and bias=2000 and read noise = 100 and QE=0.9 (quantum efficiency), from mocks.create_photon_countable_frames() + for f in dataset_err.frames: + f.data = f.data.astype(float)*7 - 2000. + # masked for half of the bright and half of the dark frames dataset_err.all_dq[30:90,3,3] = 1 dataset_err.all_dq[:, 2, 2] = 1 #masked all the way through for (2,2) @@ -127,6 +128,7 @@ def test_masking_increases_err(): f.ext_hdr['RN'] = 100 f.ext_hdr['KGAIN'] = 7 pc_dataset_err = get_pc_mean(dataset_err, detector_params) + assert np.isclose(pc_dataset_err.all_data.mean(), ill_mean - dark_mean, rtol=0.001) # the DQ for that pixel should be 1 assert pc_dataset_err[0].dq[2,2] == 1 # err for (3,3) is above the 95th percentile of error: From 3b1b96947efc3b3bc10e0d3578e00205f9525657 Mon Sep 17 00:00:00 2001 From: Kevin Ludwick Date: Sat, 23 Nov 2024 18:47:06 -0600 Subject: [PATCH 11/26] oops, forgot to actually change the number --- tests/test_photon_counting.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/test_photon_counting.py b/tests/test_photon_counting.py index aa497d9a..69f9613a 100644 --- a/tests/test_photon_counting.py +++ b/tests/test_photon_counting.py @@ -91,8 +91,8 @@ def test_expected_results(): pc_frame = fits.getdata(pc_file) pc_frame_err = fits.getdata(pc_file, 'ERR') pc_ext_hdr = fits.getheader(pc_file, 1) - # more frames (which would take longer to run) would give an even better agreement than the 5% agreement below - assert np.isclose(pc_frame.mean(), ill_mean - dark_mean, rtol=0.05) + # more frames (which would take longer to run) would give an even better agreement than the 25% agreement below + assert np.isclose(pc_frame.mean(), ill_mean - dark_mean, rtol=0.25) assert 'niter=2' in pc_ext_hdr["HISTORY"][-1] assert 'T_factor=5' in pc_ext_hdr["HISTORY"][-1] assert pc_frame_err.min() >= 0 From 5d1f521a03b43f1ea47e154a154e6721c93e4f99 Mon Sep 17 00:00:00 2001 From: Kevin Ludwick Date: Thu, 26 Dec 2024 10:29:00 -0600 Subject: [PATCH 12/26] made separate recipes for PC dark and PC illuminated, made e2e, addressed minor issues 200 and 244 --- corgidrp/caldb.py | 10 +- corgidrp/calibrate_kgain.py | 34 ++- corgidrp/calibrate_nonlin.py | 44 +-- corgidrp/combine.py | 4 +- corgidrp/darks.py | 5 +- corgidrp/data.py | 25 +- corgidrp/l2a_to_l2b.py | 19 +- corgidrp/mocks.py | 35 ++- corgidrp/photon_counting.py | 281 +++++++++++------- corgidrp/recipe_templates/l1_to_l2b_pc.json | 2 +- .../recipe_templates/l1_to_l2b_pc_dark.json | 52 ++++ corgidrp/recipe_templates/l2a_to_l2b_pc.json | 2 +- .../recipe_templates/l2a_to_l2b_pc_dark.json | 26 ++ tests/e2e_tests/flatfield_e2e.py | 24 +- tests/e2e_tests/noisemap_cal_e2e.py | 2 +- tests/e2e_tests/photon_count_e2e.py | 152 ++++++++++ tests/test_badpixelmap.py | 4 +- tests/test_build_synthesized_dark.py | 1 + tests/test_caldb.py | 8 +- tests/test_calibrate_darks_lsq.py | 2 + tests/test_combine.py | 4 +- tests/test_dark_sub.py | 8 + tests/test_data_levels.py | 3 + tests/test_dataset.py | 2 + tests/test_err_dq.py | 1 + tests/test_frame_selection.py | 2 + tests/test_photon_counting.py | 180 +++++------ tests/test_walker.py | 1 + 28 files changed, 653 insertions(+), 280 deletions(-) create mode 100644 corgidrp/recipe_templates/l1_to_l2b_pc_dark.json create mode 100644 corgidrp/recipe_templates/l2a_to_l2b_pc_dark.json create mode 100644 tests/e2e_tests/photon_count_e2e.py diff --git a/corgidrp/caldb.py b/corgidrp/caldb.py index f60b8d84..a38c7f9f 100644 --- a/corgidrp/caldb.py +++ b/corgidrp/caldb.py @@ -168,7 +168,7 @@ def _get_values_from_entry(self, entry, is_calib=True): drp_version, obsid, naxis1, - naxis2, + naxis2 ] # rest are ext_hdr keys we can copy @@ -243,7 +243,7 @@ def remove_entry(self, entry, to_disk=True): def get_calib(self, frame, dtype, to_disk=True): """ - Outputs the best calibration file of the given type for the input sciene frame. + Outputs the best calibration file of the given type for the input science frame. Args: frame (corgidrp.data.Image): an image frame to request a calibration for. If None is passed in, looks for the @@ -288,14 +288,12 @@ def get_calib(self, frame, dtype, to_disk=True): options = calibdf.loc[ ( (calibdf["EXPTIME"] == frame_dict["EXPTIME"]) - & (calibdf["NAXIS1"] == frame_dict["NAXIS1"]) - & (calibdf["NAXIS2"] == frame_dict["NAXIS2"]) ) ] if len(options) == 0: - raise ValueError("No valid Dark with EXPTIME={0} and dimension ({1},{2})" - .format(frame_dict["EXPTIME"], frame_dict["NAXIS1"], frame_dict["NAXIS2"])) + raise ValueError("No valid Dark with EXPTIME={0})" + .format(frame_dict["EXPTIME"])) # select the one closest in time result_index = np.abs(options["MJD"] - frame_dict["MJD"]).argmin() diff --git a/corgidrp/calibrate_kgain.py b/corgidrp/calibrate_kgain.py index 8a2ef1fe..3fefc0f4 100644 --- a/corgidrp/calibrate_kgain.py +++ b/corgidrp/calibrate_kgain.py @@ -34,25 +34,25 @@ def check_kgain_params( ): """ Checks integrity of kgain parameters in the dictionary kgain_params. """ if 'offset_colroi1' not in kgain_params: - raise ValueError('Missing parameter in directory pointer YAML file.') + raise ValueError('Missing parameter.') if 'offset_colroi2' not in kgain_params: - raise ValueError('Missing parameter in directory pointer YAML file.') + raise ValueError('Missing parameter.') if 'rowroi1' not in kgain_params: - raise ValueError('Missing parameter in directory pointer YAML file.') + raise ValueError('Missing parameter.') if 'rowroi2' not in kgain_params: - raise ValueError('Missing parameter in directory pointer YAML file.') + raise ValueError('Missing parameter.') if 'colroi1' not in kgain_params: - raise ValueError('Missing parameter in directory pointer YAML file.') + raise ValueError('Missing parameter.') if 'colroi2' not in kgain_params: - raise ValueError('Missing parameter in directory pointer YAML file.') + raise ValueError('Missing parameter.') if 'rn_bins1' not in kgain_params: - raise ValueError('Missing parameter in directory pointer YAML file.') + raise ValueError('Missing parameter.') if 'rn_bins2' not in kgain_params: - raise ValueError('Missing parameter in directory pointer YAML file.') + raise ValueError('Missing parameter.') if 'max_DN_val' not in kgain_params: - raise ValueError('Missing parameter in directory pointer YAML file.') + raise ValueError('Missing parameter.') if 'signal_bins_N' not in kgain_params: - raise ValueError('Missing parameter in directory pointer YAML file.') + raise ValueError('Missing parameter.') if not isinstance(kgain_params['offset_colroi1'], (float, int)): raise TypeError('offset_colroi1 is not a number') @@ -281,7 +281,7 @@ def calibrate_kgain(dataset_kgain, n_cal=10, n_mean=30, min_val=800, max_val=3000, binwidth=68, make_plot=True,plot_outdir='figures', show_plot=False, logspace_start=-1, logspace_stop=4, logspace_num=200, - verbose=False, detector_regions=None): + verbose=False, detector_regions=None, kgain_params=kgain_params): """ Given an array of frame stacks for various exposure times, each sub-stack having at least 5 illuminated pupil L1 SCI-size frames having the same @@ -343,6 +343,18 @@ def calibrate_kgain(dataset_kgain, detector_regions (dict): a dictionary of detector geometry properties. Keys should be as found in detector_areas in detector.py. Defaults to that dictionary. + kgain_params (dict): (Optional) Dictionary containing row and col specifications + for the region of interest (indicated by 'rowroi1','rowroi2','colroi1',and 'colroi2'). + The 'roi' needs one square region specified, and 'back' needs two square regions, + where a '1' ending indicates the smaller of two values, and a '2' ending indicates the larger + of two values. The coordinates of the square region are specified by matching + up as follows: (rowroi1, colroi1), (rowroi2, colroi1), etc. + Also must contain: + 'rn_bins1': lower bound of counts histogram for fitting or read noise + 'rn_bins2': upper bound of counts histogram for fitting or read noise + 'max_DN_val': maximum DN value to be included in photon transfer curve (PTC) + 'signal_bins_N': number of bins in the signal variables of PTC curve + Defaults to kgain_params included in this file. Returns: corgidrp.data.KGain: kgain estimate from the least-squares fit to the photon diff --git a/corgidrp/calibrate_nonlin.py b/corgidrp/calibrate_nonlin.py index d1a68aea..a762fab4 100644 --- a/corgidrp/calibrate_nonlin.py +++ b/corgidrp/calibrate_nonlin.py @@ -48,37 +48,37 @@ def check_nonlin_params( ): """ Checks integrity of kgain parameters in the dictionary nonlin_params. """ if 'rowroi1' not in nonlin_params: - raise ValueError('Missing parameter in directory pointer YAML file.') + raise ValueError('Missing parameter') if 'rowroi2' not in nonlin_params: - raise ValueError('Missing parameter in directory pointer YAML file.') + raise ValueError('Missing parameter') if 'colroi1' not in nonlin_params: - raise ValueError('Missing parameter in directory pointer YAML file.') + raise ValueError('Missing parameter') if 'colroi2' not in nonlin_params: - raise ValueError('Missing parameter in directory pointer YAML file.') + raise ValueError('Missing parameter') if 'rowback11' not in nonlin_params: - raise ValueError('Missing parameter in directory pointer YAML file.') + raise ValueError('Missing parameter') if 'rowback12' not in nonlin_params: - raise ValueError('Missing parameter in directory pointer YAML file.') + raise ValueError('Missing parameter') if 'rowback21' not in nonlin_params: - raise ValueError('Missing parameter in directory pointer YAML file.') + raise ValueError('Missing parameter') if 'rowback22' not in nonlin_params: - raise ValueError('Missing parameter in directory pointer YAML file.') + raise ValueError('Missing parameter') if 'colback11' not in nonlin_params: - raise ValueError('Missing parameter in directory pointer YAML file.') + raise ValueError('Missing parameter') if 'colback12' not in nonlin_params: - raise ValueError('Missing parameter in directory pointer YAML file.') + raise ValueError('Missing parameter') if 'colback21' not in nonlin_params: - raise ValueError('Missing parameter in directory pointer YAML file.') + raise ValueError('Missing parameter') if 'colback22' not in nonlin_params: - raise ValueError('Missing parameter in directory pointer YAML file.') + raise ValueError('Missing parameter') if 'min_exp_time' not in nonlin_params: - raise ValueError('Missing parameter in directory pointer YAML file.') + raise ValueError('Missing parameter') if 'num_bins' not in nonlin_params: - raise ValueError('Missing parameter in directory pointer YAML file.') + raise ValueError('Missing parameter') if 'min_bin' not in nonlin_params: - raise ValueError('Missing parameter in directory pointer YAML file.') + raise ValueError('Missing parameter') if 'min_mask_factor' not in nonlin_params: - raise ValueError('Missing parameter in directory pointer YAML file.') + raise ValueError('Missing parameter') if not isinstance(nonlin_params['rowroi1'], (float, int)): raise TypeError('rowroi1 is not a number') @@ -124,7 +124,7 @@ def calibrate_nonlin(dataset_nl, pfit_upp_cutoff1 = -2, pfit_upp_cutoff2 = -3, pfit_low_cutoff1 = 2, pfit_low_cutoff2 = 1, make_plot=True, plot_outdir='figures', show_plot=False, - verbose=False): + verbose=False, nonlin_params=nonlin_params): """ Given a large array of stacks with 1 or more EM gains, and sub-stacks of frames ranging over exposure time, each sub-stack having at least 1 illuminated @@ -211,6 +211,16 @@ def calibrate_nonlin(dataset_nl, show_plot (bool): (Optional) display the plots. Default is False. verbose (bool): (Optional) display various diagnostic print messages. Default is False. + nonlin_params (dict): (Optional) Dictionary of row and col specifications + for the region of interest (indicated by 'roi') where the frame is illuminated and for + two background regions (indicated by 'back1' and 'back2') where the frame is not illuminated. + Must contain 'rowroi1','rowroi2','colroi1','colroi2','rowback11','rowback12', + 'rowback21','rowback22','colback11','colback12','colback21',and 'colback22'. + The 'roi' needs one square region specified, and 'back' needs two square regions, + where a '1' ending indicates the smaller of two values, and a '2' ending indicates the larger + of two values. The coordinates of each square are specified by matching + up as follows: (rowroi1, colroi1), (rowroi1, colroi2), (rowback11, colback11), + (rowback11, colback12), etc. Defaults to nonlin_params specified in this file. Returns: nonlin_arr (NonLinearityCalibration): 2-D array with nonlinearity values diff --git a/corgidrp/combine.py b/corgidrp/combine.py index 7a7c1abb..ca49391a 100644 --- a/corgidrp/combine.py +++ b/corgidrp/combine.py @@ -43,7 +43,7 @@ def combine_images(data_subset, err_subset, dq_subset, collapse, num_frames_scal # dq collpase: keep all flags on dq_collapse = np.bitwise_or.reduce(dq_subset, axis=0) - # except those pixels that have been replaced + # except those pixels that have not been replaced dq_collapse[np.where((dq_collapse > 0) & (~np.isnan(data_collapse)))] = 0 return data_collapse, err_collapse, dq_collapse @@ -72,7 +72,7 @@ def combine_subexposures(input_dataset, num_frames_per_group=None, collapse="mea num_frames_per_group = len(input_dataset) if len(input_dataset) % num_frames_per_group != 0: - raise ValueError("Input dataset of lenght {0} cannot be grouped in sets of {1}".format(len(input_dataset, num_frames_per_group))) + raise ValueError("Input dataset of lenght {0} cannot be grouped in sets of {1}".format(len(input_dataset), num_frames_per_group)) if collapse.lower() not in ["mean", "median"]: raise ValueError("combine_subexposures can only collapse with mean or median") diff --git a/corgidrp/darks.py b/corgidrp/darks.py index dfb33ff4..86d5eae8 100644 --- a/corgidrp/darks.py +++ b/corgidrp/darks.py @@ -1,6 +1,7 @@ import numpy as np import warnings from astropy.io import fits +import os from corgidrp.detector import slice_section, imaging_slice, imaging_area_geom, detector_areas import corgidrp.check as check @@ -577,7 +578,7 @@ def calibrate_darks_lsq(dataset, detector_params, detector_regions=None): # input data error comes from .err arrays; could use this for error bars # in input data for weighted least squares, but we'll just simply get the # std error and add it in quadrature to least squares fit standard dev - stacks_err = np.sqrt(np.sum(mean_err_stack**2, axis=0)/len(mean_err_stack)) + stacks_err = np.sqrt(np.sum(mean_err_stack**2, axis=0)/np.sqrt(len(mean_err_stack))) # matrix to be used for least squares and covariance matrix X = np.array([np.ones([len(EMgain_arr)]), EMgain_arr, EMgain_arr*exptime_arr]).T @@ -865,5 +866,5 @@ def build_synthesized_dark(dataset, noisemaps, detector_regions=None, full_frame master_dark = Dark(md_data, prihdr, exthdr, input_data, md_noise, FDCdq, errhdr) - + return master_dark \ No newline at end of file diff --git a/corgidrp/data.py b/corgidrp/data.py index d18aa99e..f6198bfc 100644 --- a/corgidrp/data.py +++ b/corgidrp/data.py @@ -623,7 +623,7 @@ class Dark(Image): data_or_filepath (str or np.array): either the filepath to the FITS file to read in OR the 2D image data pri_hdr (astropy.io.fits.Header): the primary header (required only if raw 2D data is passed in) ext_hdr (astropy.io.fits.Header): the image extension header (required only if raw 2D data is passed in) - input_dataset (corgidrp.data.Dataset): the Image files combined together to make this noise map (required only if raw 2D data is passed in and if raw data filenames not already archived in ext_hdr) + input_dataset (corgidrp.data.Dataset): the Image files combined together to make this dark (required only if raw 2D data is passed in and if raw data filenames not already archived in ext_hdr) err (np.array): the error array (required only if raw data is passed in) err_hdr (astropy.io.fits.Header): the error header (required only if raw data is passed in) dq (np.array): the DQ array (required only if raw data is passed in) @@ -640,7 +640,8 @@ def __init__(self, data_or_filepath, pri_hdr=None, ext_hdr=None, input_dataset=N raise ValueError("This appears to be a new dark. The dataset of input files needs to be passed in to the input_dataset keyword to record history of this dark.") self.ext_hdr['DATATYPE'] = 'Dark' # corgidrp specific keyword for saving to disk self.ext_hdr['BUNIT'] = 'detected electrons' - + if 'PC_STAT' not in ext_hdr: + self.ext_hdr['PC_STAT'] = 'analog master dark' # log all the data that went into making this calibration file if 'DRPNFILE' not in ext_hdr.keys() and input_dataset is not None: self._record_parent_filenames(input_dataset) @@ -651,8 +652,15 @@ def __init__(self, data_or_filepath, pri_hdr=None, ext_hdr=None, input_dataset=N # give it a default filename using the first input file as the base # strip off everything starting at .fits if input_dataset is not None: - orig_input_filename = input_dataset[0].filename.split(".fits")[0] - self.filename = "{0}_dark.fits".format(orig_input_filename) + if self.ext_hdr['PC_STAT'] != 'photon-counted master dark': + orig_input_filename = input_dataset[0].filename.split(".fits")[0] + self.filename = "{0}_dark.fits".format(orig_input_filename) + else: + orig_input_filename = input_dataset[0].filename.split(".fits")[0] + self.filename = "{0}_pc_dark.fits".format(orig_input_filename) + + if 'PC_STAT' not in self.ext_hdr: + self.ext_hdr['PC_STAT'] = 'analog master dark' if err_hdr is not None: self.err_hdr['BUNIT'] = 'detected electrons' @@ -830,7 +838,8 @@ def __init__(self, data_or_filepath, pri_hdr=None, ext_hdr=None, raise ValueError("File that was loaded was not a NonLinearityCalibration file.") if self.ext_hdr['DATATYPE'] != 'NonLinearityCalibration': raise ValueError("File that was loaded was not a NonLinearityCalibration file.") - + self.dq_hdr['COMMENT'] = 'DQ not meaningful for this calibration; just present for class consistency' + class KGain(Image): """ Class for KGain calibration file. Until further insights it is just one float value. @@ -1024,6 +1033,8 @@ def __init__(self, data_or_filepath, pri_hdr=None, ext_hdr=None, input_dataset=N raise ValueError("File that was loaded was not a BadPixelMap file.") if self.ext_hdr['DATATYPE'] != 'BadPixelMap': raise ValueError("File that was loaded was not a BadPixelMap file.") + self.dq_hdr['COMMENT'] = 'DQ not meaningful for this calibration; just present for class consistency' + self.err_hdr['COMMENT'] = 'err not meaningful for this calibration; just present for class consistency' def copy(self, copy_data = True): """ @@ -1354,7 +1365,8 @@ def __init__(self, data_or_filepath, pri_hdr=None, ext_hdr=None, input_dataset=N # check that this is actually an AstrometricCalibration file that was read in if 'DATATYPE' not in self.ext_hdr or self.ext_hdr['DATATYPE'] != 'AstrometricCalibration': raise ValueError("File that was loaded was not an AstrometricCalibration file.") - + self.dq_hdr['COMMENT'] = 'DQ not meaningful for this calibration; just present for class consistency' + class TrapCalibration(Image): """ @@ -1401,6 +1413,7 @@ def __init__(self,data_or_filepath, pri_hdr=None,ext_hdr=None, input_dataset=Non # since if only a filepath was passed in, any file could have been read in if 'DATATYPE' not in self.ext_hdr or self.ext_hdr['DATATYPE'] != 'TrapCalibration': raise ValueError("File that was loaded was not a TrapCalibration file.") + self.dq_hdr['COMMENT'] = 'DQ not meaningful for this calibration; just present for class consistency' datatypes = { "Image" : Image, "Dark" : Dark, diff --git a/corgidrp/l2a_to_l2b.py b/corgidrp/l2a_to_l2b.py index 13d0ef5a..1a7ced19 100644 --- a/corgidrp/l2a_to_l2b.py +++ b/corgidrp/l2a_to_l2b.py @@ -81,10 +81,14 @@ def dark_subtraction(input_dataset, dark, detector_regions=None, outputdir=None) outputdir = '.' #current directory dark.save(filedir=outputdir) elif type(dark) is data.Dark: - # In this case, the Dark loaded in should already match the arry dimensions - # of input_dataset, specified by full_frame argument of build_trad_dark - # when this Dark was built - pass + if 'PC_STAT' in dark.ext_hdr: + if dark.ext_hdr['PC_STAT'] != 'analog master dark': + raise Exception('The input \'dark\' is a photon-counted dark and cannot be used in the dark_subtraction step function, which is intended for analog frames.') + # In this case, the Dark loaded in should already match the arry dimensions + # of input_dataset, specified by full_frame argument of build_trad_dark + # when this Dark was built + if (dark.ext_hdr['EXPTIME'], dark.ext_hdr['CMDGAIN'], dark.ext_hdr['KGAIN']) != unique_vals[0]: + raise Exception('Dark should have the same EXPTIME, CMDGAIN, and KGAIN as input_dataset.') else: raise Exception('dark type should be either corgidrp.data.Dark or corgidrp.data.DetectorNoiseMaps.') @@ -252,15 +256,16 @@ def convert_to_electrons(input_dataset, k_gain): kgain = k_gain.value #extract from caldb kgainerr = k_gain.error[0] # x = c*kgain, where c (counts beforehand) and kgain both have error, so do propogation of error due to the product of 2 independent sources - #[:,0:1,:,:] to get same dimensions as input_dataset.all_err - kgain_error = (np.sqrt((kgain*kgain_dataset.all_err)**2 + (kgain_cube*kgainerr)**2))[:,0:1,:,:] + kgain_error = np.zeros_like(input_dataset.all_err) + kgain_tot_err = (np.sqrt((kgain*kgain_dataset.all_err[:,0,:,:])**2 + (kgain_cube*kgainerr)**2)) + kgain_error[:,0,:,:] = kgain_tot_err kgain_cube *= kgain history_msg = "data converted to detected EM electrons by kgain {0}".format(str(kgain)) # update the output dataset with this converted data and update the history kgain_dataset.update_after_processing_step(history_msg, new_all_data=kgain_cube, new_all_err=kgain_error, header_entries = {"BUNIT":"detected EM electrons", "KGAIN":kgain, - "KGAIN_ER": k_gain.error[0], "RN":k_gain.ext_hdr['RN'], "RN_ERR":k_gain.ext_hdr["RN_ERR"]}) + "KGAIN_ER": k_gain.error[0], "RN":k_gain.ext_hdr['RN'], "RN_ERR":k_gain.ext_hdr["RN_ERR"]}) return kgain_dataset def em_gain_division(input_dataset): diff --git a/corgidrp/mocks.py b/corgidrp/mocks.py index 6ed993d5..b8b6c249 100644 --- a/corgidrp/mocks.py +++ b/corgidrp/mocks.py @@ -167,7 +167,7 @@ def create_synthesized_master_dark_calib(detector_areas): # image area, including "shielded" rows and cols: imrows, imcols, imr0c0 = imaging_area_geom('SCI', detector_areas) prerows, precols, prer0c0 = unpack_geom('SCI', 'prescan', detector_areas) - np.random.seed(4567) + frame_list = [] for i in range(len(EMgain_arr)): for l in range(N): #number of frames to produce @@ -235,7 +235,8 @@ def create_dark_calib_files(filedir=None, numfiles=10): for i in range(numfiles): prihdr, exthdr = create_default_headers() exthdr['KGAIN'] = 7 - np.random.seed(456+i); sim_data = np.random.poisson(lam=150., size=(1200, 2200)).astype(np.float64) + #np.random.seed(456+i); + sim_data = np.random.poisson(lam=150., size=(1200, 2200)).astype(np.float64) frame = data.Image(sim_data, pri_hdr=prihdr, ext_hdr=exthdr) if filedir is not None: frame.save(filedir=filedir, filename=filepattern.format(i)) @@ -264,7 +265,8 @@ def create_simflat_dataset(filedir=None, numfiles=10): for i in range(numfiles): prihdr, exthdr = create_default_headers() # generate images in normal distribution with mean 1 and std 0.01 - np.random.seed(456+i); sim_data = np.random.poisson(lam=150., size=(1024, 1024)).astype(np.float64) + #np.random.seed(456+i); + sim_data = np.random.poisson(lam=150., size=(1024, 1024)).astype(np.float64) frame = data.Image(sim_data, pri_hdr=prihdr, ext_hdr=exthdr) if filedir is not None: frame.save(filedir=filedir, filename=filepattern.format(i)) @@ -451,7 +453,8 @@ def create_flatfield_dummy(filedir=None, numfiles=2): frames=[] for i in range(numfiles): prihdr, exthdr = create_default_headers() - np.random.seed(456+i); sim_data = np.random.normal(loc=1.0, scale=0.01, size=(1024, 1024)) + #np.random.seed(456+i); + sim_data = np.random.normal(loc=1.0, scale=0.01, size=(1024, 1024)) frame = data.Image(sim_data, pri_hdr=prihdr, ext_hdr=exthdr) if filedir is not None: frame.save(filedir=filedir, filename=filepattern.format(i)) @@ -490,7 +493,8 @@ def create_nonlinear_dataset(nonlin_filepath, filedir=None, numfiles=2,em_gain=2 data_range = np.linspace(800,65536,size) # Generate data for each row, where the mean increase from 10 to 65536 for x in range(size): - np.random.seed(120+x); sim_data[:, x] = np.random.poisson(data_range[x], size).astype(np.float64) + #np.random.seed(120+x); + sim_data[:, x] = np.random.poisson(data_range[x], size).astype(np.float64) non_linearity_correction = data.NonLinearityCalibration(nonlin_filepath) @@ -543,7 +547,7 @@ def create_cr_dataset(nonlin_filepath, filedir=None, datetime=None, numfiles=2, im_width = dataset.all_data.shape[-1] # Overwrite dataset with a poisson distribution - np.random.seed(123) + #np.random.seed(123) dataset.all_data[:,:,:] = np.random.poisson(lam=150,size=dataset.all_data.shape).astype(np.float64) # Loop over images in dataset @@ -554,7 +558,7 @@ def create_cr_dataset(nonlin_filepath, filedir=None, datetime=None, numfiles=2, # Pick random locations to add a cosmic ray for x in range(numCRs): - np.random.seed(123+x) + #np.random.seed(123+x) loc = np.round(np.random.uniform(0,im_width-1, size=2)).astype(int) # Add the CR plateau @@ -1859,8 +1863,8 @@ def add_defect(sch_imgs, prob, ori, temp): else: hdul.writeto(filename, overwrite = True) -def create_photon_countable_frames(Nbrights=30, Ndarks=40, EMgain=5000, kgain=7, exptime=0.1, cosmic_rate=0, full_frame=True): - '''This creates mock Dataset containing frames with large gain and short exposure time, illuminated and dark frames. +def create_photon_countable_frames(Nbrights=30, Ndarks=40, EMgain=5000, kgain=7, exptime=0.05, cosmic_rate=0, full_frame=True): + '''This creates mock L1 Dataset containing frames with large gain and short exposure time, illuminated and dark frames. Used for unit tests for photon counting. Args: @@ -1873,11 +1877,11 @@ def create_photon_countable_frames(Nbrights=30, Ndarks=40, EMgain=5000, kgain=7, full_frame: (bool) If True, simulated frames are SCI full frames. If False, 50x50 images are simulated. Defaults to True. Returns: - dataset (corgidrp.data.Dataset): Dataset containing both the illuminated and dark frames + ill_dataset (corgidrp.data.Dataset): Dataset containing the illuminated frames + dark_dataset (corgidrp.data.Dataset): Dataset containing the dark frames ill_mean (float): mean electron count value simulated in the illuminated frames dark_mean (float): mean electron count value simulated in the dark frames ''' - np.random.seed(1234) pix_row = 1024 #number of rows and number of columns fluxmap = np.ones((pix_row,pix_row)) #photon flux map, photons/s @@ -1888,7 +1892,7 @@ def create_photon_countable_frames(Nbrights=30, Ndarks=40, EMgain=5000, kgain=7, dark_current=8.33e-4, # e-/pix/s cic=0.01, # e-/pix/frame read_noise=100., # e-/pix/frame - bias=2000, # e- + bias=20000, # e- qe=0.9, # quantum efficiency, e-/photon cr_rate=cosmic_rate, # cosmic rays incidence, hits/cm^2/s pixel_pitch=13e-6, # m @@ -1929,6 +1933,7 @@ def create_photon_countable_frames(Nbrights=30, Ndarks=40, EMgain=5000, kgain=7, frame.ext_hdr['CMDGAIN'] = EMgain frame.ext_hdr['EXPTIME'] = exptime frame.pri_hdr["VISTYPE"] = "TDEMO" + frame.filename = 'L1_for_pc_ill_{0}.fits'.format(i) frame_e_list.append(frame) for i in range(Ndarks): @@ -1941,8 +1946,10 @@ def create_photon_countable_frames(Nbrights=30, Ndarks=40, EMgain=5000, kgain=7, frame_dark.ext_hdr['CMDGAIN'] = EMgain frame_dark.ext_hdr['EXPTIME'] = exptime frame_dark.pri_hdr["VISTYPE"] = "DARK" + frame.filename = 'L1_for_pc_dark_{0}.fits'.format(i) frame_e_dark_list.append(frame_dark) - dataset = data.Dataset(frame_e_list+frame_e_dark_list) + ill_dataset = data.Dataset(frame_e_list) + dark_dataset = data.Dataset(frame_e_dark_list) - return dataset, ill_mean, dark_mean \ No newline at end of file + return ill_dataset, dark_dataset, ill_mean, dark_mean \ No newline at end of file diff --git a/corgidrp/photon_counting.py b/corgidrp/photon_counting.py index b6e2b705..be807736 100644 --- a/corgidrp/photon_counting.py +++ b/corgidrp/photon_counting.py @@ -4,6 +4,7 @@ """ import warnings import numpy as np +from astropy.io import fits import corgidrp.data as data def varL23(g, L, T, N): @@ -59,12 +60,11 @@ def photon_count(e_image, thresh): return pc_image -def get_pc_mean(input_dataset, detector_params, niter=2): - """Take a stack of images, illuminated and dark frames of the same exposure +def get_pc_mean(input_dataset, pc_master_dark=None, T_factor=None, pc_ecount_max=None, niter=2, mask_filepath=None, safemode=True): + """Take a stack of images, frames of the same exposure time, k gain, read noise, and EM gain, and return the mean expected value per pixel. The frames are taken in photon-counting mode, which means high EM - gain and very low exposure time. The number of darks can be different from - the number of illuminated frames. + gain and very low exposure time. The frames in each stack should be SCI image (1024x1024) frames that have had some of the L2b steps applied: @@ -82,24 +82,44 @@ def get_pc_mean(input_dataset, detector_params, niter=2): fixed-pattern noise, and read noise. For such low-exposure frames, though, this really should not be a concern. - This algorithm will photon count each frame in each of the 2 sets (illuminated and dark) individually, - then co-add the photon counted frames. The co-added frame is then averaged + This algorithm will photon count each frame individually, + then co-add the photon-counted frames. The co-added frame is then averaged and corrected for thresholding and coincidence loss, returning the mean - expected array in units of photoelectrons. The threshold is - determined by "T_factor" stored in DetectorParams. The dark mean expected array - is then subtracted from the illuminated mean expected array. + expected array in units of photoelectrons if dark-subtracted (detected electrons + if not dark-subtracted). The threshold is determined by the input "T_factor", + and the value stored in DetectorParams is used if this input is None. + + This function can be used for photon-counting illuminated and dark datasets + (see Args below). Args: input_dataset (corgidrp.data.Dataset): This is an instance of corgidrp.data.Dataset containing the - photon-counted illuminated frames as well as the - photon-counted dark frames, and all the frames must have the same - exposure time, EM gain, k gain, and read noise. - detector_params (corgidrp.data.DetectorParams): A calibration file storing detector calibration values. + frames to photon-count. All the frames must have the same + exposure time, EM gain, k gain, and read noise. If the input dataset's header key 'VISTYPE' is equal to 'DARK', a + photon-counted master dark calibration product will be produced. Otherwise, the input dataset is assumed to consist of illuminated + frames intended for photon-counting. + pc_master_dark (corgidrp.data.Dark): Photon-counted master dark to be used for dark subtraction. + If None, no dark subtraction is done. + T_factor (float): The number of read noise standard deviations at which to set the threshold for photon counting. If None, the value is drawn from corgidrp.data.DetectorParams. + Defaults to None. + pc_ecount_max (float): Maximum allowed electrons/pixel/frame for photon counting. If None, the value is drawn from corgidrp.data.DetectorParams. + Defaults to None. niter (int, optional): Number of Newton's method iterations (used for photometric correction). Defaults to 2. + mask_filepath (str): Filepath to a .fits file with a pixel mask in the default HDU. The mask should be an array of the same shape as that found in each frame of the input_dataset, and the mask + should be 0 where the region of interest is and 1 for pixels to be masked (not considered). + safemode (bool): If False, the function does not halt due to an exception (useful for the iterative process of digging the dark hole). + If True, the function halts with an exception if the mean intensity of the unmasked pixels (or all pixels if no mask provided) is greater than pc_ecount_max or if the minimum photon-counted pixel value is negative. + If False, the function gives a warning for these instead. + Defaults to True. Returns: - pc_dataset (corgidrp.data.Dataset): Contains mean expected photoelectron array (dark-subtracted) + corgidrp.data.Dataset or corgidrp.data.Dark: If If the input dataset's header key 'VISTYPE' is not equal to 'DARK', + corgidrp.data.Dataset is the output type, and the output is the processed illuminated set, whether + dark subtraction happened or not. Contains mean expected array (detected electrons if not dark-subtracted, + photoelectrons if dark-subtracted). + If the input dataset's header key 'VISTYPE' is equal to 'DARK', corgidrp.data.Dark is the output type, and the output + is the processed dark set. Contains mean expected array (detected electrons). References ---------- @@ -110,111 +130,164 @@ def get_pc_mean(input_dataset, detector_params, niter=2): Kevin Ludwick - UAH - 2023 """ + if not isinstance(niter, (int, np.integer)) or niter < 1: + raise PhotonCountException('niter must be an integer greater than ' + '0') test_dataset, _ = input_dataset.copy().split_dataset(exthdr_keywords=['EXPTIME', 'CMDGAIN', 'KGAIN', 'RN']) if len(test_dataset) > 1: raise PhotonCountException('All frames must have the same exposure time, ' 'commanded EM gain, and k gain.') - datasets, vals = test_dataset[0].split_dataset(prihdr_keywords=['VISTYPE']) - if len(vals) != 2: - raise PhotonCountException('There should only be 2 datasets, and one should have \'VISTYPE\' = \'DARK\'.') - dark_dataset = datasets[vals.index('DARK')] - ill_dataset = datasets[1-vals.index('DARK')] + datasets, val = test_dataset[0].split_dataset(prihdr_keywords=['VISTYPE']) + if len(val) != 1: + raise PhotonCountException('There should only be 1 \'VISTYPE\' value for the dataset.') + if val[0] == 'DARK': + if pc_master_dark is not None: + raise PhotonCountException('The input frames are \'VISTYPE\'=\'DARK\' frames, so no pc_master_dark should be provided.') + + dataset = datasets[0] pc_means = [] errs = [] dqs = [] - for dataset in [ill_dataset, dark_dataset]: - # getting number of read noise standard deviations at which to set the - # photon-counting threshold + # getting number of read noise standard deviations at which to set the + # photon-counting threshold + if T_factor is None: + detector_params = data.DetectorParams({}) T_factor = detector_params.params['T_factor'] - # now get threshold to use for photon-counting - read_noise = dataset.frames[0].ext_hdr['RN'] - thresh = T_factor*read_noise - if thresh < 0: - raise PhotonCountException('thresh must be nonnegative') - if not isinstance(niter, (int, np.integer)) or niter < 1: - raise PhotonCountException('niter must be an integer greater than ' - '0') - try: # if EM gain measured directly from frame TODO change hdr name if necessary - em_gain = dataset.frames[0].ext_hdr['EMGAIN_M'] - except: - try: # use applied EM gain if available - em_gain = dataset.frames[0].ext_hdr['EMGAIN_A'] - except: # use commanded gain otherwise - em_gain = dataset.frames[0].ext_hdr['CMDGAIN'] - if thresh >= em_gain: + # getting maximum allowed electrons/pixel/frame for photon counting + if pc_ecount_max is None: + detector_params = data.DetectorParams({}) + pc_ecount_max = detector_params.params['pc_ecount_max'] + if mask_filepath is None: + mask = np.zeros_like(dataset.frames[0].data) + else: + mask = fits.getdata(mask_filepath) + + considered_indices = (mask==0) + + # now get threshold to use for photon-counting + read_noise = test_dataset[0].frames[0].ext_hdr['RN'] + thresh = T_factor*read_noise + if thresh < 0: + raise PhotonCountException('thresh must be nonnegative') + + if 'EMGAIN_M' in dataset.frames[0].ext_hdr: # if EM gain measured directly from frame TODO change hdr name if necessary + em_gain = dataset.frames[0].ext_hdr['EMGAIN_M'] + elif 'EMGAIN_A' in dataset.frames[0].ext_hdr: # use applied EM gain if available + em_gain = dataset.frames[0].ext_hdr['EMGAIN_A'] + elif 'CMDGAIN' in dataset.frames[0].ext_hdr: # use commanded gain otherwise + em_gain = dataset.frames[0].ext_hdr['CMDGAIN'] + + if thresh >= em_gain: + if safemode: + raise Exception('thresh should be less than em_gain for effective ' + 'photon counting') + else: warnings.warn('thresh should be less than em_gain for effective ' 'photon counting') - if np.average(dataset.all_data) > 0.1: - warnings.warn('average # of photons/pixel is > 0.1. Decrease frame ' + if np.nanmean(dataset.all_data[:, considered_indices]/em_gain) > pc_ecount_max: + if safemode: + raise Exception('average # of electrons/pixel is > pc_ecount_max, which means ' + 'the average # of photons/pixel may be > pc_ecount_max (depending on the QE). Can decrease frame ' 'time to get lower average # of photons/pixel.') - if read_noise <=0: - warnings.warn('read noise should be greater than 0 for effective ' - 'photon counting') - if thresh < 4*read_noise: - warnings.warn('thresh should be at least 4 or 5 times read_noise for ' - 'accurate photon counting') - - # Photon count stack of frames - frames_pc = photon_count(dataset.all_data, thresh) - bool_map = dataset.all_dq.astype(bool).astype(float) - bool_map[bool_map > 0] = np.nan - bool_map[bool_map == 0] = 1 - nframes = np.nansum(bool_map, axis=0) - # upper and lower bounds for PC (for getting accurate err) - frames_pc_up = photon_count(dataset.all_data+dataset.all_err[:,0], thresh) - frames_pc_low = photon_count(dataset.all_data-dataset.all_err[:,0], thresh) - frames_pc_masked = frames_pc * bool_map - frames_pc_masked_up = frames_pc_up * bool_map - frames_pc_masked_low = frames_pc_low * bool_map - # Co-add frames - frame_pc_coadded = np.nansum(frames_pc_masked, axis=0) - frame_pc_coadded_up = np.nansum(frames_pc_masked_up, axis=0) - frame_pc_coadded_low = np.nansum(frames_pc_masked_low, axis=0) - - # Correct for thresholding and coincidence loss; any pixel masked all the - # way through the stack may give NaN, but it should be masked in lam_newton_fit(); - # and it doesn't matter anyways since its DQ value will be 1 (it will be masked when the - # bad pixel correction is run, which comes after this photon-counting step) - mean_expected = corr_photon_count(frame_pc_coadded, nframes, thresh, - em_gain, niter) - mean_expected_up = corr_photon_count(frame_pc_coadded_up, nframes, thresh, - em_gain, niter) - mean_expected_low = corr_photon_count(frame_pc_coadded_low, nframes, thresh, - em_gain, niter) - ##### error calculation: accounts for err coming from input dataset and - # statistical error from the photon-counting and photometric correction process. - # expected error from photon counting (biggest source from the actual values, not - # mean_expected_up or mean_expected_low): - pc_variance = varL23(em_gain,mean_expected,thresh,nframes) - up = mean_expected_up + pc_variance - low = mean_expected_low - pc_variance - errs.append(np.max([up - mean_expected, mean_expected - low], axis=0)) - dq = (nframes == 0).astype(int) - pc_means.append(mean_expected) - dqs.append(dq) + else: + warnings.warn('average # of electrons/pixel is > pc_ecount_max, which means ' + 'the average # of photons/pixel may be > pc_ecount_max (depending on the QE). Can decrease frame ' + 'time to get lower average # of photons/pixel.') + if read_noise <=0: + raise Exception('read noise should be greater than 0 for effective ' + 'photon counting') + if thresh < 4*read_noise: # leave as just a warning + warnings.warn('thresh should be at least 4 or 5 times read_noise for ' + 'accurate photon counting') + + # Photon count stack of frames + frames_pc = photon_count(dataset.all_data, thresh) + bool_map = dataset.all_dq.astype(bool).astype(float) + bool_map[bool_map > 0] = np.nan + bool_map[bool_map == 0] = 1 + nframes = np.nansum(bool_map, axis=0) + # upper and lower bounds for PC (for getting accurate err) + frames_pc_up = photon_count(dataset.all_data+dataset.all_err[:,0], thresh) + frames_pc_low = photon_count(dataset.all_data-dataset.all_err[:,0], thresh) + frames_pc_masked = frames_pc * bool_map + frames_pc_masked_up = frames_pc_up * bool_map + frames_pc_masked_low = frames_pc_low * bool_map + # Co-add frames + frame_pc_coadded = np.nansum(frames_pc_masked, axis=0) + frame_pc_coadded_up = np.nansum(frames_pc_masked_up, axis=0) + frame_pc_coadded_low = np.nansum(frames_pc_masked_low, axis=0) + + # Correct for thresholding and coincidence loss; any pixel masked all the + # way through the stack may give NaN, but it should be masked in lam_newton_fit(); + # and it doesn't matter anyways since its DQ value will be 1 (it will be masked when the + # bad pixel correction is run, which comes after this photon-counting step) + mean_expected = corr_photon_count(frame_pc_coadded, nframes, thresh, + em_gain, considered_indices, niter) + mean_expected_up = corr_photon_count(frame_pc_coadded_up, nframes, thresh, + em_gain, considered_indices, niter) + mean_expected_low = corr_photon_count(frame_pc_coadded_low, nframes, thresh, + em_gain, considered_indices, niter) + ##### error calculation: accounts for err coming from input dataset and + # statistical error from the photon-counting and photometric correction process. + # expected error from photon counting (biggest source from the actual values, not + # mean_expected_up or mean_expected_low): + pc_variance = varL23(em_gain,mean_expected,thresh,nframes) + up = mean_expected_up + pc_variance + low = mean_expected_low - pc_variance + errs.append(np.max([up - mean_expected, mean_expected - low], axis=0)) + dq = (nframes == 0).astype(int) + pc_means.append(mean_expected) + dqs.append(dq) + if pc_master_dark is not None: + if type(pc_master_dark) is not data.Dark: + raise Exception('Input type for pc_master_dark must be corgidrp.data.Dark.') + if 'PC_STAT' not in pc_master_dark.ext_hdr: + raise PhotonCountException('\'PC_STAT\' must be a key in the extension header of pc_master_dark.') + if pc_master_dark.ext_hdr['PC_STAT'] != 'photon-counted master dark': + raise PhotonCountException('pc_master_dark must be a photon-counted master dark (i.e., ' + 'the extension header key \'PC_STAT\' must be \'photon-counted master dark\').') + pc_means.append(pc_master_dark.data) + dqs.append(pc_master_dark.dq) + errs.append(pc_master_dark.err) + dark_sub = "yes" + filename_end = '' + else: + pc_means.append(np.zeros_like(pc_means[0])) + dqs.append(np.zeros_like(pc_means[0]).astype(int)) + errs.append(np.zeros_like(pc_means[0])) + dark_sub = "no" + filename_end = '_no_ds' + # now subtract the dark PC mean combined_pc_mean = pc_means[0] - pc_means[1] combined_pc_mean[combined_pc_mean<0] = 0 combined_err = np.sqrt(errs[0]**2 + errs[1]**2) combined_dq = np.bitwise_or(dqs[0], dqs[1]) - pri_hdr = ill_dataset[0].pri_hdr.copy() - ext_hdr = ill_dataset[0].ext_hdr.copy() - err_hdr = ill_dataset[0].err_hdr.copy() - dq_hdr = ill_dataset[0].dq_hdr.copy() - hdulist = ill_dataset[0].hdu_list.copy() - new_image = data.Image(combined_pc_mean, pri_hdr=pri_hdr, ext_hdr=ext_hdr, err=combined_err, dq=combined_dq, err_hdr=err_hdr, + pri_hdr = dataset[0].pri_hdr.copy() + ext_hdr = dataset[0].ext_hdr.copy() + err_hdr = dataset[0].err_hdr.copy() + dq_hdr = dataset[0].dq_hdr.copy() + hdulist = dataset[0].hdu_list.copy() + + if val[0] != "DARK": + new_image = data.Image(combined_pc_mean, pri_hdr=pri_hdr, ext_hdr=ext_hdr, err=combined_err, dq=combined_dq, err_hdr=err_hdr, dq_hdr=dq_hdr, input_hdulist=hdulist) - - new_image.filename = ill_dataset[0].filename.split('.')[0]+'_pc.fits' - new_image._record_parent_filenames(input_dataset) - pc_dataset = data.Dataset([new_image]) - pc_dataset.update_after_processing_step("Photon-counted {0} frames using T_factor={1} and niter={2}".format(len(input_dataset), T_factor, niter)) - - return pc_dataset - -def corr_photon_count(nobs, nfr, t, g, niter=2): + new_image.filename = dataset[0].filename.split('.')[0]+'_pc'+filename_end+'.fits' + new_image._record_parent_filenames(input_dataset) + pc_ill_dataset = data.Dataset([new_image]) + pc_ill_dataset.update_after_processing_step("Photon-counted {0} frames using T_factor={1} and niter={2}. Dark-subtracted: {3}.".format(len(input_dataset), T_factor, niter, dark_sub)) + return pc_ill_dataset + else: + ext_hdr['PC_STAT'] = 'photon-counted master dark' + ext_hdr['NAXIS1'] = combined_pc_mean.shape[0] + ext_hdr['NAXIS2'] = combined_pc_mean.shape[1] + pc_dark = data.Dark(combined_pc_mean, pri_hdr=pri_hdr, ext_hdr=ext_hdr, err=combined_err, dq=combined_dq, err_hdr=err_hdr, input_dataset=input_dataset) + return pc_dark + + +def corr_photon_count(nobs, nfr, t, g, mask_indices, niter=2): """Correct photon counted images. Args: @@ -222,6 +295,7 @@ def corr_photon_count(nobs, nfr, t, g, niter=2): nfr (int): Number of coadded frames, accounting for masked pixels in the frames. t (float): Photon-counting threshold. g (float): EM gain. + mask_indices (array-like): indices of pixel positions to use. niter (int, optional): Number of Newton's method iterations. Defaults to 2. Returns: @@ -232,7 +306,7 @@ def corr_photon_count(nobs, nfr, t, g, niter=2): lam0 = calc_lam_approx(nobs, nfr, t, g) # Use Newton's method to converge at a value for lambda - lam = lam_newton_fit(nobs, nfr, t, g, lam0, niter) + lam = lam_newton_fit(nobs, nfr, t, g, lam0, niter, mask_indices) return lam @@ -271,7 +345,7 @@ def calc_lam_approx(nobs, nfr, t, g): return lam1 -def lam_newton_fit(nobs, nfr, t, g, lam0, niter): +def lam_newton_fit(nobs, nfr, t, g, lam0, niter, mask_indices): """Newton fit for finding lambda. Args: @@ -281,6 +355,7 @@ def lam_newton_fit(nobs, nfr, t, g, lam0, niter): g (float): EM gain used when taking images. lam0 (array_like): Initial guess for lambda. niter (int): Number of Newton's fit iterations to take. + mask_indices (array-like): indices of pixel positions to use. Returns: lam (array_like): Mean expected (lambda). @@ -296,7 +371,7 @@ def lam_newton_fit(nobs, nfr, t, g, lam0, niter): dfunc = _calc_dfunc(nfr, t, g, lam_est_m) lam_est_m -= func / dfunc - if np.nanmin(lam_est_m.data) < 0: + if np.nanmin(lam_est_m.data[mask_indices]) < 0: raise PhotonCountException('negative number of photon counts; ' 'try decreasing the frametime') diff --git a/corgidrp/recipe_templates/l1_to_l2b_pc.json b/corgidrp/recipe_templates/l1_to_l2b_pc.json index b35144ba..a7d1b379 100644 --- a/corgidrp/recipe_templates/l1_to_l2b_pc.json +++ b/corgidrp/recipe_templates/l1_to_l2b_pc.json @@ -45,7 +45,7 @@ { "name": "get_pc_mean", "calibs" : { - "DetectorParams" : "AUTOMATIC" + "Dark" : "AUTOMATIC" } }, { diff --git a/corgidrp/recipe_templates/l1_to_l2b_pc_dark.json b/corgidrp/recipe_templates/l1_to_l2b_pc_dark.json new file mode 100644 index 00000000..3e2c41d7 --- /dev/null +++ b/corgidrp/recipe_templates/l1_to_l2b_pc_dark.json @@ -0,0 +1,52 @@ +{ + "name" : "l1_to_l2b_pc_dark", + "template" : true, + "drpconfig" : { + "track_individual_errors" : false + }, + "inputs" : [], + "outputdir" : "", + "steps" : [ + { + "name" : "prescan_biassub", + "calibs" : { + "DetectorNoiseMaps" : "AUTOMATIC" + }, + "keywords" : { + "return_full_frame" : false + } + }, + { + "name" : "detect_cosmic_rays", + "calibs" : { + "DetectorParams" : "AUTOMATIC", + "KGain" : "AUTOMATIC" + } + }, + { + "name" : "correct_nonlinearity", + "calibs" : { + "NonLinearityCalibration" : "AUTOMATIC" + } + }, + { + "name" : "update_to_l2a" + }, + + { + "name" : "frame_select" + }, + { + "name" : "convert_to_electrons", + "calibs" : { + "KGain" : "AUTOMATIC" + } + }, + { + "name": "get_pc_mean" + }, + { + "name" : "save" + } + ] +} \ No newline at end of file diff --git a/corgidrp/recipe_templates/l2a_to_l2b_pc.json b/corgidrp/recipe_templates/l2a_to_l2b_pc.json index 82550c2e..4775ae38 100644 --- a/corgidrp/recipe_templates/l2a_to_l2b_pc.json +++ b/corgidrp/recipe_templates/l2a_to_l2b_pc.json @@ -19,7 +19,7 @@ { "name": "get_pc_mean", "calibs" : { - "DetectorParams" : "AUTOMATIC" + "Dark" : "AUTOMATIC" } }, { diff --git a/corgidrp/recipe_templates/l2a_to_l2b_pc_dark.json b/corgidrp/recipe_templates/l2a_to_l2b_pc_dark.json new file mode 100644 index 00000000..c1b3df62 --- /dev/null +++ b/corgidrp/recipe_templates/l2a_to_l2b_pc_dark.json @@ -0,0 +1,26 @@ +{ + "name" : "l2a_to_l2b_pc_dark", + "template" : true, + "drpconfig" : { + "track_individual_errors" : false + }, + "inputs" : [], + "outputdir" : "", + "steps" : [ + { + "name" : "frame_select" + }, + { + "name" : "convert_to_electrons", + "calibs" : { + "KGain" : "AUTOMATIC" + } + }, + { + "name": "get_pc_mean" + }, + { + "name" : "save" + } + ] +} \ No newline at end of file diff --git a/tests/e2e_tests/flatfield_e2e.py b/tests/e2e_tests/flatfield_e2e.py index dc27edd4..3c9885c2 100644 --- a/tests/e2e_tests/flatfield_e2e.py +++ b/tests/e2e_tests/flatfield_e2e.py @@ -120,11 +120,21 @@ def test_flat_creation_neptune(tvacdata_path, e2eoutput_path): this_caldb.create_entry(nonlinear_cal) # KGain + # remove other KGain calibrations that may exist in case they don't have the added header keywords + for i in range(len(this_caldb._db['Type'])): + if this_caldb._db['Type'][i] == 'KGain': + this_caldb._db = this_caldb._db.drop(i) + elif this_caldb._db['Type'][i] == 'FlatField': + this_caldb._db = this_caldb._db.drop(i) kgain_val = 8.7 + # add in keywords not provided by create_default_headers() (since L1 headers are simulated from that function) + ext_hdr['RN'] = 100 + ext_hdr['RN_ERR'] = 0 kgain = data.KGain(np.array([[kgain_val]]), pri_hdr=pri_hdr, ext_hdr=ext_hdr, input_dataset=mock_input_dataset) kgain.save(filedir=flat_outputdir, filename="mock_kgain.fits") - this_caldb.create_entry(kgain) + this_caldb.create_entry(kgain, to_disk=False) + this_caldb.save() # NoiseMap with fits.open(fpn_path) as hdulist: @@ -290,11 +300,21 @@ def test_flat_creation_uranus(tvacdata_path, e2eoutput_path): this_caldb.create_entry(nonlinear_cal) # KGain + # remove other KGain calibrations that may exist in case they don't have the added header keywords + for i in range(len(this_caldb._db['Type'])): + if this_caldb._db['Type'][i] == 'KGain': + this_caldb._db = this_caldb._db.drop(i) + elif this_caldb._db['Type'][i] == 'FlatField': + this_caldb._db = this_caldb._db.drop(i) kgain_val = 8.7 + # add in keywords not provided by create_default_headers() (since L1 headers are simulated from that function) + ext_hdr['RN'] = 100 + ext_hdr['RN_ERR'] = 0 kgain = data.KGain(np.array([[kgain_val]]), pri_hdr=pri_hdr, ext_hdr=ext_hdr, input_dataset=mock_input_dataset) kgain.save(filedir=flat_outputdir, filename="mock_kgain.fits") - this_caldb.create_entry(kgain) + this_caldb.create_entry(kgain, to_disk=False) + this_caldb.save() # NoiseMap with fits.open(fpn_path) as hdulist: diff --git a/tests/e2e_tests/noisemap_cal_e2e.py b/tests/e2e_tests/noisemap_cal_e2e.py index 916c782c..ee68f29f 100644 --- a/tests/e2e_tests/noisemap_cal_e2e.py +++ b/tests/e2e_tests/noisemap_cal_e2e.py @@ -425,5 +425,5 @@ def test_noisemap_calibration_from_l2a(tvacdata_path, e2eoutput_path): args = ap.parse_args(args_here) tvacdata_dir = args.tvacdata_dir outputdir = args.outputdir - test_noisemap_calibration_from_l1(tvacdata_dir, outputdir) + #test_noisemap_calibration_from_l1(tvacdata_dir, outputdir) test_noisemap_calibration_from_l2a(tvacdata_dir, outputdir) diff --git a/tests/e2e_tests/photon_count_e2e.py b/tests/e2e_tests/photon_count_e2e.py new file mode 100644 index 00000000..1bab4294 --- /dev/null +++ b/tests/e2e_tests/photon_count_e2e.py @@ -0,0 +1,152 @@ +# Photon Counting E2E Test Code + +import argparse +import os +import pytest +import numpy as np +import astropy.time as time +from astropy.io import fits + +import corgidrp +import corgidrp.data as data +import corgidrp.mocks as mocks +import corgidrp.walker as walker +import corgidrp.caldb as caldb +import corgidrp.detector as detector + +@pytest.mark.e2e +def test_expected_results_e2e(file_dir): + '''Results are as expected theoretically. Also runs raw frames through pre-processing pipeline.''' + np.random.seed(1234) + ill_dataset, dark_dataset, ill_mean, dark_mean = mocks.create_photon_countable_frames(Nbrights=160, Ndarks=161, cosmic_rate=1) + output_dir = os.path.join(file_dir, 'pc_sim_test_data') + output_ill_dir = os.path.join(output_dir, 'ill_frames') + output_dark_dir = os.path.join(output_dir, 'dark_frames') + if not os.path.exists(output_dir): + os.mkdir(output_dir) + # empty out directory of any previous files + for f in os.listdir(output_dir): + if os.path.isdir(os.path.join(output_dir,f)): + continue + os.remove(os.path.join(output_dir,f)) + if not os.path.exists(output_ill_dir): + os.mkdir(output_ill_dir) + if not os.path.exists(output_dark_dir): + os.mkdir(output_dark_dir) + # empty out directory of any previous files + for f in os.listdir(output_ill_dir): + os.remove(os.path.join(output_ill_dir,f)) + for f in os.listdir(output_dark_dir): + os.remove(os.path.join(output_dark_dir,f)) + ill_dataset.save(output_ill_dir, ['pc_frame_ill_{0}.fits'.format(i) for i in range(len(ill_dataset))]) + dark_dataset.save(output_dark_dir, ['pc_frame_dark_{0}.fits'.format(i) for i in range(len(dark_dataset))]) + l1_data_ill_filelist = [] + l1_data_dark_filelist = [] + for f in os.listdir(output_ill_dir): + l1_data_ill_filelist.append(os.path.join(output_ill_dir, f)) + for f in os.listdir(output_dark_dir): + l1_data_dark_filelist.append(os.path.join(output_dark_dir, f)) + + this_caldb = caldb.CalDB() # connection to cal DB + # remove other KGain calibrations that may exist in case they don't have the added header keywords + for i in range(len(this_caldb._db['Type'])): + if this_caldb._db['Type'][i] == 'KGain': + this_caldb._db = this_caldb._db.drop(i) + elif this_caldb._db['Type'][i] == 'Dark': + this_caldb._db = this_caldb._db.drop(i) + this_caldb.save() + + # KGain + kgain_val = 7 # default value used in mocks.create_photon_countable_frames() + pri_hdr, ext_hdr = mocks.create_default_headers() + ext_hdr["DRPCTIME"] = time.Time.now().isot + ext_hdr['DRPVERSN'] = corgidrp.__version__ + mock_input_dataset = data.Dataset(l1_data_ill_filelist) + kgain = data.KGain(np.array([[kgain_val]]), pri_hdr=pri_hdr, ext_hdr=ext_hdr, + input_dataset=mock_input_dataset) + # add in keywords that didn't make it into mock_kgain.fits, using values used in mocks.create_photon_countable_frames() + kgain.ext_hdr['RN'] = 100 + kgain.ext_hdr['RN_ERR'] = 0 + kgain.save(filedir=output_dir, filename="mock_kgain.fits") + this_caldb.create_entry(kgain) + + # NoiseMap + noise_map_dat = np.zeros((3, detector.detector_areas['SCI']['frame_rows'], detector.detector_areas['SCI']['frame_cols'])) + noise_map_noise = np.zeros([1,] + list(noise_map_dat.shape)) + noise_map_dq = np.zeros(noise_map_dat.shape, dtype=int) + err_hdr = fits.Header() + err_hdr['BUNIT'] = 'detected electrons' + ext_hdr['B_O'] = 0 + ext_hdr['B_O_ERR'] = 0 + noise_map = data.DetectorNoiseMaps(noise_map_dat, pri_hdr=pri_hdr, ext_hdr=ext_hdr, + input_dataset=mock_input_dataset, err=noise_map_noise, + dq = noise_map_dq, err_hdr=err_hdr) + noise_map.save(filedir=output_dir, filename="mock_detnoisemaps.fits") + this_caldb.create_entry(noise_map) + + corgidrp_dir = os.path.split(corgidrp.__path__[0])[0] + tests_dir = os.path.join(corgidrp_dir, 'tests') + # Nonlinearity calibration + ### Create a dummy non-linearity file #### + #Create a mock dataset because it is a required input when creating a NonLinearityCalibration + dummy_dataset = mocks.create_prescan_files() + # Make a non-linearity correction calibration file + input_non_linearity_filename = "nonlin_table_TVAC.txt" + input_non_linearity_path = os.path.join(tests_dir, "test_data", input_non_linearity_filename) + test_non_linearity_filename = input_non_linearity_filename.split(".")[0] + ".fits" + nonlin_fits_filepath = os.path.join(tests_dir, "test_data", test_non_linearity_filename) + tvac_nonlin_data = np.genfromtxt(input_non_linearity_path, delimiter=",") + + pri_hdr, ext_hdr = mocks.create_default_headers() + new_nonlinearity = data.NonLinearityCalibration(tvac_nonlin_data,pri_hdr=pri_hdr,ext_hdr=ext_hdr,input_dataset = dummy_dataset) + new_nonlinearity.filename = nonlin_fits_filepath + new_nonlinearity.pri_hdr = pri_hdr + new_nonlinearity.ext_hdr = ext_hdr + this_caldb.create_entry(new_nonlinearity) + # make PC dark + walker.walk_corgidrp(l1_data_dark_filelist, '', output_dir, template="l1_to_l2b_pc_dark.json") + for f in os.listdir(output_dir): + if f.endswith('_pc_dark.fits'): + pc_dark_filename = f + # add the Dark to the Caldb for processing the illuminated + pc_dark_file = os.path.join(output_dir, pc_dark_filename) + dark_entry = data.Dark(pc_dark_file) + this_caldb.create_entry(dark_entry) + # make PC illuminated, subtracting the PC dark + walker.walk_corgidrp(l1_data_ill_filelist, '', output_dir, template="l1_to_l2b_pc.json") + # get photon-counted frame + for f in os.listdir(output_dir): + if f.endswith('_pc.fits'): + pc_filename = f + pc_file = os.path.join(output_dir, pc_filename) + pc_frame = fits.getdata(pc_file) + pc_frame_err = fits.getdata(pc_file, 'ERR') + pc_dark_frame = fits.getdata(pc_dark_file) + pc_dark_frame_err = fits.getdata(pc_dark_file, 'ERR') + + # more frames gets a better agreement; agreement to 1% for ~160 darks and illuminated + assert np.isclose(np.nanmean(pc_frame), ill_mean - dark_mean, rtol=0.01) + assert np.isclose(np.nanmean(pc_dark_frame), dark_mean, rtol=0.01) + assert pc_frame_err.min() >= 0 + assert pc_dark_frame_err.min() >= 0 + + this_caldb.remove_entry(kgain) + this_caldb.remove_entry(noise_map) + this_caldb.remove_entry(new_nonlinearity) + + +if __name__ == "__main__": + # Use arguments to run the test. Users can then write their own scripts + # that call this script with the correct arguments and they do not need + # to edit the file. The arguments use the variables in this file as their + # defaults allowing the user to edit the file if that is their preferred + # workflow. + thisfile_dir = os.path.dirname(__file__) + outputdir = thisfile_dir + + ap = argparse.ArgumentParser(description="run the l1->l2b PC end-to-end test") + ap.add_argument("-o", "--outputdir", default=outputdir, + help="directory to write results to [%(default)s]") + args = ap.parse_args() + outputdir = args.outputdir + test_expected_results_e2e(outputdir) diff --git a/tests/test_badpixelmap.py b/tests/test_badpixelmap.py index 022cd05d..645b8c07 100644 --- a/tests/test_badpixelmap.py +++ b/tests/test_badpixelmap.py @@ -7,7 +7,7 @@ from corgidrp.bad_pixel_calibration import create_bad_pixel_map from corgidrp.darks import build_trad_dark - +np.random.seed(456) def test_badpixelmap(): ''' @@ -70,7 +70,7 @@ def test_badpixelmap(): flat_frame.data[i_col, i_row] = 0.3 ###### make the badpixel map (input the flat_dataset just as a dummy): - badpixelmap = create_bad_pixel_map(flat_dataset, dark_frame,flat_frame) + badpixelmap = create_bad_pixel_map(flat_dataset, dark_frame,flat_frame, dthresh=6) # Use np.unpackbits to unpack the bits - big endien badpixelmap_bits = np.unpackbits(badpixelmap.data[:, :, np.newaxis], axis=2) diff --git a/tests/test_build_synthesized_dark.py b/tests/test_build_synthesized_dark.py index 51834677..d1550f55 100644 --- a/tests/test_build_synthesized_dark.py +++ b/tests/test_build_synthesized_dark.py @@ -26,6 +26,7 @@ C = noise_maps.CIC_map D = noise_maps.DC_map # just need some a dataset for input +np.random.seed(456) dataset = create_dark_calib_files() # values used in create_dark_calib_files() g = 1 diff --git a/tests/test_caldb.py b/tests/test_caldb.py index e8a99455..ba6556f7 100644 --- a/tests/test_caldb.py +++ b/tests/test_caldb.py @@ -1,4 +1,5 @@ import os +import numpy as np import pytest import corgidrp import corgidrp.caldb as caldb @@ -16,6 +17,7 @@ testcaldb_filepath = os.path.join(calibdir, "test_caldb.csv") # make some fake test data to use +np.random.seed(456) dark_dataset = mocks.create_dark_calib_files() master_dark = data.Dark(dark_dataset[0].data, dark_dataset[0].pri_hdr, dark_dataset[0].ext_hdr, dark_dataset) # save master dark to disk to be loaded later @@ -163,4 +165,8 @@ def test_caldb_scan(): os.remove(testcaldb_filepath) if __name__ == "__main__": - test_get_calib() \ No newline at end of file + test_get_calib() + test_caldb_create_default() + test_caldb_custom_filepath() + test_caldb_insert_and_remove() + test_caldb_scan() \ No newline at end of file diff --git a/tests/test_calibrate_darks_lsq.py b/tests/test_calibrate_darks_lsq.py index 2e85a144..8b6d8fa3 100644 --- a/tests/test_calibrate_darks_lsq.py +++ b/tests/test_calibrate_darks_lsq.py @@ -11,6 +11,8 @@ from corgidrp.mocks import detector_areas_test as dat from corgidrp.data import DetectorNoiseMaps, DetectorParams, Dataset +# make test reproducible +np.random.seed(4567) # use default parameters detector_params = DetectorParams({}) # specified parameters simulated in simulated data from diff --git a/tests/test_combine.py b/tests/test_combine.py index e6ed9faf..28f9edcf 100644 --- a/tests/test_combine.py +++ b/tests/test_combine.py @@ -125,4 +125,6 @@ def test_median_combine_subexposures(): if __name__ == "__main__": - test_mean_combine_subexposures() \ No newline at end of file + test_mean_combine_subexposures() + test_mean_combine_subexposures_with_bad() + test_median_combine_subexposures() \ No newline at end of file diff --git a/tests/test_dark_sub.py b/tests/test_dark_sub.py index 0d09546c..380c7d84 100644 --- a/tests/test_dark_sub.py +++ b/tests/test_dark_sub.py @@ -10,6 +10,8 @@ import corgidrp.l2a_to_l2b as l2a_to_l2b from corgidrp.darks import build_trad_dark +np.random.seed(456) + old_err_tracking = corgidrp.track_individual_errors # use default parameters detector_params = data.DetectorParams({}) @@ -81,6 +83,7 @@ def test_dark_sub(): # load in the dark dark_filepath = os.path.join(calibdir, dark_filename) new_dark = data.Dark(dark_filepath) + assert new_dark.ext_hdr['PC_STAT'] == 'analog master dark' # check the dark can be pickled (for CTC operations) pickled = pickle.dumps(new_dark) @@ -91,6 +94,11 @@ def test_dark_sub(): darkest_dataset = l2a_to_l2b.dark_subtraction(dark_dataset, new_dark, outputdir=calibdir) assert(dark_filename in str(darkest_dataset[0].ext_hdr["HISTORY"])) + # PC dark cannot be used for this step function + new_dark.ext_hdr['PC_STAT'] = 'photon-counted master dark' + with pytest.raises(Exception): + l2a_to_l2b.dark_subtraction(dark_dataset, new_dark, outputdir=calibdir) + # check the level of the dataset is now approximately 0, leaving off telemetry row assert np.mean(darkest_dataset.all_data[:,:-1,:]) == pytest.approx(0, abs=1e-2) diff --git a/tests/test_data_levels.py b/tests/test_data_levels.py index 588070c9..f0de7b03 100644 --- a/tests/test_data_levels.py +++ b/tests/test_data_levels.py @@ -2,12 +2,15 @@ Test the data level specification """ import pytest +import numpy as np import corgidrp.mocks as mocks import corgidrp.l1_to_l2a as l1_to_l2a import corgidrp.l2a_to_l2b as l2a_to_l2b import corgidrp.l2b_to_l3 as l2b_to_l3 import corgidrp.l3_to_l4 as l3_to_l4 +np.random.seed(456) + def test_l1_to_l4(): """ Tests a successful upgrade of L1 to L4 in bookkeeping only diff --git a/tests/test_dataset.py b/tests/test_dataset.py index 6124dc2c..8316dd16 100644 --- a/tests/test_dataset.py +++ b/tests/test_dataset.py @@ -8,6 +8,8 @@ from corgidrp.data import Image, Dataset from corgidrp.mocks import create_default_headers, create_dark_calib_files +np.random.seed(123) + data = np.ones([1024,1024]) * 2 err = np.zeros([1024,1024]) err1 = np.ones([1024,1024]) diff --git a/tests/test_err_dq.py b/tests/test_err_dq.py index 31c29828..7d2dd4a4 100644 --- a/tests/test_err_dq.py +++ b/tests/test_err_dq.py @@ -10,6 +10,7 @@ import corgidrp.caldb as caldb from corgidrp.darks import build_trad_dark +np.random.seed(123) data = np.ones([1024,1024]) * 2 err = np.zeros([1024,1024]) diff --git a/tests/test_frame_selection.py b/tests/test_frame_selection.py index b6f41b15..89e49edd 100644 --- a/tests/test_frame_selection.py +++ b/tests/test_frame_selection.py @@ -1,7 +1,9 @@ import pytest +import numpy as np import corgidrp.mocks as mocks from corgidrp.l2a_to_l2b import frame_select +np.random.seed(123) def test_no_selection(): """ diff --git a/tests/test_photon_counting.py b/tests/test_photon_counting.py index 69f9613a..2fddcadb 100644 --- a/tests/test_photon_counting.py +++ b/tests/test_photon_counting.py @@ -11,97 +11,6 @@ import corgidrp.detector as detector from corgidrp.photon_counting import get_pc_mean, photon_count, PhotonCountException -detector_params = data.DetectorParams({}) - -def test_expected_results(): - '''Results are as expected theoretically. Also runs raw frames through pre-processing pipeline.''' - dataset, ill_mean, dark_mean = mocks.create_photon_countable_frames(Nbrights=5, Ndarks=6, cosmic_rate=0) - thisfile_dir = os.path.dirname(__file__) # this file's folder - outputdir = os.path.join(thisfile_dir, 'simdata', 'pc_test_data') - if not os.path.exists(outputdir): - os.mkdir(outputdir) - # empty out directory of any previous files - for f in os.listdir(outputdir): - os.remove(os.path.join(outputdir,f)) - dataset.save(outputdir, ['pc_frame_{0}.fits'.format(i) for i in range(len(dataset))]) - l1_data_filelist = [] - for f in os.listdir(outputdir): - l1_data_filelist.append(os.path.join(outputdir, f)) - - this_caldb = caldb.CalDB() # connection to cal DB - # remove other KGain calibrations that may exist in case they don't have the added header keywords - for i in range(len(this_caldb._db['Type'])): - if this_caldb._db['Type'][i] == 'KGain': - this_caldb._db = this_caldb._db.drop(i) - this_caldb.save() - - # KGain - kgain_val = 7 # default value used in mocks.create_photon_countable_frames() - pri_hdr, ext_hdr = mocks.create_default_headers() - ext_hdr["DRPCTIME"] = time.Time.now().isot - ext_hdr['DRPVERSN'] = corgidrp.__version__ - mock_input_dataset = data.Dataset(l1_data_filelist) - kgain = data.KGain(np.array([[kgain_val]]), pri_hdr=pri_hdr, ext_hdr=ext_hdr, - input_dataset=mock_input_dataset) - # add in keywords that didn't make it into mock_kgain.fits, using values used in mocks.create_photon_countable_frames() - kgain.ext_hdr['RN'] = 100 - kgain.ext_hdr['RN_ERR'] = 0 - kgain.save(filedir=outputdir, filename="mock_kgain.fits") - this_caldb.create_entry(kgain) - - # NoiseMap - noise_map_dat = np.zeros((3, detector.detector_areas['SCI']['frame_rows'], detector.detector_areas['SCI']['frame_cols'])) - noise_map_noise = np.zeros([1,] + list(noise_map_dat.shape)) - noise_map_dq = np.zeros(noise_map_dat.shape, dtype=int) - err_hdr = fits.Header() - err_hdr['BUNIT'] = 'detected electrons' - ext_hdr['B_O'] = 0 - ext_hdr['B_O_ERR'] = 0 - noise_map = data.DetectorNoiseMaps(noise_map_dat, pri_hdr=pri_hdr, ext_hdr=ext_hdr, - input_dataset=mock_input_dataset, err=noise_map_noise, - dq = noise_map_dq, err_hdr=err_hdr) - noise_map.save(filedir=outputdir, filename="mock_detnoisemaps.fits") - this_caldb.create_entry(noise_map) - - # Nonlinearity calibration - ### Create a dummy non-linearity file #### - #Create a mock dataset because it is a required input when creating a NonLinearityCalibration - dummy_dataset = mocks.create_prescan_files() - # Make a non-linearity correction calibration file - input_non_linearity_filename = "nonlin_table_TVAC.txt" - input_non_linearity_path = os.path.join(os.path.dirname(__file__), "test_data", input_non_linearity_filename) - test_non_linearity_filename = input_non_linearity_filename.split(".")[0] + ".fits" - nonlin_fits_filepath = os.path.join(os.path.dirname(__file__), "test_data", test_non_linearity_filename) - tvac_nonlin_data = np.genfromtxt(input_non_linearity_path, delimiter=",") - - pri_hdr, ext_hdr = mocks.create_default_headers() - new_nonlinearity = data.NonLinearityCalibration(tvac_nonlin_data,pri_hdr=pri_hdr,ext_hdr=ext_hdr,input_dataset = dummy_dataset) - new_nonlinearity.filename = nonlin_fits_filepath - # fake the headers because this frame doesn't have the proper headers - new_nonlinearity.pri_hdr = pri_hdr - new_nonlinearity.ext_hdr = ext_hdr - this_caldb.create_entry(new_nonlinearity) - - walker.walk_corgidrp(l1_data_filelist, '', outputdir, template="l1_to_l2b_pc.json") - # get photon-counted frame - for f in os.listdir(outputdir): - if f.endswith('_pc.fits'): - pc_filename = f - pc_file = os.path.join(outputdir, pc_filename) - pc_frame = fits.getdata(pc_file) - pc_frame_err = fits.getdata(pc_file, 'ERR') - pc_ext_hdr = fits.getheader(pc_file, 1) - # more frames (which would take longer to run) would give an even better agreement than the 25% agreement below - assert np.isclose(pc_frame.mean(), ill_mean - dark_mean, rtol=0.25) - assert 'niter=2' in pc_ext_hdr["HISTORY"][-1] - assert 'T_factor=5' in pc_ext_hdr["HISTORY"][-1] - assert pc_frame_err.min() >= 0 - assert pc_frame_err[0][3,3] >= 0 - - this_caldb.remove_entry(kgain) - this_caldb.remove_entry(noise_map) - this_caldb.remove_entry(new_nonlinearity) - def test_negative(): """Values at or below the threshold are photon-counted as 0, including negative values.""" @@ -110,36 +19,101 @@ def test_negative(): assert (pc_fr == np.array([0,0,0])).all() -def test_masking_increases_err(): +def test_pc(): '''Test that a pixel that is masked heavily (reducing the number of usable frames for that pixel) has a bigger err than the average pixel. And if masking is all the way through for a certain pixel, that pixel will - be flagged in the DQ. Also, make sure there is failure when niter<1. Also, very good agreement with theoretically expected value since we use more frames.''' + be flagged in the DQ. Also, make sure there is failure when niter<1 and threshold<0, and test other built-in catches. + Test safemode. Test mask_filepath to ensure failure for bright pixels outside region of interest when no mask used + and success when mask used to mask those bright pixels.''' # exposure time too long to get reasonable photon-counted result (after photometric correction) - dataset_err, ill_mean, dark_mean = mocks.create_photon_countable_frames(Nbrights=60, Ndarks=60, cosmic_rate=0, full_frame=False) + np.random.seed(555) + dataset_err, dark_dataset_err, ill_mean, dark_mean = mocks.create_photon_countable_frames(Nbrights=160, Ndarks=160, cosmic_rate=0, full_frame=False) # instead of running through walker, just do the pre-processing steps simply - # using kgain=7 and bias=2000 and read noise = 100 and QE=0.9 (quantum efficiency), from mocks.create_photon_countable_frames() + # using EM gain=5000 and kgain=7 and bias=20000 and read noise = 100 and QE=0.9 (quantum efficiency), from mocks.create_photon_countable_frames() for f in dataset_err.frames: - f.data = f.data.astype(float)*7 - 2000. - + f.data = f.data.astype(float)*7 - 20000. + for f in dark_dataset_err.frames: + f.data = f.data.astype(float)*7 - 20000. + dataset_err.all_data = dataset_err.all_data.astype(float)*7 - 20000. + dark_dataset_err.all_data = dark_dataset_err.all_data.astype(float)*7 - 20000. # masked for half of the bright and half of the dark frames - dataset_err.all_dq[30:90,3,3] = 1 + dataset_err.all_dq[80:,3,3] = 1 + dark_dataset_err.all_dq[80:,3,3] = 1 dataset_err.all_dq[:, 2, 2] = 1 #masked all the way through for (2,2) for f in dataset_err.frames: f.ext_hdr['RN'] = 100 f.ext_hdr['KGAIN'] = 7 - pc_dataset_err = get_pc_mean(dataset_err, detector_params) - assert np.isclose(pc_dataset_err.all_data.mean(), ill_mean - dark_mean, rtol=0.001) + for f in dark_dataset_err.frames: + f.ext_hdr['RN'] = 100 + f.ext_hdr['KGAIN'] = 7 + # process the frames to make PC dark + pc_dark = get_pc_mean(dark_dataset_err) + assert pc_dark.ext_hdr['PC_STAT'] == 'photon-counted master dark' + # now process illuminated frames and subtract the PC dark + pc_dataset_err = get_pc_mean(dataset_err, pc_master_dark=pc_dark) + assert pc_dataset_err.frames[0].filename[-7:] == 'pc.fits' + assert pc_dataset_err.frames[0].filepath[-7:] == 'pc.fits' + history = '' + for line in pc_dataset_err.frames[0].ext_hdr["HISTORY"]: + history += line + assert "Dark-subtracted: yes" in history + assert np.isclose(pc_dataset_err.all_data.mean(), ill_mean - dark_mean, rtol=0.01) # the DQ for that pixel should be 1 assert pc_dataset_err[0].dq[2,2] == 1 # err for (3,3) is above the 95th percentile of error: assert pc_dataset_err[0].err[0][3,3]>np.nanpercentile(pc_dataset_err[0].err,95) + # also when niter<1, exception with pytest.raises(PhotonCountException): - get_pc_mean(dataset_err, detector_params, niter=0) + get_pc_mean(dataset_err, niter=0) + # when thresh<0, exception + with pytest.raises(PhotonCountException): + get_pc_mean(dataset_err, T_factor=-1, niter=2) + # when thresh>=em_gain, exception + with pytest.raises(Exception): + get_pc_mean(dataset_err, T_factor=50, niter=2) + # can't provide Dark class instance while "VISTYPE" of input dataset is 'DARK' + with pytest.raises(PhotonCountException): + get_pc_mean(dark_dataset_err, pc_master_dark=pc_dark) + # must have same 'CMDGAIN' and other header values throughout input dataset + dark_dataset_err.frames[0].ext_hdr['CMDGAIN'] = 4999 + with pytest.raises(PhotonCountException): + get_pc_mean(dark_dataset_err, pc_master_dark=pc_dark) + + # to trigger pc_ecount_max error + for f in dataset_err.frames[:-2]: #all but 2 of the frames + f.data[22:40,23:49] = np.abs(f.data.astype(float)[22:40,23:49]*3000) + dataset_err.all_data[:-2,22:40,23:49] = np.abs(dataset_err.all_data.astype(float)[:-2,22:40,23:49]*3000) + with pytest.raises(Exception): + get_pc_mean(dataset_err, pc_ecount_max=1) + # now use mask to exclude 22,23 from region of interest: + mask = np.zeros_like(dataset_err.all_data[0]) + mask[22:40,23:49] = 1 #exclude + thisfile_dir = os.path.dirname(__file__) + mask_filepath = os.path.join(thisfile_dir, 'test_data', 'pc_mask.fits') + hdr = fits.Header() + prim = fits.PrimaryHDU(header=hdr) + img = fits.ImageHDU(mask) + hdul = fits.HDUList([prim,img]) + hdul.writeto(mask_filepath, overwrite=True) + pc_dataset_err = get_pc_mean(dataset_err, pc_ecount_max=1, mask_filepath=mask_filepath) + copy_pc = pc_dataset_err.all_data.copy() + # mask out all the irrelevant pixels: + copy_pc[0,22:40,23:49] = np.nan + # didn't dark-subtract this time: + assert pc_dataset_err.frames[0].filename.endswith('pc_no_ds.fits') + assert pc_dataset_err.frames[0].filepath.endswith('pc_no_ds.fits') + history = '' + for line in pc_dataset_err.frames[0].ext_hdr["HISTORY"]: + history += line + assert "Dark-subtracted: no" in history + assert np.isclose(np.nanmean(copy_pc), ill_mean, rtol=0.01) + # test safemode=False: should get warning instead of exception + with pytest.warns(UserWarning): + get_pc_mean(dataset_err, safemode=False) if __name__ == '__main__': - test_masking_increases_err() + test_pc() test_negative() - test_expected_results() \ No newline at end of file diff --git a/tests/test_walker.py b/tests/test_walker.py index 53ba5045..daa19d06 100644 --- a/tests/test_walker.py +++ b/tests/test_walker.py @@ -14,6 +14,7 @@ import corgidrp.walker as walker import corgidrp.detector as detector +np.random.seed(456) def test_autoreducing(): """ From e3e27c45437ea3f347c2c5b1505eca6e40e7994f Mon Sep 17 00:00:00 2001 From: Kevin Ludwick Date: Thu, 26 Dec 2024 10:33:09 -0600 Subject: [PATCH 13/26] linting fix --- tests/e2e_tests/photon_count_e2e.py | 1 - 1 file changed, 1 deletion(-) diff --git a/tests/e2e_tests/photon_count_e2e.py b/tests/e2e_tests/photon_count_e2e.py index 1bab4294..fe26d4d6 100644 --- a/tests/e2e_tests/photon_count_e2e.py +++ b/tests/e2e_tests/photon_count_e2e.py @@ -16,7 +16,6 @@ @pytest.mark.e2e def test_expected_results_e2e(file_dir): - '''Results are as expected theoretically. Also runs raw frames through pre-processing pipeline.''' np.random.seed(1234) ill_dataset, dark_dataset, ill_mean, dark_mean = mocks.create_photon_countable_frames(Nbrights=160, Ndarks=161, cosmic_rate=1) output_dir = os.path.join(file_dir, 'pc_sim_test_data') From 326e8807ba6a12b1c4694169e7b5415fe5340de6 Mon Sep 17 00:00:00 2001 From: Kevin Ludwick Date: Thu, 26 Dec 2024 10:51:51 -0600 Subject: [PATCH 14/26] fixed failed tests --- tests/test_flat_div.py | 1 + tests/test_kgain_cal.py | 2 ++ 2 files changed, 3 insertions(+) diff --git a/tests/test_flat_div.py b/tests/test_flat_div.py index 48b6c80e..5a45f369 100644 --- a/tests/test_flat_div.py +++ b/tests/test_flat_div.py @@ -10,6 +10,7 @@ old_err_tracking = corgidrp.track_individual_errors +np.random.seed(9292) def test_flat_div(): """ Generate mock input data and pass into flat division function diff --git a/tests/test_kgain_cal.py b/tests/test_kgain_cal.py index 4783724b..67614052 100644 --- a/tests/test_kgain_cal.py +++ b/tests/test_kgain_cal.py @@ -19,6 +19,8 @@ from corgidrp.mocks import (create_default_headers, make_fluxmap_image, nonlin_coefs) from corgidrp.calibrate_kgain import (calibrate_kgain, CalKgainException, kgain_params) + +np.random.seed(8585) ######################## function definitions ############################### def count_contiguous_repeats(arr): From 44ccc5dbc8364a36bd08ed8dbf46c9e2fa7f1793 Mon Sep 17 00:00:00 2001 From: Kevin Ludwick Date: Thu, 26 Dec 2024 11:03:28 -0600 Subject: [PATCH 15/26] fix --- tests/test_flat_div.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/test_flat_div.py b/tests/test_flat_div.py index 5a45f369..9b94c9dc 100644 --- a/tests/test_flat_div.py +++ b/tests/test_flat_div.py @@ -62,7 +62,7 @@ def test_flat_div(): # perform checks after the flat divison assert(flat_filename in str(flatdivided_dataset[0].ext_hdr["HISTORY"])) # check the level of the dataset is now approximately 100 - assert np.mean(flatdivided_dataset.all_data) == pytest.approx(150, abs=1e-2) + assert np.mean(flatdivided_dataset.all_data) == pytest.approx(150, abs=2e-2) # check the propagated errors assert flatdivided_dataset[0].err_hdr["Layer_2"] == "FlatField_error" print("mean of all simulated data",np.mean(simflat_dataset.all_data)) From 3c17a299d2ccc2f0c8bc355a02eef971c34b88df Mon Sep 17 00:00:00 2001 From: "Maxwell A. Millar-Blanchaer" Date: Wed, 8 Jan 2025 09:55:04 -0800 Subject: [PATCH 16/26] updated the PR template --- .github/pull_request_template.md | 1 + 1 file changed, 1 insertion(+) diff --git a/.github/pull_request_template.md b/.github/pull_request_template.md index 09cdcbf7..d30cfbcd 100644 --- a/.github/pull_request_template.md +++ b/.github/pull_request_template.md @@ -16,3 +16,4 @@ Please delete options that are not relevant (and this line). - [ ] I have linted my code - [ ] I have verified that all unit tests pass in a clean environment and added new unit tests, as needed - [ ] I have verified that all docstrings are properly formatted and added new documentation, as needed +- [ ] I have filled out the Unit Test Definition Table on confluence and gotten my proposed unit tests approved by a DRP lead, a PS Team member and the VAP lead. From 062e98150ad11780c44d762bb965f619f9049628 Mon Sep 17 00:00:00 2001 From: "Maxwell A. Millar-Blanchaer" Date: Wed, 8 Jan 2025 09:55:31 -0800 Subject: [PATCH 17/26] trivial syntax fix --- .github/pull_request_template.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/pull_request_template.md b/.github/pull_request_template.md index d30cfbcd..22d47c85 100644 --- a/.github/pull_request_template.md +++ b/.github/pull_request_template.md @@ -16,4 +16,4 @@ Please delete options that are not relevant (and this line). - [ ] I have linted my code - [ ] I have verified that all unit tests pass in a clean environment and added new unit tests, as needed - [ ] I have verified that all docstrings are properly formatted and added new documentation, as needed -- [ ] I have filled out the Unit Test Definition Table on confluence and gotten my proposed unit tests approved by a DRP lead, a PS Team member and the VAP lead. +- [ ] I have filled out the Unit Test Definition Table on confluence and gotten my proposed unit tests approved by a DRP lead, a PS Team member and the VAP lead From 60444f5c19d2cbc40cb6b31afa5c00e156a27636 Mon Sep 17 00:00:00 2001 From: "Maxwell A. Millar-Blanchaer" Date: Wed, 8 Jan 2025 11:09:46 -0800 Subject: [PATCH 18/26] language update --- .github/pull_request_template.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/pull_request_template.md b/.github/pull_request_template.md index 22d47c85..8f85a7ac 100644 --- a/.github/pull_request_template.md +++ b/.github/pull_request_template.md @@ -16,4 +16,4 @@ Please delete options that are not relevant (and this line). - [ ] I have linted my code - [ ] I have verified that all unit tests pass in a clean environment and added new unit tests, as needed - [ ] I have verified that all docstrings are properly formatted and added new documentation, as needed -- [ ] I have filled out the Unit Test Definition Table on confluence and gotten my proposed unit tests approved by a DRP lead, a PS Team member and the VAP lead +- [ ] I have filled out the Unit Test Definition Table on confluence, as needed \ No newline at end of file From 7156e496d7df8cef1360b7566587e49326ccb899 Mon Sep 17 00:00:00 2001 From: Kevin Ludwick Date: Thu, 23 Jan 2025 15:45:00 -0600 Subject: [PATCH 19/26] update --- corgidrp/calibrate_kgain.py | 20 +++--- corgidrp/calibrate_nonlin.py | 32 ++++----- corgidrp/combine.py | 2 +- corgidrp/darks.py | 27 ++++---- corgidrp/data.py | 7 ++ corgidrp/mocks.py | 21 +++++- corgidrp/photon_counting.py | 26 +++++--- corgidrp/recipe_templates/l1_to_l2b_pc.json | 27 ++++++++ .../recipe_templates/l1_to_l2b_pc_dark.json | 5 +- corgidrp/recipe_templates/l2a_to_l2b_pc.json | 27 ++++++++ .../recipe_templates/l2a_to_l2b_pc_dark.json | 5 +- tests/e2e_tests/astrom_e2e.py | 2 +- tests/e2e_tests/noisemap_cal_e2e.py | 14 +++- tests/e2e_tests/photon_count_e2e.py | 48 +++++++++++--- tests/test_photon_counting.py | 65 ++++++++++++++++--- 15 files changed, 254 insertions(+), 74 deletions(-) diff --git a/corgidrp/calibrate_kgain.py b/corgidrp/calibrate_kgain.py index 3fefc0f4..6bf90094 100644 --- a/corgidrp/calibrate_kgain.py +++ b/corgidrp/calibrate_kgain.py @@ -34,25 +34,25 @@ def check_kgain_params( ): """ Checks integrity of kgain parameters in the dictionary kgain_params. """ if 'offset_colroi1' not in kgain_params: - raise ValueError('Missing parameter.') + raise ValueError('Missing parameter: offset_colroi1.') if 'offset_colroi2' not in kgain_params: - raise ValueError('Missing parameter.') + raise ValueError('Missing parameter: offset_colroi2.') if 'rowroi1' not in kgain_params: - raise ValueError('Missing parameter.') + raise ValueError('Missing parameter: rowroi1.') if 'rowroi2' not in kgain_params: - raise ValueError('Missing parameter.') + raise ValueError('Missing parameter: rowroi2.') if 'colroi1' not in kgain_params: - raise ValueError('Missing parameter.') + raise ValueError('Missing parameter: colroi1.') if 'colroi2' not in kgain_params: - raise ValueError('Missing parameter.') + raise ValueError('Missing parameter: colroi2.') if 'rn_bins1' not in kgain_params: - raise ValueError('Missing parameter.') + raise ValueError('Missing parameter: rn_bins1.') if 'rn_bins2' not in kgain_params: - raise ValueError('Missing parameter.') + raise ValueError('Missing parameter: rn_bins2.') if 'max_DN_val' not in kgain_params: - raise ValueError('Missing parameter.') + raise ValueError('Missing parameter: max_DN_val.') if 'signal_bins_N' not in kgain_params: - raise ValueError('Missing parameter.') + raise ValueError('Missing parameter: signal_bins_N.') if not isinstance(kgain_params['offset_colroi1'], (float, int)): raise TypeError('offset_colroi1 is not a number') diff --git a/corgidrp/calibrate_nonlin.py b/corgidrp/calibrate_nonlin.py index a762fab4..fb41f4a5 100644 --- a/corgidrp/calibrate_nonlin.py +++ b/corgidrp/calibrate_nonlin.py @@ -48,37 +48,37 @@ def check_nonlin_params( ): """ Checks integrity of kgain parameters in the dictionary nonlin_params. """ if 'rowroi1' not in nonlin_params: - raise ValueError('Missing parameter') + raise ValueError('Missing parameter: rowroi1.') if 'rowroi2' not in nonlin_params: - raise ValueError('Missing parameter') + raise ValueError('Missing parameter: rowroi2.') if 'colroi1' not in nonlin_params: - raise ValueError('Missing parameter') + raise ValueError('Missing parameter: colroi1.') if 'colroi2' not in nonlin_params: - raise ValueError('Missing parameter') + raise ValueError('Missing parameter: colroi2.') if 'rowback11' not in nonlin_params: - raise ValueError('Missing parameter') + raise ValueError('Missing parameter: rowback11.') if 'rowback12' not in nonlin_params: - raise ValueError('Missing parameter') + raise ValueError('Missing parameter: rowback12.') if 'rowback21' not in nonlin_params: - raise ValueError('Missing parameter') + raise ValueError('Missing parameter: rowback21.') if 'rowback22' not in nonlin_params: - raise ValueError('Missing parameter') + raise ValueError('Missing parameter: rowback22.') if 'colback11' not in nonlin_params: - raise ValueError('Missing parameter') + raise ValueError('Missing parameter: colback11.') if 'colback12' not in nonlin_params: - raise ValueError('Missing parameter') + raise ValueError('Missing parameter: colback12.') if 'colback21' not in nonlin_params: - raise ValueError('Missing parameter') + raise ValueError('Missing parameter: colback21.') if 'colback22' not in nonlin_params: - raise ValueError('Missing parameter') + raise ValueError('Missing parameter: colback22.') if 'min_exp_time' not in nonlin_params: - raise ValueError('Missing parameter') + raise ValueError('Missing parameter: min_exp_time.') if 'num_bins' not in nonlin_params: - raise ValueError('Missing parameter') + raise ValueError('Missing parameter: num_bins.') if 'min_bin' not in nonlin_params: - raise ValueError('Missing parameter') + raise ValueError('Missing parameter: min_bin.') if 'min_mask_factor' not in nonlin_params: - raise ValueError('Missing parameter') + raise ValueError('Missing parameter: min_mask_factor.') if not isinstance(nonlin_params['rowroi1'], (float, int)): raise TypeError('rowroi1 is not a number') diff --git a/corgidrp/combine.py b/corgidrp/combine.py index ca49391a..0fc914fa 100644 --- a/corgidrp/combine.py +++ b/corgidrp/combine.py @@ -43,7 +43,7 @@ def combine_images(data_subset, err_subset, dq_subset, collapse, num_frames_scal # dq collpase: keep all flags on dq_collapse = np.bitwise_or.reduce(dq_subset, axis=0) - # except those pixels that have not been replaced + # except for those pixels that have been replaced with good values dq_collapse[np.where((dq_collapse > 0) & (~np.isnan(data_collapse)))] = 0 return data_collapse, err_collapse, dq_collapse diff --git a/corgidrp/darks.py b/corgidrp/darks.py index 86d5eae8..0abec4aa 100644 --- a/corgidrp/darks.py +++ b/corgidrp/darks.py @@ -1,7 +1,6 @@ import numpy as np import warnings from astropy.io import fits -import os from corgidrp.detector import slice_section, imaging_slice, imaging_area_geom, detector_areas import corgidrp.check as check @@ -134,12 +133,13 @@ def build_trad_dark(dataset, detector_params, detector_regions=None, full_frame= - have had masks made for cosmic rays - have been corrected for nonlinearity - have been converted from DN to e- - - have been desmeared if desmearing is appropriate. Under normal - circumstances, darks should not be desmeared. The only time desmearing - would be useful is in the unexpected case that, for example, - dark current is so high that it stands far above other noise that is - not smeared upon readout, such as clock-induced charge, - fixed-pattern noise, and read noise. + - have NOT been desmeared. Darks should not be desmeared. The only component + of dark frames that would be subject to a smearing effect is dark current + since it linearly increases with time, so the extra row read time affects + the dark current per pixel. However, illuminated images + would also contain this smeared dark current, so dark subtraction should + remove this smeared dark current (and then desmearing may be applied to the + processed image if appropriate). Also, add_photon_noise() should NOT have been applied to the frames in dataset. And note that creation of the @@ -269,12 +269,13 @@ def calibrate_darks_lsq(dataset, detector_params, detector_regions=None): - have had masks made for cosmic rays - have been corrected for nonlinearity - have been converted from DN to e- - - have been desmeared if desmearing is appropriate. Under normal - circumstances, darks should not be desmeared. The only time desmearing - would be useful is in the unexpected case that, for example, - dark current is so high that it stands far above other noise that is - not smeared upon readout, such as clock-induced charge, - fixed-pattern noise, and read noise. + - have NOT been desmeared. Darks should not be desmeared. The only component + of dark frames that would be subject to a smearing effect is dark current + since it linearly increases with time, so the extra row read time affects + the dark current per pixel. However, illuminated images + would also contain this smeared dark current, so dark subtraction should + remove this smeared dark current (and then desmearing may be applied to the + processed image if appropriate). Also, add_photon_noise() should NOT have been applied to the frames in dataset. And note that creation of the diff --git a/corgidrp/data.py b/corgidrp/data.py index f6198bfc..e04f1b51 100644 --- a/corgidrp/data.py +++ b/corgidrp/data.py @@ -863,6 +863,13 @@ def __init__(self, data_or_filepath, err = None, ptc = None, pri_hdr=None, ext_h # run the image class contructor super().__init__(data_or_filepath, err=err, pri_hdr=pri_hdr, ext_hdr=ext_hdr, err_hdr=err_hdr) + # initialize these headers that have been recently added so that older calib files still contain this keyword when initialized and allow for tests that don't require + # these values to run smoothly; if these values are actually required for + # a particular process, the user would be alerted since these values below would result in an error as they aren't numerical + if 'RN' not in self.ext_hdr: + self.ext_hdr['RN'] = '' + if 'RN_ERR' not in self.ext_hdr: + self.ext_hdr['RN_ERR'] = '' # File format checks if self.data.shape != (1,1): raise ValueError('The KGain calibration data should be just one float value') diff --git a/corgidrp/mocks.py b/corgidrp/mocks.py index b8b6c249..59f7ff12 100644 --- a/corgidrp/mocks.py +++ b/corgidrp/mocks.py @@ -17,6 +17,7 @@ import corgidrp.detector as detector from corgidrp.detector import imaging_area_geom, unpack_geom from corgidrp.pump_trap_calibration import (P1, P1_P1, P1_P2, P2, P2_P2, P3, P2_P3, P3_P3, tau_temp) +from corgidrp.data import DetectorParams from emccd_detect.emccd_detect import EMCCDDetect from emccd_detect.util.read_metadata_wrapper import MetadataWrapper @@ -1863,7 +1864,7 @@ def add_defect(sch_imgs, prob, ori, temp): else: hdul.writeto(filename, overwrite = True) -def create_photon_countable_frames(Nbrights=30, Ndarks=40, EMgain=5000, kgain=7, exptime=0.05, cosmic_rate=0, full_frame=True): +def create_photon_countable_frames(Nbrights=30, Ndarks=40, EMgain=5000, kgain=7, exptime=0.05, cosmic_rate=0, full_frame=True, smear=True, flux=1): '''This creates mock L1 Dataset containing frames with large gain and short exposure time, illuminated and dark frames. Used for unit tests for photon counting. @@ -1875,6 +1876,8 @@ def create_photon_countable_frames(Nbrights=30, Ndarks=40, EMgain=5000, kgain=7, exptime (float): exposure time (in s) cosmic_rate: (float) simulated cosmic rays incidence, hits/cm^2/s full_frame: (bool) If True, simulated frames are SCI full frames. If False, 50x50 images are simulated. Defaults to True. + smear: (bool) If True, smear is simulated. Defaults to True. + flux: (float) Number of photons/s per pixel desired. Defaults to 1. Returns: ill_dataset (corgidrp.data.Dataset): Dataset containing the illuminated frames @@ -1883,7 +1886,7 @@ def create_photon_countable_frames(Nbrights=30, Ndarks=40, EMgain=5000, kgain=7, dark_mean (float): mean electron count value simulated in the dark frames ''' pix_row = 1024 #number of rows and number of columns - fluxmap = np.ones((pix_row,pix_row)) #photon flux map, photons/s + fluxmap = flux*np.ones((pix_row,pix_row)) #photon flux map, photons/s emccd = EMCCDDetect( em_gain=EMgain, @@ -1919,6 +1922,20 @@ def create_photon_countable_frames(Nbrights=30, Ndarks=40, EMgain=5000, kgain=7, ill_mean = avg_ph_flux*emccd.qe*exptime + emccd.dark_current*exptime + emccd.cic # theoretical electron flux for darks dark_mean = emccd.dark_current*exptime + emccd.cic + + if smear: + #simulate smear to fluxmap + detector_params = DetectorParams({}) + rowreadtime = detector_params.params['rowreadtime'] + smear = np.zeros_like(fluxmap) + m = len(smear) + for r in range(m): + columnsum = 0 + for i in range(r+1): + columnsum = columnsum + rowreadtime*fluxmap[i,:] + smear[r,:] = columnsum + + fluxmap = fluxmap + smear/exptime frame_e_list = [] frame_e_dark_list = [] diff --git a/corgidrp/photon_counting.py b/corgidrp/photon_counting.py index be807736..132b556e 100644 --- a/corgidrp/photon_counting.py +++ b/corgidrp/photon_counting.py @@ -60,7 +60,7 @@ def photon_count(e_image, thresh): return pc_image -def get_pc_mean(input_dataset, pc_master_dark=None, T_factor=None, pc_ecount_max=None, niter=2, mask_filepath=None, safemode=True): +def get_pc_mean(input_dataset, pc_master_dark=None, T_factor=None, pc_ecount_max=None, niter=2, mask_filepath=None, safemode=True, inputmode='illuminated'): """Take a stack of images, frames of the same exposure time, k gain, read noise, and EM gain, and return the mean expected value per pixel. The frames are taken in photon-counting mode, which means high EM @@ -74,13 +74,6 @@ def get_pc_mean(input_dataset, pc_master_dark=None, T_factor=None, pc_ecount_max - have been corrected for nonlinearity (These first 3 steps make the frames L2a.) - have been frame-selected (to weed out bad frames) - have been converted from DN to e- - - have been desmeared if desmearing is appropriate. Under normal - circumstances, these frames should not be desmeared. The only time desmearing - would be useful is in the unexpected case that, for example, - dark current is so high that it stands far above other noise that is - not smeared upon readout, such as clock-induced charge, - fixed-pattern noise, and read noise. For such low-exposure frames, though, - this really should not be a concern. This algorithm will photon count each frame individually, then co-add the photon-counted frames. The co-added frame is then averaged @@ -112,6 +105,9 @@ def get_pc_mean(input_dataset, pc_master_dark=None, T_factor=None, pc_ecount_max If True, the function halts with an exception if the mean intensity of the unmasked pixels (or all pixels if no mask provided) is greater than pc_ecount_max or if the minimum photon-counted pixel value is negative. If False, the function gives a warning for these instead. Defaults to True. + inputmode (str): If 'illuminated', the frames are assumed to be illuminated frames. If 'darks', frames are assumed to be dark frames input for creation of a photon-counted master dark. + This flag shows the user's intention with the input, and this input is checked against the file type of the dataset for compatibility (e.g., if this input is 'darks' while 'VISTYPE' is not equal + to 'DARK', an exception is raised). Returns: corgidrp.data.Dataset or corgidrp.data.Dark: If If the input dataset's header key 'VISTYPE' is not equal to 'DARK', @@ -141,8 +137,16 @@ def get_pc_mean(input_dataset, pc_master_dark=None, T_factor=None, pc_ecount_max if len(val) != 1: raise PhotonCountException('There should only be 1 \'VISTYPE\' value for the dataset.') if val[0] == 'DARK': + if inputmode != 'darks': + raise PhotonCountException('Inputmode is not \'darks\', but the input dataset has \'VISTYPE\' = \'DARK\'.') if pc_master_dark is not None: raise PhotonCountException('The input frames are \'VISTYPE\'=\'DARK\' frames, so no pc_master_dark should be provided.') + if val[0] != 'DARK': + if inputmode != 'illuminated': + raise PhotonCountException('Inputmode is not \'illuminated\', but the input dataset has \'VISTYPE\' not equal to \'DARK\'.') + if 'ISPC' in datasets[0].frames[0].ext_hdr: + if datasets[0].frames[0].ext_hdr['ISPC'] != True: + raise PhotonCountException('\'ISPC\' header value must be True if these frames are to be photon-counted.') dataset = datasets[0] @@ -248,6 +252,10 @@ def get_pc_mean(input_dataset, pc_master_dark=None, T_factor=None, pc_ecount_max if pc_master_dark.ext_hdr['PC_STAT'] != 'photon-counted master dark': raise PhotonCountException('pc_master_dark must be a photon-counted master dark (i.e., ' 'the extension header key \'PC_STAT\' must be \'photon-counted master dark\').') + if 'PCTHRESH' not in pc_master_dark.ext_hdr: + raise PhotonCountException('Threshold should be stored under the header \'PCTHRESH\'.') + if pc_master_dark.ext_hdr['PCTHRESH'] != thresh: + raise PhotonCountException('Threshold used for photon-counted master dark should match the threshold to be used for the illuminated frames.') pc_means.append(pc_master_dark.data) dqs.append(pc_master_dark.dq) errs.append(pc_master_dark.err) @@ -275,6 +283,7 @@ def get_pc_mean(input_dataset, pc_master_dark=None, T_factor=None, pc_ecount_max new_image = data.Image(combined_pc_mean, pri_hdr=pri_hdr, ext_hdr=ext_hdr, err=combined_err, dq=combined_dq, err_hdr=err_hdr, dq_hdr=dq_hdr, input_hdulist=hdulist) new_image.filename = dataset[0].filename.split('.')[0]+'_pc'+filename_end+'.fits' + new_image.ext_hdr['PCTHRESH'] = thresh new_image._record_parent_filenames(input_dataset) pc_ill_dataset = data.Dataset([new_image]) pc_ill_dataset.update_after_processing_step("Photon-counted {0} frames using T_factor={1} and niter={2}. Dark-subtracted: {3}.".format(len(input_dataset), T_factor, niter, dark_sub)) @@ -283,6 +292,7 @@ def get_pc_mean(input_dataset, pc_master_dark=None, T_factor=None, pc_ecount_max ext_hdr['PC_STAT'] = 'photon-counted master dark' ext_hdr['NAXIS1'] = combined_pc_mean.shape[0] ext_hdr['NAXIS2'] = combined_pc_mean.shape[1] + ext_hdr['PCTHRESH'] = thresh pc_dark = data.Dark(combined_pc_mean, pri_hdr=pri_hdr, ext_hdr=ext_hdr, err=combined_err, dq=combined_dq, err_hdr=err_hdr, input_dataset=input_dataset) return pc_dark diff --git a/corgidrp/recipe_templates/l1_to_l2b_pc.json b/corgidrp/recipe_templates/l1_to_l2b_pc.json index a7d1b379..7067704b 100644 --- a/corgidrp/recipe_templates/l1_to_l2b_pc.json +++ b/corgidrp/recipe_templates/l1_to_l2b_pc.json @@ -48,6 +48,33 @@ "Dark" : "AUTOMATIC" } }, + { + "name" : "desmear", + "calibs" : { + "DetectorParams" : "AUTOMATIC" + } + }, + { + "name" : "cti_correction", + "calibs" : { + "TrapCalibration" : "AUTOMATIC,OPTIONAL" + } + }, + { + "name" : "flat_division", + "calibs" : { + "FlatField" : "AUTOMATIC" + } + }, + { + "name" : "correct_bad_pixels", + "calibs" : { + "BadPixelMap" : "AUTOMATIC" + } + }, + { + "name" : "update_to_l2b" + }, { "name" : "save" } diff --git a/corgidrp/recipe_templates/l1_to_l2b_pc_dark.json b/corgidrp/recipe_templates/l1_to_l2b_pc_dark.json index 3e2c41d7..02e97622 100644 --- a/corgidrp/recipe_templates/l1_to_l2b_pc_dark.json +++ b/corgidrp/recipe_templates/l1_to_l2b_pc_dark.json @@ -43,7 +43,10 @@ } }, { - "name": "get_pc_mean" + "name": "get_pc_mean", + "keywords" : { + "inputmode" : "darks" + } }, { "name" : "save" diff --git a/corgidrp/recipe_templates/l2a_to_l2b_pc.json b/corgidrp/recipe_templates/l2a_to_l2b_pc.json index 4775ae38..0d69f670 100644 --- a/corgidrp/recipe_templates/l2a_to_l2b_pc.json +++ b/corgidrp/recipe_templates/l2a_to_l2b_pc.json @@ -22,6 +22,33 @@ "Dark" : "AUTOMATIC" } }, + { + "name" : "desmear", + "calibs" : { + "DetectorParams" : "AUTOMATIC" + } + }, + { + "name" : "cti_correction", + "calibs" : { + "TrapCalibration" : "AUTOMATIC,OPTIONAL" + } + }, + { + "name" : "flat_division", + "calibs" : { + "FlatField" : "AUTOMATIC" + } + }, + { + "name" : "correct_bad_pixels", + "calibs" : { + "BadPixelMap" : "AUTOMATIC" + } + }, + { + "name" : "update_to_l2b" + }, { "name" : "save" } diff --git a/corgidrp/recipe_templates/l2a_to_l2b_pc_dark.json b/corgidrp/recipe_templates/l2a_to_l2b_pc_dark.json index c1b3df62..e53ac9d6 100644 --- a/corgidrp/recipe_templates/l2a_to_l2b_pc_dark.json +++ b/corgidrp/recipe_templates/l2a_to_l2b_pc_dark.json @@ -17,7 +17,10 @@ } }, { - "name": "get_pc_mean" + "name": "get_pc_mean", + "keywords" : { + "inputmode" : "darks" + } }, { "name" : "save" diff --git a/tests/e2e_tests/astrom_e2e.py b/tests/e2e_tests/astrom_e2e.py index ff491ebd..facf1856 100644 --- a/tests/e2e_tests/astrom_e2e.py +++ b/tests/e2e_tests/astrom_e2e.py @@ -178,7 +178,7 @@ def test_astrom_e2e(tvacdata_path, e2eoutput_path): assert dec == pytest.approx(target[1], abs=8.333e-7) if __name__ == "__main__": - tvacdata_dir = "/Users/macuser/Roman/corgi_contributions/Callibration_Notebooks/TVAC" + tvacdata_dir = "/Users/kevinludwick/Library/CloudStorage/Box-Box/CGI_TVAC_Data/Working_Folder/" #"/Users/macuser/Roman/corgi_contributions/Callibration_Notebooks/TVAC" outputdir = thisfile_dir ap = argparse.ArgumentParser(description="run the l1->l2b->boresight end-to-end test") diff --git a/tests/e2e_tests/noisemap_cal_e2e.py b/tests/e2e_tests/noisemap_cal_e2e.py index ee68f29f..c0ddc1ed 100644 --- a/tests/e2e_tests/noisemap_cal_e2e.py +++ b/tests/e2e_tests/noisemap_cal_e2e.py @@ -145,6 +145,9 @@ def test_noisemap_calibration_from_l1(tvacdata_path, e2eoutput_path): kgain_val = 8.7 # From TVAC-20 noise characterization measurements kgain = data.KGain(np.array([[kgain_val]]), pri_hdr=pri_hdr, ext_hdr=ext_hdr, input_dataset=mock_input_dataset) + # add in keywords that didn't make it into mock_kgain.fits, using values used in mocks.create_photon_countable_frames() + kgain.ext_hdr['RN'] = 100 + kgain.ext_hdr['RN_ERR'] = 0 kgain.save(filedir=noisemap_outputdir, filename="mock_kgain.fits") this_caldb.create_entry(kgain) @@ -325,11 +328,20 @@ def test_noisemap_calibration_from_l2a(tvacdata_path, e2eoutput_path): mock_input_dataset = data.Dataset(mock_cal_filelist) this_caldb = caldb.CalDB() # connection to cal DB - + # remove other KGain calibrations that may exist in case they don't have the added header keywords + for i in range(len(this_caldb._db['Type'])): + if this_caldb._db['Type'][i] == 'KGain': + this_caldb._db = this_caldb._db.drop(i) + elif this_caldb._db['Type'][i] == 'Dark': + this_caldb._db = this_caldb._db.drop(i) + this_caldb.save() # KGain calibration kgain_val = 8.7 # From TVAC-20 noise characterization measurements kgain = data.KGain(np.array([[kgain_val]]), pri_hdr=pri_hdr, ext_hdr=ext_hdr, input_dataset=mock_input_dataset) + # add in keywords that didn't make it into mock_kgain.fits, using values used in mocks.create_photon_countable_frames() + kgain.ext_hdr['RN'] = 100 + kgain.ext_hdr['RN_ERR'] = 0 kgain.save(filedir=noisemap_outputdir, filename="mock_kgain.fits") this_caldb.create_entry(kgain) diff --git a/tests/e2e_tests/photon_count_e2e.py b/tests/e2e_tests/photon_count_e2e.py index fe26d4d6..9c2c2454 100644 --- a/tests/e2e_tests/photon_count_e2e.py +++ b/tests/e2e_tests/photon_count_e2e.py @@ -15,9 +15,13 @@ import corgidrp.detector as detector @pytest.mark.e2e -def test_expected_results_e2e(file_dir): +def test_expected_results_e2e(tvacdata_dir, file_dir): + processed_cal_path = os.path.join(tvacdata_dir, "TV-36_Coronagraphic_Data", "Cals") + flat_path = os.path.join(processed_cal_path, "flat.fits") + bp_path = os.path.join(processed_cal_path, "bad_pix.fits") + np.random.seed(1234) - ill_dataset, dark_dataset, ill_mean, dark_mean = mocks.create_photon_countable_frames(Nbrights=160, Ndarks=161, cosmic_rate=1) + ill_dataset, dark_dataset, ill_mean, dark_mean = mocks.create_photon_countable_frames(Nbrights=160, Ndarks=161, cosmic_rate=1, flux=0.5) output_dir = os.path.join(file_dir, 'pc_sim_test_data') output_ill_dir = os.path.join(output_dir, 'ill_frames') output_dark_dir = os.path.join(output_dir, 'dark_frames') @@ -102,15 +106,29 @@ def test_expected_results_e2e(file_dir): new_nonlinearity.pri_hdr = pri_hdr new_nonlinearity.ext_hdr = ext_hdr this_caldb.create_entry(new_nonlinearity) + + ## Flat field + with fits.open(flat_path) as hdulist: + flat_dat = hdulist[0].data + flat = data.FlatField(flat_dat, pri_hdr=pri_hdr, ext_hdr=ext_hdr, input_dataset=mock_input_dataset) + flat.save(filedir=output_dir, filename="mock_flat.fits") + this_caldb.create_entry(flat) + + # bad pixel map + with fits.open(bp_path) as hdulist: + bp_dat = hdulist[0].data + bp_map = data.BadPixelMap(bp_dat, pri_hdr=pri_hdr, ext_hdr=ext_hdr, input_dataset=mock_input_dataset) + bp_map.save(filedir=output_dir, filename="mock_bpmap.fits") + this_caldb.create_entry(bp_map) + # make PC dark walker.walk_corgidrp(l1_data_dark_filelist, '', output_dir, template="l1_to_l2b_pc_dark.json") for f in os.listdir(output_dir): if f.endswith('_pc_dark.fits'): pc_dark_filename = f - # add the Dark to the Caldb for processing the illuminated + # calDB was just updated with the PC Dark that was created with the walker above pc_dark_file = os.path.join(output_dir, pc_dark_filename) - dark_entry = data.Dark(pc_dark_file) - this_caldb.create_entry(dark_entry) + # make PC illuminated, subtracting the PC dark walker.walk_corgidrp(l1_data_ill_filelist, '', output_dir, template="l1_to_l2b_pc.json") # get photon-counted frame @@ -129,9 +147,16 @@ def test_expected_results_e2e(file_dir): assert pc_frame_err.min() >= 0 assert pc_dark_frame_err.min() >= 0 - this_caldb.remove_entry(kgain) - this_caldb.remove_entry(noise_map) - this_caldb.remove_entry(new_nonlinearity) + # load in CalDB again to reflect the PC Dark that was implicitly added in (but not found in this_caldb, which was loaded before the Dark was created) + post_caldb = caldb.CalDB() + post_caldb.remove_entry(kgain) + post_caldb.remove_entry(noise_map) + post_caldb.remove_entry(new_nonlinearity) + post_caldb.remove_entry(flat) + post_caldb.remove_entry(bp_map) + for i in range(len(post_caldb._db['Type'])): + if post_caldb._db['Type'][i] == 'Dark': + post_caldb._db = post_caldb._db.drop(i) if __name__ == "__main__": @@ -142,10 +167,13 @@ def test_expected_results_e2e(file_dir): # workflow. thisfile_dir = os.path.dirname(__file__) outputdir = thisfile_dir + tvacdata_dir = "/Users/kevinludwick/Library/CloudStorage/Box-Box/CGI_TVAC_Data/Working_Folder/" - ap = argparse.ArgumentParser(description="run the l1->l2b PC end-to-end test") + ap = argparse.ArgumentParser(description="run the l1->l2a end-to-end test") + ap.add_argument("-tvac", "--tvacdata_dir", default=tvacdata_dir, + help="Path to CGI_TVAC_Data Folder [%(default)s]") ap.add_argument("-o", "--outputdir", default=outputdir, help="directory to write results to [%(default)s]") args = ap.parse_args() outputdir = args.outputdir - test_expected_results_e2e(outputdir) + test_expected_results_e2e(tvacdata_dir, outputdir) diff --git a/tests/test_photon_counting.py b/tests/test_photon_counting.py index 2fddcadb..ac63e819 100644 --- a/tests/test_photon_counting.py +++ b/tests/test_photon_counting.py @@ -20,14 +20,36 @@ def test_negative(): def test_pc(): - '''Test that a pixel that is masked heavily (reducing the number of usable frames for that pixel) has a bigger err than the average pixel. - And if masking is all the way through for a certain pixel, that pixel will - be flagged in the DQ. Also, make sure there is failure when niter<1 and threshold<0, and test other built-in catches. - Test safemode. Test mask_filepath to ensure failure for bright pixels outside region of interest when no mask used - and success when mask used to mask those bright pixels.''' + ''' + Tests that a pixel that is masked heavily (reducing the number of usable frames for that pixel) has a bigger err than the average pixel. + + Tests that if masking is all the way through for a certain pixel, that pixel will + be flagged in the DQ. + + Make sure there is failure when the number of iterations niter<1. + + Tests that the use of a mask to mask out bright pixels when photon-counting, to ensure failure for bright pixels outside region of interest when no mask used + and success when mask used to mask those bright pixels. + + Tests safemode, which makes the function issue warnings instead of crashing (useful for HOWFSC's iterative digging of dark holes). + + Checks various inputs: threshold<0 gives exception, thresh>=em_gain gives exception, providing dataset of VISTYPE='DARK' while also inputting a photon-counted master dark raises exception, exception raised when 'CMDGAIN' not the same for all frames. + + Checks various metadata changes: the expected filename and history edit for the output is done (for case of dark subtraction and the case of no dark subtraction), the master dark is indicated as such via a header keyword 'PC_STAT'. + + Tests that output average value of the photon-counted frames agrees with what was simulated (for the case of dark subtraction and no dark subtraction). + + Tests that if 'ISPC' is in image header and if it is False, get_pc_mean() raises an exception. + + Tests that an exception occurs if the user input "inputmode" is inconsistent with the input dataset type. + + Tests that an exception occurs if the photon-counted master dark does not have the header 'PCTHRESH'. + + Tests that an exception occurs if the photon-counted master dark has a 'PCTHRESH' value than the one to be used on the illuminated dataset. + ''' # exposure time too long to get reasonable photon-counted result (after photometric correction) np.random.seed(555) - dataset_err, dark_dataset_err, ill_mean, dark_mean = mocks.create_photon_countable_frames(Nbrights=160, Ndarks=160, cosmic_rate=0, full_frame=False) + dataset_err, dark_dataset_err, ill_mean, dark_mean = mocks.create_photon_countable_frames(Nbrights=160, Ndarks=160, cosmic_rate=0, full_frame=False, smear=False) # instead of running through walker, just do the pre-processing steps simply # using EM gain=5000 and kgain=7 and bias=20000 and read noise = 100 and QE=0.9 (quantum efficiency), from mocks.create_photon_countable_frames() for f in dataset_err.frames: @@ -47,7 +69,7 @@ def test_pc(): f.ext_hdr['RN'] = 100 f.ext_hdr['KGAIN'] = 7 # process the frames to make PC dark - pc_dark = get_pc_mean(dark_dataset_err) + pc_dark = get_pc_mean(dark_dataset_err, inputmode='darks') assert pc_dark.ext_hdr['PC_STAT'] == 'photon-counted master dark' # now process illuminated frames and subtract the PC dark pc_dataset_err = get_pc_mean(dataset_err, pc_master_dark=pc_dark) @@ -72,14 +94,26 @@ def test_pc(): # when thresh>=em_gain, exception with pytest.raises(Exception): get_pc_mean(dataset_err, T_factor=50, niter=2) - # can't provide Dark class instance while "VISTYPE" of input dataset is 'DARK' + # can't provide master dark while "VISTYPE" of input dataset is 'DARK' with pytest.raises(PhotonCountException): - get_pc_mean(dark_dataset_err, pc_master_dark=pc_dark) + get_pc_mean(dark_dataset_err, pc_master_dark=pc_dark, inputmode='darks') # must have same 'CMDGAIN' and other header values throughout input dataset dark_dataset_err.frames[0].ext_hdr['CMDGAIN'] = 4999 with pytest.raises(PhotonCountException): - get_pc_mean(dark_dataset_err, pc_master_dark=pc_dark) + get_pc_mean(dark_dataset_err, inputmode='dark') + # change back: + dark_dataset_err.frames[0].ext_hdr['CMDGAIN'] = 5000 + # test to make sure PC dark's threshold matches the one used for illuminated frames + with pytest.raises(PhotonCountException): + pc_dark.ext_hdr['PCTHRESH'] = 499 + get_pc_mean(dataset_err, pc_master_dark=pc_dark, inputmode='illuminated') + # PC dark should have header 'PCTHRESH' + with pytest.raises(PhotonCountException): + pc_dark.ext_hdr.pop("PCTHRESH") + get_pc_mean(dataset_err, pc_master_dark=pc_dark, inputmode='illuminated') + # set it back: + pc_dark.ext_hdr['PCTHRESH'] = 500 # to trigger pc_ecount_max error for f in dataset_err.frames[:-2]: #all but 2 of the frames f.data[22:40,23:49] = np.abs(f.data.astype(float)[22:40,23:49]*3000) @@ -111,6 +145,17 @@ def test_pc(): # test safemode=False: should get warning instead of exception with pytest.warns(UserWarning): get_pc_mean(dataset_err, safemode=False) + # test ISPC header value + with pytest.raises(PhotonCountException): + dataset_err.frames[0].ext_hdr['ISPC'] = False + get_pc_mean(dataset_err, pc_master_dark=pc_dark) + # set to True now + dataset_err.frames[0].ext_hdr['ISPC'] = True + # test inputmode's compatibility with dataset type + with pytest.raises(PhotonCountException): + get_pc_mean(dataset_err, pc_master_dark=pc_dark, inputmode='darks') + with pytest.raises(PhotonCountException): + get_pc_mean(dark_dataset_err, inputmode='illuminated') if __name__ == '__main__': test_pc() From 535ff290bcb5f4009065bf24bae32a582601c76c Mon Sep 17 00:00:00 2001 From: Kevin Ludwick Date: Mon, 27 Jan 2025 08:55:41 -0600 Subject: [PATCH 20/26] using remove_entry now for PC dark --- tests/e2e_tests/photon_count_e2e.py | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/tests/e2e_tests/photon_count_e2e.py b/tests/e2e_tests/photon_count_e2e.py index 9c2c2454..7b1df530 100644 --- a/tests/e2e_tests/photon_count_e2e.py +++ b/tests/e2e_tests/photon_count_e2e.py @@ -154,9 +154,8 @@ def test_expected_results_e2e(tvacdata_dir, file_dir): post_caldb.remove_entry(new_nonlinearity) post_caldb.remove_entry(flat) post_caldb.remove_entry(bp_map) - for i in range(len(post_caldb._db['Type'])): - if post_caldb._db['Type'][i] == 'Dark': - post_caldb._db = post_caldb._db.drop(i) + pc_dark = data.Dark(pc_dark_file) + post_caldb.remove_entry(pc_dark) if __name__ == "__main__": From 752a1ae20fc637b4e6c691f2f626b2f16d906034 Mon Sep 17 00:00:00 2001 From: Kevin Ludwick Date: Mon, 10 Feb 2025 10:13:28 -0600 Subject: [PATCH 21/26] walker can guess PC templates --- .../l2a_build_trad_dark_image.json | 26 +++++++++++++++++++ corgidrp/walker.py | 19 ++++++++++++-- tests/e2e_tests/noisemap_cal_e2e.py | 16 +++++++++--- tests/e2e_tests/photon_count_e2e.py | 6 ++--- 4 files changed, 58 insertions(+), 9 deletions(-) create mode 100644 corgidrp/recipe_templates/l2a_build_trad_dark_image.json diff --git a/corgidrp/recipe_templates/l2a_build_trad_dark_image.json b/corgidrp/recipe_templates/l2a_build_trad_dark_image.json new file mode 100644 index 00000000..7ba7146b --- /dev/null +++ b/corgidrp/recipe_templates/l2a_build_trad_dark_image.json @@ -0,0 +1,26 @@ +{ + "name" : "build_trad_dark_l2a", + "template" : true, + "drpconfig" : { + "track_individual_errors" : false + }, + "inputs" : [], + "outputdir" : "", + "steps" : [ + { + "name" : "convert_to_electrons", + "calibs" : { + "KGain" : "AUTOMATIC" + } + }, + { + "name" : "build_trad_dark", + "calibs" : { + "DetectorParams" : "AUTOMATIC" + } + }, + { + "name" : "save" + } + ] +} \ No newline at end of file diff --git a/corgidrp/walker.py b/corgidrp/walker.py index d7a1774d..81d5c981 100644 --- a/corgidrp/walker.py +++ b/corgidrp/walker.py @@ -224,6 +224,7 @@ def guess_template(dataset): str or list: the best template filename or a list of multiple template filenames """ image = dataset[0] # first image for convenience + _, unique_vals = dataset.split_dataset(exthdr_keywords=['EXPTIME', 'CMDGAIN', 'KGAIN']) if image.ext_hdr['DATA_LEVEL'] == "L1": if 'VISTYPE' not in image.pri_hdr: # this is probably IIT test data. Do generic processing @@ -237,14 +238,28 @@ def guess_template(dataset): elif image.pri_hdr['VISTYPE'] == "FFIELD": recipe_filename = "l1_flat_and_bp.json" elif image.pri_hdr['VISTYPE'] == "DARK": - recipe_filename = "l1_to_l2a_noisemap.json" + if image.ext_hdr['ISPC']: + recipe_filename = "l1_to_l2b_pc_dark.json" + elif len(unique_vals) > 1: # darks for noisemap creation + recipe_filename = "l1_to_l2a_noisemap.json" + elif len(unique_vals) == 1: # traditional darks + recipe_filename = "build_trad_dark_image.json" elif image.pri_hdr['VISTYPE'] == "PUPILIMG": recipe_filename = ["l1_to_l2a_nonlin.json", "l1_to_kgain.json"] else: recipe_filename = "l1_to_l2b.json" + if image.ext_hdr['ISPC'] and image.pri_hdr['VISTYPE'] != "DARK": + recipe_filename = "l1_to_l2b_pc.json" elif image.ext_hdr['DATA_LEVEL'] == "L2a": if image.pri_hdr['VISTYPE'] == "DARK": - recipe_filename = "l2a_to_l2a_noisemap.json" + if image.ext_hdr['ISPC']: + recipe_filename = "l2a_to_l2b_pc_dark.json" + elif len(unique_vals) > 1: # darks for noisemap creation + recipe_filename = "l2a_to_l2a_noisemap.json" + elif len(unique_vals) == 1: # traditional darks + recipe_filename = "l2a_build_trad_dark_image.json" + if image.ext_hdr['ISPC'] and image.pri_hdr['VISTYPE'] != "DARK": + recipe_filename = "l2a_to_l2b_pc.json" else: raise NotImplementedError() else: diff --git a/tests/e2e_tests/noisemap_cal_e2e.py b/tests/e2e_tests/noisemap_cal_e2e.py index c0ddc1ed..e75682aa 100644 --- a/tests/e2e_tests/noisemap_cal_e2e.py +++ b/tests/e2e_tests/noisemap_cal_e2e.py @@ -129,6 +129,15 @@ def test_noisemap_calibration_from_l1(tvacdata_path, e2eoutput_path): mock_input_dataset = data.Dataset(mock_cal_filelist) this_caldb = caldb.CalDB() # connection to cal DB + # remove other KGain calibrations that may exist in case they don't have the added header keywords + for i in range(len(this_caldb._db['Type'])): + if this_caldb._db['Type'][i] == 'KGain': + this_caldb._db = this_caldb._db.drop(i) + elif this_caldb._db['Type'][i] == 'Dark': + this_caldb._db = this_caldb._db.drop(i) + elif this_caldb._db['Type'][i] == 'NonLinearityCalibration': + this_caldb._db = this_caldb._db.drop(i) + this_caldb.save() pri_hdr, ext_hdr = mocks.create_default_headers() ext_hdr["DRPCTIME"] = time.Time.now().isot @@ -170,16 +179,15 @@ def test_noisemap_calibration_from_l1(tvacdata_path, e2eoutput_path): ##### Check against II&T ("TVAC") data corgidrp_noisemap_fname = os.path.join(noisemap_outputdir,output_filename) # iit_noisemap_fname = os.path.join(iit_noisemap_datadir,"iit_test_noisemaps.fits") - corgidrp_noisemap = data.autoload(corgidrp_noisemap_fname) - this_caldb.remove_entry(corgidrp_noisemap) - + assert(np.nanmax(np.abs(corgidrp_noisemap.data[0]- F_map)) < 1e-11) assert(np.nanmax(np.abs(corgidrp_noisemap.data[1]- C_map)) < 1e-11) assert(np.nanmax(np.abs(corgidrp_noisemap.data[2]- D_map)) < 1e-11) assert(np.abs(corgidrp_noisemap.ext_hdr['B_O']- bias_offset) < 1e-11) pass + this_caldb.remove_entry(corgidrp_noisemap) # for noise_ext in ["FPN_map","CIC_map","DC_map"]: # corgi_dat = detector.imaging_slice('SCI', corgidrp_noisemap.__dict__[noise_ext]) # iit_dat = detector.imaging_slice('SCI', iit_noisemap.__dict__[noise_ext]) @@ -437,5 +445,5 @@ def test_noisemap_calibration_from_l2a(tvacdata_path, e2eoutput_path): args = ap.parse_args(args_here) tvacdata_dir = args.tvacdata_dir outputdir = args.outputdir - #test_noisemap_calibration_from_l1(tvacdata_dir, outputdir) + test_noisemap_calibration_from_l1(tvacdata_dir, outputdir) test_noisemap_calibration_from_l2a(tvacdata_dir, outputdir) diff --git a/tests/e2e_tests/photon_count_e2e.py b/tests/e2e_tests/photon_count_e2e.py index 7b1df530..83e04a4e 100644 --- a/tests/e2e_tests/photon_count_e2e.py +++ b/tests/e2e_tests/photon_count_e2e.py @@ -15,14 +15,14 @@ import corgidrp.detector as detector @pytest.mark.e2e -def test_expected_results_e2e(tvacdata_dir, file_dir): - processed_cal_path = os.path.join(tvacdata_dir, "TV-36_Coronagraphic_Data", "Cals") +def test_expected_results_e2e(tvacdata_path, e2eoutput_path): + processed_cal_path = os.path.join(tvacdata_path, "TV-36_Coronagraphic_Data", "Cals") flat_path = os.path.join(processed_cal_path, "flat.fits") bp_path = os.path.join(processed_cal_path, "bad_pix.fits") np.random.seed(1234) ill_dataset, dark_dataset, ill_mean, dark_mean = mocks.create_photon_countable_frames(Nbrights=160, Ndarks=161, cosmic_rate=1, flux=0.5) - output_dir = os.path.join(file_dir, 'pc_sim_test_data') + output_dir = os.path.join(e2eoutput_path, 'pc_sim_test_data') output_ill_dir = os.path.join(output_dir, 'ill_frames') output_dark_dir = os.path.join(output_dir, 'dark_frames') if not os.path.exists(output_dir): From 8ce95862759a880cb478b30386e3f69eb6c9e206 Mon Sep 17 00:00:00 2001 From: Kevin Ludwick Date: Mon, 10 Feb 2025 11:35:17 -0600 Subject: [PATCH 22/26] fixed walker PC checks 'ISPC' is not currently in any files; it will be added, though --- corgidrp/calibrate_kgain.py | 25 +++++++++--------- corgidrp/calibrate_nonlin.py | 16 ++++++++---- corgidrp/walker.py | 50 ++++++++++++++++++++++++------------ tests/test_kgain_cal.py | 11 ++++++-- tests/test_nonlin_cal.py | 13 +++++++++- 5 files changed, 77 insertions(+), 38 deletions(-) diff --git a/corgidrp/calibrate_kgain.py b/corgidrp/calibrate_kgain.py index 9bc6c622..ae4e7a5b 100644 --- a/corgidrp/calibrate_kgain.py +++ b/corgidrp/calibrate_kgain.py @@ -12,7 +12,7 @@ from corgidrp.detector import slice_section, detector_areas # Dictionary with constant kgain calibration parameters -kgain_params= { +kgain_params_default= { # ROI constants 'rowroi1': 9, 'rowroi2': 1000, @@ -30,13 +30,13 @@ 'signal_bins_N': 400, } -def check_kgain_params( - ): - """ Checks integrity of kgain parameters in the dictionary kgain_params. """ - if 'offset_colroi1' not in kgain_params: - raise ValueError('Missing parameter: offset_colroi1.') - if 'offset_colroi2' not in kgain_params: - raise ValueError('Missing parameter: offset_colroi2.') +def check_kgain_params(kgain_params): + """ Checks integrity of kgain parameters in the dictionary kgain_params. + + Args: + kgain_params (dict): Dictionary of parameters used for calibrating the k gain. + """ + if 'rowroi1' not in kgain_params: raise ValueError('Missing parameter: rowroi1.') if 'rowroi2' not in kgain_params: @@ -54,10 +54,6 @@ def check_kgain_params( if 'signal_bins_N' not in kgain_params: raise ValueError('Missing parameter: signal_bins_N.') - if not isinstance(kgain_params['offset_colroi1'], (float, int)): - raise TypeError('offset_colroi1 is not a number') - if not isinstance(kgain_params['offset_colroi2'], (float, int)): - raise TypeError('offset_colroi2 is not a number') if not isinstance(kgain_params['rowroi1'], (float, int)): raise TypeError('rowroi1 is not a number') if not isinstance(kgain_params['rowroi2'], (float, int)): @@ -281,7 +277,7 @@ def calibrate_kgain(dataset_kgain, n_cal=10, n_mean=30, min_val=800, max_val=3000, binwidth=68, make_plot=True,plot_outdir='figures', show_plot=False, logspace_start=-1, logspace_stop=4, logspace_num=200, - verbose=False, detector_regions=None, kgain_params=kgain_params): + verbose=False, detector_regions=None, kgain_params=kgain_params_default): """ kgain (e-/DN) is calculated from the means and variances within the defined minimum and maximum mean values. A photon transfer curve @@ -362,6 +358,9 @@ def calibrate_kgain(dataset_kgain, transfer curve (in e-/DN). The expected value of kgain for EXCAM with flight readout sequence should be between 8 and 9 e-/DN """ + + check_kgain_params(kgain_params) + if detector_regions is None: detector_regions = detector_areas diff --git a/corgidrp/calibrate_nonlin.py b/corgidrp/calibrate_nonlin.py index d2f3894a..020189f8 100644 --- a/corgidrp/calibrate_nonlin.py +++ b/corgidrp/calibrate_nonlin.py @@ -16,7 +16,7 @@ from corgidrp.calibrate_kgain import CalKgainException # Dictionary with constant non-linearity calibration parameters -nonlin_params = { +nonlin_params_default = { # ROI constants 'rowroi1': 305, 'rowroi2': 736, @@ -44,9 +44,12 @@ 'min_mask_factor': 1.1, } -def check_nonlin_params( - ): - """ Checks integrity of kgain parameters in the dictionary nonlin_params. """ +def check_nonlin_params(nonlin_params): + """ Checks integrity of kgain parameters in the dictionary nonlin_params. + + Args: + nonlin_params (dict): Dictionary of parameters used for calibrating nonlinearity. + """ if 'rowroi1' not in nonlin_params: raise ValueError('Missing parameter: rowroi1.') if 'rowroi2' not in nonlin_params: @@ -124,7 +127,7 @@ def calibrate_nonlin(dataset_nl, pfit_upp_cutoff1 = -2, pfit_upp_cutoff2 = -3, pfit_low_cutoff1 = 2, pfit_low_cutoff2 = 1, make_plot=True, plot_outdir='figures', show_plot=False, - verbose=False, nonlin_params=nonlin_params): + verbose=False, nonlin_params=nonlin_params_default): """ Function that derives the non-linearity calibration table for a set of DN and EM values. @@ -227,6 +230,9 @@ def calibrate_nonlin(dataset_nl, input signal in DN is the first column. Signal values start with min_write and run through max_write in steps of 20 DN. """ + + check_nonlin_params(nonlin_params) + # dataset_nl.all_data must be 3-D if np.ndim(dataset_nl.all_data) != 3: raise Exception('dataset_nl.all_data must be 3-D') diff --git a/corgidrp/walker.py b/corgidrp/walker.py index 81d5c981..11d9a649 100644 --- a/corgidrp/walker.py +++ b/corgidrp/walker.py @@ -224,7 +224,6 @@ def guess_template(dataset): str or list: the best template filename or a list of multiple template filenames """ image = dataset[0] # first image for convenience - _, unique_vals = dataset.split_dataset(exthdr_keywords=['EXPTIME', 'CMDGAIN', 'KGAIN']) if image.ext_hdr['DATA_LEVEL'] == "L1": if 'VISTYPE' not in image.pri_hdr: # this is probably IIT test data. Do generic processing @@ -238,28 +237,45 @@ def guess_template(dataset): elif image.pri_hdr['VISTYPE'] == "FFIELD": recipe_filename = "l1_flat_and_bp.json" elif image.pri_hdr['VISTYPE'] == "DARK": - if image.ext_hdr['ISPC']: - recipe_filename = "l1_to_l2b_pc_dark.json" - elif len(unique_vals) > 1: # darks for noisemap creation - recipe_filename = "l1_to_l2a_noisemap.json" - elif len(unique_vals) == 1: # traditional darks - recipe_filename = "build_trad_dark_image.json" + _, unique_vals = dataset.split_dataset(exthdr_keywords=['EXPTIME', 'CMDGAIN', 'KGAIN']) + if 'ISPC' in image.ext_hdr: # to be added eventually, so check if it's present first + if image.ext_hdr['ISPC']: + recipe_filename = "l1_to_l2b_pc_dark.json" + elif len(unique_vals) > 1: # darks for noisemap creation + recipe_filename = "l1_to_l2a_noisemap.json" + elif len(unique_vals) == 1: # traditional darks + recipe_filename = "build_trad_dark_image.json" + else: # can't distinguish PC or analog then; default to analog + if len(unique_vals) > 1: # darks for noisemap creation + recipe_filename = "l1_to_l2a_noisemap.json" + elif len(unique_vals) == 1: # traditional darks + recipe_filename = "build_trad_dark_image.json" + elif image.pri_hdr['VISTYPE'] == "PUPILIMG": recipe_filename = ["l1_to_l2a_nonlin.json", "l1_to_kgain.json"] else: recipe_filename = "l1_to_l2b.json" - if image.ext_hdr['ISPC'] and image.pri_hdr['VISTYPE'] != "DARK": - recipe_filename = "l1_to_l2b_pc.json" + if 'ISPC' in image.ext_hdr: + if image.ext_hdr['ISPC'] and image.pri_hdr['VISTYPE'] != "DARK": + recipe_filename = "l1_to_l2b_pc.json" elif image.ext_hdr['DATA_LEVEL'] == "L2a": if image.pri_hdr['VISTYPE'] == "DARK": - if image.ext_hdr['ISPC']: - recipe_filename = "l2a_to_l2b_pc_dark.json" - elif len(unique_vals) > 1: # darks for noisemap creation - recipe_filename = "l2a_to_l2a_noisemap.json" - elif len(unique_vals) == 1: # traditional darks - recipe_filename = "l2a_build_trad_dark_image.json" - if image.ext_hdr['ISPC'] and image.pri_hdr['VISTYPE'] != "DARK": - recipe_filename = "l2a_to_l2b_pc.json" + _, unique_vals = dataset.split_dataset(exthdr_keywords=['EXPTIME', 'CMDGAIN', 'KGAIN']) + if 'ISPC' in image.ext_hdr: + if image.ext_hdr['ISPC']: + recipe_filename = "l2a_to_l2b_pc_dark.json" + elif len(unique_vals) > 1: # darks for noisemap creation + recipe_filename = "l2a_to_l2a_noisemap.json" + elif len(unique_vals) == 1: # traditional darks + recipe_filename = "l2a_build_trad_dark_image.json" + else: # can't distinguish PC or analog then; default to analog + if len(unique_vals) > 1: # darks for noisemap creation + recipe_filename = "l2a_to_l2a_noisemap.json" + elif len(unique_vals) == 1: # traditional darks + recipe_filename = "l2a_build_trad_dark_image.json" + if 'ISPC' in image.ext_hdr: + if image.ext_hdr['ISPC'] and image.pri_hdr['VISTYPE'] != "DARK": + recipe_filename = "l2a_to_l2b_pc.json" else: raise NotImplementedError() else: diff --git a/tests/test_kgain_cal.py b/tests/test_kgain_cal.py index 67614052..23ebfa8c 100644 --- a/tests/test_kgain_cal.py +++ b/tests/test_kgain_cal.py @@ -17,7 +17,7 @@ from corgidrp import check from corgidrp.data import Image, Dataset from corgidrp.mocks import (create_default_headers, make_fluxmap_image, nonlin_coefs) -from corgidrp.calibrate_kgain import (calibrate_kgain, CalKgainException, kgain_params) +from corgidrp.calibrate_kgain import (calibrate_kgain, CalKgainException, kgain_params_default) np.random.seed(8585) @@ -137,11 +137,18 @@ def test_expected_results_sub(): """Outputs are as expected, for imported frames.""" kgain = calibrate_kgain(dataset_kg, n_cal, n_mean, min_val, max_val, binwidth) - signal_bins_N = kgain_params['signal_bins_N'] + signal_bins_N = kgain_params_default['signal_bins_N'] # kgain - should be close to the assumed value assert np.isclose(round(kgain.value,1), kgain_in, atol=0.5) assert np.all(np.equal(kgain.ptc.shape, (signal_bins_N,2))) + # test bad input for kgain_params + kgain_params_bad = kgain_params_default.copy() + kgain_params_bad['colroi2'] = 'foo' + with pytest.raises(TypeError): + calibrate_kgain(dataset_kg, n_cal, n_mean, min_val, max_val, binwidth, + kgain_params=kgain_params_bad) + def test_psi(): """These three below must be positive scalar integers.""" check_list = test_check.psilist diff --git a/tests/test_nonlin_cal.py b/tests/test_nonlin_cal.py index 9a5122c8..a1461a16 100644 --- a/tests/test_nonlin_cal.py +++ b/tests/test_nonlin_cal.py @@ -15,7 +15,7 @@ import test_check from corgidrp import check from corgidrp.data import Image, Dataset -from corgidrp.calibrate_nonlin import (calibrate_nonlin, CalNonlinException) +from corgidrp.calibrate_nonlin import (calibrate_nonlin, CalNonlinException, nonlin_params_default) from corgidrp.mocks import (create_default_headers, make_fluxmap_image, nonlin_coefs) ############################# prepare simulated frames ####################### @@ -248,8 +248,19 @@ def test_psi(): for mmerr in check_list: with pytest.raises(TypeError): calibrate_nonlin(dataset_nl, n_cal, n_mean, mmerr, min_write, max_write) + +def test_nonlin_params(): + '''Check that the check function for nonlin_params works as expected.''' + # test bad input for nonlin_params + nonlin_params_bad = nonlin_params_default.copy() + nonlin_params_bad['colroi2'] = 'foo' + with pytest.raises(TypeError): + calibrate_nonlin(dataset_nl, n_cal, n_mean, norm_val, min_write, max_write, + nonlin_params=nonlin_params_bad) if __name__ == '__main__': + print('Running test_nonlin_params') + test_nonlin_params() print('Running test_expected_results_nom_sub') test_expected_results_nom_sub() print('Running test_expected_results_time_sub') From 2641d5d6140d63a47f32f3d9777a58c73f673e0a Mon Sep 17 00:00:00 2001 From: Kevin Ludwick Date: Mon, 10 Feb 2025 18:47:55 -0600 Subject: [PATCH 23/26] fixed unit test With the addition of simulation of smearing, the tolerance needed to be increased slightly --- tests/e2e_tests/photon_count_e2e.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/e2e_tests/photon_count_e2e.py b/tests/e2e_tests/photon_count_e2e.py index 83e04a4e..09e82842 100644 --- a/tests/e2e_tests/photon_count_e2e.py +++ b/tests/e2e_tests/photon_count_e2e.py @@ -142,7 +142,7 @@ def test_expected_results_e2e(tvacdata_path, e2eoutput_path): pc_dark_frame_err = fits.getdata(pc_dark_file, 'ERR') # more frames gets a better agreement; agreement to 1% for ~160 darks and illuminated - assert np.isclose(np.nanmean(pc_frame), ill_mean - dark_mean, rtol=0.01) + assert np.isclose(np.nanmean(pc_frame), ill_mean - dark_mean, rtol=0.02) assert np.isclose(np.nanmean(pc_dark_frame), dark_mean, rtol=0.01) assert pc_frame_err.min() >= 0 assert pc_dark_frame_err.min() >= 0 From 7a81d1c872487caa512c92510d3db8123c6b151d Mon Sep 17 00:00:00 2001 From: Kevin Ludwick Date: Fri, 14 Feb 2025 08:57:02 -0600 Subject: [PATCH 24/26] addressed Jason's comments --- corgidrp/calibrate_kgain.py | 8 +++++--- corgidrp/calibrate_nonlin.py | 8 +++++--- corgidrp/mocks.py | 4 ++++ corgidrp/walker.py | 20 +++++++++++--------- tests/e2e_tests/photon_count_e2e.py | 8 +++++--- 5 files changed, 30 insertions(+), 18 deletions(-) diff --git a/corgidrp/calibrate_kgain.py b/corgidrp/calibrate_kgain.py index ae4e7a5b..d4212164 100644 --- a/corgidrp/calibrate_kgain.py +++ b/corgidrp/calibrate_kgain.py @@ -277,7 +277,7 @@ def calibrate_kgain(dataset_kgain, n_cal=10, n_mean=30, min_val=800, max_val=3000, binwidth=68, make_plot=True,plot_outdir='figures', show_plot=False, logspace_start=-1, logspace_stop=4, logspace_num=200, - verbose=False, detector_regions=None, kgain_params=kgain_params_default): + verbose=False, detector_regions=None, kgain_params=None): """ kgain (e-/DN) is calculated from the means and variances within the defined minimum and maximum mean values. A photon transfer curve @@ -351,14 +351,16 @@ def calibrate_kgain(dataset_kgain, 'rn_bins2': upper bound of counts histogram for fitting or read noise 'max_DN_val': maximum DN value to be included in photon transfer curve (PTC) 'signal_bins_N': number of bins in the signal variables of PTC curve - Defaults to kgain_params included in this file. + Defaults to kgain_params_default included in this file. Returns: corgidrp.data.KGain: kgain estimate from the least-squares fit to the photon transfer curve (in e-/DN). The expected value of kgain for EXCAM with flight readout sequence should be between 8 and 9 e-/DN """ - + if kgain_params is None: + kgain_params = kgain_params_default + check_kgain_params(kgain_params) if detector_regions is None: diff --git a/corgidrp/calibrate_nonlin.py b/corgidrp/calibrate_nonlin.py index 020189f8..c5a6b2dc 100644 --- a/corgidrp/calibrate_nonlin.py +++ b/corgidrp/calibrate_nonlin.py @@ -127,7 +127,7 @@ def calibrate_nonlin(dataset_nl, pfit_upp_cutoff1 = -2, pfit_upp_cutoff2 = -3, pfit_low_cutoff1 = 2, pfit_low_cutoff2 = 1, make_plot=True, plot_outdir='figures', show_plot=False, - verbose=False, nonlin_params=nonlin_params_default): + verbose=False, nonlin_params=None): """ Function that derives the non-linearity calibration table for a set of DN and EM values. @@ -222,7 +222,7 @@ def calibrate_nonlin(dataset_nl, where a '1' ending indicates the smaller of two values, and a '2' ending indicates the larger of two values. The coordinates of each square are specified by matching up as follows: (rowroi1, colroi1), (rowroi1, colroi2), (rowback11, colback11), - (rowback11, colback12), etc. Defaults to nonlin_params specified in this file. + (rowback11, colback12), etc. Defaults to nonlin_params_default specified in this file. Returns: nonlin_arr (NonLinearityCalibration): 2-D array with nonlinearity values @@ -230,7 +230,9 @@ def calibrate_nonlin(dataset_nl, input signal in DN is the first column. Signal values start with min_write and run through max_write in steps of 20 DN. """ - + if nonlin_params is None: + nonlin_params = nonlin_params_default + check_nonlin_params(nonlin_params) # dataset_nl.all_data must be 3-D diff --git a/corgidrp/mocks.py b/corgidrp/mocks.py index f4532d93..546e2f6c 100644 --- a/corgidrp/mocks.py +++ b/corgidrp/mocks.py @@ -1950,6 +1950,8 @@ def create_photon_countable_frames(Nbrights=30, Ndarks=40, EMgain=5000, kgain=7, frame = data.Image(frame_dn, pri_hdr=prihdr, ext_hdr=exthdr) frame.ext_hdr['CMDGAIN'] = EMgain frame.ext_hdr['EXPTIME'] = exptime + frame.ext_hdr['KGAIN'] = kgain + frame.ext_hdr['ISPC'] = True frame.pri_hdr["VISTYPE"] = "TDEMO" frame.filename = 'L1_for_pc_ill_{0}.fits'.format(i) frame_e_list.append(frame) @@ -1963,6 +1965,8 @@ def create_photon_countable_frames(Nbrights=30, Ndarks=40, EMgain=5000, kgain=7, frame_dark = data.Image(frame_dn_dark, pri_hdr=prihdr.copy(), ext_hdr=exthdr.copy()) frame_dark.ext_hdr['CMDGAIN'] = EMgain frame_dark.ext_hdr['EXPTIME'] = exptime + frame_dark.ext_hdr['KGAIN'] = kgain + frame_dark.ext_hdr['ISPC'] = True frame_dark.pri_hdr["VISTYPE"] = "DARK" frame.filename = 'L1_for_pc_dark_{0}.fits'.format(i) frame_e_dark_list.append(frame_dark) diff --git a/corgidrp/walker.py b/corgidrp/walker.py index 11d9a649..ffb1da40 100644 --- a/corgidrp/walker.py +++ b/corgidrp/walker.py @@ -250,14 +250,15 @@ def guess_template(dataset): recipe_filename = "l1_to_l2a_noisemap.json" elif len(unique_vals) == 1: # traditional darks recipe_filename = "build_trad_dark_image.json" - elif image.pri_hdr['VISTYPE'] == "PUPILIMG": recipe_filename = ["l1_to_l2a_nonlin.json", "l1_to_kgain.json"] else: - recipe_filename = "l1_to_l2b.json" - if 'ISPC' in image.ext_hdr: - if image.ext_hdr['ISPC'] and image.pri_hdr['VISTYPE'] != "DARK": - recipe_filename = "l1_to_l2b_pc.json" + if 'ISPC' in image.ext_hdr: + if image.ext_hdr['ISPC']: + recipe_filename = "l1_to_l2b_pc.json" + else: + recipe_filename = "l1_to_l2b.json" + elif image.ext_hdr['DATA_LEVEL'] == "L2a": if image.pri_hdr['VISTYPE'] == "DARK": _, unique_vals = dataset.split_dataset(exthdr_keywords=['EXPTIME', 'CMDGAIN', 'KGAIN']) @@ -273,11 +274,12 @@ def guess_template(dataset): recipe_filename = "l2a_to_l2a_noisemap.json" elif len(unique_vals) == 1: # traditional darks recipe_filename = "l2a_build_trad_dark_image.json" - if 'ISPC' in image.ext_hdr: - if image.ext_hdr['ISPC'] and image.pri_hdr['VISTYPE'] != "DARK": - recipe_filename = "l2a_to_l2b_pc.json" else: - raise NotImplementedError() + if 'ISPC' in image.ext_hdr: + if image.ext_hdr['ISPC']: + recipe_filename = "l2a_to_l2b_pc.json" + else: + raise NotImplementedError() else: raise NotImplementedError() diff --git a/tests/e2e_tests/photon_count_e2e.py b/tests/e2e_tests/photon_count_e2e.py index 09e82842..622ea811 100644 --- a/tests/e2e_tests/photon_count_e2e.py +++ b/tests/e2e_tests/photon_count_e2e.py @@ -21,7 +21,7 @@ def test_expected_results_e2e(tvacdata_path, e2eoutput_path): bp_path = os.path.join(processed_cal_path, "bad_pix.fits") np.random.seed(1234) - ill_dataset, dark_dataset, ill_mean, dark_mean = mocks.create_photon_countable_frames(Nbrights=160, Ndarks=161, cosmic_rate=1, flux=0.5) + ill_dataset, dark_dataset, ill_mean, dark_mean = mocks.create_photon_countable_frames(Nbrights=2, Ndarks=2, cosmic_rate=1, flux=0.5)#(Nbrights=160, Ndarks=161, cosmic_rate=1, flux=0.5) output_dir = os.path.join(e2eoutput_path, 'pc_sim_test_data') output_ill_dir = os.path.join(output_dir, 'ill_frames') output_dark_dir = os.path.join(output_dir, 'dark_frames') @@ -122,7 +122,8 @@ def test_expected_results_e2e(tvacdata_path, e2eoutput_path): this_caldb.create_entry(bp_map) # make PC dark - walker.walk_corgidrp(l1_data_dark_filelist, '', output_dir, template="l1_to_l2b_pc_dark.json") + # below I leave out the template specification to check that the walker recipe guesser works as expected + walker.walk_corgidrp(l1_data_dark_filelist, '', output_dir)#, template="l1_to_l2b_pc_dark.json") for f in os.listdir(output_dir): if f.endswith('_pc_dark.fits'): pc_dark_filename = f @@ -130,7 +131,8 @@ def test_expected_results_e2e(tvacdata_path, e2eoutput_path): pc_dark_file = os.path.join(output_dir, pc_dark_filename) # make PC illuminated, subtracting the PC dark - walker.walk_corgidrp(l1_data_ill_filelist, '', output_dir, template="l1_to_l2b_pc.json") + # below I leave out the template specification to check that the walker recipe guesser works as expected + walker.walk_corgidrp(l1_data_ill_filelist, '', output_dir)#, template="l1_to_l2b_pc.json") # get photon-counted frame for f in os.listdir(output_dir): if f.endswith('_pc.fits'): From 9ed4ade6f00cf32faa24f479c229033a789ac61d Mon Sep 17 00:00:00 2001 From: Kevin Ludwick Date: Fri, 14 Feb 2025 10:16:26 -0600 Subject: [PATCH 25/26] addressed comments also set the number of frames in photon_count_e2e.py back to what it should be. (I set the number of frames to 2 for quick testing.) --- corgidrp/walker.py | 21 +++++++++------------ tests/e2e_tests/photon_count_e2e.py | 2 +- 2 files changed, 10 insertions(+), 13 deletions(-) diff --git a/corgidrp/walker.py b/corgidrp/walker.py index ffb1da40..6da487b3 100644 --- a/corgidrp/walker.py +++ b/corgidrp/walker.py @@ -243,22 +243,20 @@ def guess_template(dataset): recipe_filename = "l1_to_l2b_pc_dark.json" elif len(unique_vals) > 1: # darks for noisemap creation recipe_filename = "l1_to_l2a_noisemap.json" - elif len(unique_vals) == 1: # traditional darks + else: # then len(unique_vals) is 1 and not PC: traditional darks recipe_filename = "build_trad_dark_image.json" else: # can't distinguish PC or analog then; default to analog if len(unique_vals) > 1: # darks for noisemap creation recipe_filename = "l1_to_l2a_noisemap.json" - elif len(unique_vals) == 1: # traditional darks + else: # then len(unique_vals) is 1: traditional darks recipe_filename = "build_trad_dark_image.json" elif image.pri_hdr['VISTYPE'] == "PUPILIMG": recipe_filename = ["l1_to_l2a_nonlin.json", "l1_to_kgain.json"] else: - if 'ISPC' in image.ext_hdr: - if image.ext_hdr['ISPC']: - recipe_filename = "l1_to_l2b_pc.json" - else: - recipe_filename = "l1_to_l2b.json" - + if ('ISPC' in image.ext_hdr) & (image.ext_hdr['ISPC']): + recipe_filename = "l1_to_l2b_pc.json" + else: + recipe_filename = "l1_to_l2b.json" elif image.ext_hdr['DATA_LEVEL'] == "L2a": if image.pri_hdr['VISTYPE'] == "DARK": _, unique_vals = dataset.split_dataset(exthdr_keywords=['EXPTIME', 'CMDGAIN', 'KGAIN']) @@ -267,16 +265,15 @@ def guess_template(dataset): recipe_filename = "l2a_to_l2b_pc_dark.json" elif len(unique_vals) > 1: # darks for noisemap creation recipe_filename = "l2a_to_l2a_noisemap.json" - elif len(unique_vals) == 1: # traditional darks + else: # then len(unique_vals) is 1 and not PC: traditional darks recipe_filename = "l2a_build_trad_dark_image.json" else: # can't distinguish PC or analog then; default to analog if len(unique_vals) > 1: # darks for noisemap creation recipe_filename = "l2a_to_l2a_noisemap.json" - elif len(unique_vals) == 1: # traditional darks + else: # then len(unique_vals) is 1: traditional darks recipe_filename = "l2a_build_trad_dark_image.json" else: - if 'ISPC' in image.ext_hdr: - if image.ext_hdr['ISPC']: + if ('ISPC' in image.ext_hdr) & (image.ext_hdr['ISPC']): recipe_filename = "l2a_to_l2b_pc.json" else: raise NotImplementedError() diff --git a/tests/e2e_tests/photon_count_e2e.py b/tests/e2e_tests/photon_count_e2e.py index 622ea811..f507f695 100644 --- a/tests/e2e_tests/photon_count_e2e.py +++ b/tests/e2e_tests/photon_count_e2e.py @@ -21,7 +21,7 @@ def test_expected_results_e2e(tvacdata_path, e2eoutput_path): bp_path = os.path.join(processed_cal_path, "bad_pix.fits") np.random.seed(1234) - ill_dataset, dark_dataset, ill_mean, dark_mean = mocks.create_photon_countable_frames(Nbrights=2, Ndarks=2, cosmic_rate=1, flux=0.5)#(Nbrights=160, Ndarks=161, cosmic_rate=1, flux=0.5) + ill_dataset, dark_dataset, ill_mean, dark_mean = mocks.create_photon_countable_frames(Nbrights=160, Ndarks=161, cosmic_rate=1, flux=0.5) output_dir = os.path.join(e2eoutput_path, 'pc_sim_test_data') output_ill_dir = os.path.join(output_dir, 'ill_frames') output_dark_dir = os.path.join(output_dir, 'dark_frames') From af86a28cbac67bb385c501df9ee5426cebd57747 Mon Sep 17 00:00:00 2001 From: Kevin Ludwick Date: Fri, 14 Feb 2025 10:35:00 -0600 Subject: [PATCH 26/26] changed logic of walker back --- corgidrp/walker.py | 12 +++++++++--- 1 file changed, 9 insertions(+), 3 deletions(-) diff --git a/corgidrp/walker.py b/corgidrp/walker.py index 6da487b3..6ca02469 100644 --- a/corgidrp/walker.py +++ b/corgidrp/walker.py @@ -253,8 +253,11 @@ def guess_template(dataset): elif image.pri_hdr['VISTYPE'] == "PUPILIMG": recipe_filename = ["l1_to_l2a_nonlin.json", "l1_to_kgain.json"] else: - if ('ISPC' in image.ext_hdr) & (image.ext_hdr['ISPC']): - recipe_filename = "l1_to_l2b_pc.json" + if 'ISPC' in image.ext_hdr: + if image.ext_hdr['ISPC']: + recipe_filename = "l1_to_l2b_pc.json" + else: + recipe_filename = "l1_to_l2b.json" else: recipe_filename = "l1_to_l2b.json" elif image.ext_hdr['DATA_LEVEL'] == "L2a": @@ -273,8 +276,11 @@ def guess_template(dataset): else: # then len(unique_vals) is 1: traditional darks recipe_filename = "l2a_build_trad_dark_image.json" else: - if ('ISPC' in image.ext_hdr) & (image.ext_hdr['ISPC']): + if 'ISPC' in image.ext_hdr: + if image.ext_hdr['ISPC']: recipe_filename = "l2a_to_l2b_pc.json" + else: + raise NotImplementedError() else: raise NotImplementedError() else: