diff --git a/CHANGES.rst b/CHANGES.rst index b22ce3e2..cf7bc1f8 100644 --- a/CHANGES.rst +++ b/CHANGES.rst @@ -1,7 +1,64 @@ -1.5.3 (unreleased) +1.6.2 (unreleased) ================== -- +Changes to API +-------------- + +Bug Fixes +--------- + +Other +----- + +1.6.1 (2024-02-29) +================== + +Changes to API +-------------- + +ramp_fitting +~~~~~~~~~~~~ + +- Add ``average_dark_current`` to calculations of poisson variance. [#243] + +1.6.0 (2024-02-15) +================== + +Changes to API +-------------- + +jump +~~~~ + +- Add in the flagging of groups in the integration after a snowball + occurs. The saturated core of the snowball gets flagged as jump + for a number of groups passed in as a parameter [#238] + +Bug Fixes +--------- + +jump +~~~~ + +- Fixed the computation of the number of rows per slice for multiprocessing, which + was causing different results when running the step with multiprocess [#239] + +- Fix the code to at least always flag the group with the shower and the requested + groups after the primary shower. [#237] + +Other +----- + +jump +~~~~ + +- Reorganize jump docs between the jwst and stcal repos. [#240] + +ramp_fitting +~~~~~~~~~~~~ + +- Reorganize ramp_fitting docs between the jwst and stcal repos. [#240] + 1.5.2 (2023-12-13) ================== @@ -24,7 +81,7 @@ Other - Enable automatic linting and code style checks [#187] ramp_fitting ------------- +~~~~~~~~~~~~ - Refactor Casertano, et.al, 2022 uneven ramp fitting and incorporate the matching jump detection algorithm into it. [#215] @@ -84,13 +141,16 @@ jump within a group. [#207] - Added more allowable selections for the number of cores to use for - multiprocessing [#183]. + multiprocessing [#183] + +- Fixed the computation of the number of rows per slice for multiprocessing, + which caused different results when running the step with multiprocess [#239] ramp_fitting ~~~~~~~~~~~~ - Added more allowable selections for the number of cores to use for - multiprocessing [#183]. + multiprocessing [#183] - Updating variance computation for invalid integrations, as well as updating the median rate computation by excluding groups marked as diff --git a/docs/stcal/jump/description.rst b/docs/stcal/jump/description.rst index c2d4298f..ac33721e 100644 --- a/docs/stcal/jump/description.rst +++ b/docs/stcal/jump/description.rst @@ -1,6 +1,8 @@ +.. _jump_algorithm: + Algorithm --------- -This routine detects jumps in an exposure by looking for outliers +This routine detects jumps by looking for outliers in the up-the-ramp signal for each pixel in each integration within an input exposure. On output, the GROUPDQ array is updated with the DQ flag "JUMP_DET" to indicate the location of each jump that was found. @@ -10,34 +12,39 @@ output PIXELDQ array. The SCI and ERR arrays of the input data are not modified. The current implementation uses the two-point difference method described -in Anderson&Gordon2011_. +in `Anderson & Gordon (2011) `_. Two-Point Difference Method ^^^^^^^^^^^^^^^^^^^^^^^^^^^ The two-point difference method is applied to each integration as follows: -* Compute the first differences for each pixel (the difference between - adjacent groups) -* Compute the clipped (dropping the largest difference) median of the first differences for each pixel. -* Use the median to estimate the Poisson noise for each group and combine it - with the read noise to arrive at an estimate of the total expected noise for - each difference. -* Compute the "difference ratio" as the difference between the first differences - of each group and the median, divided by the expected noise. -* If the largest "difference ratio" is greater than the rejection threshold, - flag the group corresponding to that ratio as having a jump. -* If a jump is found in a given pixel, iterate the above steps with the - jump-impacted group excluded, looking for additional lower-level jumps - that still exceed the rejection threshold. -* Stop iterating on a given pixel when no new jumps are found or only one - difference remains. -* If the there are only three differences (four groups), the standard median - is used rather than the clipped median. -* If there are only two differences (three groups), the smallest one is compared to the larger - one and if the larger one is above a threshold, it is flagged as a jump. +#. Compute the first differences for each pixel (the difference between + adjacent groups) +#. Compute the clipped median (dropping the largest difference) of the first differences for each pixel. + If there are only three first difference values (four groups), no clipping is + performed when computing the median. +#. Use the median to estimate the Poisson noise for each group and combine it + with the read noise to arrive at an estimate of the total expected noise for + each difference. +#. Compute the "difference ratio" as the difference between the first differences + of each group and the median, divided by the expected noise. +#. If the largest "difference ratio" is greater than the rejection threshold, + flag the group corresponding to that ratio as having a jump. +#. If a jump is found in a given pixel, iterate the above steps with the + jump-impacted group excluded, looking for additional lower-level jumps + that still exceed the rejection threshold. +#. Stop iterating on a given pixel when no new jumps are found or only one + difference remains. +#. If there are only two differences (three groups), the smallest one is compared to the larger + one and if the larger one is above a threshold, it is flagged as a jump. +#. If flagging of the 4 neighbors is requested, then the 4 adjacent pixels will + have ramp jumps flagged in the same group as the central pixel as long as it has + a jump between the min and max requested levels for this option. +#. If flagging of groups after a ramp jump is requested, then the groups in the + requested time since a detected ramp jump will be flagged as ramp jumps if + the ramp jump is above the requested threshold. Two thresholds and times are + possible for this option. -Note that any ramp values flagged as SATURATED in the input GROUPDQ array +Note that any ramp groups flagged as SATURATED in the input GROUPDQ array are not used in any of the above calculations and hence will never be marked as containing a jump. - -.. _Anderson&Gordon2011: https://ui.adsabs.harvard.edu/abs/2011PASP..123.1237A diff --git a/docs/stcal/ramp_fitting/description.rst b/docs/stcal/ramp_fitting/description.rst index d0d12c88..9ea6830d 100644 --- a/docs/stcal/ramp_fitting/description.rst +++ b/docs/stcal/ramp_fitting/description.rst @@ -1,22 +1,17 @@ Description -============ +=========== This step determines the mean count rate, in units of counts per second, for each pixel by performing a linear fit to the data in the input file. The fit is done using the "ordinary least squares" method. -The fit is performed independently for each pixel. There can be up to three -output files created by the step. The primary output file ("rate") contains the -slope at each pixel averaged over all integrations. -Slope images from each integration are stored as a data cube in a second output -data product ("rateints"). -A third, optional output product is also available, containing detailed fit -information for each pixel. The three types of output files are described in -more detail below. +The fit is performed independently for each pixel. The count rate for each pixel is determined by a linear fit to the cosmic-ray-free and saturation-free ramp intervals for each pixel; hereafter this interval will be referred to as a "segment." The fitting algorithm uses an -'optimal' weighting scheme, as described by Fixsen et al, PASP, 112, 1350. +'optimal' weighting scheme, as described by +`Fixsen et al. (2011) `_. + Segments are determined using the 4-D GROUPDQ array of the input data set, under the assumption that the jump step will have already flagged CR's. Segments are terminated where @@ -24,113 +19,107 @@ saturation flags are found. Pixels are processed simultaneously in blocks using the array-based functionality of numpy. The size of the block depends on the image size and the number of groups. -Multiprocessing -=============== -This step has the option of running in multiprocessing mode. In that mode it will -split the input data cube into a number of row slices based on the number of available -cores on the host computer and the value of the max_cores input parameter. By -default the step runs on a single processor. At the other extreme if max_cores is -set to 'all', it will use all available cores (real and virtual). Testing has shown -a reduction in the elapsed time for the step proportional to the number of real -cores used. Using the virtual cores also reduces the elapsed time but at a slightly -lower rate than the real cores. Since the data is sliced based on the number -of rows, if the number of cores requested for multiprocessing is greater than -the number of rows, the number of cores actually used will be no more than the -number of rows. This prevents any additional cores from operating on empty -datasets, which would cause errors during ramp fitting. +.. _ramp_output_products: -Special Cases -+++++++++++++ +Output Products +--------------- -If the input dataset has only a single group in each integration, the count rate -for all unsaturated pixels in that integration will be calculated as the -value of the science data in that group divided by the group time. If the -input dataset has only two groups per integration, the count rate for all -unsaturated pixels in each integration will be calculated using the differences -between the two valid groups of the science data. +There are two output products created by default, with a third optional +product also available: -For datasets having more than a single group in each integration, a ramp having -a segment with only a single group is processed differently depending on the -number and size of the other segments in the ramp. If a ramp has only one -segment and that segment contains a single group, the count rate will be calculated -to be the value of the science data in that group divided by the group time. If a ramp -has a segment having a single group, and at least one other segment having more -than one good group, only data from the segment(s) having more than a single -good group will be used to calculate the count rate. +#. The primary output file ("rate") contains slope and variance/error + estimates for each pixel that are the result of averaging over all + integrations in the exposure. This is a product with 2-D data arrays. +#. The secondary product ("rateints") contains slope and variance/error + estimates for each pixel on a per-integration basis, stored as 3-D + data cubes. +#. The third, optional, output product contains detailed + fit information for every ramp segment for each pixel. -The data are checked for ramps in which there is good data in the first group, -but all first differences for the ramp are undefined because the remainder of -the groups are either saturated or affected by cosmic rays. For such ramps, -the first differences will be set to equal the data in the first group. The -first difference is used to estimate the slope of the ramp, as explained in the -'segment-specific computations' section below. - -If any input dataset contains ramps saturated in their second group, the count -rates for those pixels in that integration will be calculated as the value -of the science data in the first group divided by the group time. - -The MIRI first frame correction step flags all pixels in the first group of -each integration, so that those data do not get used in either the jump detection -or ramp fitting steps. -Similarly, the MIRI last frame correction step flags all pixels in the last -group of each integration. -The ramp fitting will only fit data if there are at least 2 good groups -of data and will log a warning otherwise. - -All Cases -+++++++++ -For all input datasets, including the special cases described above, arrays for -the primary output (rate) product are computed as follows. - -After computing the slopes for all segments for a given pixel, the final slope is +RATE Product +++++++++++++ +After computing the slopes and variances for all segments for a given pixel, the final slope is determined as a weighted average from all segments in all integrations, and is -written as the primary output product. In this output product, the +written to the "rate" output product. In this output product, the 4-D GROUPDQ from all integrations is collapsed into 2-D, merged (using a bitwise OR) with the input 2-D PIXELDQ, and stored as a 2-D DQ array. The 3-D VAR_POISSON and VAR_RNOISE arrays from all integrations are averaged -into corresponding 2-D output arrays. There is a case where the median rate -for a pixel can be computed as negative. This value is used in the numerator -when computing the VAR_POISSON. If the median rate is negative, the VAR_POISSON -is computed as negative, which is nonsnse. In this case, the VAR_POISSON is -set to zero for all output products. - -The slope images for each integration are stored as a data cube in a second output data -product (rateints). Each plane of the 3-D SCI, ERR, DQ, VAR_POISSON, and VAR_RNOISE +into corresponding 2-D output arrays. In cases where the median rate +for a pixel is negative, the VAR_POISSON is set to zero, in order to avoid the +unphysical situation of having a negative variance. + +RATEINTS Product +++++++++++++++++ +The slope images for each integration are stored as a data cube in "rateints" output data +product. Each plane of the 3-D SCI, ERR, DQ, VAR_POISSON, and VAR_RNOISE arrays in this product corresponds to the result for a given integration. In this output -product, the GROUPDQ data for a given integration is collapsed into 2-D, which -is then merged with the input 2-D PIXELDQ to create the output DQ array for each +product, the GROUPDQ data for a given integration is collapsed into 2-D and then +merged with the input 2-D PIXELDQ array to create the output DQ array for each integration. The 3-D VAR_POISSON and VAR_RNOISE arrays are calculated by averaging over the fit segments in the corresponding 4-D variance arrays. +FITOPT Product +++++++++++++++ A third, optional output product is also available and is produced only when -the step parameter 'save_opt' is True (the default is False). This optional +the step parameter ``save_opt`` is True (the default is False). This optional product contains 4-D arrays called SLOPE, SIGSLOPE, YINT, SIGYINT, WEIGHTS, -VAR_POISSON, and VAR_RNOISE that contain the slopes, uncertainties in the -slopes, y-intercept, uncertainty in the y-intercept, fitting weights, the -variance of the slope due to poisson noise only, and the variance of the slope -due to read noise only for each segment of each pixel, respectively. The y-intercept refers +VAR_POISSON, and VAR_RNOISE, which contain the slopes, uncertainties in the +slopes, y-intercept, uncertainty in the y-intercept, fitting weights, +variance of the slope due to poisson noise, and the variance of the slope +due to read noise for each segment of each pixel, respectively. The y-intercept refers to the result of the fit at an effective exposure time of zero. This product also contains a 3-D array called PEDESTAL, which gives the signal at zero exposure time for each pixel, and the 4-D CRMAG array, which contains the magnitude of -each group that was flagged as having a CR hit. By default, the name of this -output file will have the suffix "_fitopt". +each group that was flagged as having a CR hit. + +By default, the name of this +output file will have the product type suffix "_fitopt". In this optional output product, the pedestal array is calculated for each integration by extrapolating the final slope (the weighted average of the slopes of all ramp segments in the integration) for each pixel from its value at the first group to an exposure time of zero. Any pixel that is saturated on the first group is given a pedestal value of 0. Before compression, -the cosmic ray magnitude array is equivalent to the input SCI array but with the +the cosmic-ray magnitude array is equivalent to the input SCI array but with the only nonzero values being those whose pixel locations are flagged in the input GROUPDQ as cosmic ray hits. The array is compressed, removing all groups in which all the values are 0 for pixels having at least one group with a non-zero magnitude. The order of the cosmic rays within the ramp is preserved. +.. _ramp_special_cases: + +Special Cases +------------- +If the input dataset has only one group in each integration (NGROUPS=1), the count rate +for all unsaturated pixels in each integration will be calculated as the +value of the science data in the one group divided by the group time. If the +input dataset has only two groups per integration (NGROUPS=2), the count rate for all +unsaturated pixels in each integration will be calculated using the differences +between the two valid groups of the science data divided by the group time. + +For datasets having more than one group in each integration (NGROUPS>1), a ramp having +a segment with only one good group is processed differently depending on the +number and size of the other segments in the ramp. If a ramp has only one +segment and that segment contains a single group, the count rate will be calculated +to be the value of the science data in that group divided by the group time. If a ramp +has a segment with only one good group, and at least one other segment having more +than one good group, only data from the segment(s) having more than one +good group will be used to calculate the count rate. + +For ramps in a given integration that are saturated beginning in their second group, +the count rate for that integration will be calculated as the value of the science data +in the first group divided by the group time, but only if the step parameter +``suppress_one_group`` is set to ``False``. If set to ``True``, the computation of +slopes for pixels that have only one good group will be suppressed and the slope +for that integration will be set to zero. + +.. _ramp_slopes_and_variances: + Slope and Variance Calculations -+++++++++++++++++++++++++++++++ +------------------------------- Slopes and their variances are calculated for each segment, for each integration, and for the entire exposure. As defined above, a segment is a set of contiguous -groups where none of the groups are saturated or cosmic ray-affected. The +groups where none of the groups is saturated or cosmic ray-affected. The appropriate slopes and variances are output to the primary output product, the integration-specific output product, and the optional output product. The following is a description of these computations. The notation in the equations @@ -147,10 +136,12 @@ the rate product will be set to NaN. An example of invalid data would be a fully saturated integration for a pixel. Optimal Weighting Algorithm ---------------------------- ++++++++++++++++++++++++++++ The slope of each segment is calculated using the least-squares method with optimal -weighting, as described by Fixsen et al. 2000, PASP, 112, 1350; Regan 2007, -JWST-STScI-001212. Optimal weighting determines the relative weighting of each sample +weighting, as described by +`Fixsen et al 2000 `_ +and Regan 2007, JWST-STScI-001212. +Optimal weighting determines the relative weighting of each sample when calculating the least-squares fit to the ramp. When the data have low signal-to-noise ratio :math:`S`, the data are read noise dominated and equal weighting of samples is the best approach. In the high signal-to-noise regime, data are Poisson-noise dominated and @@ -190,8 +181,8 @@ sufficient; they are given as: | 100 | | 10 | +-------------------+------------------------+----------+ -Segment-specific Computations: ------------------------------- +Segment-specific Computations ++++++++++++++++++++++++++++++ The variance of the slope of a segment due to read noise is: .. math:: @@ -204,14 +195,16 @@ time in seconds (from the keyword TGROUP). The variance of the slope in a segment due to Poisson noise is: .. math:: - var^P_{s} = \frac{ slope_{est} }{ tgroup \times gain\ (ngroups_{s} -1)} \,, + var^P_{s} = \frac{ slope_{est} + darkcurrent}{ tgroup \times gain\ (ngroups_{s} -1)} \,, where :math:`gain` is the gain for the pixel (from the GAIN reference file), in e/DN. The :math:`slope_{est}` is an overall estimated slope of the pixel, calculated by taking the median of the first differences of the groups that are unaffected by saturation and cosmic rays, in all integrations. This is a more robust estimate of the slope than the segment-specific slope, which may be noisy -for short segments. +for short segments. The contributions from the dark current are added when present; +the value can be provided by the user during the `jwst.dark_current.DarkCurrentStep`, +or it can be specified in scalar or 2D array form by the dark reference file. The combined variance of the slope of a segment is the sum of the variances: @@ -219,8 +212,8 @@ The combined variance of the slope of a segment is the sum of the variances: var^C_{s} = var^R_{s} + var^P_{s} -Integration-specific computations: ----------------------------------- +Integration-specific computations ++++++++++++++++++++++++++++++++++ The variance of the slope for an integration due to read noise is: .. math:: @@ -244,8 +237,8 @@ The slope for an integration depends on the slope and the combined variance of e .. math:: slope_{i} = \frac{ \sum_{s}{ \frac{slope_{s}} {var^C_{s}}}} { \sum_{s}{ \frac{1} {var^C_{s}}}} -Exposure-level computations: ----------------------------- +Exposure-level computations ++++++++++++++++++++++++++++ The variance of the slope due to read noise depends on a sum over all integrations: @@ -262,61 +255,57 @@ The combined variance of the slope is the sum of the variances: .. math:: var^C_{o} = var^R_{o} + var^P_{o} -The square root of the combined variance is stored in the ERR array of the primary output. +The square-root of the combined variance is stored in the ERR array of the output product. The overall slope depends on the slope and the combined variance of the slope of each integration's -segments, so is a sum over integrations and segments: +segments, and hence is a sum over integrations and segments: .. math:: slope_{o} = \frac{ \sum_{i,s}{ \frac{slope_{i,s}} {var^C_{i,s}}}} { \sum_{i,s}{ \frac{1} {var^C_{i,s}}}} -Upon successful completion of this step, the status keyword S_RAMP will be set -to "COMPLETE". +.. _ramp_error_propagation: Error Propagation -================= - -Error propagation in the ramp fitting step is implemented by storing the -square-root of the exposure-level combined variance in the ERR array of the primary -output product. This combined variance of the exposure-level slope is the sum -of the variance of the slope due to the Poisson noise and the variance of the -slope due to the read noise. These two variances are also separately written -to the extensions VAR_POISSON and VAR_RNOISE in the primary output. - -At the integration-level, the variance of the per-integration slope due to -Poisson noise is written to the VAR_POISSON extension in the -integration-specific product, and the variance of the per-integration slope -due to read noise is written to the VAR_RNOISE extension. The square-root of -the combined variance of the slope due to both Poisson and read noise -is written to the ERR extension. - -For the optional output product, the variance of the slope due to the Poisson -noise of the segment-specific slope is written to the VAR_POISSON extension. -Similarly, the variance of the slope due to the read noise of the -segment-specific slope is written to the VAR_RNOISE extension. +----------------- +Error propagation in the ``ramp_fitting`` step is implemented by carrying along +the individual variances in the slope due to Poisson noise and read noise at all +levels of calculations. The total error estimate at each level is computed as +the square-root of the sum of the two variance estimates. + +In each type of output product generated by the step, the variance in the slope +due to Poisson noise is stored in the "VAR_POISSON" extension, the variance in +the slope due to read noise is stored in the "VAR_RNOISE" extension, and the +total error is stored in the "ERR" extension. In the optional output product, +these arrays contain information for every segment used in the fitting for each +pixel. In the "rateints" product they contain values for each integration, and +in the "rate" product they contain values for the exposure as a whole. + +.. _ramp_dq_propagation: Data Quality Propagation -======================== +------------------------ For a given pixel, if all groups in an integration are flagged as DO_NOT_USE or SATURATED, then that pixel will be flagged as DO_NOT_USE in the corresponding -integration in the rateints product. Note this does NOT mean that all groups +integration in the "rateints" product. Note this does NOT mean that all groups are flagged as SATURATED, nor that all groups are flagged as DO_NOT_USE. For -example, suppressed one ramp groups will be flagged as DO_NOT_USE in the -zeroeth group, but not necessarily any other group, while only groups one and -on are flagged as SATURATED. Further, only if all integrations in the rateints -product are marked as DO_NOT_USE, then the pixel will be flagged as DO_NOT_USE -in the rate product. +example, slope calculations that are suppressed due to a ramp containing only +one good group will be flagged as DO_NOT_USE in the +first group, but not necessarily any other group, while only groups two and +beyond are flagged as SATURATED. Further, only if all integrations in the "rateints" +product are flagged as DO_NOT_USE, then the pixel will be flagged as DO_NOT_USE +in the "rate" product. For a given pixel, if all groups in an integration are flagged as SATURATED, then that pixel will be flagged as SATURATED and DO_NOT_USE in the corresponding -integration in the rateints product. This is different from the above case in +integration in the "rateints" product. This is different from the above case in that this is only for all groups flagged as SATURATED, not for some combination -of DO_NOT_USE and SATURATED. Further, only if all integrations in the rateints -product are marked as SATURATED, then the pixel will be flagged as SATURATED -and DO_NOT_USE in the rate product. +of DO_NOT_USE and SATURATED. Further, only if all integrations in the "rateints" +product are flagged as SATURATED, then the pixel will be flagged as SATURATED +and DO_NOT_USE in the "rate" product. For a given pixel, if any group in an integration is flagged as JUMP_DET, then that pixel will be flagged as JUMP_DET in the corresponding integration in the -rateints product. Also, that pixel will be flagged as JUMP_DET in the rate +"rateints" product. That pixel will also be flagged as JUMP_DET in the "rate" product. + diff --git a/src/stcal/jump/jump.py b/src/stcal/jump/jump.py index d4b109d3..c2103e7c 100644 --- a/src/stcal/jump/jump.py +++ b/src/stcal/jump/jump.py @@ -2,10 +2,13 @@ import multiprocessing import time import warnings -import cv2 as cv import numpy as np -from astropy import stats -from astropy.convolution import Ring2DKernel, convolve +import cv2 as cv +import astropy.stats as stats + +from astropy.convolution import Ring2DKernel +from astropy.convolution import convolve + from . import constants from . import twopoint_difference as twopt @@ -54,6 +57,8 @@ def detect_jumps( minimum_groups=3, minimum_sigclip_groups=100, only_use_ints=True, + mask_persist_grps_next_int=True, + persist_grps_flagged=25, ): """ This is the high-level controlling routine for the jump detection process. @@ -236,8 +241,8 @@ def detect_jumps( err *= gain_2d readnoise_2d *= gain_2d # also apply to the after_jump thresholds - after_jump_flag_e1 = after_jump_flag_dn1 * gain_2d - after_jump_flag_e2 = after_jump_flag_dn2 * gain_2d + after_jump_flag_e1 = after_jump_flag_dn1 * np.nanmedian(gain_2d) + after_jump_flag_e2 = after_jump_flag_dn2 * np.nanmedian(gain_2d) # Apply the 2-point difference method as a first pass log.info("Executing two-point difference method") @@ -278,9 +283,15 @@ def detect_jumps( minimum_sigclip_groups=minimum_sigclip_groups, only_use_ints=only_use_ints, ) + # remove redundant bits in pixels that have jump flagged but were + # already flagged as do_not_use or saturated. + gdq[gdq == np.bitwise_or(dqflags['DO_NOT_USE'], dqflags['JUMP_DET'])] = \ + dqflags['DO_NOT_USE'] + gdq[gdq == np.bitwise_or(dqflags['SATURATED'], dqflags['JUMP_DET'])] = \ + dqflags['SATURATED'] # This is the flag that controls the flagging of snowballs. if expand_large_events: - total_snowballs = flag_large_events( + gdq, total_snowballs = flag_large_events( gdq, jump_flag, sat_flag, @@ -292,6 +303,8 @@ def detect_jumps( edge_size=edge_size, sat_expand=sat_expand, max_extended_radius=max_extended_radius, + mask_persist_grps_next_int=mask_persist_grps_next_int, + persist_grps_flagged=persist_grps_flagged, ) log.info("Total snowballs = %i", total_snowballs) number_extended_events = total_snowballs @@ -326,11 +339,10 @@ def detect_jumps( # must copy arrays here, find_crs will make copies but if slices # are being passed in for multiprocessing then the original gdq will be - # modified unless copied beforehand + # modified unless copied beforehand. gdq = gdq.copy() data = data.copy() copy_arrs = False # we don't need to copy arrays again in find_crs - for i in range(n_slices - 1): slices.insert( i, @@ -384,6 +396,8 @@ def detect_jumps( ) log.info("Creating %d processes for jump detection ", n_slices) pool = multiprocessing.Pool(processes=n_slices) +######### JUST FOR DEBUGGING ######################### +# pool = multiprocessing.Pool(processes=1) # Starts each slice in its own process. Starmap allows more than one # parameter to be passed. real_result = pool.starmap(twopt.find_crs, slices) @@ -439,7 +453,7 @@ def detect_jumps( # This is the flag that controls the flagging of snowballs. if expand_large_events: - total_snowballs = flag_large_events( + gdq, total_snowballs = flag_large_events( gdq, jump_flag, sat_flag, @@ -451,6 +465,8 @@ def detect_jumps( edge_size=edge_size, sat_expand=sat_expand, max_extended_radius=max_extended_radius, + mask_persist_grps_next_int=mask_persist_grps_next_int, + persist_grps_flagged=persist_grps_flagged, ) log.info("Total snowballs = %i", total_snowballs) number_extended_events = total_snowballs @@ -483,7 +499,6 @@ def detect_jumps( data /= gain_2d err /= gain_2d readnoise_2d /= gain_2d - # Return the updated data quality arrays return gdq, pdq, total_primary_crs, number_extended_events, stddev @@ -500,6 +515,8 @@ def flag_large_events( sat_expand=2, edge_size=25, max_extended_radius=200, + mask_persist_grps_next_int=True, + persist_grps_flagged=5, ): """ This routine controls the creation of expanded regions that are flagged as @@ -537,7 +554,12 @@ def flag_large_events( required for a snowball to be created max_extended_radius : int The largest radius that a snowball or shower can be extended - + mask_persist_grps_next_int : bool + The flag to turn on the extension of the flagging of the saturated cores of + snowballs. + persist_grps_flagged : int + How many groups to be flagged when the saturated cores are extended into + subsequent integrations Returns ------- total Snowballs @@ -547,8 +569,8 @@ def flag_large_events( n_showers_grp = [] total_snowballs = 0 - nints = gdq.shape[0] - ngrps = gdq.shape[1] + nints, ngrps, nrows, ncols = gdq.shape + persist_jumps = np.zeros(shape=(nints, gdq.shape[2], gdq.shape[3]), dtype=np.uint8) for integration in range(nints): for group in range(1, ngrps): current_gdq = gdq[integration, group, :, :] @@ -562,10 +584,8 @@ def flag_large_events( jump_ellipses = find_ellipses(gdq[integration, group, :, :], jump_flag, min_jump_area) if sat_required_snowball: low_threshold = edge_size - nrows = gdq.shape[2] high_threshold = max(0, nrows - edge_size) - - gdq, snowballs = make_snowballs( + gdq, snowballs, persist_jumps = make_snowballs( gdq, integration, group, @@ -576,7 +596,9 @@ def flag_large_events( min_sat_radius_extend, sat_expand, sat_flag, + jump_flag, max_extended_radius, + persist_jumps, ) else: snowballs = jump_ellipses @@ -590,17 +612,29 @@ def flag_large_events( sat_flag, jump_flag, expansion=expand_factor, + num_grps_masked=0, max_extended_radius=max_extended_radius, ) - return total_snowballs + # Test to see if the flagging of the saturated cores will be extended into the + # subsequent integrations. Persist_jumps contains all the pixels that were saturated + # in the cores of snowballs. + if mask_persist_grps_next_int: + for intg in range(1, nints): + if persist_grps_flagged >= 1: + last_grp_flagged = min(persist_grps_flagged, ngrps) + gdq[intg, 1:last_grp_flagged, :, :] = np.bitwise_or(gdq[intg, 1:last_grp_flagged, :, :], + np.repeat(persist_jumps[intg - 1, np.newaxis, :, :], + last_grp_flagged - 1, axis=0)) + return gdq, total_snowballs def extend_saturation( - cube, grp, sat_ellipses, sat_flag, min_sat_radius_extend, expansion=2, max_extended_radius=200 + cube, grp, sat_ellipses, sat_flag, jump_flag, min_sat_radius_extend, persist_jumps, + expansion=2, max_extended_radius=200 ): - ncols = cube.shape[2] - nrows = cube.shape[1] + ngroups, nrows, ncols = cube.shape image = np.zeros(shape=(nrows, ncols, 3), dtype=np.uint8) + persist_image = np.zeros(shape=(nrows, ncols, 3), dtype=np.uint8) outcube = cube.copy() for ellipse in sat_ellipses: ceny = ellipse[0][0] @@ -619,13 +653,29 @@ def extend_saturation( alpha, 0, 360, - (0, 0, 22), + (0, 0, 22), # in the RGB cube, set blue plane pixels of the ellipse to 22 -1, ) - sat_ellipse = image[:, :, 2] - saty, satx = np.where(sat_ellipse == 22) + # Create another non-extended ellipse that is used to create the + # persist_jumps for this integration. This will be used to mask groups + # in subsequent integrations. + sat_ellipse = image[:, :, 2] # extract the Blue plane of the image + saty, satx = np.where(sat_ellipse == 22) # find all the ellipse pixels in the ellipse outcube[grp:, saty, satx] = sat_flag - return outcube + persist_image = cv.ellipse( + persist_image, + (round(ceny), round(cenx)), + (round(ellipse[1][0] / 2), round(ellipse[1][1] / 2)), + alpha, + 0, + 360, + (0, 0, 22), + -1, + ) + persist_ellipse = persist_image[:, :, 2] + persist_saty, persist_satx = np.where(persist_ellipse == 22) + persist_jumps[persist_saty, persist_satx] = jump_flag + return outcube, persist_jumps def extend_ellipses( @@ -643,7 +693,8 @@ def extend_ellipses( # For a given DQ plane it will use the list of ellipses to create # expanded ellipses of pixels with # the jump flag set. - plane = gdq_cube[intg, grp, :, :] + out_gdq_cube = gdq_cube.copy() + plane = gdq_cube[intg, grp, :, :].copy() ncols = plane.shape[1] nrows = plane.shape[0] image = np.zeros(shape=(nrows, ncols, 3), dtype=np.uint8) @@ -689,12 +740,31 @@ def extend_ellipses( sat_pix = np.bitwise_and(gdq_cube[intg, flg_grp, :, :], sat_flag) saty, satx = np.where(sat_pix == sat_flag) jump_ellipse[saty, satx] = 0 - gdq_cube[intg, flg_grp, :, :] = np.bitwise_or(gdq_cube[intg, flg_grp, :, :], jump_ellipse) - return gdq_cube, num_ellipses + out_gdq_cube[intg, flg_grp, :, :] = np.bitwise_or(gdq_cube[intg, flg_grp, :, :], jump_ellipse) + diff_cube = out_gdq_cube - gdq_cube + return out_gdq_cube, num_ellipses + +def find_last_grp(grp, ngrps, num_grps_masked): + """ + Parameters + ---------- + grp : int + The location of the shower + + ngrps : int + The number of groups in the integration + + num_grps_masked : int + The requested number of groups to be flagged after the shower -def find_last_grp(grp, ngrps, num_grps_masked_after): - num_grps_masked_after += 1 - last_grp = min(grp + num_grps_masked_after, ngrps) + Returns + ------- + last_grp : int + The index of the last group to flag for the shower + + """ + num_grps_masked += 1 + last_grp = min(grp + num_grps_masked, ngrps) return last_grp def find_circles(dqplane, bitmask, min_area): @@ -730,7 +800,9 @@ def make_snowballs( min_sat_radius, expansion, sat_flag, + jump_flag, max_extended_radius, + persist_jumps, ): # This routine will create a list of snowballs (ellipses) that have the # center @@ -762,17 +834,19 @@ def make_snowballs( ): snowballs.append(jump) # extend the saturated ellipses that are larger than the min_sat_radius - gdq[integration, :, :, :] = extend_saturation( + gdq[integration, :, :, :], persist_jumps[integration, :, :] = extend_saturation( gdq[integration, :, :, :], group, sat_ellipses, sat_flag, + jump_flag, min_sat_radius, + persist_jumps[integration, :, :], expansion=expansion, max_extended_radius=max_extended_radius, ) - return gdq, snowballs + return gdq, snowballs, persist_jumps def point_inside_ellipse(point, ellipse): @@ -945,6 +1019,7 @@ def find_faint_extended( extended_emission = np.zeros(shape=(nrows, ncols), dtype=np.uint8) exty, extx = np.where(masked_smoothed_ratio > snr_threshold) extended_emission[exty, extx] = 1 + # find the contours of the extended emission contours, hierarchy = cv.findContours(extended_emission, cv.RETR_EXTERNAL, cv.CHAIN_APPROX_SIMPLE) # get the contours that are above the minimum size @@ -952,7 +1027,6 @@ def find_faint_extended( # get the minimum enclosing rectangle which is the same as the # minimum enclosing ellipse ellipses = [cv.minAreaRect(con) for con in bigcontours] - expand_by_ratio = True expansion = 1.0 plane = gdq[intg, grp, :, :] @@ -999,6 +1073,7 @@ def find_faint_extended( all_ellipses.append([intg, grp, ellipses]) # Reset the warnings filter to its original state warnings.resetwarnings() + total_showers = 0 if all_ellipses: # Now we actually do the flagging of the pixels inside showers. @@ -1009,6 +1084,7 @@ def find_faint_extended( intg = showers[0] grp = showers[1] ellipses = showers[2] + total_showers += len(ellipses) gdq, num = extend_ellipses( gdq, intg, @@ -1019,9 +1095,9 @@ def find_faint_extended( expansion=ellipse_expand, expand_by_ratio=True, num_grps_masked=num_grps_masked, - max_extended_radius=max_extended_radius, + max_extended_radius=max_extended_radius ) - return gdq, len(all_ellipses) + return gdq, total_showers def find_first_good_group(int_gdq, do_not_use): ngrps = int_gdq.shape[0] diff --git a/src/stcal/jump/twopoint_difference.py b/src/stcal/jump/twopoint_difference.py index b6d92b59..b9111870 100644 --- a/src/stcal/jump/twopoint_difference.py +++ b/src/stcal/jump/twopoint_difference.py @@ -402,19 +402,16 @@ def find_crs( flag_e_threshold = [after_jump_flag_e1, after_jump_flag_e2] flag_groups = [after_jump_flag_n1, after_jump_flag_n2] for cthres, cgroup in zip(flag_e_threshold, flag_groups): - if cgroup > 0: + if cgroup > 0 and cthres > 0: cr_intg, cr_group, cr_row, cr_col = np.where(np.bitwise_and(gdq, jump_flag)) for j in range(len(cr_group)): intg = cr_intg[j] group = cr_group[j] row = cr_row[j] col = cr_col[j] - if e_jump_4d[intg, group - 1, row, col] >= cthres[row, col]: - for kk in range(group, min(group + cgroup + 1, ngroups)): - if (gdq[intg, kk, row, col] & sat_flag) == 0 and ( - gdq[intg, kk, row, col] & dnu_flag - ) == 0: - gdq[intg, kk, row, col] = np.bitwise_or(gdq[integ, kk, row, col], jump_flag) + if e_jump_4d[intg, group - 1, row, col] >= cthres: + for kk in range(group + 1, min(group + cgroup + 1, ngroups)): + gdq[intg, kk, row, col] = np.bitwise_or(gdq[intg, kk, row, col], jump_flag) if "stddev" in locals(): return gdq, row_below_gdq, row_above_gdq, num_primary_crs, stddev diff --git a/src/stcal/ramp_fitting/ols_fit.py b/src/stcal/ramp_fitting/ols_fit.py index f2c6d5a2..a52e2fba 100644 --- a/src/stcal/ramp_fitting/ols_fit.py +++ b/src/stcal/ramp_fitting/ols_fit.py @@ -542,8 +542,9 @@ def slice_ramp_data(ramp_data, start_row, nrows): err = ramp_data.err[:, :, start_row : start_row + nrows, :].copy() groupdq = ramp_data.groupdq[:, :, start_row : start_row + nrows, :].copy() pixeldq = ramp_data.pixeldq[start_row : start_row + nrows, :].copy() + average_dark_current = ramp_data.average_dark_current[start_row : start_row + nrows, :].copy() - ramp_data_slice.set_arrays(data, err, groupdq, pixeldq) + ramp_data_slice.set_arrays(data, err, groupdq, pixeldq, average_dark_current) if ramp_data.zeroframe is not None: ramp_data_slice.zeroframe = ramp_data.zeroframe[:, start_row : start_row + nrows, :].copy() @@ -1138,7 +1139,8 @@ def ramp_fit_compute_variances(ramp_data, gain_2d, readnoise_2d, fit_slopes_ans) # Suppress harmless arithmetic warnings for now warnings.filterwarnings("ignore", ".*invalid value.*", RuntimeWarning) warnings.filterwarnings("ignore", ".*divide by zero.*", RuntimeWarning) - var_p4[num_int, :, rlo:rhi, :] = den_p3 * med_rates[rlo:rhi, :] + var_p4[num_int, :, rlo:rhi, :] = (den_p3 * med_rates[rlo:rhi, :]) + \ + ramp_data.average_dark_current[rlo:rhi, :] # Find the segment variance due to read noise and convert back to DN var_r4[num_int, :, rlo:rhi, :] = num_r3 * den_r3 / gain_sect**2 @@ -1175,10 +1177,11 @@ def ramp_fit_compute_variances(ramp_data, gain_2d, readnoise_2d, fit_slopes_ans) # Huge variances correspond to non-existing segments, so are reset to 0 # to nullify their contribution. var_p3[var_p3 > utils.LARGE_VARIANCE_THRESHOLD] = 0.0 - var_p3[:, med_rates <= 0.0] = 0.0 + med_rate_mask = med_rates <= 0.0 + var_p3[:, med_rate_mask] = 0.0 warnings.resetwarnings() - var_p4[num_int, :, med_rates <= 0.0] = 0.0 + var_p4[num_int, :, med_rate_mask] = ramp_data.average_dark_current[med_rate_mask][..., np.newaxis] var_both4[num_int, :, :, :] = var_r4[num_int, :, :, :] + var_p4[num_int, :, :, :] inv_var_both4[num_int, :, :, :] = 1.0 / var_both4[num_int, :, :, :] diff --git a/src/stcal/ramp_fitting/ramp_fit.py b/src/stcal/ramp_fitting/ramp_fit.py index 514429f1..11ba5b9e 100755 --- a/src/stcal/ramp_fitting/ramp_fit.py +++ b/src/stcal/ramp_fitting/ramp_fit.py @@ -53,10 +53,16 @@ def create_ramp_fit_class(model, dqflags=None, suppress_one_group=False): """ ramp_data = ramp_fit_class.RampData() + if not hasattr(model, 'average_dark_current'): + dark_current_array = np.zeros_like(model.pixeldq) + else: + dark_current_array = model.average_dark_current + if isinstance(model.data, u.Quantity): - ramp_data.set_arrays(model.data.value, model.err.value, model.groupdq, model.pixeldq) + ramp_data.set_arrays(model.data.value, model.err.value, model.groupdq, + model.pixeldq, dark_current_array) else: - ramp_data.set_arrays(model.data, model.err, model.groupdq, model.pixeldq) + ramp_data.set_arrays(model.data, model.err, model.groupdq, model.pixeldq, dark_current_array) # Attribute may not be supported by all pipelines. Default is NoneType. drop_frames1 = model.meta.exposure.drop_frames1 if hasattr(model, "drop_frames1") else None diff --git a/src/stcal/ramp_fitting/ramp_fit_class.py b/src/stcal/ramp_fitting/ramp_fit_class.py index f8a78efd..ef2df6bc 100644 --- a/src/stcal/ramp_fitting/ramp_fit_class.py +++ b/src/stcal/ramp_fitting/ramp_fit_class.py @@ -6,6 +6,7 @@ def __init__(self): self.err = None self.groupdq = None self.pixeldq = None + self.average_dark_current = None # Meta information self.instrument_name = None @@ -41,7 +42,7 @@ def __init__(self): self.current_integ = -1 - def set_arrays(self, data, err, groupdq, pixeldq): + def set_arrays(self, data, err, groupdq, pixeldq, average_dark_current): """ Set the arrays needed for ramp fitting. @@ -62,12 +63,17 @@ def set_arrays(self, data, err, groupdq, pixeldq): pixeldq : ndarray (uint32) 4-D array containing the pixel data quality information. It has dimensions (nintegrations, ngroups, nrows, ncols) + + average_dark_current : ndarray (float32) + 2-D array containing the average dark current. It has + dimensions (nrows, ncols) """ # Get arrays from the data model self.data = data self.err = err self.groupdq = groupdq self.pixeldq = pixeldq + self.average_dark_current = average_dark_current def set_meta(self, name, frame_time, group_time, groupgap, nframes, drop_frames1=None): """ diff --git a/tests/test_jump.py b/tests/test_jump.py index 5bca3f07..0ac49cf2 100644 --- a/tests/test_jump.py +++ b/tests/test_jump.py @@ -9,6 +9,8 @@ flag_large_events, point_inside_ellipse, find_first_good_group, + detect_jumps, + find_last_grp ) DQFLAGS = {"JUMP_DET": 4, "SATURATED": 2, "DO_NOT_USE": 1, "GOOD": 0, "NO_GAIN_VALUE": 8} @@ -31,6 +33,99 @@ def _cube(ngroups, readnoise=10): return _cube +def test_multiprocessing(): + nints = 1 + nrows = 13 + ncols = 2 + ngroups = 13 + readnoise = 10 + frames_per_group = 1 + + data = np.zeros(shape=(nints, ngroups, nrows, ncols), dtype=np.float32) + readnoise_2d = np.ones((nrows, ncols), dtype=np.float32) * readnoise + gain_2d = np.ones((nrows, ncols), dtype=np.float32) * 4 + gdq = np.zeros(shape=(nints, ngroups, nrows, ncols), dtype=np.uint32) + pdq = np.zeros(shape=(nrows, ncols), dtype=np.uint32) + err = np.zeros(shape=(nrows, ncols), dtype=np.float32) + num_cores = "1" + data[0, 4:, 5, 1] = 2000 + gdq[0, 4:, 6, 1] = DQFLAGS['DO_NOT_USE'] + gdq, pdq, total_primary_crs, number_extended_events, stddev = detect_jumps( + frames_per_group, data, gdq, pdq, err, gain_2d, readnoise_2d, rejection_thresh=5, three_grp_thresh=6, + four_grp_thresh=7, max_cores=num_cores, max_jump_to_flag_neighbors=10000, min_jump_to_flag_neighbors=100, + flag_4_neighbors=True, dqflags=DQFLAGS) + print(data[0, 4, :, :]) + print(gdq[0, 4, :, :]) + assert gdq[0, 4, 5, 1] == DQFLAGS['JUMP_DET'] + assert gdq[0, 4, 6, 1] == DQFLAGS['DO_NOT_USE'] + + # This section of code will fail without the fixes for PR #239 that prevent + # the double flagging pixels with jump which already have do_not_use or saturation set. + num_cores = "5" + data = np.zeros(shape=(nints, ngroups, nrows, ncols), dtype=np.float32) + gdq = np.zeros(shape=(nints, ngroups, nrows, ncols), dtype=np.uint32) + pdq = np.zeros(shape=(nrows, ncols), dtype=np.uint32) + readnoise_2d = np.ones((nrows, ncols), dtype=np.float32) * readnoise + gain_2d = np.ones((nrows, ncols), dtype=np.float32) * 3 + err = np.zeros(shape=(nrows, ncols), dtype=np.float32) + data[0, 4:, 5, 1] = 2000 + gdq[0, 4:, 6, 1] = DQFLAGS['DO_NOT_USE'] + gdq, pdq, total_primary_crs, number_extended_events, stddev = detect_jumps( + frames_per_group, data, gdq, pdq, err, gain_2d, readnoise_2d, rejection_thresh=5, three_grp_thresh=6, + four_grp_thresh=7, max_cores=num_cores, max_jump_to_flag_neighbors=10000, min_jump_to_flag_neighbors=100, + flag_4_neighbors=True, dqflags=DQFLAGS) + assert gdq[0, 4, 5, 1] == DQFLAGS['JUMP_DET'] + assert gdq[0, 4, 6, 1] == DQFLAGS['DO_NOT_USE'] #This value would have been 5 without the fix. + + +def test_multiprocessing_big(): + nints = 1 + nrows = 2048 + ncols = 7 + ngroups = 13 + readnoise = 10 + frames_per_group = 1 + + data = np.zeros(shape=(nints, ngroups, nrows, ncols), dtype=np.float32) + readnoise_2d = np.ones((nrows, ncols), dtype=np.float32) * readnoise + gain_2d = np.ones((nrows, ncols), dtype=np.float32) * 4 + gdq = np.zeros(shape=(nints, ngroups, nrows, ncols), dtype=np.uint32) + pdq = np.zeros(shape=(nrows, ncols), dtype=np.uint32) + err = np.zeros(shape=(nrows, ncols), dtype=np.float32) + num_cores = "1" + data[0, 4:, 204, 5] = 2000 + gdq[0, 4:, 204, 6] = DQFLAGS['DO_NOT_USE'] + gdq, pdq, total_primary_crs, number_extended_events, stddev = detect_jumps( + frames_per_group, data, gdq, pdq, err, gain_2d, readnoise_2d, rejection_thresh=5, three_grp_thresh=6, + four_grp_thresh=7, max_cores=num_cores, max_jump_to_flag_neighbors=10000, min_jump_to_flag_neighbors=100, + flag_4_neighbors=True, dqflags=DQFLAGS) + print(data[0, 4, :, :]) + print(gdq[0, 4, :, :]) + assert gdq[0, 4, 204, 5] == DQFLAGS['JUMP_DET'] + assert gdq[0, 4, 205, 5] == DQFLAGS['JUMP_DET'] + assert gdq[0, 4, 204, 6] == DQFLAGS['DO_NOT_USE'] + + # This section of code will fail without the fixes for PR #239 that prevent + # the double flagging pixels with jump which already have do_not_use or saturation set. + num_cores = "10" + data = np.zeros(shape=(nints, ngroups, nrows, ncols), dtype=np.float32) + gdq = np.zeros(shape=(nints, ngroups, nrows, ncols), dtype=np.uint32) + pdq = np.zeros(shape=(nrows, ncols), dtype=np.uint32) + readnoise_2d = np.ones((nrows, ncols), dtype=np.float32) * readnoise + gain_2d = np.ones((nrows, ncols), dtype=np.float32) * 3 + err = np.zeros(shape=(nrows, ncols), dtype=np.float32) + data[0, 4:, 204, 5] = 2000 + gdq[0, 4:, 204, 6] = DQFLAGS['DO_NOT_USE'] + gdq, pdq, total_primary_crs, number_extended_events, stddev = detect_jumps( + frames_per_group, data, gdq, pdq, err, gain_2d, readnoise_2d, rejection_thresh=5, three_grp_thresh=6, + four_grp_thresh=7, max_cores=num_cores, max_jump_to_flag_neighbors=10000, min_jump_to_flag_neighbors=100, + flag_4_neighbors=True, dqflags=DQFLAGS) + assert gdq[0, 4, 204, 5] == DQFLAGS['JUMP_DET'] + assert gdq[0, 4, 205, 5] == DQFLAGS['JUMP_DET'] + assert gdq[0, 4, 204, 6] == DQFLAGS['DO_NOT_USE'] #This value would have been 5 without the fix. + + + def test_find_simple_ellipse(): plane = np.zeros(shape=(5, 5), dtype=np.uint8) plane[2, 2] = DQFLAGS["JUMP_DET"] @@ -62,6 +157,7 @@ def test_find_ellipse2(): def test_extend_saturation_simple(): cube = np.zeros(shape=(5, 7, 7), dtype=np.uint8) + persist_jumps = np.zeros(shape=(7, 7), dtype=np.uint8) grp = 1 min_sat_radius_extend = 1 cube[1, 3, 3] = DQFLAGS["SATURATED"] @@ -71,8 +167,9 @@ def test_extend_saturation_simple(): cube[1, 3, 2] = DQFLAGS["SATURATED"] cube[1, 2, 2] = DQFLAGS["JUMP_DET"] sat_circles = find_ellipses(cube[grp, :, :], DQFLAGS["SATURATED"], 1) - new_cube = extend_saturation( - cube, grp, sat_circles, DQFLAGS["SATURATED"], min_sat_radius_extend, expansion=1.1 + new_cube, persist_jumps = extend_saturation( + cube, grp, sat_circles, DQFLAGS["SATURATED"], DQFLAGS["JUMP_DET"], + 1.1, persist_jumps, ) assert new_cube[grp, 2, 2] == DQFLAGS["SATURATED"] @@ -127,7 +224,7 @@ def test_flag_large_events_withsnowball(): cube[0, 2, 5, 1:6] = DQFLAGS["JUMP_DET"] cube[0, 2, 1:6, 1] = DQFLAGS["JUMP_DET"] cube[0, 2, 1:6, 5] = DQFLAGS["JUMP_DET"] - flag_large_events( + cube, total_snowballs = flag_large_events( cube, DQFLAGS["JUMP_DET"], DQFLAGS["SATURATED"], @@ -188,7 +285,7 @@ def test_flag_large_events_withsnowball_noextension(): cube[0, 2, 5, 1:6] = DQFLAGS["JUMP_DET"] cube[0, 2, 1:6, 1] = DQFLAGS["JUMP_DET"] cube[0, 2, 1:6, 5] = DQFLAGS["JUMP_DET"] - flag_large_events( + cube, num_snowballs = flag_large_events( cube, DQFLAGS["JUMP_DET"], DQFLAGS["SATURATED"], @@ -358,7 +455,26 @@ def test_inside_ellipes5(): result = point_inside_ellipse(point, ellipse) assert result - +@pytest.mark.skip(" used for local testing") +def test_flag_persist_groups(): +# gdq = fits.getdata("persistgdq.fits") + gdq = np.zeros(shape=(2,2,2,2)) + print(gdq.shape[0]) + gdq = gdq[:, 0:10, :, :] + total_snowballs = flag_large_events( + gdq, + DQFLAGS["JUMP_DET"], + DQFLAGS["SATURATED"], + min_sat_area=1, + min_jump_area=6, + expand_factor=1.9, + edge_size=0, + sat_required_snowball=True, + min_sat_radius_extend=2.5, + sat_expand=1.1, + mask_persist_grps_next_int=True, + persist_grps_flagged=0) +# fits.writeto("persitflaggedgdq.fits", gdq, overwrite=True) def test_calc_num_slices(): n_rows = 20 max_available_cores = 10 @@ -376,15 +492,14 @@ def test_calc_num_slices(): n_rows = 9 assert calc_num_slices(n_rows, "21", max_available_cores) == 9 -def test_find_first_good_grp(): - ngrps = 5 - ncols = 2 - nrows = 2 - intg_gdq = np.zeros(shape=(ngrps, ncols, nrows), dtype=np.uint32) - assert find_first_good_group(intg_gdq, DQFLAGS['DO_NOT_USE']) == 0 - intg_gdq[0, :, :] = 5 - assert find_first_good_group(intg_gdq, DQFLAGS['DO_NOT_USE']) == 1 - intg_gdq[1, :, :] = 5 - assert find_first_good_group(intg_gdq, DQFLAGS['DO_NOT_USE']) == 2 - intg_gdq[0, 0, 1] = 4 - assert find_first_good_group(intg_gdq, DQFLAGS['DO_NOT_USE']) == 0 + +def test_find_last_grp(): + assert (find_last_grp(grp=5, ngrps=7, num_grps_masked=0) == 6) + assert (find_last_grp(grp=5, ngrps=7, num_grps_masked=2) == 7) + assert (find_last_grp(grp=5, ngrps=7, num_grps_masked=3) == 7) + assert (find_last_grp(grp=5, ngrps=6, num_grps_masked=1) == 6) + assert (find_last_grp(grp=5, ngrps=6, num_grps_masked=0) == 6) + assert (find_last_grp(grp=5, ngrps=6, num_grps_masked=2) == 6) + assert (find_last_grp(grp=5, ngrps=8, num_grps_masked=0) == 6) + assert (find_last_grp(grp=5, ngrps=8, num_grps_masked=1) == 7) + assert (find_last_grp(grp=5, ngrps=8, num_grps_masked=2) == 8) diff --git a/tests/test_ramp_fitting.py b/tests/test_ramp_fitting.py index d8e90610..aab91a12 100644 --- a/tests/test_ramp_fitting.py +++ b/tests/test_ramp_fitting.py @@ -322,12 +322,13 @@ def jp_2326_test_setup(): gdq = np.zeros((nints, ngroups, nrows, ncols), dtype=np.uint8) err = np.zeros((nints, ngroups, nrows, ncols)) pdq = np.zeros((nrows, ncols), dtype=np.uint32) + dark_current = np.zeros((nrows, ncols)) data[0, :, 0, 0] = ramp.copy() gdq[0, :, 0, 0] = dq.copy() ramp_data = RampData() - ramp_data.set_arrays(data=data, err=err, groupdq=gdq, pixeldq=pdq) + ramp_data.set_arrays(data=data, err=err, groupdq=gdq, pixeldq=pdq, average_dark_current=dark_current) ramp_data.set_meta( name="MIRI", frame_time=2.77504, group_time=2.77504, groupgap=0, nframes=1, drop_frames1=None ) @@ -417,6 +418,7 @@ def test_2_group_cases(): rnoise = np.ones((1, npix)) * rnoise_val gain = np.ones((1, npix)) * gain_val pixeldq = np.zeros((1, npix), dtype=np.uint32) + dark_current = np.zeros((nrows, ncols), dtype=np.float32) data = np.zeros(dims, dtype=np.float32) # Science data for k in range(npix): @@ -433,7 +435,7 @@ def test_2_group_cases(): # Setup the RampData class to run ramp fitting on. ramp_data = RampData() - ramp_data.set_arrays(data, err, groupdq, pixeldq) + ramp_data.set_arrays(data, err, groupdq, pixeldq, average_dark_current=dark_current) ramp_data.set_meta( name="NIRSPEC", frame_time=14.58889, group_time=14.58889, groupgap=0, nframes=1, drop_frames1=None @@ -730,6 +732,8 @@ def create_zero_frame_data(): pixdq = np.zeros(shape=(nrows, ncols), dtype=np.uint32) gdq = np.zeros(shape=(nints, ngroups, nrows, ncols), dtype=np.uint8) zframe = np.ones(shape=(nints, nrows, ncols), dtype=np.float32) + dark_current = np.zeros((nrows, ncols), dtype=np.float32) + # Create base ramps for each pixel in each integration. base_slope = 2000.0 @@ -756,7 +760,7 @@ def create_zero_frame_data(): # Create RampData for testing. ramp_data = RampData() - ramp_data.set_arrays(data=data, err=err, groupdq=gdq, pixeldq=pixdq) + ramp_data.set_arrays(data=data, err=err, groupdq=gdq, pixeldq=pixdq, average_dark_current=dark_current) ramp_data.set_meta( name="NIRCam", frame_time=frame_time, @@ -855,6 +859,8 @@ def create_only_good_0th_group_data(): err = np.ones(shape=(nints, ngroups, nrows, ncols), dtype=np.float32) pixdq = np.zeros(shape=(nrows, ncols), dtype=np.uint32) gdq = np.zeros(shape=(nints, ngroups, nrows, ncols), dtype=np.uint8) + dark_current = np.zeros((nrows, ncols), dtype=np.float32) + # Create base ramps for each pixel in each integration. base_slope = 2000.0 @@ -877,7 +883,7 @@ def create_only_good_0th_group_data(): # Create RampData for testing. ramp_data = RampData() - ramp_data.set_arrays(data=data, err=err, groupdq=gdq, pixeldq=pixdq) + ramp_data.set_arrays(data=data, err=err, groupdq=gdq, pixeldq=pixdq, average_dark_current=dark_current) ramp_data.set_meta( name="NIRCam", frame_time=frame_time, @@ -1422,9 +1428,10 @@ def create_blank_ramp_data(dims, var, tm): err = np.ones(shape=(nints, ngroups, nrows, ncols), dtype=np.float32) pixdq = np.zeros(shape=(nrows, ncols), dtype=np.uint32) gdq = np.zeros(shape=(nints, ngroups, nrows, ncols), dtype=np.uint8) + dark_current = np.zeros(shape=(nrows, ncols), dtype = np.float32) ramp_data = RampData() - ramp_data.set_arrays(data=data, err=err, groupdq=gdq, pixeldq=pixdq) + ramp_data.set_arrays(data=data, err=err, groupdq=gdq, pixeldq=pixdq, average_dark_current=dark_current) ramp_data.set_meta( name="NIRSpec", frame_time=frame_time, @@ -1476,6 +1483,7 @@ def setup_inputs(dims, var, tm): err = np.ones(shape=(nints, ngroups, nrows, ncols), dtype=np.float32) pixdq = np.zeros(shape=(nrows, ncols), dtype=np.uint32) gdq = np.zeros(shape=(nints, ngroups, nrows, ncols), dtype=np.uint8) + dark_current = np.zeros(shape=(nrows, ncols), dtype=np.float32) base_array = np.array([k + 1 for k in range(ngroups)]) base, inc = 1.5, 1.5 @@ -1488,7 +1496,7 @@ def setup_inputs(dims, var, tm): data[c_int, :, :, :] = data[0, :, :, :].copy() ramp_data = RampData() - ramp_data.set_arrays(data=data, err=err, groupdq=gdq, pixeldq=pixdq) + ramp_data.set_arrays(data=data, err=err, groupdq=gdq, pixeldq=pixdq, average_dark_current=dark_current) ramp_data.set_meta( name="MIRI", frame_time=dtime, group_time=gtime, groupgap=0, nframes=nframes, drop_frames1=None ) diff --git a/tests/test_ramp_fitting_cases.py b/tests/test_ramp_fitting_cases.py index 675e6a75..9e32813b 100644 --- a/tests/test_ramp_fitting_cases.py +++ b/tests/test_ramp_fitting_cases.py @@ -786,9 +786,10 @@ def create_blank_ramp_data(dims, var, timing, ts_name="NIRSpec"): err = np.ones(shape=(nints, ngroups, nrows, ncols), dtype=np.float32) pixdq = np.zeros(shape=(nrows, ncols), dtype=np.uint32) gdq = np.zeros(shape=(nints, ngroups, nrows, ncols), dtype=np.uint8) + dark_current = np.zeros(shape=(nrows, ncols), dtype=np.float32) ramp_data = RampData() - ramp_data.set_arrays(data=data, err=err, groupdq=gdq, pixeldq=pixdq) + ramp_data.set_arrays(data=data, err=err, groupdq=gdq, pixeldq=pixdq, average_dark_current=dark_current) ramp_data.set_meta( name=ts_name, frame_time=frame_time, diff --git a/tests/test_ramp_fitting_gls_fit.py b/tests/test_ramp_fitting_gls_fit.py index e940bd94..fef7c0a1 100644 --- a/tests/test_ramp_fitting_gls_fit.py +++ b/tests/test_ramp_fitting_gls_fit.py @@ -63,9 +63,11 @@ def setup_inputs(dims, gain, rnoise, group_time, frame_time): err = np.ones(shape=(nints, ngroups, nrows, ncols), dtype=np.float32) groupdq = np.zeros(shape=(nints, ngroups, nrows, ncols), dtype=np.uint8) pixeldq = np.zeros(shape=(nrows, ncols), dtype=np.uint32) + dark_current = np.zeros(shape=(nrows, ncols), dtype=np.float32) + # Set clas arrays - ramp_class.set_arrays(data, err, groupdq, pixeldq) + ramp_class.set_arrays(data, err, groupdq, pixeldq, average_dark_current=dark_current) # Set class meta ramp_class.set_meta( diff --git a/tests/test_twopoint_difference.py b/tests/test_twopoint_difference.py index c6443bc7..f066b4b1 100644 --- a/tests/test_twopoint_difference.py +++ b/tests/test_twopoint_difference.py @@ -855,7 +855,7 @@ def test_10grps_1cr_afterjump(setup_cube): data[0, 8, 100, 100] = 1190 data[0, 9, 100, 100] = 1209 - after_jump_flag_e1 = np.full(data.shape[2:4], 1.0) * 0.0 + after_jump_flag_e1 = 1.0 out_gdq, row_below_gdq, rows_above_gdq, total_crs, stddev = find_crs( data, gdq, @@ -891,7 +891,7 @@ def test_10grps_1cr_afterjump_2group(setup_cube): data[0, 8, 100, 100] = 1190 data[0, 9, 100, 100] = 1209 - after_jump_flag_e1 = np.full(data.shape[2:4], 1.0) * 0.0 + after_jump_flag_e1 = 1.0 out_gdq, row_below_gdq, rows_above_gdq, total_crs, stddev = find_crs( data, gdq, @@ -932,7 +932,7 @@ def test_10grps_1cr_afterjump_toosmall(setup_cube): data[0, 8, 100, 100] = 1190 data[0, 9, 100, 100] = 1209 - after_jump_flag_e1 = np.full(data.shape[2:4], 1.0) * 10000.0 + after_jump_flag_e1 = 10000.0 out_gdq, row_below_gdq, rows_above_gdq, total_crs, stddev = find_crs( data, gdq, @@ -968,8 +968,8 @@ def test_10grps_1cr_afterjump_twothresholds(setup_cube): data[0, 8, 100, 100] = 1190 data[0, 9, 100, 100] = 1209 - after_jump_flag_e1 = np.full(data.shape[2:4], 1.0) * 500.0 - after_jump_flag_e2 = np.full(data.shape[2:4], 1.0) * 10.0 + after_jump_flag_e1 = 500.0 + after_jump_flag_e2 = 10.0 out_gdq, row_below_gdq, rows_above_gdq, total_crs, stddev = find_crs( data, gdq,