Fitting/Modeling/Sigma Clipping Crashes with NANs, whereas before it just gave a warning

[orifox]

I'm working on PSF Fitting using Photutils. My script previously worked (astropy == 5.0.4 ) worked. I upgraded astropy, and now it fails. Although I'm not exactly sure why. Error message below. Nothing has changed.

Now (astropy == 5.2.dev75+ga3c27114b):

---------------------------------------------------------------------------
NonFiniteValueError                       Traceback (most recent call last)
Input In [26], in <cell line: 22>()
     15 tic = time.perf_counter()
     17 phot = IterativelySubtractedPSFPhotometry(finder=daofind, group_maker=daogroup,
     18                                             bkg_estimator=mmm_bkg, psf_model=psf_model,
     19                                             fitter=fitter,
     20                                             niters=10, fitshape=(7, 7), aperture_radius=ap_radius, 
     21                                             extra_output_cols=('sharpness', 'roundness2'))
---> 22 psf_phot_results = phot(data)
     23 toc = time.perf_counter()
     25 print('Time needed to perform photometry:', '%.2f' % ((toc - tic) / 3600), 'hours')

File ~/miniconda3/envs/hackphot/lib/python3.8/site-packages/photutils/psf/photometry.py:228, in BasicPSFPhotometry.__call__(self, image, init_guesses)
    223 def __call__(self, image, init_guesses=None):
    224     """
    225     Perform PSF photometry. See `do_photometry` for more details
    226     including the `__call__` signature.
    227     """
--> 228     return self.do_photometry(image, init_guesses)

File ~/miniconda3/envs/hackphot/lib/python3.8/site-packages/photutils/psf/photometry.py:758, in IterativelySubtractedPSFPhotometry.do_photometry(self, image, init_guesses)
    755     if self.aperture_radius is None:
    756         self.set_aperture_radius()
--> 758     output_table = self._do_photometry()
    760 output_table.meta = {'version': _get_version_info()}
    762 return QTable(output_table)

File ~/miniconda3/envs/hackphot/lib/python3.8/site-packages/photutils/psf/photometry.py:817, in IterativelySubtractedPSFPhotometry._do_photometry(self, n_start)
    810         init_guess_tab.add_column(
    811             Column(name=param_tab_name,
    812                    data=(getattr(self.psf_model,
    813                                  param_name) *
    814                          np.ones(len(sources)))))
    816 star_groups = self.group_maker(init_guess_tab)
--> 817 table, self._residual_image = super().nstar(
    818     self._residual_image, star_groups)
    820 star_groups = star_groups.group_by('group_id')
    821 table = hstack([star_groups, table])

File ~/miniconda3/envs/hackphot/lib/python3.8/site-packages/photutils/psf/photometry.py:425, in BasicPSFPhotometry.nstar(self, image, star_groups)
    419 for row in star_groups.groups[n]:
    420     usepixel[overlap_slices(large_array_shape=image.shape,
    421                             small_array_shape=self.fitshape,
    422                             position=(row['y_0'], row['x_0']),
    423                             mode='trim')[0]] = True
--> 425 fit_model = self.fitter(group_psf, x[usepixel], y[usepixel],
    426                         image[usepixel])
    427 param_table = self._model_params2table(fit_model,
    428                                        star_groups.groups[n])
    429 result_tab = vstack([result_tab, param_table])

File ~/miniconda3/envs/hackphot/lib/python3.8/site-packages/astropy/modeling/fitting.py:269, in fitter_unit_support.<locals>.wrapper(self, model, x, y, z, **kwargs)
    264         raise NotImplementedError("This model does not support being "
    265                                   "fit to data with units.")
    267 else:
--> 269     return func(self, model, x, y, z=z, **kwargs)

File ~/miniconda3/envs/hackphot/lib/python3.8/site-packages/astropy/modeling/fitting.py:1222, in _NonLinearLSQFitter.__call__(self, model, x, y, z, weights, maxiter, acc, epsilon, estimate_jacobian)
   1219 model_copy.sync_constraints = False
   1220 farg = (model_copy, weights, ) + _convert_input(x, y, z)
-> 1222 init_values, fitparams, cov_x = self._run_fitter(model_copy, farg,
   1223                                                  maxiter, acc, epsilon, estimate_jacobian)
   1225 self._compute_param_cov(model_copy, y, init_values, cov_x, fitparams, farg)
   1227 model.sync_constraints = True

File ~/miniconda3/envs/hackphot/lib/python3.8/site-packages/astropy/modeling/fitting.py:1283, in LevMarLSQFitter._run_fitter(self, model, farg, maxiter, acc, epsilon, estimate_jacobian)
   1281     dfunc = self._wrap_deriv
   1282 init_values, _, _ = model_to_fit_params(model)
-> 1283 fitparams, cov_x, dinfo, mess, ierr = optimize.leastsq(
   1284     self.objective_function, init_values, args=farg, Dfun=dfunc,
   1285     col_deriv=model.col_fit_deriv, maxfev=maxiter, epsfcn=epsilon,
   1286     xtol=acc, full_output=True)
   1287 fitter_to_model_params(model, fitparams)
   1288 self.fit_info.update(dinfo)

File ~/miniconda3/envs/hackphot/lib/python3.8/site-packages/scipy/optimize/_minpack_py.py:410, in leastsq(func, x0, args, Dfun, full_output, col_deriv, ftol, xtol, gtol, maxfev, epsfcn, factor, diag)
    408 if not isinstance(args, tuple):
    409     args = (args,)
--> 410 shape, dtype = _check_func('leastsq', 'func', func, x0, args, n)
    411 m = shape[0]
    413 if n > m:

File ~/miniconda3/envs/hackphot/lib/python3.8/site-packages/scipy/optimize/_minpack_py.py:24, in _check_func(checker, argname, thefunc, x0, args, numinputs, output_shape)
     22 def _check_func(checker, argname, thefunc, x0, args, numinputs,
     23                 output_shape=None):
---> 24     res = atleast_1d(thefunc(*((x0[:numinputs],) + args)))
     25     if (output_shape is not None) and (shape(res) != output_shape):
     26         if (output_shape[0] != 1):

File ~/miniconda3/envs/hackphot/lib/python3.8/site-packages/astropy/modeling/fitting.py:1077, in _NonLinearLSQFitter.objective_function(self, fps, *args)
   1074     value = np.ravel(weights * (model(*args[2: -1]) - meas))
   1076 if not np.all(np.isfinite(value)):
-> 1077     raise NonFiniteValueError("Objective function has encountered a non-finite value, "
   1078                               "this will cause the fit to fail!")
   1080 return value

NonFiniteValueError: Objective function has encountered a non-finite value, this will cause the fit to fail!`

[github-actions[bot]]

Welcome to Astropy 👋 and thank you for your first issue!

A project member will respond to you as soon as possible; in the meantime, please double-check the guidelines for submitting issues and make sure you've provided the requested details.

GitHub issues in the Astropy repository are used to track bug reports and feature requests; If your issue poses a question about how to use Astropy, please instead raise your question in the Astropy Discourse user forum and close this issue.

If you feel that this issue has not been responded to in a timely manner, please leave a comment mentioning our software support engineer @embray, or send a message directly to the development mailing list. If the issue is urgent or sensitive in nature (e.g., a security vulnerability) please send an e-mail directly to the private e-mail feedback@astropy.org.

[pllim]

@orifox , what version of photutils?

@larrybradley and/or @eteq , I transferred this issue from astropy to photutils in case this is something that downstream (photutils) need to fix. Also cc @WilliamJamieson since he has been refactoring modeling.

[larrybradley]

@pllim Nothing has changed in IterativelySubtractedPSFPhotometry recently. I suspect the change was in in astropy modeling.

[larrybradley]

This looks like the change in astropy that is causing the error: https://github.com/astropy/astropy/pull/12811

[pllim]

Does this mean photutils need to adapt to new fitting algorithm or is this astropy bug?

[WilliamJamieson]

This was my fear with about going with an error rather than a warning in astropy/astropy#12811, technically as soon as non-finite values are encountered during fitting the results of the fit become unreliable.

The short answer this is when NaN values are present, scipy.optimize.leastsq fails to actually do any optimization leading to fitting "stalling" at some non-optimized value. For example if a NaN is present in the initial step no actual fitting takes place, such as what was reported in astropy/astropy#12809.

[WilliamJamieson]

My suggestion is that the model being fit have a filter which removes NaN (or other non-finite values) from the output of the model.

[WilliamJamieson]

This maybe a good place to revisit the issue of whether the error being raised should instead be a warning.

[larrybradley]

If fitting with NaN inputs effectively does nothing, why not have the modeling tools automatically remove/mask the NaN values (with a warning that such values were excluded)? This is what we currently do in sigma_clip / SigmaClip.

[WilliamJamieson]

If fitting with NaN inputs effectively does nothing, why not have the modeling tools automatically remove/mask the NaN values (with a warning that such values were excluded)? This is what we currently do in sigma_clip / SigmaClip.

I proposed that as an optional feature (better to be optional in case someone does not want it?), and it was rejected with the comment that persons could just wrap their models in such a feature.

[pllim]

Does this block v5.1 release? If so, please update photutils entry at https://github.com/astropy/astropy/wiki/v5.1-RC-testing

[hamogu]

@orifox when you said "your script previously worked" do you just mean "It run without and error" or "it gave sensible fits"? I was one of the people that argued that it's better to fail explicitly, than to give silently wrong results. My understanding is that with an older version of astropy, your script passed, but in fact it would not have executed the fit, but just returned your input values or some random value the fitter happened to try when it first encountered a nan or inf.

I imagine that could be bad, let's say you end up with a catalog of fitted sources and believe their magnitudes, when in fact, the fit never fully run.

Do you know if the numbers you got with the previous astropy are good (maybe I'm not understanding correctly how the fitting code works)? If not, would be agree that an error is the best way to alert a user (as it's not in the new astropy version) or what would you prefer to happen?

[orifox]

Hi @hamogu Good question. By "previously worked", I meant it "ran without error." I am not particularly in the business (yet) of science validation of photutils. But you bring up a great point. If I was getting bad photometry, I didn't know it. Of course, I'd prefer to know, but perhaps that output can be as a NAN or upper limit or huge error bar rather than crashing the code. My data sets are so large, there is no way I can fix every pixel/star so that the models work every time. In other words, I need the code to complete the fitting in some capacity so that I can view the results rather than just crash. I'm happy to discuss other feasible solutions.

[larrybradley]

My suggestion is that the model being fit have a filter which removes NaN (or other non-finite values) from the output of the model.

@WilliamJamieson How would that work exactly? For example, if I'm trying to fit my data (with NaNs) with the astropy Gaussian2D model, then I would need to rewrite the Gaussian2D model to add a filter? Wouldn't every model in astropy need to be rewritten to filter non-finite values? If that's the case, why not do it for every model (even user-defined custom models) in a single place at the fitting stage?

[larrybradley]

IMHO, I think the sensible thing to do here is have the fitting code automatically remove the non-finite values (with a warning). I see that @astrofrog proposed exactly that 7 years ago (https://github.com/astropy/astropy/issues/3575), but apparently it never went anywhere.

[hamogu]

As I said 7 years ago, I still think that removing nans silently in general is not a good idea since usually a nan or inf means that something went wrong. However, in the specific case of photutils it's probably safe to do that. In images, we know that nans will comes from cosmics, masked overexposed regions etc. and the rest of the image is fine. So, I suggest that photutils, not astropy, does something like (pseudo-code):

index = ~np.isnan(image)
m_best = fit(m_init, xy[index], image[index])

or, more precisely, that we insert one new line before https://github.com/astropy/photutils/blob/main/photutils/psf/photometry.py#L425 :

            usepixel = usepixel & np.isfinite(image)
            fit_model = self.fitter(group_psf, x[usepixel], y[usepixel],
                                    image[usepixel])

(This code was written 6 years ago by @mirca as part of a GSOC project.)

[WilliamJamieson]

@WilliamJamieson How would that work exactly? For example, if I'm trying to fit my data (with NaNs) with the astropy Gaussian2D model, then I would need to rewrite the Gaussian2D model to add a filter? Wouldn't every model in astropy need to be rewritten to filter non-finite values? If that's the case, why not do it for every model (even user-defined custom models) in a single place at the fitting stage?

I think the simplest place to put this is in the fitter itself, as I originally proposed in astropy/astropy#12811, that is one can add an option to the fitter __call__, say filter_non_finite=True, which turns on the fitter. I do not think it should be enabled by default; however, because filtering out some data points changes the fitting problem itself meaning users should be aware of this.

[pllim]

astropy/astropy#13257 is merged. Is there a PR here that should go in first before RC2 or actual v5.1 release?

[larrybradley]

I think @WilliamJamieson is planning to submit a PR on adding the filter_non_finite option, but that can wait until after the 5.1 release.

[pllim]

I see that #1350 is merged. Does that mean we no longer need https://github.com/astropy/astropy/pull/13259 ?

[larrybradley]

No, I think https://github.com/astropy/astropy/pull/13259 is a useful addition!