Resomat clipping produces warnings

[p-slash]

While running picca_delta_extraction, the following warning pops up:
/global/u1/n/naimgk/picca/py/picca/delta_extraction/utils_pk1d.py:248: RuntimeWarning: divide by zero encountered in true_divide
np.sqrt(4.0 / 2.0 / np.log(
/global/u1/n/naimgk/picca/py/picca/delta_extraction/utils_pk1d.py:242: RuntimeWarning: invalid value encountered in sqrt
(np.sqrt(1.0 / 2.0 / np.log(
which is here:

picca/py/picca/delta_extraction/utils_pk1d.py

Lines 236 to 250 in 54a94a8

	reso = np.clip(reso_matrix, 1.0e-6, 1.0e6)
	#assume reso = A*exp(-(x-central_pixel_pos)**2 / 2 / sigma**2)
	#=> sigma = sqrt((x-central_pixel_pos)/2)**2 / log(A/reso)
	# A = reso(central_pixel_pos)
	# the following averages over estimates for four symmetric values of x
	rms_in_pixel = (
	(np.sqrt(1.0 / 2.0 / np.log(
	reso[reso.shape[0] // 2, :] / reso[reso.shape[0] // 2 - 1, :])) +
	np.sqrt(4.0 / 2.0 / np.log(
	reso[reso.shape[0] // 2, :] / reso[reso.shape[0] // 2 - 2, :])) +
	np.sqrt(1.0 / 2.0 / np.log(
	reso[reso.shape[0] // 2, :] / reso[reso.shape[0] // 2 + 1, :])) +
	np.sqrt(4.0 / 2.0 / np.log(
	reso[reso.shape[0] // 2, :] / reso[reso.shape[0] // 2 + 2, :])))
	/ 4.0) #this is rms

I think this is because two clipped values coincide, so the log becomes 0. I have been using an additive shift instead of clipping in my oversampling method. Here's how I do it:
shift = reso_matrix.min() - np.abs(reso_matrix[reso_matrix.nonzero()]).min()
reso = reso_matrix - shift

This works better then an arbitrary epsilon in my case. reso_matrix may have to be ravelled first.

Would this work here and yield the correct average rms in the end? Thoughts?

[Waelthus]

I'm not sure if this would help, currently we're not using these rough resolution estimates within the analysis anyway.
But agreed that the warnings are annoying.

Note that if we're ever actually impacted by these estimates becoming nan/inf (i.e. if we have the reso matrix diagonal match the first/second off-diagonal in a pixel after clipping resulting in division by zero ) we should be getting a nan/inf mean_reso in the outputs as well (we're using mean there and not nanmean), so if we don't get those (I didn't check for this on the whole dataset) I'd think we have the zero-divisions only for pixels that are removed later on.
Note that the calculation happens in the read-in phase before spectra are even masked/rebinned/truncated to the forest region.

But agreed that if matrix values close to zero are causing this, your solution should work better than the simple clipping.

[p-slash]

There are also surprising amounts of negative elements in resolution matrix, so this is probably quite common.
Yeah, the final results are numerical, so I don't think this is a crucial bug. But some pixels might be spuriously masked by np.isnan checks or something down the line.

[Waelthus]

I did search the codebase for isnan, isfinite, nanmean etc, and none of those appears in combination with resolution. So I think we are currently fine. Nevertheless could be good to fix this...

[p-slash]

Branch resomat-nan-fix addresses this issue. It changes MEAN_RESO values slightly. Need to test on data to check against crashing behaviors.