Questions from the jwst/stcal/romancal pipeline team

[emolter]

Hi there, I'm part of the jwst pipeline team. We are in the process of establishing a public API for stcal, and I noticed that your repository is using some of our functionality. I saw, though, that you aren't importing directly from stcal. Is there something we could do that would make your life easier when using our code, or any issues that are preventing you from importing stcal?

I ask because:

• We want to make sure we retain as public any methods that are actually being used by other projects like yours and
• We are continuously making updates to stcal functionality, including bugfixes and runtime improvements, that we hope will benefit external users in addition to JWST and Roman

Disclaimer: I'm reaching out on my own to better understand current usages of our repositories. We do not have any official direction from on-high regarding supporting other observatories. However it should be relatively straightforward for us to help you use our code, as it'll likely also help us organize our own code, find bugs, and reduce technical debt.

[jehturner]

Hello! Are you referring to the NIR ramp fitting (& CR jump detection) code here?

ScorpioDR/scorpiodr/scorpio/primitives_scorpio_nearIR.py

Lines 93 to 554 in 33a4476

	def calculateSignalByRegression(self, adinputs=None, **params):
	"""
	Iteratively fit a slope, intercept, and cosmic rays to a ramp.
	This function fits a ramp, possibly with discontinuities (cosmic-ray
	hits), to a 3-D data "cube" with shape (number of groups, number of
	pixels high, number of pixels wide). The fit will be done multiple
	times, with the previous fit being used to assign weights (via the
	covariance matrix) for the current fit. The iterations stop either
	when the maximum number of iterations has been reached or when the
	maximum difference between the previous fit and the current fit is
	below a cutoff. This function calls compute_slope and evaluate_fit.
	compute_slope creates arrays for the slope, intercept, and cosmic-ray
	amplitudes (the arrays that will be returned by determine_slope). Then
	it loops over the number of cosmic rays, from 0 to max_num_cr
	inclusive. Within this loop, compute_slope copies to temporary arrays
	the ramp data for all the pixels that have the current number of cosmic
	ray hits, calls gls_fit to compute the fit, then copies the results
	of the fit (slope, etc.) to the output arrays for just those pixels.
	The input to gls_fit is ramp data for a subset of pixels (nz in number)
	that all have the same number of cosmic-ray hits. gls_fit solves
	matrix equations (one for each of the nz pixels) of the form:
	y = x * p
	where y is a column vector containing the observed data values in
	electrons for each group (the length of y is ngroups, the number of
	groups); x is a matrix with ngroups rows and 2 + num_cr columns, where
	num_cr is the number of cosmic rays being included in the fit; and p
	is the solution, a column vector containing the intercept, slope, and
	the amplitude of each of the num_cr cosmic rays. The first column of
	x is all ones, for fitting to the intercept. The second column of x
	is the time (seconds) at the beginning of each group. The remaining
	num_cr columns (if num_cr > 0) are Heaviside functions, 0 for the
	first rows and 1 for all rows at and following the group containing a
	cosmic-ray hit (each row corresponds to a group). There will be one
	such column for each cosmic ray, so that the cosmic rays will be fit
	independently of each other. Whether a cosmic ray hit the detector
	during a particular group was determined by a previous step, and the
	affected groups are flagged in the group data quality array. In order
	to account for the variance of each observed value and the covariance
	between them (since they're measurements along a ramp), the solution
	is computed in this form (the @ sign represents matrix multiplication):
	(xT @ C^-1 @ x)^-1 @ [xT @ C^-1 @ y]
	where C is the ngroups by ngroups covariance matrix, ^-1 means matrix
	inverse, and xT is the transpose of x.
	Summary of the notation:
	data_sect is 3-D, (ngroups, ny, nx); this is the ramp of science data.
	cr_flagged is 3-D, (ngroups, ny, nx); 1 indicates a cosmic ray, e.g.:
	cr_flagged = np.where(np.bitwise_and(gdq_sect, jump_flag), 1, 0)
	cr_flagged_2d is 2-D, (ngroups, nz); this gives the location within
	the ramp of each cosmic ray, for the subset of pixels (nz of them)
	that have a total of num_cr cosmic ray hits at each pixel. This
	is passed to gls_fit(), which fits a slope to the ramp.
	ramp_data has shape (ngroups, nz); this will be a ramp with a 1-D
	array of pixels copied out of data_sect. The pixels will be those
	that have a particular number of cosmic-ray hits, somewhere within
	the ramp.
	Sum cr_flagged over groups to get an (ny, nx) image of the number of
	cosmic rays (i.e. accumulated over the ramp) in each pixel.
	sum_flagged = cr_flagged.sum(axis=0, ...)
	sum_flagged is used to extract the nz pixels from (ny, nx) that have a
	specified number of cosmic ray hits, e.g.:
	for num_cr in range(max_num_cr + 1):
	ncr_mask = (sum_flagged == num_cr)
	nz = ncr_mask.sum(dtype=np.int32)
	for k in range(ngroups):
	ramp_data[k] = data_sect[k][ncr_mask]
	cr_flagged_2d[k] = cr_flagged[k][ncr_mask]
	gls_fit is called for the subset of pixels (nz of them) that have
	num_cr cosmic ray hits within the ramp, the same number for every
	pixel.
	License
	-------
	Copyright (C) 2020 Association of Universities for Research in Astronomy (AURA)
	Redistribution and use in source and binary forms, with or without
	modification, are permitted provided that the following conditions are met:
	1. Redistributions of source code must retain the above copyright
	notice, this list of conditions and the following disclaimer.
	2. Redistributions in binary form must reproduce the above
	copyright notice, this list of conditions and the following
	disclaimer in the documentation and/or other materials provided
	with the distribution.
	3. The name of AURA and its representatives may not be used to
	endorse or promote products derived from this software without
	specific prior written permission.
	THIS SOFTWARE IS PROVIDED BY AURA ``AS IS'' AND ANY EXPRESS OR IMPLIED
	WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
	MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
	DISCLAIMED. IN NO EVENT SHALL AURA BE LIABLE FOR ANY DIRECT, INDIRECT,
	INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
	BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS
	OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
	ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR
	TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE
	USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH
	DAMAGE.
	Parameters
	----------
	suffix: str
	suffix to be added to output files
	"""
	log = self.log
	log.debug(gt.log_message("primitive", self.myself(), "starting"))
	sfx = params["suffix"]
	slope_diff_cutoff = 1.e-5
	# Lower and upper limits of the number of iterations that will be done
	# by determineSlope.
	min_iter = 1
	max_iter = 1
	adoutputs = []
	for ad in adinputs:
	new_ad = astrodata.create(ad.phu)
	for e, ext in enumerate(ad):
	ints, groups, rows, cols = ext.data.shape
	new_ext_data = np.empty((ints, rows, cols), dtype=ext.data.dtype)
	new_ext_mask = np.empty((ints, rows, cols), dtype=ext.mask.dtype)
	new_ext_variance = np.empty((ints, rows, cols), dtype=ext.variance.dtype)
	new_ad.append(new_ext_data, header=ext.hdr)
	for intn in range(len(ext.data)):
	# Get data from this extension.
	input_var_sect_2 = ext.variance[intn] ** 2
	gdq_sect = ext.mask[intn]
	readnoise_sect = gt.array_from_descriptor_value(ext, 'read_noise')[0][0]
	gain_sect = gt.array_from_descriptor_value(ext, 'gain')[0][0]
	# todo: KL: why is this hardcoded?
	# From the simulator files, this is the frame_time
	# INTTOTT - INTTIME (seconds) could be use. It's the readout time.
	frame_time = 0.38
	group_time = frame_time * (ext.hdr['UTRFRAME'] + ext.hdr['UTRSKIP'])
	nframes = ext.hdr['UTRFRAME']
	max_num_cr = self._get_max_num_cr(gdq_sect, DQ.cosmic_ray)
	data = deepcopy(ext.data[intn])
	# These will be updated in the loop.
	prev_slope_sect = (data[1, :, :] - data[0, :, :]) / group_time
	prev_fit = deepcopy(data)
	iter = 0
	done = False
	while not done:
	(intercept_sect, int_var_sect, slope_sect, slope_var_sect,
	cr_sect, cr_var_sect) = \
	self._compute_slope(data, input_var_sect_2, gdq_sect,
	readnoise_sect, gain_sect,
	prev_fit, prev_slope_sect,
	frame_time, group_time, nframes,
	max_num_cr, DQ.saturated, DQ.cosmic_ray)
	iter += 1
	if iter >= max_iter:
	done = True
	else:
	# If there are pixels with zero or negative variance,
	# ignore them when taking the difference between the
	# slopes computed in the current and previous
	# iterations.
	slope_diff = np.where(slope_var_sect > 0., prev_slope_sect - slope_sect, 0.)
	max_slope_diff = np.abs(slope_diff).max()
	if iter >= min_iter and max_slope_diff < slope_diff_cutoff:
	done = True
	current_fit = self._evaluate_fit(intercept_sect, slope_sect, cr_sect,
	frame_time, group_time, gdq_sect, DQ.cosmic_ray)
	prev_fit = self._positive_fit(current_fit) # use for next iteration
	del current_fit
	prev_slope_sect = slope_sect.copy()
	# Create a basis for the 2D data quality array
	temp_dq = np.zeros((gdq_sect.shape[1], gdq_sect.shape[2]), dtype=DQ.datatype)
	# Replace zero or negative variances with this:
	LARGE_VARIANCE = 1.e8
	v_mask = (slope_var_sect <= 0.)
	if v_mask.any():
	# Replace negative or zero variances with a large value.
	slope_var_sect[v_mask] = LARGE_VARIANCE
	# Also set a flag in the pixel dq array.
	temp_dq[v_mask] = DQ.bad_pixel
	del v_mask
	# If a pixel was flagged (by an earlier step) as saturated in
	# the first group, flag the pixel as bad.
	s_mask = (gdq_sect[0] == DQ.saturated)
	if s_mask.any():
	temp_dq[s_mask] = DQ.bad_pixel
	# Compress the DQ array 3D -> 2D
	pixel_dq = self._dq_compress_sect(gdq_sect, temp_dq)
	new_ext_data[intn] = slope_sect
	new_ext_mask[intn] = pixel_dq
	new_ext_variance[intn] = np.sqrt(slope_var_sect)
	# Both science data & cals now need the NDR group axis removing
	# and the function below can only remove one World axis at a
	# time, so do that one first for cals:
	new_ad[e].reset(new_ext_data, mask=new_ext_mask,
	variance=new_ext_variance)
	new_ad[e].wcs = ext.nddata.window[:, 0].wcs
	astrodata.wcs.remove_unused_world_axis(new_ad[e])
	# Remove the extra integration axis for CAL images.
	# Eventually this should be moved elsewhere in ScorpioDR.
	# Return a 2D array per extension for CAL images.
	# Otherwise return a 3D array for SCI images.
	if "CAL" in ad.tags:
	try:
	assert len(new_ext_data) == 1
	assert len(new_ext_mask) == 1
	assert len(new_ext_variance) == 1
	except AssertionError:
	log.error("CAL integration axis dimension exceeds 1")
	raise ValueError("CAL integration axis dimension exceeds 1")
	new_ad[e].reset(new_ad[e].nddata[0])
	astrodata.wcs.remove_unused_world_axis(new_ad[e])
	# Update the filename.
	ad.update_filename(suffix=sfx, strip=True)
	new_ad.filename = ad.filename
	adoutputs.append(new_ad)
	return adoutputs
	def flagCosmicRaysFromNDRs(self, adinputs=None, **params):
	"""
	Two-Point difference method for finding outliers in a 3-D data array.
	The scheme used in this variation of the method uses numpy array methods
	to compute first differences and find the max outlier in each pixel
	while still working in the full 3-D data array. This makes detection of
	the first outlier very fast. We then iterate pixel-by-pixel over only
	those pixels that are already known to contain an outlier, to look for
	any additional outliers and set the appropriate DQ mask for all outliers
	in the pixel.
	This function was copied and modified from stcal. The code was accessed
	between June and July 2021.
	stcal:
	https://github.com/spacetelescope/stcal
	The two-point difference method used in this function is based on the
	method by Anderson & Gordon, 2011:
	https://iopscience.iop.org/article/10.1086/662593
	TODO:
	- add rej_threshold input
	License
	-------
	Copyright (C) 2020 Association of Universities for Research in Astronomy (AURA)
	Redistribution and use in source and binary forms, with or without
	modification, are permitted provided that the following conditions are met:
	1. Redistributions of source code must retain the above copyright
	notice, this list of conditions and the following disclaimer.
	2. Redistributions in binary form must reproduce the above
	copyright notice, this list of conditions and the following
	disclaimer in the documentation and/or other materials provided
	with the distribution.
	3. The name of AURA and its representatives may not be used to
	endorse or promote products derived from this software without
	specific prior written permission.
	THIS SOFTWARE IS PROVIDED BY AURA ``AS IS'' AND ANY EXPRESS OR IMPLIED
	WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
	MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
	DISCLAIMED. IN NO EVENT SHALL AURA BE LIABLE FOR ANY DIRECT, INDIRECT,
	INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
	BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS
	OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
	ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR
	TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE
	USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH
	DAMAGE.
	Parameters
	----------
	suffix: str
	suffix to be added to output files
	"""
	log = self.log
	log.debug(gt.log_message("primitive", self.myself(), "starting"))
	sfx = params["suffix"]
	# TODO - add input to the primitive to set this value
	# This value yields ~0.04% pixels whose largest non-saturated ratio is above the threshold, using test data from PhoSim
	rej_threshold = 50
	for ad in adinputs:
	for e, ext in enumerate(ad):
	for intn in range(len(ext.data)):
	# Make a copy of the data so we don't make changes to it
	data = deepcopy(ext.data[intn])
	# Get the number of frames averaged per group
	nframes = ext.hdr["UTRFRAME"]
	# Get data characteristics
	ngroups, nrows, ncols = data.shape
	ndiffs = ngroups - 1
	# Get the read noise array and square it, for use later
	read_noise = gt.array_from_descriptor_value(ext, "read_noise")[0][0]
	read_noise_2 = read_noise ** 2
	# Set saturated values and pixels flagged as bad pixels in the
	# input data array to NaN, so they don't get used in any of the
	# subsequent calculations.
	data[np.where(np.bitwise_and(ext.mask[intn], DQ.saturated))] = np.nan
	data[np.where(np.bitwise_and(ext.mask[intn], DQ.bad_pixel))] = np.nan
	# Compute first differences of adjacent groups up the ramp.
	# Note: roll the ngroups axis of data array to the end, to make
	# memory access to the values for a given pixel faster.
	# New form of the array has dimensions [nrows, ncols, ngroups].
	first_diffs = np.diff(np.rollaxis(data, axis=0, start=3), axis=2)
	positive_first_diffs = np.abs(first_diffs)
	# sat_groups is a 3D array that is true when the group is
	# saturated.
	sat_groups = np.isnan(positive_first_diffs)
	# number_sat_groups is a 2D array with the count of saturated
	# groups for each pixel.
	number_sat_groups = sat_groups.sum(axis=2)
	# Make all the first diffs for saturated groups be equal to
	# 100,000 to put them above the good values in the sorted index.
	first_diffs[np.isnan(first_diffs)] = 100000.
	# Here we sort the 3D array along the last axis, which is the
	# group axis. np.argsort returns a 3D array with the last axis
	# containing the indices that would yield the groups/diffs in
	# order.
	sort_index = np.argsort(positive_first_diffs)
	# median_diffs is a 2D array with the clipped median of each pixel
	median_diffs = self._get_clipped_median(ndiffs, number_sat_groups, first_diffs, sort_index)
	# Compute uncertainties as the quadrature sum of the poisson
	# noise in the first difference signal and read noise. Because
	# the first differences can be biased by CRs/jumps, we use the
	# median signal for computing the poisson noise. Here we lower
	# the read noise by the square root of number of frames in
	# the group. Sigma is a 3D array.
	poisson_noise = np.sqrt(np.abs(median_diffs))
	poisson_noise_2 = poisson_noise ** 2
	sigma = np.sqrt(poisson_noise_2 + read_noise_2 / nframes)
	# Reset sigma to exclude pixels with both readnoise and signal=0
	sigma_0_pixels = np.where(sigma == 0.)
	if len(sigma_0_pixels[0] > 0):
	#log.debug(f'Twopt found {len(sigma_0_pixels[0])} pixels with sigma=0')
	#log.debug('which will be reset so that no jump will be detected.')
	huge_num = np.finfo(np.float32).max
	sigma[sigma_0_pixels] = huge_num
	# Compute distance of each sample from the median in units of
	# sigma; note that the use of "abs" means we'll detect positive
	# and negative outliers. ratio is a 3D array with the units of
	# sigma deviation of the difference from the median.
	ratio = np.abs(first_diffs - median_diffs[:, :, np.newaxis]) / sigma[:, :, np.newaxis]
	# Get the group index for each pixel of the largest
	# non-saturated group, assuming the indices are sorted. 2 is
	# subtracted from ngroups because we are using differences
	# and there is one less difference than the number of groups.
	# This is a 2D array.
	max_value_index = ngroups - 2 - number_sat_groups
	# Extract from the sorted group indices the index of the largest
	# non-saturated group.
	row, col = np.where(number_sat_groups >= 0)
	max_index1d = sort_index[row, col, max_value_index[row, col]]
	max_index1 = np.reshape(max_index1d, (nrows, ncols)) # reshape to a 2-D array
	# Get the row and column indices of pixels whose largest
	# non-saturated ratio is above the threshold.
	r, c = np.indices(max_index1.shape)
	row1, col1 = np.where(ratio[r, c, max_index1] > rej_threshold)
	#log.info(f'From highest outlier Two-point found {len(row1)} pixels with at least one CR')
	# Loop over all pixels that we found the first CR in
	number_pixels_with_cr = len(row1)
	for j in range(number_pixels_with_cr):
	# Extract the first diffs for this pixel with at least one
	# CR, yielding a 1-D array.
	pixel_masked_diffs = first_diffs[row1[j], col1[j]]
	# Get the scalar readnoise^2 and number of saturated groups
	# for this pixel
	pixel_rn2 = read_noise_2[row1[j], col1[j]]
	pixel_sat_groups = number_sat_groups[row1[j], col1[j]]
	# Create a CR mask and set 1st CR to be found
	# cr_mask=0 designates a CR
	pixel_cr_mask = np.ones(pixel_masked_diffs.shape, dtype=bool)
	number_CRs_found = 1
	pixel_sorted_index = sort_index[row1[j], col1[j], :]
	pixel_cr_mask[pixel_sorted_index[ndiffs - pixel_sat_groups - 1]] = 0 # setting largest diff to be a CR
	new_CR_found = True
	# Loop and see if there is more than one CR, setting the
	# mask as you go
	while new_CR_found and ((ndiffs - number_CRs_found - pixel_sat_groups) > 1):
	new_CR_found = False
	# For this pixel get a new median difference excluding
	# the number of CRs found and the number of saturated
	# groups
	pixel_med_diff = self._get_clipped_median(ndiffs, number_CRs_found + pixel_sat_groups,
	pixel_masked_diffs, pixel_sorted_index)
	# Recalculate the noise and ratio for this pixel now
	# that we have rejected a CR
	pixel_poisson_noise = np.sqrt(np.abs(pixel_med_diff))
	pixel_poisson_noise_2 = pixel_poisson_noise ** 2
	pixel_sigma = np.sqrt(pixel_poisson_noise_2 + pixel_rn2 / nframes)
	pixel_ratio = np.abs(pixel_masked_diffs - pixel_med_diff) / pixel_sigma
	# Check if largest remaining difference is above threshold
	if pixel_ratio[pixel_sorted_index[ndiffs - number_CRs_found - pixel_sat_groups - 1]] > rej_threshold:
	new_CR_found = True
	pixel_cr_mask[pixel_sorted_index[ndiffs - number_CRs_found - pixel_sat_groups - 1]] = 0
	number_CRs_found += 1
	# Found all CRs for this pixel. Set CR flags in input DQ array for this pixel
	ext.mask[intn, 1:, row1[j], col1[j]] = \
	np.bitwise_or(ext.mask[intn, 1:, row1[j], col1[j]], DQ.cosmic_ray * np.invert(pixel_cr_mask))
	# Update the filename.
	ad.update_filename(suffix=sfx, strip=True)
	return adinputs

If not, I'd appreciate you pointing me to the lines you're looking at (we picked this up unfinished from a contractor last year).

That ramp fitting code is already due for replacement. It was added by a student, based on an old version of the JWST code, and has significant memory issues etc. About 9 months ago, I had a quick look at your then-latest implementation with C optimizations and it seemed too deeply tied into the JWST pipeline for us to be able to use directly, but the recently-published algorithm from Timothy Brandt (also affiliated with STScI) was looking promising (except that the jump detection won't work as-is with ground-based data). I came across some discussion suggesting that you were also considering switching to Brandt's method instead? Since then, our instrument/project scientist has also looked at the implementations used by some other ground-based instruments, eg. at ESO.

To make things even more moot, there has been a recent, high-level decision to change the data format for this upcoming instrument, such that the ramp fitting will get pushed out of our DRAGONS pipeline and into the instrument software. Thanks for your kind offer to help facilitate this, but I think we won't continue using the STCAL code here (whatever we do have will live in another repo, outside the standard user pipeline, and be maintained by other colleagues). I'd be interested to hear briefly about your plans though, eg. have you been testing Brandt's method?

We do collaborate a little bit with your colleagues on other things, so I'll probably see you around!

James.

[emolter]

Hi James, thanks for your detailed response.

I saw two mentions of stcal, at

ScorpioDR/scorpiodr/scorpio/primitives_scorpio_nearIR.py

Line 341 in 33a4476

def flagCosmicRaysFromNDRs(self, adinputs=None, **params):

and

ScorpioDR/scorpiodr/scorpio/primitives_scorpio_nearIR.py

Line 1298 in 33a4476

def _get_clipped_median(self, num_differences, diffs_to_ignore, differences, sorted_index):

I didn't look too closely, but it seemed to me it might be easier to import directly from stcal.

Stcal is supposed to be mostly observatory-agnostic, so if you have suggestions for how we could make it generic enough for you to use easily, we would love to hear them. That's one of the main reasons I reached out.

With regards to ramp fitting, I think you mean the likelihood algorithm. Yes, we did test and implement this (see spacetelescope/stcal#278), and AFAIK it's now our default for jwst and romancal.

Jump detection is also available separately in stcal, but it's too bad that it sounds like it isn't doing what you need?

If you're moving in a different direction then absolutely no worries, of course. We will just continue with business as usual.