From 62a928b67b06a8eea296f2a84e9f5e24e488b8f7 Mon Sep 17 00:00:00 2001 From: Kevin Ludwick Date: Wed, 20 Nov 2024 23:26:42 -0600 Subject: [PATCH 01/11] 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/11] 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/11] 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/11] 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/11] 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/11] 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/11] 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/11] 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/11] 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/11] 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/11] 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