[Error] during model making correction model "RuntimeWarning division by zero error"

[sumanchapai]

I have been getting this error when making array correction model:

[REDACTED]\Python310\site-packages\regularizepsf\corrector.py:239: RuntimeWarning: invalid value encountered in divide
  [v / v.sum() for v in self._evaluations.values()], dtype=float)

but what's been surprising is that I am only getting this when making model some images and not others.

I am using this code:

from typing import List
from pathlib import Path
import regularizepsf as rpsf
import numpy as np

COMA_PSF_SIZE = 16 # size of the PSF model to use in pixels
COMA_PATCH_SIZE = 128  # square side dimension PSF will be applied over
COMA_ALPHA = 3  # see paper
COMA_EPSILON = 0.3  # see paper

fwhm_x_target, fwhm_y_target = 2.855, 2.8514
IMAGES_FOLDER = Path("G:/Raw Data/Summer 2003/September 8, 2003/m23")

def make_coma_correction_model(
    images: List[str] | List[Path], xfwhm_target: float, yfwhm_target: float
) -> rpsf.ArrayCorrector:
    # Define the target PSF
    @rpsf.simple_psf
    def target(
        x,
        y,
        x0=COMA_PATCH_SIZE / 2,
        y0=COMA_PATCH_SIZE / 2,
        sigma_x=xfwhm_target / 2.355,
        sigma_y=yfwhm_target / 2.355,
    ):
        return np.exp(
            -(
                np.square(x - x0) / (2 * np.square(sigma_x))
                + np.square(y - y0) / (2 * np.square(sigma_y))
            )
        )

    target_evaluation = target(
        *np.meshgrid(np.arange(COMA_PATCH_SIZE), np.arange(COMA_PATCH_SIZE))
    )

    # Extract all the stars from that image and create a PSF model with a target PSF
    image_paths = [str(p) for p in images]  # No support for pathlib.Path yet
    cpc = rpsf.CoordinatePatchCollection.find_stars_and_average(
        image_paths, COMA_PSF_SIZE, COMA_PATCH_SIZE
    )

    return cpc.to_array_corrector(target_evaluation)

def ac1():
    ac_1_images = [ IMAGES_FOLDER / f"m23_3.5-{x}.fit" for x in range (101, 121)]
    return make_coma_correction_model(
        ac_1_images, fwhm_x_target, fwhm_y_target
    )

def ac2():
    ac_2_images = [ IMAGES_FOLDER / f"m23_3.5-{x}.fit" for x in range (521, 541)]
    return make_coma_correction_model(ac_2_images, fwhm_x_target, fwhm_y_target)

if __name__ == "__main__":
    # Note that when generating array corrector model using ac1, we don't get any warning
    # While when using ac2, we do. See this by just ac1 and commenting ac2 and vice versa.
    ac1()
    ac2()

Note that I am getting this error when running the function ac1 but not ac2. The only difference in these two are the images used. The images used are available at: https://drive.google.com/drive/folders/1XuYReGwjFqnEm1GV0O_uTO31aczy8BHT?usp=sharing

The impact of this is that when applying correction using model from ac2, my images have a value of nan in certain boxes. These are the black rectangles in these images.

Are these regions in the original images where it wasn't possible to generate a psf model?

[sumanchapai]

Additionally, it also looks like some pixel values become < 0 after applying the correction to an image. I am also confused as to how to make sense of that.

[jmbhughes]

I believe you're right. That warning arises probably because you divide by zero, indicating you likely didn't have enough stars in that region to build a model. I can look into it a bit more though. We might want to warn the user to expect this if their model is missing patches instead of just springing NaN regions on them.

As to the negative values, it may be because the model is built using a background subtracted stellar PSF. Another good thing for me to look into.

Give me a few days and I'll get back to you. Thanks!

[sumanchapai]

In the case of NaN values, I think that we'll internally just replace those values with zeros which mean that there's no data in this region. This is equivalent of saying that in this region, we cannot do coma correction reasonably well thus we will not extract star ADU values from that region at all. I don't know if we can do the same when the ADU is negative. By can I mean whether it's scientifically appropriate. I would be grateful for your suggestion.

[sumanchapai]

Hi @jmbhughes , I think that we are in the final phase of integrating coma correction into our software. We have a couple of questions, and it wanted to clarify with you:

1. I see that most bright stars have a black square surrounding them after coma correction. Is this because of background subtraction? I see that not all stars that I can identify with the eye have this feature. Is this because those stars weren't identified as stars by the correction model? What things are identified as stars and what not? Does this mean that the stars where I don't see a black square around aren't coma corrected?

[jmbhughes]

The correction is applied to every pixel, not just on a star by star basis.

The model is only built from stars (although conceivably you could do it a different way). There is a finicky way to determine which stars are being used in each image; I'll see if I can find time to make that easier. We also have this visualization tool that shows how many stars are in each region when building the model.

As to the black squares, those stars are being coma corrected but an artifact is being introduced in the procedure. My guess, without thinking deeply about it, is that your target PSF is too small maybe.

Action items for Marcus (although you're welcome to take a stab at them too):

• [x] improve docs to clarify more simply how the technique works
• [x] investigate background subtraction and impacts on negative values in corrected images
• [x] open an issue to make a warning when building a model that a patch is empty... also would be good to remind the user this is happening when applying the PSF model
• [x] open an issue to improve star detection feedback, i.e. a way to know which stars are used when building a model