Repository: gammapy/gammapy
Pull request #5974 (2025-07-09 08:42:50)

Fix links label for CompoundSpectralModel
As discussed in the help channel the copy (and likely serialisation) breaks the link between CompoundSpectralModel. This PR resolves it by collecting correctly the  compound spectral models in `_set_models_link`


Patch:
diff --git a/examples/tutorials/analysis-1d/ebl.py b/examples/tutorials/analysis-1d/ebl.py
index a8b408835..87ab9b38d 100644
--- a/examples/tutorials/analysis-1d/ebl.py
+++ b/examples/tutorials/analysis-1d/ebl.py
@@ -39,9 +39,9 @@ from gammapy.modeling.models import (
 #
 # We will use 6 observations of the blazars PKS 2155-304 taken in 2008 by
 # H.E.S.S. when it was in a steady state. The data have already been
-# reduced to OGIP format `~gammapy.datasets.SpectrumDatasetOnOff` following the procedure
+# reduced to OGIP format `SpectrumDatasetOnOff` following the procedure
 # :doc:`/tutorials/analysis-1d/spectral_analysis` tutorial using a
-# reflected regions background estimation. The spectra and IRFs from the
+# `ReflectedRegions` background estimation. The spectra and IRFs from the
 # 6 observations have been stacked together.
 #
 # We will load this dataset as a `~gammapy.datasets.SpectrumDatasetOnOff` and proceed with
@@ -61,8 +61,8 @@ print(dataset)
 #
 # The observed spectrum is already attenuated due to the EBL. Assuming
 # that the intrinsic spectrum is a power law, the observed spectrum is a
-# `~gammapy.modeling.models.CompoundSpectralModel` given by the product of an
-# EBL model with the intrinsic model.
+# `gammapy.modeling.models.CompoundSpectralModel` given by the product of an EBL model with the
+# intrinsic model.
 #
 
 
@@ -146,11 +146,9 @@ flux_points_intrinsic = flux_points_obs.copy(
     reference_model=SkyModel(spectral_model=pwl)
 )
 
-######################################################################
 #
 print(flux_points_obs.reference_model)
 
-######################################################################
 #
 print(flux_points_intrinsic.reference_model)
 
diff --git a/examples/tutorials/analysis-1d/spectral_analysis.py b/examples/tutorials/analysis-1d/spectral_analysis.py
index 5771f7970..ad73e9180 100644
--- a/examples/tutorials/analysis-1d/spectral_analysis.py
+++ b/examples/tutorials/analysis-1d/spectral_analysis.py
@@ -75,8 +75,8 @@ In practice, we have to:
 
 - Create the necessary makers :
 
-  - the spectrum dataset maker : `~gammapy.makers.SpectrumDatasetMaker`
-  - the OFF background maker, here a `~gammapy.makers.ReflectedRegionsBackgroundMaker`
+  - the spectrum dataset maker : `~gammapy.makers.SpectrumDatasetMaker` -
+    the OFF background maker, here a `~gammapy.makers.ReflectedRegionsBackgroundMaker`
   - and the safe range maker : `~gammapy.makers.SafeMaskMaker`
 
 - Perform the data reduction loop. And for every observation:
@@ -126,12 +126,12 @@ from gammapy.modeling.models import (
     SkyModel,
     create_crab_spectral_model,
 )
-from gammapy.visualization import plot_spectrum_datasets_off_regions
 
 ######################################################################
 # Check setup
 # -----------
 from gammapy.utils.check import check_tutorials_setup
+from gammapy.visualization import plot_spectrum_datasets_off_regions
 
 check_tutorials_setup()
 
@@ -157,7 +157,7 @@ observations = datastore.get_observations(obs_ids)
 #
 # The next step is to define a signal extraction region, also known as on
 # region. In the simplest case this is just a
-# `~regions.CircleSkyRegion`.
+# `CircleSkyRegion <http://astropy-regions.readthedocs.io/en/latest/api/regions.CircleSkyRegion.html>`__.
 #
 
 target_position = SkyCoord(ra=83.63, dec=22.01, unit="deg", frame="icrs")
diff --git a/examples/tutorials/analysis-1d/spectral_analysis_hli.py b/examples/tutorials/analysis-1d/spectral_analysis_hli.py
index befbf53c3..051b6f367 100644
--- a/examples/tutorials/analysis-1d/spectral_analysis_hli.py
+++ b/examples/tutorials/analysis-1d/spectral_analysis_hli.py
@@ -34,7 +34,7 @@ Proposed approach
 -----------------
 
 This notebook uses the high level `~gammapy.analysis.Analysis` class to orchestrate data
-reduction and run the data fits. In its current state, `~gammapy.analysis.Analysis`
+reduction and run the data fits. In its current state, `Analysis`
 supports the standard analysis cases of joint or stacked 3D and 1D
 analyses. It is instantiated with an `~gammapy.analysis.AnalysisConfig` object that
 gives access to analysis parameters either directly or via a YAML config
@@ -130,7 +130,7 @@ print(config)
 
 ######################################################################
 # Note that you can save this string into a yaml file and load it as
-# follows:
+# follow:
 #
 
 # config = AnalysisConfig.read("config-1d.yaml")
@@ -222,8 +222,8 @@ analysis.get_observations()
 
 
 ######################################################################
-# The observations are now available on the `~gammapy.analysis.Analysis` object.
-# The selection corresponds to the following ids:
+# The observations are now available on the `Analysis` object. The
+# selection corresponds to the following ids:
 #
 
 print(analysis.observations.ids)
@@ -361,7 +361,7 @@ analysis.set_models(model_1d)
 # Setting fitting parameters
 # ~~~~~~~~~~~~~~~~~~~~~~~~~~
 #
-# `~gammapy.analysis.Analysis` can perform a few modeling and fitting tasks besides data
+# `Analysis` can perform a few modeling and fitting tasks besides data
 # reduction. Parameters have then to be passed to the configuration
 # object.
 #
@@ -475,7 +475,7 @@ plt.show()
 # ------------
 #
 # You can look at the same analysis without the high level interface in
-# :doc:`/tutorials/analysis-1d/spectral_analysis`.
+# :doc:`/tutorials/analysis-1d/spectral_analysis`
 #
 # As we can store the best model fit, you can overlay the fit results of
 # both methods on a unique plot.
diff --git a/examples/tutorials/analysis-1d/spectral_analysis_rad_max.py b/examples/tutorials/analysis-1d/spectral_analysis_rad_max.py
index 81314dece..01dc5994c 100644
--- a/examples/tutorials/analysis-1d/spectral_analysis_rad_max.py
+++ b/examples/tutorials/analysis-1d/spectral_analysis_rad_max.py
@@ -21,7 +21,7 @@ Context
 
 As already explained in the :doc:`/tutorials/analysis-1d/spectral_analysis`
 tutorial, the background is estimated from the field of view of the observation.
-In particular, the source and background events are counted within a circular
+In particular, the source and background events are counted within a circular 
 ON region enclosing the source. The background to be subtracted is then estimated
 from one or more OFF regions with an expected background rate similar to the one
 in the ON region (i.e. from regions with similar acceptance).
@@ -51,7 +51,7 @@ specifications offer two different ways to store this information:
   values are stored in additional `FITS` data unit, named
   `RAD_MAX_2D <https://gamma-astro-data-formats.readthedocs.io/en/latest/irfs/point_like/#rad-max-2d>`__.
 
-Gammapy provides a class to automatically read these values,
+`Gammapy` provides a class to automatically read these values,
 `~gammapy.irf.RadMax2D`, for both cases (fixed or energy-dependent
 :math:`\\theta` cut). In this notebook we will focus on how to perform a
 spectral extraction with a point-like IRF with an energy-dependent
@@ -79,8 +79,8 @@ We load two MAGIC observations in the
 IRF component with a :math:`\\theta` cut.
 
 We define the ON region, this time as a `~regions.PointSkyRegion` instead of a
-`~regions.CircleSkyRegion`, i.e. we specify only the center of our ON region.
-We create a `~gammapy.maps.RegionGeom` adding to the region the estimated energy
+`CircleSkyRegion`, i.e. we specify only the center of our ON region.
+We create a `RegionGeom` adding to the region the estimated energy
 axis of the `~gammapy.datasets.SpectrumDataset` object we want to
 produce. The corresponding dataset maker will automatically use the
 :math:`\\theta` values in `~gammapy.irf.RadMax2D` to set the
@@ -129,12 +129,12 @@ from gammapy.modeling.models import (
     SkyModel,
     create_crab_spectral_model,
 )
-from gammapy.visualization import plot_spectrum_datasets_off_regions
 
 ######################################################################
 # Check setup
 # -----------
 from gammapy.utils.check import check_tutorials_setup
+from gammapy.visualization import plot_spectrum_datasets_off_regions
 
 check_tutorials_setup()
 
@@ -206,14 +206,14 @@ dataset_empty = SpectrumDataset.create(geom=geom, energy_axis_true=energy_axis_t
 
 
 ######################################################################
-# The `~gammapy.datasets.SpectrumDataset` is now based on a geometry consisting of a
+# The `SpectrumDataset` is now based on a geometry consisting of a
 # single coordinate and an estimated energy axis. The
-# `~gammapy.datasets.SpectrumDatasetMaker` and `~gammapy.datasets.ReflectedRegionsBackgroundMaker` will
+# `SpectrumDatasetMaker` and `ReflectedRegionsBackgroundMaker` will
 # take care of producing ON and OFF with the proper sizes, automatically
-# adopting the :math:`\theta` values in `~gammapy.data.Observation.rad_max`.
+# adopting the :math:`\theta` values in `Observation.rad_max`.
 #
-# As explained in the introduction, we use a `~gammapy.makers.WobbleRegionsFinder`, to
-# determine the OFF positions. The parameter ``n_off_positions`` specifies
+# As explained in the introduction, we use a `WobbleRegionsFinder`, to
+# determine the OFF positions. The parameter `n_off_positions` specifies
 # the number of OFF regions to be considered.
 #
 
@@ -259,8 +259,9 @@ plt.show()
 # Fit spectrum
 # ------------
 #
-# We perform a joint likelihood fit of the two datasets. For this particular datasets
-# we select a fit range between :math:`80\,{\rm GeV}` and :math:`20\,{\rm TeV}`.
+# | We perform a joint likelihood fit of the two datasets.
+# | For this particular datasets we select a fit range between
+#   :math:`80\,{\rm GeV}` and :math:`20\,{\rm TeV}`.
 #
 
 e_min = 80 * u.GeV
diff --git a/examples/tutorials/analysis-1d/spectrum_simulation.py b/examples/tutorials/analysis-1d/spectrum_simulation.py
index 9ec763af1..f83585688 100644
--- a/examples/tutorials/analysis-1d/spectrum_simulation.py
+++ b/examples/tutorials/analysis-1d/spectrum_simulation.py
@@ -169,8 +169,8 @@ print(dataset)
 # ~~~~~~~~~~~~~~~
 #
 # To do an on off spectral analysis, which is the usual science case, the
-# standard would be to use `~gammapy.datasets.SpectrumDatasetOnOff`, which uses the
-# acceptance to fake off-counts. Please also refer to the ``Dataset simulations``
+# standard would be to use `SpectrumDatasetOnOff`, which uses the
+# acceptance to fake off-counts. Please also refer to the `Dataset simulations`
 # section in the :doc:`/tutorials/analysis-1d/spectral_analysis_rad_max` tutorial,
 # dealing with simulations based on observations of real off counts.
 #
diff --git a/examples/tutorials/analysis-2d/detect.py b/examples/tutorials/analysis-2d/detect.py
index 4a21c5129..52b8948f3 100644
--- a/examples/tutorials/analysis-2d/detect.py
+++ b/examples/tutorials/analysis-2d/detect.py
@@ -112,8 +112,8 @@ dataset = MapDataset(
 # `Ebeling et
 # al. (2006) <https://ui.adsabs.harvard.edu/abs/2006MNRAS.368...65E/abstract>`__.
 #
-# In the following example the `~gammapy.estimators.ASmoothMapEstimator.threshold`
-# argument gives the minimum significance expected, values below are clipped.
+# In the following example the `ASmoothMapEstimator.threshold` argument gives the minimum
+# significance expected, values below are clipped.
 #
 
 # %%time
@@ -185,7 +185,7 @@ maps = estimator.run(dataset)
 #
 # Below we print the result of the `~gammapy.estimators.TSMapEstimator`. We have access to a number of
 # different quantities, as shown below. We can also access the quantities names
-# through ``maps.available_quantities``.
+# through ``map_result.available_quantities``.
 #
 
 print(maps)
diff --git a/examples/tutorials/analysis-2d/modeling_2D.py b/examples/tutorials/analysis-2d/modeling_2D.py
index 431447620..f1de504a3 100644
--- a/examples/tutorials/analysis-2d/modeling_2D.py
+++ b/examples/tutorials/analysis-2d/modeling_2D.py
@@ -28,7 +28,7 @@ Proposed approach
 -----------------
 
 The first step here, as in most analysis with DL3 data, is to create
-reduced datasets. For this, we will use the `~gammapy.analysis.Analysis` class to create
+reduced datasets. For this, we will use the `Analysis` class to create
 a single set of stacked maps with a single bin in energy (thus, an
 *image* which behaves as a *cube*). This, we will then model with a
 spatial model of our choice, while keeping the spectral model fixed to
@@ -112,7 +112,7 @@ print(config)
 
 
 ######################################################################
-# We now use the config file and create a single `~gammapy.datasets.MapDataset` containing
+# We now use the config file and create a single `MapDataset` containing
 # `counts`, `background`, `exposure`, `psf` and `edisp` maps.
 #
 
diff --git a/examples/tutorials/data/cta.py b/examples/tutorials/data/cta.py
index f7f85eb9b..b1782aea0 100644
--- a/examples/tutorials/data/cta.py
+++ b/examples/tutorials/data/cta.py
@@ -100,7 +100,8 @@ check_tutorials_setup()
 # **For convenience**, since the 1DC data files are large and not publicly
 # available to anyone, we have taken a tiny subset of the CTA 1DC data,
 # four observations with the southern array from the GPS survey, pointing
-# near the Galactic center, and **included them at** `$GAMMAPY_DATA/cta-1dc`.
+# near the Galactic center, and **included them at `$GAMMAPY_DATA/cta-1dc`**
+# which you get via `gammapy download datasets`.
 #
 # Files
 # ~~~~~
@@ -142,7 +143,7 @@ print(data_store)
 
 ######################################################################
 # If you can’t download the index files, or got errors related to the data
-# access using them, you can generate the `~gammapy.data.DataStore` directly from the
+# access using them, you can generate the `DataStore` directly from the
 # event files.
 #
 
@@ -388,7 +389,7 @@ plt.show()
 #
 # -  Load the EVENTS file for `obs_id=111159` as a
 #    `~gammapy.data.EventList` object.
-# -  Use ``EventList.table`` to find the energy, sky coordinate and time of
+# -  Use `~gammapy.data.EventList.table` to find the energy, sky coordinate and time of
 #    the highest-energy event.
 # -  Use `~gammapy.data.EventList.pointing_radec` to find the pointing position of this
 #    observation, and use `astropy.coordinates.SkyCoord` methods to find
@@ -403,6 +404,8 @@ plt.show()
 #    `South_z20_50h` configurations.
 #
 
+# start typing here ...
+
 
 ######################################################################
 # Next steps
diff --git a/examples/tutorials/data/fermi_lat.py b/examples/tutorials/data/fermi_lat.py
index 3e7cfe70f..256ba6a66 100644
--- a/examples/tutorials/data/fermi_lat.py
+++ b/examples/tutorials/data/fermi_lat.py
@@ -177,7 +177,7 @@ plt.show()
 # in energy on the original exposure cube. Probably log interpolation
 # would be better, but it doesn’t matter much here, because the energy
 # binning is fine. Finally, we just copy the counts map geometry, which
-# contains an energy axis with ``node_type="edges"``. This is non-ideal
+# contains an energy axis with `node_type="edges"`. This is non-ideal
 # for exposure cubes, but again, acceptable because exposure doesn’t vary
 # much from bin to bin, so the exact way interpolation occurs in later use
 # of that exposure cube doesn’t matter a lot. Of course you could define
@@ -192,7 +192,7 @@ exposure_hpx.plot()
 plt.show()
 
 ######################################################################
-# For exposure, we choose a geometry with ``node_type='center'``,
+# For exposure, we choose a geometry with node_type='center',
 axis = MapAxis.from_energy_bounds(
     "10 GeV",
     "2 TeV",
diff --git a/examples/tutorials/starting/analysis_1.py b/examples/tutorials/starting/analysis_1.py
index 1d6d8bcba..fd19184f5 100644
--- a/examples/tutorials/starting/analysis_1.py
+++ b/examples/tutorials/starting/analysis_1.py
@@ -192,7 +192,7 @@ config.datasets.background.exclusion = (
 # Setting modeling and fitting parameters
 # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 #
-# `~gammapy.analysis.Analysis` can perform a few modeling and fitting tasks besides data
+# ``Analysis`` can perform a few modeling and fitting tasks besides data
 # reduction. Parameters have then to be passed to the configuration
 # object.
 #
@@ -209,7 +209,7 @@ config.excess_map.correlation_radius = 0.1 * u.deg
 
 ######################################################################
 # We’re all set. But before we go on let’s see how to save or import
-# `~gammapy.analysis.AnalysisConfig` objects though YAML files.
+# ``AnalysisConfig`` objects though YAML files.
 #
 
 
@@ -217,7 +217,7 @@ config.excess_map.correlation_radius = 0.1 * u.deg
 # Using YAML configuration files
 # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 #
-# One can export/import the `~gammapy.analysis.AnalysisConfig` to/from a YAML file.
+# One can export/import the `~gammapy.modeling.AnalysisConfig` to/from a YAML file.
 #
 
 config.write("config.yaml", overwrite=True)
@@ -249,7 +249,7 @@ analysis.get_observations()
 
 
 ######################################################################
-# The observations are now available on the `~gammapy.analysis.Analysis` object. The
+# The observations are now available on the `~gammapy.modeling.Analysis` object. The
 # selection corresponds to the following ids:
 #
 
@@ -258,7 +258,7 @@ print(analysis.observations.ids)
 
 ######################################################################
 # To see how to explore observations, please refer to the following
-# notebook: :doc:`CTAO with Gammapy </tutorials/data/cta>` or :doc:`H.E.S.S. with
+# notebook: :doc:`CTA with Gammapy </tutorials/data/cta>` or :doc:`H.E.S.S. with
 # Gammapy </tutorials/data/hess>`
 #
 
@@ -269,7 +269,7 @@ print(analysis.observations.ids)
 #
 # Now we proceed to the data reduction. In the config file we have chosen
 # a WCS map geometry, energy axis and decided to stack the maps. We can
-# run the reduction using `~gammapy.analysis.Analysis.get_datasets()`:
+# run the reduction using `~gammapy.modeling.Analysis.get_datasets()`:
 #
 
 # %%time
@@ -443,8 +443,8 @@ plt.show()
 # -----------
 #
 # You can look at the same analysis without the high level interface in
-# :doc:`/tutorials/starting/analysis_2`.
+# :doc:`/tutorials/starting/analysis_2`
 #
 # You can see how to perform a 1D spectral analysis of the same data in
-# :doc:`/tutorials/analysis-1d/spectral_analysis`.
+# :doc:`/tutorials/analysis-1d/spectral_analysis`
 #
diff --git a/examples/tutorials/starting/analysis_2.py b/examples/tutorials/starting/analysis_2.py
index 3299d7099..08e9fe7ef 100644
--- a/examples/tutorials/starting/analysis_2.py
+++ b/examples/tutorials/starting/analysis_2.py
@@ -116,7 +116,7 @@ data_store = DataStore.from_dir("$GAMMAPY_DATA/hess-dl3-dr1")
 # observations. Here we use a cone search which we define with a python
 # dict.
 #
-# We then filter the `~gammapy.data.ObservationTable` with
+# We then filter the `ObservationTable` with
 # `~gammapy.data.ObservationTable.select_observations`.
 #
 
@@ -316,7 +316,7 @@ print(result)
 
 ######################################################################
 # The fitted parameters are visible from the
-# `~gammapy.modeling.models.Models` object.
+# `~astropy.modeling.models.Models` object.
 #
 
 print(stacked.models.to_parameters_table())
diff --git a/examples/tutorials/starting/overview.py b/examples/tutorials/starting/overview.py
index 8d8640762..ea9be3257 100644
--- a/examples/tutorials/starting/overview.py
+++ b/examples/tutorials/starting/overview.py
@@ -66,6 +66,7 @@ We will cover the following topics:
 # setup is correct you can execute the following cell:
 #
 
+
 import astropy.units as u
 from astropy.coordinates import SkyCoord
 import matplotlib.pyplot as plt
@@ -120,7 +121,7 @@ print(gc_3fhl.geom)
 
 
 ######################################################################
-# Let’s take a closer look at the ``.data`` attribute:
+# Let’s take a closer look a the ``.data`` attribute:
 #
 
 print(gc_3fhl.data)
@@ -161,7 +162,7 @@ plt.show()
 # The smoothed plot already looks much nicer, but still the image is
 # rather large. As we are mostly interested in the inner part of the
 # image, we will cut out a quadratic region of the size 9 deg x 9 deg
-# around Vela. Therefore, we use `~gammapy.maps.Map.cutout` to make a
+# around Vela. Therefore we use `~gammapy.maps.Map.cutout` to make a
 # cutout map:
 #
 
@@ -173,13 +174,13 @@ plt.show()
 
 
 ######################################################################
-# For a more detailed introduction to `~gammapy.maps`, take a look at the
+# For a more detailed introduction to `~gammapy.maps`, take a look a the
 # :doc:`/tutorials/api/maps` notebook.
 #
 # Exercises
 # ~~~~~~~~~
 #
-# -  Add a marker and circle at the position of Sgr A* (you can find
+# -  Add a marker and circle at the position of ``Sgr A*`` (you can find
 #    examples in
 #    `astropy.visualization.wcsaxes <https://docs.astropy.org/en/stable/visualization/wcsaxes/>`__).
 #
@@ -196,7 +197,7 @@ plt.show()
 # In this section we will learn how to:
 #
 # -  Read event lists from FITS files
-# -  Access and work with the `~gammapy.data.EventList` attributes such as ``.table``
+# -  Access and work with the ``EventList`` attributes such as ``.table``
 #    and ``.energy``
 # -  Filter events lists using convenience methods
 #
@@ -214,7 +215,8 @@ events_3fhl = EventList.read("$GAMMAPY_DATA/fermi-3fhl-gc/fermi-3fhl-gc-events.f
 
 
 ######################################################################
-# This time the actual data is stored as an `~astropy.table.Table`
+# This time the actual data is stored as an
+# `~astropy.table.Table `
 # object. It can be accessed with ``.table`` attribute:
 #
 
@@ -222,7 +224,7 @@ print(events_3fhl.table)
 
 
 ######################################################################
-# You can do ``len`` over ``event_3fhl.table`` to find the total number of
+# You can do *len* over event_3fhl.table to find the total number of
 # events.
 #
 
@@ -230,7 +232,7 @@ print(len(events_3fhl.table))
 
 
 ######################################################################
-# And we can access any other attribute of the `~astropy.table.Table` object as well:
+# And we can access any other attribute of the ``Table`` object as well:
 #
 
 print(events_3fhl.table.colnames)
@@ -238,10 +240,13 @@ print(events_3fhl.table.colnames)
 
 ######################################################################
 # For convenience we can access the most important event parameters as
-# properties on the `~gammapy.data.EventList` objects. The attributes will return
+# properties on the ``EventList`` objects. The attributes will return
 # corresponding Astropy objects to represent the data, such as
-# `~astropy.units.Quantity`, `~astropy.coordinates.SkyCoord` or
-# `~astropy.time.Time` objects:
+# `~astropy.units.Quantity`,
+# `~astropy.coordinates.SkyCoord`
+# or
+# `~astropy.time.Time`
+# objects:
 #
 
 print(events_3fhl.energy.to("GeV"))
@@ -255,7 +260,7 @@ print(events_3fhl.time)
 
 
 ######################################################################
-# In addition `~gammapy.data.EventList` provides convenience methods to filter the
+# In addition ``EventList`` provides convenience methods to filter the
 # event lists. One possible use case is to find the highest energy event
 # within a radius of 0.5 deg around the vela position:
 #
@@ -313,7 +318,7 @@ print(fermi_3fhl.table)
 ######################################################################
 # This looks very familiar again. The data is just stored as an
 # `~astropy.table.Table`
-# object. We have all the methods and attributes of the `~astropy.table.Table` object
+# object. We have all the methods and attributes of the ``Table`` object
 # available. E.g. we can sort the underlying table by ``Signif_Avg`` to
 # find the top 5 most significant sources:
 #
diff --git a/gammapy/modeling/models/core.py b/gammapy/modeling/models/core.py
index b9d95c75b..881d9df79 100644
--- a/gammapy/modeling/models/core.py
+++ b/gammapy/modeling/models/core.py
@@ -55,22 +55,21 @@ def _set_model_link(shared_register, model):
     return shared_register
 
 
-def _set_models_link(models):
-    from . import SkyModel
-
-    shared_register = {}
+def _flatten_models(models):
+    flat_models = []
     for model in models:
-        if isinstance(model, SkyModel):
-            submodels = [
-                model.spectral_model,
-                model.spatial_model,
-                model.temporal_model,
-            ]
-            for submodel in submodels:
-                if submodel is not None:
-                    shared_register = _set_model_link(shared_register, submodel)
+        if hasattr(model, "_models"):
+            flat_models.extend(_flatten_models(model._models))
         else:
-            shared_register = _set_model_link(shared_register, model)
+            flat_models.append(model)
+    return flat_models
+
+
+def _set_models_link(models):
+    shared_register = {}
+    flat_models = _flatten_models(models)
+    for model in flat_models:
+        shared_register = _set_model_link(shared_register, model)
 
 
 def _get_model_class_from_dict(data):


Test patch:
diff --git a/gammapy/modeling/models/tests/test_io.py b/gammapy/modeling/models/tests/test_io.py
index 783e5f836..c09e7d744 100644
--- a/gammapy/modeling/models/tests/test_io.py
+++ b/gammapy/modeling/models/tests/test_io.py
@@ -19,6 +19,8 @@ from gammapy.modeling.models import (
     Model,
     Models,
     PiecewiseNormSpectralModel,
+    PointSpatialModel,
+    PowerLawNormSpectralModel,
     PowerLawSpectralModel,
     PowerLawTemporalModel,
     SineTemporalModel,
@@ -496,6 +498,38 @@ def test_link_label(models):
     assert n_link == 2
 
 
+def test_link_label_compound():
+    shared_pl = PowerLawSpectralModel(
+        index=2.2, amplitude="1e-12 cm-2 s-1 TeV-1", reference="1 TeV"
+    )
+    shared_point = PointSpatialModel()
+    norm1 = PowerLawNormSpectralModel(norm=1.0, tilt=0.0)
+    norm2 = PowerLawNormSpectralModel(norm=0.5, tilt=0.0)
+    model1 = SkyModel(
+        name="model-1", spatial_model=shared_point, spectral_model=shared_pl * norm1
+    )
+    model2 = SkyModel(
+        name="model-2",
+        spatial_model=shared_point,
+        spectral_model=shared_pl * norm2 * norm1,
+    )
+    models = Models([model1, model2])
+
+    assert len(models.parameters.free_unique_parameters.to_table()) == 6
+    assert len(models.parameters.free_unique_parameters) == 6
+
+    models_copy = models.copy()
+    assert len(models_copy.parameters.free_unique_parameters.to_table()) == 6
+    assert (
+        models_copy[0].spectral_model.model1.index
+        == models_copy[1].spectral_model.model1.model1.index
+    )
+    assert (
+        models_copy[0].spectral_model.model1.amplitude
+        == models_copy[1].spectral_model.model1.model1.amplitude
+    )
+
+
 def test_to_dict_not_default():
     model = PowerLawSpectralModel()
     model.index.min = -1