Allow fitting of a wider range of compound models - #10415

@karllark The fitted scale parameter in BlackBody has different values on Linux and MacOSX. The tests are failing because they are running on Linux while I entered fitted values from a run on MacOSX. Weird!
When I run the example on Linux I get the same value for scale as the tests.

I've run into this problem before. And never found the solution though I did do a fairly deep dive. What are the differences you are seeing? I remember they were not large, but larger than expected for numerical precision.

I've solved this problem for tests by changing the tolerance allowed for the test comparison.

I've tried my problem and am now getting a different error. Likely this error is due to my not setting up my attenuation model correctly.

AttributeError: 'S07_attenuation' object has no attribute '_parameter_units_for_data_units

And here is the S07_attenuation model.

class S07_attenuation(Fittable1DModel):
    """
    Smith, Draine, et al. (2007) kvt attenuation model calculation.
    Calculation is for a fully mixed geometrically model.
    Uses an extinction curve based on the silicate profiles from
    Kemper, Vreind, & Tielens (2004, apJ, 609, 826).
    Constructed as a weighted sum of two components: silicate profiles,
    and an exponent 1.7 power-law.

    Attenuation curve for a mixed case calculated from
    .. math::

        Att(x) = \\frac{1 - e^{-\\tau_{x}}}{\\tau_{x}}

    Parameters
    ----------
    kvt_amp : float
      amplitude
    """

    # Attenuation tau
    tau_sil = Parameter(description="kvt term: amplitude",
                        default=1.0, min=0.0, max=10)

    @staticmethod
    def kvt(in_x):
        """
        Output the kvt extinction curve
        """
        kvt_wav = np.array([8.0, 8.2, 8.4, 8.6, 8.8, 9.0, 9.2, 9.4, 9.6, 9.7,
                            9.75, 9.8, 10.0, 10.2, 10.4, 10.6, 10.8, 11.0,
                            11.2, 11.4, 11.6, 11.8, 12.0, 12.2, 12.4, 12.6,
                            12.7])

        kvt_int = np.array([.06, .09, .16, .275, .415, .575, .755, .895, .98,
                            .99, 1.0, .99, .94, .83, .745, .655, .58, .525,
                            .43, .35, .27, .20, .13, .09, .06, .045, .04314])

        # Extend kvt profile to shorter wavelengths
        kvt_wav_short = in_x[in_x < min(kvt_wav)]
        kvt_int_short_tmp = (min(kvt_int)
                             * np.exp(2.03*(kvt_wav_short-min(kvt_wav))))
        # Since kvt_int_shoft_tmp does not reach min(kvt_int),
        # we scale it to stitch it.
        kvt_int_short = kvt_int_short_tmp * (kvt_int[0]/max(kvt_int_short_tmp))

        # Extend kvt profile to longer wavelengths
        kvt_wav_long = in_x[in_x > max(kvt_wav)]
        # Need to input the value instead of using zeros
        kvt_int_long = np.zeros(len(kvt_wav_long))

        spline_x = np.concatenate([kvt_wav_short, kvt_wav, kvt_wav_long])
        spline_y = np.concatenate([kvt_int_short, kvt_int, kvt_int_long])

        spline_rep = interpolate.splrep(spline_x, spline_y)
        new_spline_y = interpolate.splev(in_x, spline_rep, der=0)

        nf = Drude1D(amplitude=0.4,
                     x_0=18.,
                     fwhm=0.247*18.)
        ext = nf(in_x) + new_spline_y

        # Extend to ~2 um
        # assuing beta is 0.1
        beta = 0.1
        y = (1.0-beta)*ext + beta*(9.7/in_x)**1.7

        return y

    def evaluate(self, in_x, tau_si):
        if tau_si == 0.0:
            return np.full((len(in_x)), 0.0)
        else:
            tau_x = tau_si*self.kvt(in_x)
            return (1.0 - np.exp(-1.0*tau_x))/tau_x

`

The values of the scale parameter on Linux vs MacOSX are 5.821110 vs 3.795289, so too large.

In order to fit a model with units you need to write a method called _parameter_units_for_data_units which takes two dict - inputs_unit1 and outputs_unit` and returns a dict of {output: unit}.

In order to fit a model with units you need to write a method called _parameter_units_for_data_units which takes two dict - inputs_unit1 and outputs_unit` and returns a dict of {output: unit}.

Can I be cheeky and ask if there are docs for this or a good example you could point me to?

And using the SimplexLSQFitter on Linux and Mac produces the same result for scale but yet a different value - 2.75257143.

I wonder if the problem is in the numerical derivative implementation that's used by the LevMar algorithm. Because BlackBody does not provide a derivative a numerical one will be used.

The differences in results is a bit scary. Admittedly, the data is very simple and the model may not be well constrained. But it surprising that it is not repeatable between types of computers given the same fitter.

Let me know if you think more explanation is needed in the documentation than what's on the bottom of this page

https://docs.astropy.org/en/latest/modeling/add-units.html#fitting

Results of fitting a BlackBody model only agree on the two platforms. Const1D or * operator messing things up.

Let me know if you think more explanation is needed in the documentation than what's on the bottom of this page

https://docs.astropy.org/en/latest/modeling/add-units.html#fitting

How do I specify that tau_sil in the S07_attenuation model is unitless? That is the case for optical depths which is what tau_sil is.

I did this and the fitting works for my case! This is awesome!

    def _parameter_units_for_data_units(self, inputs_unit, outputs_unit):
        return {'tau_sil': None}

Are the results correct?

Are the results correct?

Not sure. Checking right now. Converting to units has changed things.

Check done. The fitting is returning results that don't make sense. Almost like it didn't update any of the input values. And I get the following for fit_info['message']

The relative error between two consecutive iterates is at most 0.000000

What happens if you use the Simplex fitter - SimplexLSQFitter()? It's a bit slower.

Doesn't work with SimplexLSQFitter. Different error.

TypeError: only dimensionless scalar quantities can be converted to Python scalars

Let's move this to slack next week.

The difference in the results on the two platforms may be due to an ill-conditioned model. I need to do some more verification with that model specifically, however if that's the case there's no guarantee we should get identical results. Regardless the principle issue here with handling the fitting of compound models is solved and I think this should be merged if approved.

Possibly a new issue was identified in this comment. Will follow up on it separately.

Did anyone do a plot vs the data for the two different fits on the two platforms? And do the corresponding plots for the solution on the other platform on the one that disagrees? In other words, such plots should definitively show whether there is an error in evaluation vs being poorly conditioned (the latter being confirmed by both sets of parameters showing good visual fits). Or is that not possible?