Crash when extracting flux points with a CompoundSpectralModel

[registerrier]

Gammapy version
v0.19

Bug description
If a compound model is used for the FluxPointsEstimator, an error occurs in FluxEstimator.get_scale_model() because a CompoundSpectralModel has no norm or amplitude attribute.

Expected behavior
Using a CompoundSpectralModel should work.

The issue is caused by the mechanism introduced to scale the norm according to the input parameter. It does not work in general for a CompoundSpectralModel which might have different amplitude and norm parameters. So that determining the actual range a priori might be ambiguous.

Note that we miss a test for flux points extraction with a compound model.

To Reproduce

from gammapy.modeling.models import EBLAbsorptionNormSpectralModel, PowerLawSpectralModel, SkyModel
from gammapy.estimators.flux import FluxEstimator

pl = PowerLawSpectralModel()
ebl = EBLAbsorptionNormSpectralModel.read_builtin(reference='finke')
model=SkyModel(spectral_model=pl*ebl)

fe = FluxEstimator()
fe.get_scale_model([model])

Other information
Any other information you think will be useful for us to fix the issue can go here.

[registerrier]

Note that a workaround for v0.19 is shown in #3705

[luca-giunti]

I think this still fails if one defines a custom model with a random name for the parameter to be scaled (e.g. Flux, or whatever). Also nothing prevents one to define a custom model with two parameters called amplitude and a norm, in which case we have a problem. I think our implicit convention on the names of parameters that can be scaled is not obvious, and should be better handled or documented.

Possible solutions that I see:

1. Adding an optional argument to FluxPointEstimator.__init__() to eventually specify the parameter to be scaled. But this can become problematic if one specifies e.g. index
2. Documentation. Specify that a custom model needs either an amplitude or a norm parameter. But this can become problematic e.g. if one defines a wrapper around a model coming from an external library, for which there is no freedom in the choice of parameter names. Fortunately the naima models have an amplitude, but this may not be the case in general
3. Define a flag on Parameter such as is_scalable_parameter, and adapt SkyModel to check that there is one and only one such parameter defined

I tend to prefer the option 3. Opinions @registerrier @adonath @AtreyeeS @QRemy ?

[QRemy]

In principle the scaled parameter does not have to be one of the existing parameters as it is already defined by the ScaleModel so in the get_scale_model method you pointed we could have :

        scaled_parameter = None
        if "amplitude" in ref_model.parameters.names:
            scaled_parameter = ref_model.parameters["amplitude"]
        elif "norm" in ref_model.parameters.names:
            scaled_parameter = ref_model.parameters["norm"]
        if scaled_parameter is not None:
            scale_model.norm = self._set_norm_parameter(scale_model.norm, scaled_parameter)

[AtreyeeS]

Thanks for bringing this up again @luca-giunti ! I think we orally discussed during the dev calls that we need to clarify in the documentation that either norm or amplitude must be present, and then completely forgot to document it!

I would go for (2) with the solution proposed by @QRemy . For external models without amplitude, you have to multiply it by a scale model.

[QRemy]

modifying _set_norm_parameter to deal correctly with the case where scaled_parameter=None:

   def _set_norm_parameter(self, norm=None, scaled_parameter=None):
        """Define properties of the norm spectral parameter."""
        if norm is None:
            norm = Parameter("norm", 1, unit="", interp="log")

        norm.value = 1.0
        norm.frozen = False

        if scaled_parameter is not None :
            norm.min = scaled_parameter.min / scaled_parameter.value
            norm.max = scaled_parameter.max / scaled_parameter.value
            norm.interp = scaled_parameter.interp

        norm.scan_values = self.norm_values
        norm.scan_min = self.norm_min
        norm.scan_max = self.norm_max
        norm.scan_n_values = self.norm_n_values
        return norm

    def get_scale_model(self, models):
        """Set scale model
        Parameters
        ----------
        models : `Models`
            Models
        Returns
        -------
        model : `ScaleSpectralModel`
            Scale spectral model
        """
        ref_model = models[self.source].spectral_model
        scale_model = ScaleSpectralModel(ref_model)

        scaled_parameter = None
        if "amplitude" in ref_model.parameters.names:
            scaled_parameter = ref_model.parameters["amplitude"]
        elif "norm" in ref_model.parameters.names:
            scaled_parameter = ref_model.parameters["norm"]

        scale_model.norm = self._set_norm_parameter(scale_model.norm, scaled_parameter)
        return scale_model

[luca-giunti]

I would go for (2) with the solution proposed by @QRemy . For external models without amplitude, you have to multiply it by a scale model.

@AtreyeeS Yes this also makes sense. In that case this could be documented here I guess?

modifying _set_norm_parameter to deal correctly with the case where scaled_parameter=None

@QRemy I just opened a PR on this exact point: #3845. Perhaps you can also have a look

[AtreyeeS]

Maybe it is better not to silently add scalable parameters but make the user do it ?

[adonath]

Users now have to define the is_norm property, when a SpectralModel is used with a SkyModel. See https://github.com/gammapy/gammapy/blob/master/gammapy/modeling/models/cube.py#L75 The example for the custom model was adapted accordingly: https://docs.gammapy.org/dev/tutorials/api/models.html#Implementing-a-custom-model