diff --git a/autolens/imaging/model/visualizer.py b/autolens/imaging/model/visualizer.py index 46d4ce9b86..2bdf1cb167 100644 --- a/autolens/imaging/model/visualizer.py +++ b/autolens/imaging/model/visualizer.py @@ -140,7 +140,7 @@ def visualize( try: fit.inversion.reconstruction except exc.InversionException: - logger( + logger.warning( ag.exc.invalid_linear_algebra_for_visualization_message() ) return diff --git a/autolens/interferometer/fit_interferometer.py b/autolens/interferometer/fit_interferometer.py index 2cb50cca99..cbc9152b21 100644 --- a/autolens/interferometer/fit_interferometer.py +++ b/autolens/interferometer/fit_interferometer.py @@ -342,6 +342,67 @@ def galaxy_image_dict(self) -> Dict[ag.Galaxy, np.ndarray]: return {**galaxy_image_dict, **galaxy_linear_obj_image_dict} + @property + def model_image_natural(self) -> aa.Array2D: + """ + The real-space (image-plane) model image `m` of the fit, on the dataset's `real_space_mask`: the lensed + image of every ordinary (non-linear) light profile (`profile_image`) plus, when the fit has an + inversion, the solved linear objects' reconstruction mapped to the image plane + (`inversion.mapped_reconstructed_data`, linear light profiles and pixelized sources) -- the real-space + image whose visibilities are `model_data`. + + It is built from these two terms rather than from `galaxy_image_dict`, whose entry for a galaxy with + both ordinary and linear light holds only the linear reconstruction. + + It needs neither visibilities nor a transformer, so it is available on an array-free dataset (built by + `Interferometer.from_stream` / `from_sparse_terms`), where it is the image the natural-weighted dirty + model image `dirty_model_image_natural` is formed from. (There `profile_image` is all zeros, because a + fit with ordinary light on an array-free dataset raises before it gets here.) + """ + image = np.asarray( + getattr(self.profile_image, "array", self.profile_image), dtype=np.float64 + ) + + if self.inversion is not None: + reconstruction = self.inversion.mapped_reconstructed_data + image = image + np.asarray( + getattr(reconstruction, "array", reconstruction), dtype=np.float64 + ) + + return aa.Array2D( + values=image, + mask=self.dataset.real_space_mask, + ) + + @property + def dirty_model_image_natural(self) -> aa.Array2D: + """ + The naturally weighted, normalised dirty image of the model visibilities, `W~ m / sum(w)`, formed from + `model_image_natural` with the dataset's `sparse_operator` (see + `autoarray.fit.fit_interferometer.dirty_model_image_natural_from`). + + It is the model counterpart of the dataset's `dirty_image_natural` and needs no visibilities, so it is + how a fit on an array-free dataset is visualized. It is available on any dataset carrying a + `sparse_operator` (array-free, or in-memory after `apply_sparse_operator()`); otherwise it raises an + `aa.exc.DatasetException`. + """ + return aa.fit.fit_interferometer.dirty_model_image_natural_from( + dataset=self.dataset, image=self.model_image_natural + ) + + @property + def dirty_residual_map_natural(self) -> aa.Array2D: + """ + The naturally weighted dirty residual map, `dirty_image_natural - dirty_model_image_natural`, which is + `Re(F^H W (d - F m)) / sum(w)`: the natural dirty image of the visibility residuals, computed without + them. + """ + return aa.Array2D( + values=np.asarray(self.dataset.dirty_image_natural.array) + - np.asarray(self.dirty_model_image_natural.array), + mask=self.dataset.real_space_mask, + ) + @property def galaxy_signal_to_noise_map_dict(self) -> Dict[ag.Galaxy, np.ndarray]: """ diff --git a/autolens/interferometer/model/visualizer.py b/autolens/interferometer/model/visualizer.py index 893c9d06b6..5bdc344088 100644 --- a/autolens/interferometer/model/visualizer.py +++ b/autolens/interferometer/model/visualizer.py @@ -49,9 +49,14 @@ def visualize_before_fit( positions = ag.Grid2DIrregular(positions_list) - plotter.image_with_positions( - image=analysis.dataset.dirty_image, positions=positions - ) + # An array-free dataset (from_stream / from_sparse_terms) has no visibilities to + # form the unweighted dirty image from, so its natural-weighted one is shown. + if analysis.dataset.is_array_free: + image = analysis.dataset.dirty_image_natural + else: + image = analysis.dataset.dirty_image + + plotter.image_with_positions(image=image, positions=positions) if analysis.adapt_images is not None: plotter.adapt_images(adapt_images=analysis.adapt_images) @@ -116,7 +121,7 @@ def visualize( source_plane_line_colors=sp_colors, ) except exc.InversionException: - logger(ag.exc.invalid_linear_algebra_for_visualization_message()) + logger.warning(ag.exc.invalid_linear_algebra_for_visualization_message()) return if quick_update: diff --git a/autolens/interferometer/plot/fit_interferometer_plots.py b/autolens/interferometer/plot/fit_interferometer_plots.py index 01162898ec..97b91f7ac6 100644 --- a/autolens/interferometer/plot/fit_interferometer_plots.py +++ b/autolens/interferometer/plot/fit_interferometer_plots.py @@ -142,6 +142,10 @@ def subplot_fit( * Dirty chi-squared map * Source plane image (full extent) + On an array-free dataset (``fit.dataset.is_array_free``) there are no visibilities, so a + 2 × 3 grid of the natural-weighted dirty image, dirty model image and dirty residual map + plus the zoomed and unzoomed source plane is written to the same filename instead. + Parameters ---------- fit : FitInterferometer @@ -164,6 +168,37 @@ def subplot_fit( ) _pf = (lambda t: f"{title_prefix.rstrip()} {t}") if title_prefix else (lambda t: t) + + if fit.dataset.is_array_free: + # An array-free dataset (from_stream / from_sparse_terms) has no visibilities, so the + # visibility-space and unweighted dirty panels are replaced by the natural-weighted + # dirty image / model / residual (from the sparse terms) and the source plane. + fig, axes = subplots(2, 3, figsize=conf_subplot_figsize(2, 3)) + axes_flat = list(axes.flatten()) + + plot_array(array=fit.dataset.dirty_image_natural, ax=axes_flat[0], + title=_pf("Dirty Image (Natural)"), colormap=colormap) + plot_array(array=fit.dirty_model_image_natural, ax=axes_flat[1], + title=_pf("Dirty Model Image (Natural)"), colormap=colormap, + lines=image_plane_lines, line_colors=image_plane_line_colors) + plot_array(array=fit.dirty_residual_map_natural, ax=axes_flat[2], + title=_pf("Dirty Residual Map (Natural)"), colormap=colormap) + _plot_source_plane(fit, axes_flat[3], final_plane_index, + zoom_to_brightest=True, colormap=colormap, + title=_pf("Source Plane (Zoomed)"), + lines=source_plane_lines, + line_colors=source_plane_line_colors) + _plot_source_plane(fit, axes_flat[4], final_plane_index, + zoom_to_brightest=False, colormap=colormap, + title=_pf("Source Plane (No Zoom)"), + lines=source_plane_lines, + line_colors=source_plane_line_colors) + axes_flat[5].axis("off") + + tight_layout() + save_figure(fig, path=output_path, filename="fit", format=output_format) + return + fig, axes = subplots(3, 4, figsize=conf_subplot_figsize(3, 4)) axes_flat = list(axes.flatten()) @@ -265,6 +300,9 @@ def subplot_fit_dirty_images( Dirty Image | Dirty Signal-To-Noise Map | Dirty Model Image (critical curves) Dirty Residual Map | Dirty Norm Residual Map | Dirty Chi-Squared Map + On an array-free dataset (``fit.dataset.is_array_free``) a 1 × 3 subplot of the + natural-weighted dirty image, dirty model image and dirty residual map is written instead. + Parameters ---------- fit : FitInterferometer @@ -286,6 +324,25 @@ def subplot_fit_dirty_images( ) _pf = (lambda t: f"{title_prefix.rstrip()} {t}") if title_prefix else (lambda t: t) + + if fit.dataset.is_array_free: + fig, axes = subplots(1, 3, figsize=conf_subplot_figsize(1, 3)) + axes_flat = list(axes.flatten()) + + plot_array(array=fit.dataset.dirty_image_natural, ax=axes_flat[0], + title=_pf("Dirty Image (Natural)"), colormap=colormap, + use_log10=use_log10) + plot_array(array=fit.dirty_model_image_natural, ax=axes_flat[1], + title=_pf("Dirty Model Image (Natural)"), colormap=colormap, + use_log10=use_log10, lines=image_plane_lines, + line_colors=image_plane_line_colors) + plot_array(array=fit.dirty_residual_map_natural, ax=axes_flat[2], + title=_pf("Dirty Residual Map (Natural)"), colormap=colormap) + + tight_layout() + save_figure(fig, path=output_path, filename="fit_dirty_images", format=output_format) + return + fig, axes = subplots(2, 3, figsize=conf_subplot_figsize(2, 3)) axes_flat = list(axes.flatten()) @@ -329,6 +386,9 @@ def subplot_fit_interferometer_combined( different panel choice because interferometer fits are most informatively visualised in dirty-image space. + A fit on an array-free dataset (``fit.dataset.is_array_free``) has no visibilities, so its row + shows the natural-weighted dirty image, dirty model image, source plane and dirty residual map. + Parameters ---------- fit_list : list of FitInterferometer @@ -362,17 +422,23 @@ def subplot_fit_interferometer_combined( tracer, cc_grid ) + array_free = fit.dataset.is_array_free + plot_array( - array=fit.dirty_image, + array=fit.dataset.dirty_image_natural if array_free else fit.dirty_image, ax=row_axes[0], - title=_pf(f"Dirty Image (ch {row})"), + title=_pf( + f"Dirty Image (Natural) (ch {row})" + if array_free + else f"Dirty Image (ch {row})" + ), colormap=colormap, ) plot_array( - array=fit.dirty_model_image, + array=fit.dirty_model_image_natural if array_free else fit.dirty_model_image, ax=row_axes[1], - title=_pf("Dirty Model Image"), + title=_pf("Dirty Model Image (Natural)" if array_free else "Dirty Model Image"), colormap=colormap, lines=ip_lines, line_colors=ip_colors, @@ -391,13 +457,21 @@ def subplot_fit_interferometer_combined( except Exception: row_axes[2].axis("off") - plot_array( - array=fit.dirty_normalized_residual_map, - ax=row_axes[3], - title=_pf("Dirty Norm Residual"), - colormap=colormap, - cb_unit=r"$\sigma$", - ) + if array_free: + plot_array( + array=fit.dirty_residual_map_natural, + ax=row_axes[3], + title=_pf("Dirty Residual Map (Natural)"), + colormap=colormap, + ) + else: + plot_array( + array=fit.dirty_normalized_residual_map, + ax=row_axes[3], + title=_pf("Dirty Norm Residual"), + colormap=colormap, + cb_unit=r"$\sigma$", + ) tight_layout() save_figure(fig, path=output_path, filename="fit_combined", format=output_format) @@ -456,8 +530,18 @@ def subplot_fit_real_space( lines=source_plane_lines, line_colors=source_plane_line_colors) else: # Pixelized source: dirty model image + source reconstruction - plot_array(array=fit.dirty_model_image, ax=axes_flat[0], - title=_pf("Reconstructed Image"), colormap=colormap) + # An array-free dataset has no transformer: its dirty model image is the + # natural-weighted one formed from the sparse terms. + plot_array( + array=( + fit.dirty_model_image_natural + if fit.dataset.is_array_free + else fit.dirty_model_image + ), + ax=axes_flat[0], + title=_pf("Reconstructed Image"), + colormap=colormap, + ) _plot_source_plane(fit, axes_flat[1], final_plane_index, zoom_to_brightest=True, colormap=colormap, title=_pf("Source Reconstruction"), @@ -532,9 +616,15 @@ def subplot_tracer_from_fit( axes_flat = list(axes.flatten()) # Panel 0: Dirty Model Image - plot_array(array=fit.dirty_model_image, ax=axes_flat[0], title=_pf("Dirty Model Image"), - lines=image_plane_lines, line_colors=image_plane_line_colors, - colormap=colormap) + if fit.dataset.is_array_free: + plot_array(array=fit.dirty_model_image_natural, ax=axes_flat[0], + title=_pf("Dirty Model Image (Natural)"), + lines=image_plane_lines, line_colors=image_plane_line_colors, + colormap=colormap) + else: + plot_array(array=fit.dirty_model_image, ax=axes_flat[0], title=_pf("Dirty Model Image"), + lines=image_plane_lines, line_colors=image_plane_line_colors, + colormap=colormap) # Panel 1: Lensed source image (image-plane projection). # Use galaxy_image_dict so that pixelized (inversion) sources are included. diff --git a/test_autolens/interferometer/model/files/dirty_images.fits b/test_autolens/interferometer/model/files/dirty_images.fits deleted file mode 100644 index 9aa1b9c8aa..0000000000 Binary files a/test_autolens/interferometer/model/files/dirty_images.fits and /dev/null differ diff --git a/test_autolens/interferometer/model/test_plotter_interferometer.py b/test_autolens/interferometer/model/test_plotter_interferometer.py index 8be31f7a9e..22f7496a2b 100644 --- a/test_autolens/interferometer/model/test_plotter_interferometer.py +++ b/test_autolens/interferometer/model/test_plotter_interferometer.py @@ -54,7 +54,163 @@ def test__fit_interferometer( assert image.shape == (5, 5) image = al.ndarray_via_fits_from( - file_path=plot_path / "dirty_images.fits", hdu=0 + file_path=plot_path / "fit_dirty_images.fits", hdu=0 ) assert image.shape == (5, 5) + + +def _array_free_dataset_from(dataset): + import autoarray as aa + + return aa.Interferometer.from_stream( + [(dataset.uv_wavelengths, dataset.data, dataset.noise_map)], + real_space_mask=dataset.real_space_mask, + transformer_class=type(dataset.transformer), + ) + + +def _pixelized_source_model(): + import autofit as af + + pixelization = al.Pixelization( + mesh=al.mesh.RectangularUniform(shape=(3, 3)), + regularization=al.reg.Constant(coefficient=1.0), + ) + + return af.Collection( + galaxies=af.Collection( + lens=af.Model( + al.Galaxy, + redshift=0.5, + mass=al.mp.Isothermal(centre=(0.0, 0.0), einstein_radius=1.0), + ), + source=af.Model(al.Galaxy, redshift=1.0, pixelization=pixelization), + ) + ) + + +def _visualize(dataset, image_path): + """ + Run the interferometer visualizer's `visualize_before_fit` and `visualize` for a lens + with a pixelized source (and positions) on `dataset`, as a non-linear search would. + """ + from types import SimpleNamespace + + from autolens.interferometer.model.visualizer import VisualizerInterferometer + + model = _pixelized_source_model() + instance = model.instance_from_prior_medians() + + analysis = al.AnalysisInterferometer( + dataset=dataset, + positions_likelihood_list=[ + al.PositionsLH( + positions=al.Grid2DIrregular([(1.0, 0.0), (-1.0, 0.0)]), + threshold=1.0, + ) + ], + use_jax=False, + ) + paths = SimpleNamespace(image_path=image_path, output_path=image_path) + + VisualizerInterferometer.visualize_before_fit( + analysis=analysis, paths=paths, model=model + ) + VisualizerInterferometer.visualize( + analysis=analysis, paths=paths, instance=instance, during_analysis=False + ) + + return analysis.fit_from(instance=instance) + + +def _ext_names_from(file_path): + from astropy.io import fits + + with fits.open(file_path) as hdu_list: + return [hdu.name for hdu in hdu_list] + + +_PNGS = ( + "dataset", + "image_with_positions", + "fit", + "fit_dirty_images", + "fit_real_space", + "tracer", + "inversion_0_0", +) + + +def test__visualizer__array_free_dataset(interferometer_7, tmp_path, plot_patch): + import numpy as np + + dataset = _array_free_dataset_from(interferometer_7) + + fit = _visualize(dataset=dataset, image_path=tmp_path) + + for filename in _PNGS: + assert str(tmp_path / f"{filename}.png") in plot_patch.paths, filename + + assert (tmp_path / "galaxy_images.fits").exists() + assert _ext_names_from(tmp_path / "fit_dirty_images.fits") == [ + "MASK", + "DIRTY_IMAGE_NATURAL", + "DIRTY_BEAM", + "DIRTY_MODEL_IMAGE_NATURAL", + "DIRTY_RESIDUAL_MAP_NATURAL", + ] + + for hdu, array in ( + (1, dataset.dirty_image_natural), + (2, dataset.dirty_beam), + (3, fit.dirty_model_image_natural), + (4, fit.dirty_residual_map_natural), + ): + np.testing.assert_allclose( + al.ndarray_via_fits_from( + file_path=tmp_path / "fit_dirty_images.fits", hdu=hdu + ), + array.native_for_fits, + rtol=1.0e-6, + atol=1.0e-12, + ) + + +def test__visualizer__in_memory_dataset__fit_dirty_images_unchanged( + interferometer_7, tmp_path, plot_patch +): + import numpy as np + + fit = _visualize(dataset=interferometer_7, image_path=tmp_path) + + for filename in _PNGS: + assert str(tmp_path / f"{filename}.png") in plot_patch.paths, filename + + assert _ext_names_from(tmp_path / "fit_dirty_images.fits") == [ + "MASK", + "DIRTY_IMAGE", + "DIRTY_NOISE_MAP", + "DIRTY_MODEL_IMAGE", + "DIRTY_RESIDUAL_MAP", + "DIRTY_NORMALIZED_RESIDUAL_MAP", + "DIRTY_CHI_SQUARED_MAP", + ] + + for hdu, array in enumerate( + ( + fit.dirty_image, + fit.dirty_noise_map, + fit.dirty_model_image, + fit.dirty_residual_map, + fit.dirty_normalized_residual_map, + fit.dirty_chi_squared_map, + ), + start=1, + ): + np.testing.assert_array_equal( + al.ndarray_via_fits_from( + file_path=tmp_path / "fit_dirty_images.fits", hdu=hdu + ), + np.asarray(array.native_for_fits), + ) diff --git a/test_autolens/interferometer/plot/test_fit_interferometer_plots.py b/test_autolens/interferometer/plot/test_fit_interferometer_plots.py index 1c20a439cf..a2684cd215 100644 --- a/test_autolens/interferometer/plot/test_fit_interferometer_plots.py +++ b/test_autolens/interferometer/plot/test_fit_interferometer_plots.py @@ -35,3 +35,79 @@ def test__subplot_fit_real_space( output_format="png", ) assert str(plot_path / "fit_real_space.png") in plot_patch.paths + + +def test__subplots__array_free_dataset( + interferometer_7, plot_path, plot_patch, monkeypatch +): + import autoarray as aa + import autolens as al + from autolens.interferometer.plot import fit_interferometer_plots + from autolens.interferometer.plot.fit_interferometer_plots import ( + subplot_fit_dirty_images, + subplot_fit_interferometer_combined, + subplot_tracer_from_fit, + ) + + dataset = aa.Interferometer.from_stream( + [(interferometer_7.uv_wavelengths, interferometer_7.data, interferometer_7.noise_map)], + real_space_mask=interferometer_7.real_space_mask, + transformer_class=type(interferometer_7.transformer), + ) + + lens = al.Galaxy( + redshift=0.5, mass=al.mp.Isothermal(centre=(0.0, 0.0), einstein_radius=1.0) + ) + source = al.Galaxy( + redshift=1.0, + pixelization=al.Pixelization( + mesh=al.mesh.RectangularUniform(shape=(3, 3)), + regularization=al.reg.Constant(coefficient=1.0), + ), + ) + + fit = al.FitInterferometer( + dataset=dataset, tracer=al.Tracer(galaxies=[lens, source]) + ) + + titles = [] + plot_array = fit_interferometer_plots.plot_array + + def _plot_array(*args, title=None, **kwargs): + titles.append(title) + return plot_array(*args, title=title, **kwargs) + + monkeypatch.setattr(fit_interferometer_plots, "plot_array", _plot_array) + + natural_titles = [ + "Dirty Image (Natural)", + "Dirty Model Image (Natural)", + "Dirty Residual Map (Natural)", + ] + + subplot_fit(fit=fit, output_path=plot_path, output_format="png") + assert str(plot_path / "fit.png") in plot_patch.paths + assert titles[:3] == natural_titles + + titles.clear() + subplot_fit_dirty_images(fit=fit, output_path=plot_path, output_format="png") + assert str(plot_path / "fit_dirty_images.png") in plot_patch.paths + assert titles == natural_titles + + titles.clear() + subplot_fit_interferometer_combined( + fit_list=[fit], output_path=plot_path, output_format="png" + ) + assert str(plot_path / "fit_combined.png") in plot_patch.paths + assert titles[0] == "Dirty Image (Natural) (ch 0)" + assert "Dirty Model Image (Natural)" in titles + assert "Dirty Residual Map (Natural)" in titles + + titles.clear() + subplot_fit_real_space(fit=fit, output_path=plot_path, output_format="png") + assert str(plot_path / "fit_real_space.png") in plot_patch.paths + assert titles[0] == "Reconstructed Image" + + titles.clear() + subplot_tracer_from_fit(fit=fit, output_path=plot_path, output_format="png") + assert titles[0] == "Dirty Model Image (Natural)" diff --git a/test_autolens/interferometer/test_fit_interferometer.py b/test_autolens/interferometer/test_fit_interferometer.py index 130fa1b188..b3bb6097a9 100644 --- a/test_autolens/interferometer/test_fit_interferometer.py +++ b/test_autolens/interferometer/test_fit_interferometer.py @@ -829,3 +829,135 @@ def test__fit_figure_of_merit__array_free_dataset__lens_light_profile__raises( with pytest.raises(aa.exc.DatasetException): al.FitInterferometer(dataset=dataset_array_free, tracer=tracer).figure_of_merit + + +def test__natural_dirty_images__array_free_pixelized_source_fit(interferometer_7): + dataset = _array_free_dataset_from(interferometer_7) + + fit = al.FitInterferometer(dataset=dataset, tracer=_pixelized_source_tracer()) + + assert fit.inversion.transformer is None + + model_image = fit.model_image_natural + + assert isinstance(model_image, aa.Array2D) + np.testing.assert_allclose( + model_image.array, + fit.inversion.mapped_reconstructed_data.array, + rtol=1.0e-12, + atol=1.0e-12 * np.abs(model_image.array).max(), + ) + + np.testing.assert_allclose( + fit.dirty_residual_map_natural.array, + dataset.dirty_image_natural.array - fit.dirty_model_image_natural.array, + rtol=1.0e-12, + ) + + # The in-memory sparse fit's natural dirty model image equals the transformer's + # natural dirty image of the model visibilities `F m`, and the array-free one. + dataset_memory = interferometer_7.apply_sparse_operator(use_jax=False) + + fit_memory = al.FitInterferometer( + dataset=dataset_memory, tracer=_pixelized_source_tracer() + ) + + expected = _natural_dirty_image_of_model_data( + dataset=dataset_memory, fit=fit_memory + ) + + np.testing.assert_allclose( + fit_memory.dirty_model_image_natural.array, + expected.array, + rtol=1.0e-8, + atol=1.0e-8 * np.abs(expected.array).max(), + ) + np.testing.assert_allclose( + fit.dirty_model_image_natural.array, + fit_memory.dirty_model_image_natural.array, + rtol=1.0e-6, + atol=1.0e-6 * np.abs(expected.array).max(), + ) + + +def _natural_dirty_image_of_model_data(dataset, fit): + """ + The independent reference for `fit.dirty_model_image_natural`: the transformer's natural-weighted dirty + image of the fit's actual model visibilities `fit.model_data`, `Re(F^H (w m_vis)) / sum(w)` with + `w = 1 / sigma^2` per component (the weighting of `Interferometer.dirty_image_natural`). + """ + noise_map = dataset.noise_map.array + visibilities = np.asarray(fit.model_data.array) + weighted = aa.Visibilities( + visibilities=visibilities.real * noise_map.real**-2.0 + + 1j * visibilities.imag * noise_map.imag**-2.0 + ) + return dataset.transformer.image_from(visibilities=weighted) / float( + np.sum(noise_map.real**-2.0) + ) + + +@pytest.mark.parametrize("linear_component", ["linear_light_profile", "pixelization"]) +def test__natural_dirty_images__galaxy_with_ordinary_and_linear_light__matches_model_data( + interferometer_7, linear_component +): + """ + A galaxy with both an ordinary light profile and a linear component (a lens with a linear light profile, + or a source with a pixelization) has only its linear reconstruction in `galaxy_image_dict`; + `model_image_natural` must still include the ordinary (lensed) light, so `dirty_model_image_natural` is + the natural dirty image of the fit's model visibilities `model_data`. + """ + dataset = interferometer_7.apply_sparse_operator(use_jax=False) + + pixelization = al.Pixelization( + mesh=al.mesh.RectangularUniform(shape=(3, 3)), + regularization=al.reg.Constant(coefficient=1.0e-3), + ) + + if linear_component == "linear_light_profile": + lens = al.Galaxy( + redshift=0.5, + bulge=al.lp.Sersic(intensity=0.001, centre=(0.05, 0.05)), + disk=al.lp_linear.Sersic(centre=(0.0, 0.1)), + mass=al.mp.Isothermal(centre=(0.0, 0.0), einstein_radius=1.0), + ) + source = al.Galaxy(redshift=1.0, bulge=al.lp.Sersic(intensity=0.001)) + else: + lens = al.Galaxy( + redshift=0.5, + mass=al.mp.Isothermal(centre=(0.0, 0.0), einstein_radius=1.0), + ) + source = al.Galaxy( + redshift=1.0, + bulge=al.lp.Sersic(intensity=1.0e-4), + pixelization=pixelization, + ) + + fit = al.FitInterferometer( + dataset=dataset, tracer=al.Tracer(galaxies=[lens, source]) + ) + + # Both the ordinary and the linear light contribute, so omitting either is detectable. + assert np.abs(fit.profile_image.array).max() > 0.0 + assert np.abs(fit.inversion.mapped_reconstructed_data.array).max() > 0.0 + + np.testing.assert_allclose( + fit.model_image_natural.array, + fit.profile_image.array + fit.inversion.mapped_reconstructed_data.array, + rtol=1.0e-12, + ) + + expected = _natural_dirty_image_of_model_data(dataset=dataset, fit=fit) + + np.testing.assert_allclose( + fit.dirty_model_image_natural.array, + expected.array, + rtol=1.0e-8, + atol=1.0e-8 * np.abs(expected.array).max(), + ) + np.testing.assert_allclose( + fit.dirty_residual_map_natural.array, + dataset.dirty_image_natural.array - expected.array, + rtol=1.0e-8, + atol=1.0e-8 * np.abs(expected.array).max(), + )