Fitting different datasets simultaneously, using models with parameters tied across models

[hamogu]

Given two models called a with parameters a1, a2 and b with parameters b1, b2, I want to do a joint fit, of a evaluated on data1 (in my example that data is a two-dimensional image) and b evaluated on data2 (in my example a different image which happens to have the same dimensionality, but in the general case, that could be different), such that parameters a1 and b1 are tied through some arbitrary python function (in my example: the identity function, a1 = b1).

I'll describe my specific use case in mode detail. In this case, I deal with PSF fitting, but note that this is outside of the scope of current photutils and the modeling problem is more general. I can give examples for fitting spectra or SEDs, too, I just happened to come upon this problem by working with images right now, that's why I describe this example here: I have two images. They show the same object on the same WCS, but that object has different fluxes (duh, it's in different bands). Unfortunately, there are a lot of unususable pixles in both images (saturated from a near-by object). Thus, I'd like to fit a psf (let's call the models psf1 and psf2) to both images at the same time and require that x, y are the same, but the flux is different. So, I have two models, evaluated on two different datasets (psf_621_copy on image1 with parameters x, y, amplitude_1 and psf_845_copy evaluated on image2 with parameters x, y, flux2). I want to minimize the combined fit statistic. Is there a reasonable way to make that happen with astropy now or do I have to wait for #8769?
Maybe something like psf1 | psf2 evaluated on data that is a list with two elements (image 1 and amplitude2). But then, how do I couple x and y?
Seems that a problem like this is not unique to me, but I don't find a way to express it well in the current astropy modelling paradigm.

I'm not asking for help with a specific problem, I'm describing the use case here for a feature request. If that's already possible in the current version, an example in the docs might be enough.

For reference, I'll paste below how I address this problem with a custom model, but I feel that there should be (and maybe is already?) a way to make that work by using operations on models as opposed to coding up a user model:

@custom_model
def twobandpsf(xy, x_0=0., y_0=0., amplitude1=1, amplitude2=1.):
    psf_621_copy.x_0_2 = x_0
    psf_621_copy.y_0_2 = y_0
    psf_621_copy.amplitude_3 = amplitude1
    psf_845_copy.x_0_2 = x_0
    psf_845_copy.y_0_2 = y_0
    psf_845_copy.amplitude_3 = amplitude2
    return np.stack([psf_621_copy(xy[0], xy[1]), psf_845_copy(xy[0], xy[1])])
@custom_model

[hamogu]

Here is an example with 1D datasets: I'm fitting lines in a time series of spectra. I know that the position of the lines does not change (because they are formed several au form the star in the outer disk where changes of the disk structure take long), however, the flux might change much faster as the irradiation from the star changes. So, I want to fit several spectra at the same time, for simplicity, let's assume that I use a Gaussian as the model. The positions of all Gaussians is the same, but the values is unknown a priory (it's a free parameter of the fit), the fluxes are all different.

[hamogu]

Last, a model where data1 and data2 are different things: I have a physical MHD model that predicts the spectrum and the position of a shock in an outflow from a young star. So, I want to fit, at the same time, the observed spectrum (1 D dataset) and the position (just a number, you can see that as a 0D dataset, or a 1D dataset, with just one y value). Again, I'm using parameters (mass outflow rate, velocity) that are the same in model_for_spectrum (which takes a spectrum as data and velocity and mass outflow rate as parameters) and model_for_position (which takes a position as data and velocity ans mass outflow rate as parameters). In https://ui.adsabs.harvard.edu/abs/2014ApJ...795...51G/abstract I used sherpa for fitting, but there is no reason why astropy should not be able to do that.

[hamogu]

I just discovered https://docs.astropy.org/en/latest/api/astropy.modeling.fitting.JointFitter.html#astropy.modeling.fitting.JointFitter
which solves this problem, if the two parameters across models are meant to match exactly. So, improve docs by adding example. However, this does not solve the (most extended) problem that parameters might just be tied by some function between models.

[nden]

I think something along the lines of JointFitter should work (not currently). I am hoping work on modeling will start again next month.

[karllark]

The JointFitter does work and likely can do some of the fitting described. Just submitted a PR to provide docs and an email. See #12720.

[pllim]

tl;dr -- Does #12720 completely resolve this issue?

[karllark]

Not fully. Tying parameters between models using a python function is not supported by JointFitter at this time.