Fitting compound Polynominal2D models no longer works

Description

I have code where I fit a compound model of a Gaussian2D and a Polynomial2D that stopped working as of astropy 7.x. It can be made to work if I switch to SimplexLSQFitter, but the results aren't as reliable. DogBoxLSQFitter, LevMarLSQFitter, SLSQPLSQFitter, and LMLSQFitter all fail in a similar manner. The failures happen within fit_deriv, for example:

-> [1287](https://file+.vscode-resource.vscode-cdn.net/Users/tim/MMT/mmtwfs/notebooks/~/conda/envs/mmtwfs/lib/python3.12/site-packages/astropy/modeling/fitting.py:1287)     for _ in weights * np.array(model.fit_deriv(x, y, *params))
   [1288](https://file+.vscode-resource.vscode-cdn.net/Users/tim/MMT/mmtwfs/notebooks/~/conda/envs/mmtwfs/lib/python3.12/site-packages/astropy/modeling/fitting.py:1288) ]

ValueError: setting an array element with a sequence. The requested array has an inhomogeneous shape after 1 dimensions. The detected shape was (12,)

Caused by

Implement Analytical Derivatives on CompoundModel #17034 (v7.0.0)

Expected behavior

I would expect the compound model fits to work correctly like they did in previous versions of astropy.

How to Reproduce

I built an example using this example in the astropy.modeling documentation as a starting point.

import warnings
import numpy as np
import matplotlib.pyplot as plt
from astropy.modeling import models, fitting
from astropy.utils.exceptions import AstropyUserWarning

# Generate fake data
rng = np.random.default_rng(0)
y, x = np.mgrid[:128, :128]
z = 2. * x ** 2 - 0.5 * x ** 2 + 1.5 * x * y - 1.
z += rng.normal(0., 0.1, z.shape) * 50000.
z += models.Gaussian2D(
    amplitude=50000,
    x_mean=60, 
    y_mean=60, 
    x_stddev=5,
    y_stddev=5
)(x, y)

# Fit the data using astropy.modeling
p_init = models.Polynomial2D(degree=2) + models.Gaussian2D(amplitude=50000, x_mean=60, y_mean=60)
fit_p = fitting.DogBoxLSQFitter()

with warnings.catch_warnings():
    # Ignore model linearity warning from the fitter
    warnings.filterwarnings('ignore', message='Model is linear in parameters',
                            category=AstropyUserWarning)
    p = fit_p(p_init, x, y, z)

# Plot the data with the best-fit model
plt.figure(figsize=(8, 2.5))
plt.subplot(1, 3, 1)
plt.imshow(z, origin='lower', interpolation='nearest', vmin=-1e4, vmax=5e4)
plt.title("Data")
plt.subplot(1, 3, 2)
plt.imshow(p(x, y), origin='lower', interpolation='nearest', vmin=-1e4,
           vmax=5e4)
plt.title("Model")
plt.subplot(1, 3, 3)
plt.imshow(z - p(x, y), origin='lower', interpolation='nearest', vmin=-1e4,
           vmax=5e4)
plt.title("Residual")

Versions

import astropy
try:
    astropy.system_info()
except AttributeError:
    import platform; print(platform.platform())
    import sys; print("Python", sys.version)
    import astropy; print("astropy", astropy.__version__)
    import numpy; print("Numpy", numpy.__version__)
    import erfa; print("pyerfa", erfa.__version__)
    try:
        import scipy
        print("Scipy", scipy.__version__)
    except ImportError:
        print("Scipy not installed")
    try:
        import matplotlib
        print("Matplotlib", matplotlib.__version__)
    except ImportError:
        print("Matplotlib not installed")

# Paste the result here

platform
--------
platform.platform() = 'macOS-15.1.1-arm64-arm-64bit'
platform.version() = 'Darwin Kernel Version 24.1.0: Thu Oct 10 21:03:15 PDT 2024; root:xnu-11215.41.3~2/RELEASE_ARM64_T6000'
platform.python_version() = '3.12.7'

packages
--------
astropy              7.0.0
numpy                2.1.3
scipy                1.14.1
matplotlib           3.9.3
pandas               2.2.3
pyerfa               2.0.1.5

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 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 [email protected].

Interestingly, I get a different error with the env I happened to have setup (oldest-deps oriented)

/Users/clm/dev/astropy-project/coordinated/astropy/astropy/modeling/fitting.py:1273: VisibleDeprecationWarning: Creating an ndarray from ragged nested sequences (which is a list-or-tuple of lists-or-tuples-or ndarrays with different lengths or shapes) is deprecated. If you meant to do this, you must specify 'dtype=object' when creating the ndarray.
  for _ in weights * np.array(model.fit_deriv(x, y, *params))
Traceback (most recent call last):
  File "/Users/clm/dev/astropy-project/coordinated/astropy/issues/17509/t.py", line 28, in <module>
    p = fit_p(p_init, x, y, z)
        ^^^^^^^^^^^^^^^^^^^^^^
  File "/Users/clm/dev/astropy-project/coordinated/astropy/astropy/modeling/fitting.py", line 287, in wrapper
    return func(self, model, x, y, z=z, **kwargs)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/Users/clm/dev/astropy-project/coordinated/astropy/astropy/modeling/fitting.py", line 1420, in __call__
    init_values, fitparams, cov_x = self._run_fitter(
                                    ^^^^^^^^^^^^^^^^^
  File "/Users/clm/dev/astropy-project/coordinated/astropy/astropy/modeling/fitting.py", line 1594, in _run_fitter
    self.fit_info = optimize.least_squares(
                    ^^^^^^^^^^^^^^^^^^^^^^^
  File "/Users/clm/dev/astropy-project/coordinated/astropy/.venv/lib/python3.11/site-packages/scipy/optimize/_lsq/least_squares.py", line 849, in least_squares
    J0 = jac(x0, *args, **kwargs)
         ^^^^^^^^^^^^^^^^^^^^^^^^
  File "/Users/clm/dev/astropy-project/coordinated/astropy/astropy/modeling/fitting.py", line 1573, in _dfunc
    self._wrap_deriv(
  File "/Users/clm/dev/astropy-project/coordinated/astropy/astropy/modeling/fitting.py", line 1273, in _wrap_deriv
    for _ in weights * np.array(model.fit_deriv(x, y, *params))
             ~~~~~~~~^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
TypeError: can't multiply sequence by non-int of type 'float'

platform
--------
platform.platform() = 'macOS-15.1.1-arm64-arm-64bit'
platform.version() = 'Darwin Kernel Version 24.1.0: Thu Oct 10 21:03:11 PDT 2024; root:xnu-11215.41.3~2/RELEASE_ARM64_T6020'
platform.python_version() = '3.11.10'

packages
--------
astropy              7.1.dev253+g4a55b10002
numpy                1.23.2
scipy                1.9.2
matplotlib           3.6.0
pandas               2.0.0
pyerfa               2.0.1.1

Upgrading numpy (and dependees for compatibility) to 2.1.3, I can reproduce the error from the report. I'll try to fix both !

I can also confirm that this is a regression in astropy 7.0.0; no error is raised with 6.1.7
I suspect that #17034 might be the source of the regression, but I'll run an actual bisection to check that.

@Cadair

confirmed that #17034 is the culprit.
Technically, the first breaking commit is 898f340, though it leads to a different error: NotImplementedError: SAUSUAGE which I assume is a left over from @Cadair's development process (and not squashing is on me).

I could have sworn I encountered and fixed this issue during the development.

clearly you have a head start. Let me know if you'd like a hand but otherwise I'm stepping back from this issue :)

Ok, I don't understand wtaf the analytical derivatives inside Polynominal are doing well enough to try and solve this at 5pm on a Friday.

This method on Polynominal2D:

astropy/astropy/modeling/polynomial.py

Line 1223 in 7aa40e8

def fit_deriv(self, x, y, *params):

returns an array of shape (144, 6) where 6 is the number of parameters, while the _fit_deriv method of Gaussian2D returns a 6 element list of (12, 12) arrays.

I have no idea why Polynominal2D is doing what it's doing so I would only be blindly fixing it.

@nden

@nden / @perrygreenfield / @WilliamJamieson You got any ideas?

As a start we should add a new test.

If no reply from modeling maintainers, do we have to revert #17034 ?

I am pretty sure this is limited to Polynominal2D, I don't think we should revert it we should fix the bug. I will try and take another look this week.

I am going to try and have a quick look at this

@tepickering

@tepickering - could you try out #17618?

Navigation Menu

Repository navigation

Fitting compound Polynominal2D models no longer works #17509

Description

Caused by

Expected behavior

How to Reproduce

Versions

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Fitting compound Polynominal2D models no longer works #17509

Description

Description

Caused by

Expected behavior

How to Reproduce

Versions

Activity

github-actions commented on Dec 6, 2024

neutrinoceros commented on Dec 6, 2024

neutrinoceros commented on Dec 6, 2024

neutrinoceros commented on Dec 6, 2024

neutrinoceros commented on Dec 6, 2024

Cadair commented on Dec 6, 2024

neutrinoceros commented on Dec 6, 2024

Cadair commented on Dec 6, 2024

Cadair commented on Dec 6, 2024

pllim commented on Dec 9, 2024

Cadair commented on Dec 10, 2024

marscode-agent commented on Dec 11, 2024

astrofrog commented on Jan 9, 2025

astrofrog commented on Jan 10, 2025

Metadata

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Issue actions