diff --git a/test_autoarray/fit/test_fit_dataset.py b/test_autoarray/fit/test_fit_dataset.py
index 674f25555..a9f24f490 100644
--- a/test_autoarray/fit/test_fit_dataset.py
+++ b/test_autoarray/fit/test_fit_dataset.py
@@ -3,7 +3,9 @@
 import autoarray as aa
 
 
-def test__figure_of_merit__with_inversion(masked_imaging_7x7, model_image_7x7):
+def test__figure_of_merit__with_inversion_and_regularization__equals_log_evidence(
+    masked_imaging_7x7, model_image_7x7
+):
     inversion = aa.m.MockInversion(
         linear_obj_list=[aa.m.MockMapper(regularization=aa.m.MockRegularization())],
         data_vector=1,
@@ -21,6 +23,10 @@ def test__figure_of_merit__with_inversion(masked_imaging_7x7, model_image_7x7):
 
     assert fit.figure_of_merit == fit.log_evidence
 
+
+def test__figure_of_merit__with_inversion_and_no_regularization__equals_log_likelihood(
+    masked_imaging_7x7, model_image_7x7
+):
     inversion = aa.m.MockInversion(
         linear_obj_list=[aa.m.MockLinearObj(regularization=None)], data_vector=1
     )
@@ -35,8 +41,8 @@ def test__figure_of_merit__with_inversion(masked_imaging_7x7, model_image_7x7):
     assert fit.figure_of_merit == fit.log_likelihood
 
 
-def test__figure_of_merit__with_noise_covariance_matrix_in_dataset(
-    masked_imaging_covariance_7x7, model_image_7x7, masked_imaging_7x7
+def test__chi_squared__with_noise_covariance_matrix__uses_covariance_weighted_formula(
+    masked_imaging_covariance_7x7, model_image_7x7
 ):
     fit = aa.m.MockFitImaging(
         dataset=masked_imaging_covariance_7x7,
@@ -51,21 +57,49 @@ def test__figure_of_merit__with_noise_covariance_matrix_in_dataset(
 
     assert fit.chi_squared == chi_squared
 
+
+def test__figure_of_merit__with_noise_covariance_matrix__equals_negative_half_chi_squared_plus_noise_normalization(
+    masked_imaging_covariance_7x7, model_image_7x7
+):
+    fit = aa.m.MockFitImaging(
+        dataset=masked_imaging_covariance_7x7,
+        use_mask_in_fit=False,
+        model_data=model_image_7x7,
+    )
+
     assert fit.figure_of_merit == pytest.approx(
         -0.5 * (fit.chi_squared + fit.noise_normalization), 1.0e-4
     )
 
-    fit = aa.m.MockFitImaging(
-        dataset=masked_imaging_7x7,
+
+def test__chi_squared__without_noise_covariance_matrix__differs_from_covariance_weighted_result(
+    masked_imaging_covariance_7x7, masked_imaging_7x7, model_image_7x7
+):
+    fit_with_covariance = aa.m.MockFitImaging(
+        dataset=masked_imaging_covariance_7x7,
         use_mask_in_fit=False,
         model_data=model_image_7x7,
     )
 
-    assert fit.chi_squared != pytest.approx(chi_squared, 1.0e-4)
+    chi_squared_covariance = aa.util.fit.chi_squared_with_noise_covariance_from(
+        residual_map=fit_with_covariance.residual_map,
+        noise_covariance_matrix_inv=masked_imaging_covariance_7x7.noise_covariance_matrix_inv,
+    )
 
+    fit_without_covariance = aa.m.MockFitImaging(
+        dataset=masked_imaging_7x7,
+        use_mask_in_fit=False,
+        model_data=model_image_7x7,
+    )
+
+    assert fit_without_covariance.chi_squared != pytest.approx(
+        chi_squared_covariance, 1.0e-4
+    )
 
-def test__grid_offset_via_data_model(imaging_7x7, mask_2d_7x7, model_image_7x7):
 
+def test__grids__with_dataset_model_grid_offset__lp_and_pixelization_grids_offset_correctly(
+    imaging_7x7, mask_2d_7x7, model_image_7x7
+):
     masked_imaging_7x7 = imaging_7x7.apply_mask(mask=mask_2d_7x7)
 
     fit = aa.m.MockFitImaging(
@@ -76,6 +110,5 @@ def test__grid_offset_via_data_model(imaging_7x7, mask_2d_7x7, model_image_7x7):
     )
 
     assert fit.dataset_model.grid_offset == (1.0, 2.0)
-
     assert fit.grids.lp[0] == pytest.approx((0.0, -3.0), 1.0e-4)
     assert fit.grids.pixelization[0] == pytest.approx((0.0, -3.0), 1.0e-4)
diff --git a/test_autoarray/fit/test_fit_imaging.py b/test_autoarray/fit/test_fit_imaging.py
index c90a453ac..411330d46 100644
--- a/test_autoarray/fit/test_fit_imaging.py
+++ b/test_autoarray/fit/test_fit_imaging.py
@@ -3,102 +3,205 @@
 import autoarray as aa
 
 
-def test__data_and_model_are_identical__no_masking__check_values_are_correct():
-    mask = aa.Mask2D(mask=[[False, False], [False, False]], pixel_scales=(1.0, 1.0))
+# ---------------------------------------------------------------------------
+# Helper: build the "identical model, no masking" fit used by multiple tests
+# ---------------------------------------------------------------------------
 
+def _make_identical_fit_no_mask():
+    mask = aa.Mask2D(mask=[[False, False], [False, False]], pixel_scales=(1.0, 1.0))
     data = aa.Array2D(values=[1.0, 2.0, 3.0, 4.0], mask=mask)
     noise_map = aa.Array2D(values=[2.0, 2.0, 2.0, 2.0], mask=mask)
-
     dataset = aa.Imaging(data=data, noise_map=noise_map)
-
     dataset = dataset.apply_mask(mask=mask)
-
     model_data = aa.Array2D(values=[1.0, 2.0, 3.0, 4.0], mask=mask)
+    fit = aa.m.MockFitImaging(
+        dataset=dataset, use_mask_in_fit=False, model_data=model_data
+    )
+    return fit, noise_map
 
+
+def _make_different_fit_with_mask():
+    mask = aa.Mask2D(mask=[[False, False], [True, False]], pixel_scales=(1.0, 1.0))
+    data = aa.Array2D(values=[1.0, 2.0, 4.0], mask=mask)
+    noise_map = aa.Array2D(values=[2.0, 2.0, 2.0], mask=mask)
+    dataset = aa.Imaging(data=data, noise_map=noise_map)
+    model_data = aa.Array2D(values=[1.0, 2.0, 3.0], mask=mask)
     fit = aa.m.MockFitImaging(
         dataset=dataset, use_mask_in_fit=False, model_data=model_data
     )
+    return fit, noise_map
+
+
+def _make_identical_fit_with_inversion():
+    mask = aa.Mask2D(mask=[[False, False], [False, False]], pixel_scales=(1.0, 1.0))
+    data = aa.Array2D(values=[1.0, 2.0, 3.0, 4.0], mask=mask)
+    noise_map = aa.Array2D(values=[2.0, 2.0, 2.0, 2.0], mask=mask)
+    dataset = aa.Imaging(data=data, noise_map=noise_map)
+    dataset = dataset.apply_mask(mask=mask)
+    model_data = aa.Array2D(values=[1.0, 2.0, 3.0, 4.0], mask=mask)
+    inversion = aa.m.MockInversion(
+        linear_obj_list=[aa.m.MockMapper()],
+        data_vector=1,
+        regularization_term=2.0,
+        log_det_curvature_reg_matrix_term=3.0,
+        log_det_regularization_matrix_term=4.0,
+    )
+    fit = aa.m.MockFitImaging(
+        dataset=dataset,
+        use_mask_in_fit=False,
+        model_data=model_data,
+        inversion=inversion,
+    )
+    return fit, noise_map
+
+
+# ---------------------------------------------------------------------------
+# Tests: identical data and model, no masking
+# ---------------------------------------------------------------------------
+
+
+def test__mask__no_masking__returns_2x2_all_false_mask():
+    fit, _ = _make_identical_fit_no_mask()
 
     assert (fit.mask == np.array([[False, False], [False, False]])).all()
+
+
+def test__data__no_masking__returns_correct_data_values():
+    fit, _ = _make_identical_fit_no_mask()
+
     assert (fit.data == np.array([1.0, 2.0, 3.0, 4.0])).all()
+
+
+def test__noise_map__no_masking__returns_correct_noise_map_values():
+    fit, _ = _make_identical_fit_no_mask()
+
     assert (fit.noise_map == np.array([2.0, 2.0, 2.0, 2.0])).all()
+
+
+def test__signal_to_noise_map__no_masking__returns_correct_signal_to_noise_values():
+    fit, _ = _make_identical_fit_no_mask()
+
     assert (fit.signal_to_noise_map == np.array([0.5, 1.0, 1.5, 2.0])).all()
-    assert (fit.model_data == np.array([1.0, 2.0, 3.0, 4.0])).all()
+
+
+def test__residual_map__identical_data_and_model__all_zero_residuals():
+    fit, _ = _make_identical_fit_no_mask()
+
     assert (fit.residual_map == np.array([0.0, 0.0, 0.0, 0.0])).all()
+
+
+def test__normalized_residual_map__identical_data_and_model__all_zero_normalized_residuals():
+    fit, _ = _make_identical_fit_no_mask()
+
     assert (fit.normalized_residual_map == np.array([0.0, 0.0, 0.0, 0.0])).all()
+
+
+def test__chi_squared_map__identical_data_and_model__all_zero_chi_squared_map():
+    fit, _ = _make_identical_fit_no_mask()
+
     assert (fit.chi_squared_map == np.array([0.0, 0.0, 0.0, 0.0])).all()
 
+
+def test__chi_squared__identical_data_and_model__is_zero():
+    fit, _ = _make_identical_fit_no_mask()
+
     assert fit.chi_squared == 0.0
+
+
+def test__reduced_chi_squared__identical_data_and_model__is_zero():
+    fit, _ = _make_identical_fit_no_mask()
+
     assert fit.reduced_chi_squared == 0.0
+
+
+def test__noise_normalization__uniform_noise_map__correct_log_sum_formula():
+    fit, noise_map = _make_identical_fit_no_mask()
+
     assert fit.noise_normalization == np.sum(np.log(2 * np.pi * noise_map.array**2.0))
-    assert fit.log_likelihood == -0.5 * (fit.chi_squared + fit.noise_normalization)
 
 
-def test__data_and_model_are_different__include_masking__check_values_are_correct():
-    mask = aa.Mask2D(mask=[[False, False], [True, False]], pixel_scales=(1.0, 1.0))
+def test__log_likelihood__identical_data_and_model__negative_half_noise_normalization():
+    fit, _ = _make_identical_fit_no_mask()
 
-    data = aa.Array2D(values=[1.0, 2.0, 4.0], mask=mask)
-    noise_map = aa.Array2D(values=[2.0, 2.0, 2.0], mask=mask)
+    assert fit.log_likelihood == -0.5 * (fit.chi_squared + fit.noise_normalization)
 
-    dataset = aa.Imaging(data=data, noise_map=noise_map)
 
-    model_data = aa.Array2D(values=[1.0, 2.0, 3.0], mask=mask)
+# ---------------------------------------------------------------------------
+# Tests: different data and model with partial mask
+# ---------------------------------------------------------------------------
 
-    fit = aa.m.MockFitImaging(
-        dataset=dataset, use_mask_in_fit=False, model_data=model_data
-    )
 
-    assert (fit.mask == np.array([[False, False], [True, False]])).all()
-    assert (fit.data.slim == np.array([1.0, 2.0, 4.0])).all()
-    assert (fit.noise_map.slim == np.array([2.0, 2.0, 2.0])).all()
-    assert (fit.signal_to_noise_map.slim == np.array([0.5, 1.0, 2.0])).all()
-    assert (fit.model_data.slim == np.array([1.0, 2.0, 3.0])).all()
+def test__residual_map__different_data_and_model_with_partial_mask__correct_slim_residuals():
+    fit, _ = _make_different_fit_with_mask()
+
     assert (fit.residual_map.slim == np.array([0.0, 0.0, 1.0])).all()
+
+
+def test__normalized_residual_map__different_data_and_model_with_partial_mask__correct_slim_values():
+    fit, _ = _make_different_fit_with_mask()
+
     assert (fit.normalized_residual_map.slim == np.array([0.0, 0.0, 0.5])).all()
+
+
+def test__chi_squared_map__different_data_and_model_with_partial_mask__correct_slim_chi_squared():
+    fit, _ = _make_different_fit_with_mask()
+
     assert (fit.chi_squared_map.slim == np.array([0.0, 0.0, 0.25])).all()
 
+
+def test__chi_squared__different_model_with_masked_data__correct_value():
+    fit, _ = _make_different_fit_with_mask()
+
     assert fit.chi_squared == 0.25
+
+
+def test__reduced_chi_squared__different_model_with_masked_data__divided_by_unmasked_pixel_count():
+    fit, _ = _make_different_fit_with_mask()
+
     assert fit.reduced_chi_squared == 0.25 / 3.0
+
+
+def test__log_likelihood__different_model_with_masked_data__correct_value():
+    fit, noise_map = _make_different_fit_with_mask()
+
     assert fit.noise_normalization == np.sum(np.log(2 * np.pi * noise_map.array**2.0))
     assert fit.log_likelihood == -0.5 * (fit.chi_squared + fit.noise_normalization)
 
 
-def test__data_and_model_are_identical__inversion_included__changes_certain_properties():
-    mask = aa.Mask2D(mask=[[False, False], [False, False]], pixel_scales=(1.0, 1.0))
+# ---------------------------------------------------------------------------
+# Tests: identical data and model with inversion
+# ---------------------------------------------------------------------------
 
-    data = aa.Array2D(values=[1.0, 2.0, 3.0, 4.0], mask=mask)
-    noise_map = aa.Array2D(values=[2.0, 2.0, 2.0, 2.0], mask=mask)
 
-    dataset = aa.Imaging(data=data, noise_map=noise_map)
+def test__chi_squared__identical_data_and_model_with_inversion__is_zero():
+    fit, _ = _make_identical_fit_with_inversion()
 
-    dataset = dataset.apply_mask(mask=mask)
+    assert fit.chi_squared == 0.0
 
-    model_data = aa.Array2D(values=[1.0, 2.0, 3.0, 4.0], mask=mask)
 
-    inversion = aa.m.MockInversion(
-        linear_obj_list=[aa.m.MockMapper()],
-        data_vector=1,
-        regularization_term=2.0,
-        log_det_curvature_reg_matrix_term=3.0,
-        log_det_regularization_matrix_term=4.0,
-    )
-
-    fit = aa.m.MockFitImaging(
-        dataset=dataset,
-        use_mask_in_fit=False,
-        model_data=model_data,
-        inversion=inversion,
-    )
+def test__reduced_chi_squared__identical_data_and_model_with_inversion__is_zero():
+    fit, _ = _make_identical_fit_with_inversion()
 
-    assert fit.chi_squared == 0.0
     assert fit.reduced_chi_squared == 0.0
-    assert fit.noise_normalization == np.sum(np.log(2 * np.pi * noise_map.array**2.0))
-    assert fit.log_likelihood == -0.5 * (fit.chi_squared + fit.noise_normalization)
+
+
+def test__log_likelihood_with_regularization__with_inversion__adds_regularization_term():
+    fit, noise_map = _make_identical_fit_with_inversion()
 
     assert fit.log_likelihood_with_regularization == -0.5 * (
         fit.chi_squared + 2.0 + fit.noise_normalization
     )
+
+
+def test__log_evidence__with_inversion__uses_chi_squared_reg_and_determinant_terms():
+    fit, noise_map = _make_identical_fit_with_inversion()
+
     assert fit.log_evidence == -0.5 * (
         fit.chi_squared + 2.0 + 3.0 - 4.0 + fit.noise_normalization
     )
+
+
+def test__figure_of_merit__with_inversion__equals_log_evidence():
+    fit, _ = _make_identical_fit_with_inversion()
+
     assert fit.figure_of_merit == fit.log_evidence
diff --git a/test_autoarray/fit/test_fit_interferometer.py b/test_autoarray/fit/test_fit_interferometer.py
index e22d9de19..97e7edf09 100644
--- a/test_autoarray/fit/test_fit_interferometer.py
+++ b/test_autoarray/fit/test_fit_interferometer.py
@@ -4,177 +4,306 @@
 import autoarray as aa
 
 
-def test__data_and_model_are_identical__no_masking__check_values_are_correct():
+# ---------------------------------------------------------------------------
+# Helpers
+# ---------------------------------------------------------------------------
+
+def _make_dataset():
     real_space_mask = aa.Mask2D(
         mask=[[False, False], [False, False]], pixel_scales=(1.0, 1.0)
     )
-
     data = aa.Visibilities(visibilities=[1.0 + 2.0j, 3.0 + 4.0j])
     noise_map = aa.VisibilitiesNoiseMap(visibilities=[2.0 + 2.0j, 2.0 + 2.0j])
-
     dataset = aa.Interferometer(
         data=data,
         noise_map=noise_map,
         uv_wavelengths=np.ones(shape=(2, 2)),
         real_space_mask=real_space_mask,
     )
+    return dataset, noise_map
 
+
+def _make_identical_fit():
+    dataset, noise_map = _make_dataset()
     model_data = aa.Visibilities(visibilities=[1.0 + 2.0j, 3.0 + 4.0j])
+    fit = aa.m.MockFitInterferometer(
+        dataset=dataset, use_mask_in_fit=False, model_data=model_data
+    )
+    return fit, noise_map
+
 
+def _make_different_fit():
+    dataset, noise_map = _make_dataset()
+    model_data = aa.Visibilities(visibilities=[1.0 + 2.0j, 3.0 + 3.0j])
     fit = aa.m.MockFitInterferometer(
         dataset=dataset, use_mask_in_fit=False, model_data=model_data
     )
+    return fit, noise_map
+
+
+def _make_identical_fit_with_inversion():
+    dataset, noise_map = _make_dataset()
+    data = dataset.data
+    model_data = aa.Visibilities(visibilities=[1.0 + 2.0j, 3.0 + 4.0j])
+    chi_squared = data - model_data
+    chi_squared = np.sum(
+        (chi_squared.real**2.0 / dataset.noise_map.real**2.0)
+        + (chi_squared.imag**2.0 / dataset.noise_map.imag**2.0)
+    )
+    inversion = aa.m.MockInversion(
+        linear_obj_list=[aa.m.MockMapper()],
+        data_vector=1,
+        regularization_term=2.0,
+        log_det_curvature_reg_matrix_term=3.0,
+        log_det_regularization_matrix_term=4.0,
+        fast_chi_squared=chi_squared,
+    )
+    fit = aa.m.MockFitInterferometer(
+        dataset=dataset,
+        use_mask_in_fit=False,
+        model_data=model_data,
+        inversion=inversion,
+    )
+    return fit, noise_map
+
+
+# ---------------------------------------------------------------------------
+# Tests: identical visibilities, no masking
+# ---------------------------------------------------------------------------
+
+
+def test__data__identical_visibilities__returns_correct_complex_data():
+    fit, _ = _make_identical_fit()
 
     assert (fit.data == np.array([1.0 + 2.0j, 3.0 + 4.0j])).all()
 
+
+def test__noise_map__uniform_complex_noise__returns_correct_noise_map():
+    fit, _ = _make_identical_fit()
+
     assert (fit.noise_map == np.array([2.0 + 2.0j, 2.0 + 2.0j])).all()
 
+
+def test__signal_to_noise_map__identical_visibilities__returns_correct_complex_snr():
+    fit, _ = _make_identical_fit()
+
     assert (fit.signal_to_noise_map == np.array([0.5 + 1.0j, 1.5 + 2.0j])).all()
 
-    assert (fit.model_data == np.array([1.0 + 2.0j, 3.0 + 4.0j])).all()
+
+def test__residual_map__identical_visibilities__all_zero_residuals():
+    fit, _ = _make_identical_fit()
 
     assert (fit.residual_map == np.array([0.0 + 0.0j, 0.0 + 0.0j])).all()
 
+
+def test__normalized_residual_map__identical_visibilities__all_zero_normalized_residuals():
+    fit, _ = _make_identical_fit()
+
     assert (fit.normalized_residual_map == np.array([0.0 + 0.0j, 0.0 + 0.0j])).all()
 
+
+def test__chi_squared_map__identical_visibilities__all_zero_chi_squared():
+    fit, _ = _make_identical_fit()
+
     assert (fit.chi_squared_map == np.array([0.0 + 0.0j, 0.0 + 0.0j])).all()
 
+
+def test__chi_squared__identical_visibilities__is_zero():
+    fit, _ = _make_identical_fit()
+
     assert fit.chi_squared == 0.0
+
+
+def test__reduced_chi_squared__identical_visibilities__is_zero():
+    fit, _ = _make_identical_fit()
+
     assert fit.reduced_chi_squared == 0.0
+
+
+def test__noise_normalization__uniform_complex_noise__correct_log_formula():
+    fit, _ = _make_identical_fit()
+
     assert fit.noise_normalization == pytest.approx(
         4.0 * np.log(2 * np.pi * 2.0**2.0), 1.0e-4
     )
-    assert fit.log_likelihood == -0.5 * (fit.chi_squared + fit.noise_normalization)
 
 
-def test__data_and_model_are_different__no_masking__check_values_are_correct():
-    real_space_mask = aa.Mask2D(
-        mask=[[False, False], [False, False]], pixel_scales=(1.0, 1.0)
-    )
+def test__log_likelihood__identical_visibilities__negative_half_noise_normalization():
+    fit, _ = _make_identical_fit()
 
-    data = aa.Visibilities(visibilities=[1.0 + 2.0j, 3.0 + 4.0j])
-    noise_map = aa.VisibilitiesNoiseMap(visibilities=[2.0 + 2.0j, 2.0 + 2.0j])
+    assert fit.log_likelihood == -0.5 * (fit.chi_squared + fit.noise_normalization)
 
-    dataset = aa.Interferometer(
-        data=data,
-        noise_map=noise_map,
-        uv_wavelengths=np.ones(shape=(2, 2)),
-        real_space_mask=real_space_mask,
-    )
 
-    model_data = aa.Visibilities(visibilities=[1.0 + 2.0j, 3.0 + 3.0j])
+# ---------------------------------------------------------------------------
+# Tests: different visibilities, no masking
+# ---------------------------------------------------------------------------
 
-    fit = aa.m.MockFitInterferometer(
-        dataset=dataset, use_mask_in_fit=False, model_data=model_data
-    )
 
-    assert (fit.data == np.array([1.0 + 2.0j, 3.0 + 4.0j])).all()
-    assert (fit.noise_map == np.array([2.0 + 2.0j, 2.0 + 2.0j])).all()
-    assert (fit.signal_to_noise_map == np.array([0.5 + 1.0j, 1.5 + 2.0j])).all()
-    assert (fit.model_data == np.array([1.0 + 2.0j, 3.0 + 3.0j])).all()
+def test__residual_map__different_visibilities__correct_complex_residuals():
+    fit, _ = _make_different_fit()
+
     assert (fit.residual_map == np.array([0.0 + 0.0j, 0.0 + 1.0j])).all()
+
+
+def test__normalized_residual_map__different_visibilities__correct_complex_normalized_residuals():
+    fit, _ = _make_different_fit()
+
     assert (fit.normalized_residual_map == np.array([0.0 + 0.0j, 0.0 + 0.5j])).all()
+
+
+def test__chi_squared_map__different_visibilities__correct_complex_chi_squared_map():
+    fit, _ = _make_different_fit()
+
     assert (fit.chi_squared_map == np.array([0.0 + 0.0j, 0.0 + 0.25j])).all()
+
+
+def test__chi_squared__different_visibilities__correct_value():
+    fit, _ = _make_different_fit()
+
     assert fit.chi_squared == 0.25
+
+
+def test__reduced_chi_squared__different_visibilities__divided_by_visibility_count():
+    fit, _ = _make_different_fit()
+
     assert fit.reduced_chi_squared == 0.25 / 2.0
+
+
+def test__log_likelihood__different_visibilities__correct_value():
+    fit, _ = _make_different_fit()
+
     assert fit.noise_normalization == pytest.approx(
         4.0 * np.log(2 * np.pi * 2.0**2.0), 1.0e-4
     )
     assert fit.log_likelihood == -0.5 * (fit.chi_squared + fit.noise_normalization)
 
 
-def test__data_and_model_are_identical__inversion_included__changes_certain_properties():
-    real_space_mask = aa.Mask2D(
-        mask=[[False, False], [False, False]], pixel_scales=(1.0, 1.0)
-    )
-
-    data = aa.Visibilities(visibilities=[1.0 + 2.0j, 3.0 + 4.0j])
-    noise_map = aa.VisibilitiesNoiseMap(visibilities=[2.0 + 2.0j, 2.0 + 2.0j])
+# ---------------------------------------------------------------------------
+# Tests: identical visibilities with inversion
+# ---------------------------------------------------------------------------
 
-    dataset = aa.Interferometer(
-        data=data,
-        noise_map=noise_map,
-        uv_wavelengths=np.ones(shape=(2, 2)),
-        real_space_mask=real_space_mask,
-    )
 
-    model_data = aa.Visibilities(visibilities=[1.0 + 2.0j, 3.0 + 4.0j])
+def test__chi_squared__identical_visibilities_with_inversion__is_zero():
+    fit, _ = _make_identical_fit_with_inversion()
 
-    chi_squared = data - model_data
-    chi_squared = np.sum(
-        (chi_squared.real**2.0 / noise_map.real**2.0)
-        + (chi_squared.imag**2.0 / noise_map.imag**2.0)
-    )
+    assert fit.chi_squared == 0.0
 
-    inversion = aa.m.MockInversion(
-        linear_obj_list=[aa.m.MockMapper()],
-        data_vector=1,
-        regularization_term=2.0,
-        log_det_curvature_reg_matrix_term=3.0,
-        log_det_regularization_matrix_term=4.0,
-        fast_chi_squared=chi_squared,
-    )
 
-    fit = aa.m.MockFitInterferometer(
-        dataset=dataset,
-        use_mask_in_fit=False,
-        model_data=model_data,
-        inversion=inversion,
-    )
+def test__reduced_chi_squared__identical_visibilities_with_inversion__is_zero():
+    fit, _ = _make_identical_fit_with_inversion()
 
-    assert fit.chi_squared == 0.0
     assert fit.reduced_chi_squared == 0.0
-    assert fit.noise_normalization == pytest.approx(
-        4.0 * np.log(2 * np.pi * 2.0**2.0), 1.0e-4
-    )
-    assert fit.log_likelihood == -0.5 * (fit.chi_squared + fit.noise_normalization)
+
+
+def test__log_likelihood_with_regularization__interferometer_with_inversion__adds_regularization_term():
+    fit, _ = _make_identical_fit_with_inversion()
 
     assert fit.log_likelihood_with_regularization == -0.5 * (
         fit.chi_squared + 2.0 + fit.noise_normalization
     )
+
+
+def test__log_evidence__interferometer_with_inversion__uses_chi_squared_reg_and_determinant_terms():
+    fit, _ = _make_identical_fit_with_inversion()
+
     assert fit.log_evidence == -0.5 * (
         fit.chi_squared + 2.0 + 3.0 - 4.0 + fit.noise_normalization
     )
+
+
+def test__figure_of_merit__interferometer_with_inversion__equals_log_evidence():
+    fit, _ = _make_identical_fit_with_inversion()
+
     assert fit.figure_of_merit == fit.log_evidence
 
 
-def test__dirty_quantities(transformer_7x7_7, interferometer_7, fit_interferometer_7):
+# ---------------------------------------------------------------------------
+# Tests: dirty image quantities (transformer applied to fit quantities)
+# ---------------------------------------------------------------------------
+
+
+def test__dirty_image__equals_transformer_image_from_interferometer_data(
+    transformer_7x7_7, interferometer_7, fit_interferometer_7
+):
     fit_interferometer_7.dataset.transformer = transformer_7x7_7
 
     dirty_image = transformer_7x7_7.image_from(visibilities=interferometer_7.data)
+
     assert (fit_interferometer_7.dirty_image == dirty_image).all()
 
+
+def test__dirty_noise_map__equals_transformer_image_from_noise_map(
+    transformer_7x7_7, interferometer_7, fit_interferometer_7
+):
+    fit_interferometer_7.dataset.transformer = transformer_7x7_7
+
     dirty_noise_map = transformer_7x7_7.image_from(
         visibilities=interferometer_7.noise_map
     )
+
     assert (fit_interferometer_7.dirty_noise_map == dirty_noise_map).all()
 
+
+def test__dirty_signal_to_noise_map__equals_transformer_image_from_signal_to_noise_map(
+    transformer_7x7_7, interferometer_7, fit_interferometer_7
+):
+    fit_interferometer_7.dataset.transformer = transformer_7x7_7
+
     dirty_signal_to_noise_map = transformer_7x7_7.image_from(
         visibilities=interferometer_7.signal_to_noise_map
     )
+
     assert (
         fit_interferometer_7.dirty_signal_to_noise_map == dirty_signal_to_noise_map
     ).all()
 
+
+def test__dirty_model_image__equals_transformer_image_from_model_data(
+    transformer_7x7_7, interferometer_7, fit_interferometer_7
+):
+    fit_interferometer_7.dataset.transformer = transformer_7x7_7
+
     dirty_model_image = transformer_7x7_7.image_from(
         visibilities=fit_interferometer_7.model_data
     )
+
     assert (fit_interferometer_7.dirty_model_image == dirty_model_image).all()
 
+
+def test__dirty_residual_map__equals_transformer_image_from_residual_map(
+    transformer_7x7_7, interferometer_7, fit_interferometer_7
+):
+    fit_interferometer_7.dataset.transformer = transformer_7x7_7
+
     dirty_residual_map = transformer_7x7_7.image_from(
         visibilities=fit_interferometer_7.residual_map
     )
+
     assert (fit_interferometer_7.dirty_residual_map == dirty_residual_map).all()
 
+
+def test__dirty_normalized_residual_map__equals_transformer_image_from_normalized_residual_map(
+    transformer_7x7_7, interferometer_7, fit_interferometer_7
+):
+    fit_interferometer_7.dataset.transformer = transformer_7x7_7
+
     dirty_normalized_residual_map = transformer_7x7_7.image_from(
         visibilities=fit_interferometer_7.normalized_residual_map
     )
+
     assert (
         fit_interferometer_7.dirty_normalized_residual_map
         == dirty_normalized_residual_map
     ).all()
 
+
+def test__dirty_chi_squared_map__equals_transformer_image_from_chi_squared_map(
+    transformer_7x7_7, interferometer_7, fit_interferometer_7
+):
+    fit_interferometer_7.dataset.transformer = transformer_7x7_7
+
     dirty_chi_squared_map = transformer_7x7_7.image_from(
         visibilities=fit_interferometer_7.chi_squared_map
     )
+
     assert (fit_interferometer_7.dirty_chi_squared_map == dirty_chi_squared_map).all()
diff --git a/test_autoarray/fit/test_fit_util.py b/test_autoarray/fit/test_fit_util.py
index 0702dc95a..c38c572f2 100644
--- a/test_autoarray/fit/test_fit_util.py
+++ b/test_autoarray/fit/test_fit_util.py
@@ -4,7 +4,12 @@
 import autoarray as aa
 
 
-def test__residual_map_from():
+# ---------------------------------------------------------------------------
+# residual_map_from
+# ---------------------------------------------------------------------------
+
+
+def test__residual_map_from__identical_data_and_model__all_zero_residuals():
     data = np.array([10.0, 10.0, 10.0, 10.0])
     model_data = np.array([10.0, 10.0, 10.0, 10.0])
 
@@ -12,6 +17,8 @@ def test__residual_map_from():
 
     assert (residual_map == np.array([0.0, 0.0, 0.0, 0.0])).all()
 
+
+def test__residual_map_from__different_data_and_model__correct_signed_residuals():
     data = np.array([10.0, 10.0, 10.0, 10.0])
     model_data = np.array([11.0, 10.0, 9.0, 8.0])
 
@@ -20,7 +27,12 @@ def test__residual_map_from():
     assert (residual_map == np.array([-1.0, 0.0, 1.0, 2.0])).all()
 
 
-def test__residual_map_with_mask_from():
+# ---------------------------------------------------------------------------
+# residual_map_with_mask_from
+# ---------------------------------------------------------------------------
+
+
+def test__residual_map_with_mask_from__masked_pixels__zero_at_masked_locations():
     data = np.array([10.0, 10.0, 10.0, 10.0])
     mask = np.array([True, False, False, True])
     model_data = np.array([11.0, 10.0, 9.0, 8.0])
@@ -32,13 +44,17 @@ def test__residual_map_with_mask_from():
     assert (residual_map == np.array([0.0, 0.0, 1.0, 0.0])).all()
 
 
-def test__normalized_residual_map_from():
+# ---------------------------------------------------------------------------
+# normalized_residual_map_from
+# ---------------------------------------------------------------------------
+
+
+def test__normalized_residual_map_from__identical_data_and_model__all_zero():
     data = np.array([10.0, 10.0, 10.0, 10.0])
     noise_map = np.array([2.0, 2.0, 2.0, 2.0])
     model_data = np.array([10.0, 10.0, 10.0, 10.0])
 
     residual_map = aa.util.fit.residual_map_from(data=data, model_data=model_data)
-
     normalized_residual_map = aa.util.fit.normalized_residual_map_from(
         residual_map=residual_map, noise_map=noise_map
     )
@@ -47,10 +63,13 @@ def test__normalized_residual_map_from():
         np.array([0.0, 0.0, 0.0, 0.0]), 1.0e-4
     )
 
+
+def test__normalized_residual_map_from__different_data_and_model__divided_by_noise_map():
+    data = np.array([10.0, 10.0, 10.0, 10.0])
+    noise_map = np.array([2.0, 2.0, 2.0, 2.0])
     model_data = np.array([11.0, 10.0, 9.0, 8.0])
 
     residual_map = aa.util.fit.residual_map_from(data=data, model_data=model_data)
-
     normalized_residual_map = aa.util.fit.normalized_residual_map_from(
         residual_map=residual_map, noise_map=noise_map
     )
@@ -60,7 +79,12 @@ def test__normalized_residual_map_from():
     )
 
 
-def test__normalized_residual_map_with_mask_from():
+# ---------------------------------------------------------------------------
+# normalized_residual_map_with_mask_from
+# ---------------------------------------------------------------------------
+
+
+def test__normalized_residual_map_with_mask_from__masked_pixels__zero_at_masked_locations():
     data = np.array([10.0, 10.0, 10.0, 10.0])
     mask = np.array([True, False, False, True])
     noise_map = np.array([2.0, 2.0, 2.0, 2.0])
@@ -69,7 +93,6 @@ def test__normalized_residual_map_with_mask_from():
     residual_map = aa.util.fit.residual_map_with_mask_from(
         data=data, mask=mask, model_data=model_data
     )
-
     normalized_residual_map = aa.util.fit.normalized_residual_map_with_mask_from(
         residual_map=residual_map, mask=mask, noise_map=noise_map
     )
@@ -79,13 +102,17 @@ def test__normalized_residual_map_with_mask_from():
     )
 
 
-def test__normalized_residual_map_complex_from():
+# ---------------------------------------------------------------------------
+# normalized_residual_map_complex_from
+# ---------------------------------------------------------------------------
+
+
+def test__normalized_residual_map_complex_from__complex_data__correct_real_and_imaginary_parts():
     data = np.array([10.0 + 10.0j, 10.0 + 10.0j])
     noise_map = np.array([2.0 + 2.0j, 2.0 + 2.0j])
     model_data = np.array([9.0 + 12.0j, 9.0 + 12.0j])
 
     residual_map = aa.util.fit.residual_map_from(data=data, model_data=model_data)
-
     normalized_residual_map = aa.util.fit.normalized_residual_map_complex_from(
         residual_map=residual_map, noise_map=noise_map
     )
@@ -93,23 +120,30 @@ def test__normalized_residual_map_complex_from():
     assert (normalized_residual_map == np.array([0.5 - 1.0j, 0.5 - 1.0j])).all()
 
 
-def test__chi_squared_map_from():
+# ---------------------------------------------------------------------------
+# chi_squared_map_from
+# ---------------------------------------------------------------------------
+
+
+def test__chi_squared_map_from__identical_data_and_model__all_zero():
     data = np.array([10.0, 10.0, 10.0, 10.0])
     noise_map = np.array([2.0, 2.0, 2.0, 2.0])
     model_data = np.array([10.0, 10.0, 10.0, 10.0])
 
     residual_map = aa.util.fit.residual_map_from(data=data, model_data=model_data)
-
     chi_squared_map = aa.util.fit.chi_squared_map_from(
         residual_map=residual_map, noise_map=noise_map
     )
 
     assert (chi_squared_map == np.array([0.0, 0.0, 0.0, 0.0])).all()
 
+
+def test__chi_squared_map_from__different_data_and_model__squared_normalized_residuals():
+    data = np.array([10.0, 10.0, 10.0, 10.0])
+    noise_map = np.array([2.0, 2.0, 2.0, 2.0])
     model_data = np.array([11.0, 10.0, 9.0, 8.0])
 
     residual_map = aa.util.fit.residual_map_from(data=data, model_data=model_data)
-
     chi_squared_map = aa.util.fit.chi_squared_map_from(
         residual_map=residual_map, noise_map=noise_map
     )
@@ -120,7 +154,12 @@ def test__chi_squared_map_from():
     ).all()
 
 
-def test__chi_squared_map_with_mask_from():
+# ---------------------------------------------------------------------------
+# chi_squared_map_with_mask_from
+# ---------------------------------------------------------------------------
+
+
+def test__chi_squared_map_with_mask_from__masked_pixels__zero_at_masked_locations():
     data = np.array([10.0, 10.0, 10.0, 10.0])
     mask = np.array([True, False, False, True])
     noise_map = np.array([2.0, 2.0, 2.0, 2.0])
@@ -129,33 +168,24 @@ def test__chi_squared_map_with_mask_from():
     residual_map = aa.util.fit.residual_map_with_mask_from(
         data=data, mask=mask, model_data=model_data
     )
-
     chi_squared_map = aa.util.fit.chi_squared_map_with_mask_from(
         residual_map=residual_map, mask=mask, noise_map=noise_map
     )
 
     assert (chi_squared_map == np.array([0.0, 0.0, (1.0 / 2.0) ** 2.0, 0.0])).all()
 
-    model_data = np.array([11.0, 10.0, 9.0, 8.0])
-
-    residual_map = aa.util.fit.residual_map_with_mask_from(
-        data=data, mask=mask, model_data=model_data
-    )
-
-    chi_squared_map = aa.util.fit.chi_squared_map_with_mask_from(
-        residual_map=residual_map, mask=mask, noise_map=noise_map
-    )
 
-    assert (chi_squared_map == np.array([0.0, 0.0, (1.0 / 2.0) ** 2.0, 0.0])).all()
+# ---------------------------------------------------------------------------
+# chi_squared_map_complex_from
+# ---------------------------------------------------------------------------
 
 
-def test__chi_squared_map_complex_from():
+def test__chi_squared_map_complex_from__complex_data__correct_real_and_imaginary_chi_squared():
     data = np.array([10.0 + 10.0j, 10.0 + 10.0j])
     noise_map = np.array([2.0 + 2.0j, 2.0 + 2.0j])
     model_data = np.array([9.0 + 12.0j, 9.0 + 12.0j])
 
     residual_map = aa.util.fit.residual_map_from(data=data, model_data=model_data)
-
     chi_squared_map = aa.util.fit.chi_squared_map_complex_from(
         residual_map=residual_map, noise_map=noise_map
     )
@@ -163,8 +193,13 @@ def test__chi_squared_map_complex_from():
     assert (chi_squared_map == np.array([0.25 + 1.0j, 0.25 + 1.0j])).all()
 
 
-def test__chi_squared_with_noise_covariance_from():
-    resdiual_map = aa.Array2D.no_mask([[1.0, 1.0], [2.0, 2.0]], pixel_scales=1.0)
+# ---------------------------------------------------------------------------
+# chi_squared_with_noise_covariance_from
+# ---------------------------------------------------------------------------
+
+
+def test__chi_squared_with_noise_covariance_from__known_residual_and_inverse__correct_chi_squared():
+    residual_map = aa.Array2D.no_mask([[1.0, 1.0], [2.0, 2.0]], pixel_scales=1.0)
 
     noise_covariance_matrix_inv = np.array(
         [
@@ -176,14 +211,19 @@ def test__chi_squared_with_noise_covariance_from():
     )
 
     chi_squared = aa.util.fit.chi_squared_with_noise_covariance_from(
-        residual_map=resdiual_map,
+        residual_map=residual_map,
         noise_covariance_matrix_inv=noise_covariance_matrix_inv,
     )
 
     assert chi_squared == 43
 
 
-def test__chi_squared_with_mask_fast_from():
+# ---------------------------------------------------------------------------
+# chi_squared_with_mask_fast_from
+# ---------------------------------------------------------------------------
+
+
+def test__chi_squared_with_mask_fast_from__1d_data_with_mask__equals_standard_computation():
     data = np.array([10.0, 10.0, 10.0, 10.0])
     mask = np.array([True, False, False, True])
     noise_map = np.array([1.0, 2.0, 3.0, 4.0])
@@ -192,15 +232,12 @@ def test__chi_squared_with_mask_fast_from():
     residual_map = aa.util.fit.residual_map_with_mask_from(
         data=data, mask=mask, model_data=model_data
     )
-
     chi_squared_map = aa.util.fit.chi_squared_map_with_mask_from(
         residual_map=residual_map, mask=mask, noise_map=noise_map
     )
-
     chi_squared = aa.util.fit.chi_squared_with_mask_from(
         mask=mask, chi_squared_map=chi_squared_map
     )
-
     chi_squared_fast = aa.util.fit.chi_squared_with_mask_fast_from(
         data=data,
         noise_map=noise_map,
@@ -210,6 +247,8 @@ def test__chi_squared_with_mask_fast_from():
 
     assert chi_squared == pytest.approx(chi_squared_fast, 1.0e-4)
 
+
+def test__chi_squared_with_mask_fast_from__2d_data_with_mask__equals_standard_computation():
     data = np.array([[10.0, 10.0], [10.0, 10.0]])
     mask = np.array([[True, False], [False, True]])
     noise_map = np.array([[1.0, 2.0], [3.0, 4.0]])
@@ -218,15 +257,12 @@ def test__chi_squared_with_mask_fast_from():
     residual_map = aa.util.fit.residual_map_with_mask_from(
         data=data, mask=mask, model_data=model_data
     )
-
     chi_squared_map = aa.util.fit.chi_squared_map_with_mask_from(
         residual_map=residual_map, mask=mask, noise_map=noise_map
     )
-
     chi_squared = aa.util.fit.chi_squared_with_mask_from(
         mask=mask, chi_squared_map=chi_squared_map
     )
-
     chi_squared_fast = aa.util.fit.chi_squared_with_mask_fast_from(
         data=data,
         noise_map=noise_map,
@@ -237,90 +273,78 @@ def test__chi_squared_with_mask_fast_from():
     assert chi_squared == pytest.approx(chi_squared_fast, 1.0e-4)
 
 
-def test__log_likelihood_from():
+# ---------------------------------------------------------------------------
+# log_likelihood_from
+# ---------------------------------------------------------------------------
+
+
+def test__log_likelihood_from__identical_data_and_model__correct_value_from_noise_normalization():
     data = np.array([10.0, 10.0, 10.0, 10.0])
     noise_map = np.array([2.0, 2.0, 2.0, 2.0])
     model_data = np.array([10.0, 10.0, 10.0, 10.0])
 
     residual_map = aa.util.fit.residual_map_from(data=data, model_data=model_data)
-
     chi_squared_map = aa.util.fit.chi_squared_map_from(
         residual_map=residual_map, noise_map=noise_map
     )
-
     chi_squared = aa.util.fit.chi_squared_from(chi_squared_map=chi_squared_map)
-
     noise_normalization = aa.util.fit.noise_normalization_from(noise_map=noise_map)
-
     log_likelihood = aa.util.fit.log_likelihood_from(
         chi_squared=chi_squared, noise_normalization=noise_normalization
     )
 
-    chi_squared = 0.0
-    noise_normalization = (
-        np.log(2.0 * np.pi * (2.0**2.0))
-        + np.log(2.0 * np.pi * (2.0**2.0))
-        + np.log(2.0 * np.pi * (2.0**2.0))
-        + np.log(2.0 * np.pi * (2.0**2.0))
-    )
+    expected_chi_squared = 0.0
+    expected_noise_normalization = 4 * np.log(2.0 * np.pi * (2.0**2.0))
 
     assert log_likelihood == pytest.approx(
-        -0.5 * (chi_squared + noise_normalization), 1.0e-4
+        -0.5 * (expected_chi_squared + expected_noise_normalization), 1.0e-4
     )
 
+
+def test__log_likelihood_from__different_data_and_model_uniform_noise__correct_value():
+    data = np.array([10.0, 10.0, 10.0, 10.0])
+    noise_map = np.array([2.0, 2.0, 2.0, 2.0])
     model_data = np.array([11.0, 10.0, 9.0, 8.0])
 
     residual_map = aa.util.fit.residual_map_from(data=data, model_data=model_data)
-
     chi_squared_map = aa.util.fit.chi_squared_map_from(
         residual_map=residual_map, noise_map=noise_map
     )
-
     chi_squared = aa.util.fit.chi_squared_from(chi_squared_map=chi_squared_map)
-
     noise_normalization = aa.util.fit.noise_normalization_from(noise_map=noise_map)
-
     log_likelihood = aa.util.fit.log_likelihood_from(
         chi_squared=chi_squared, noise_normalization=noise_normalization
     )
 
     # chi squared = 0.25, 0, 0.25, 1.0
-    # log_likelihood = -0.5*(0.25+0+0.25+1.0)
-
-    chi_squared = (
+    expected_chi_squared = (
         ((1.0 / 2.0) ** 2.0) + 0.0 + ((1.0 / 2.0) ** 2.0) + ((2.0 / 2.0) ** 2.0)
     )
-    noise_normalization = (
-        np.log(2.0 * np.pi * (2.0**2.0))
-        + np.log(2.0 * np.pi * (2.0**2.0))
-        + np.log(2.0 * np.pi * (2.0**2.0))
-        + np.log(2.0 * np.pi * (2.0**2.0))
-    )
+    expected_noise_normalization = 4 * np.log(2.0 * np.pi * (2.0**2.0))
 
     assert log_likelihood == pytest.approx(
-        -0.5 * (chi_squared + noise_normalization), 1.0e-4
+        -0.5 * (expected_chi_squared + expected_noise_normalization), 1.0e-4
     )
 
+
+def test__log_likelihood_from__different_data_and_model_varied_noise__correct_value():
+    data = np.array([10.0, 10.0, 10.0, 10.0])
     noise_map = np.array([1.0, 2.0, 3.0, 4.0])
+    model_data = np.array([11.0, 10.0, 9.0, 8.0])
 
     residual_map = aa.util.fit.residual_map_from(data=data, model_data=model_data)
-
     chi_squared_map = aa.util.fit.chi_squared_map_from(
         residual_map=residual_map, noise_map=noise_map
     )
-
     chi_squared = aa.util.fit.chi_squared_from(chi_squared_map=chi_squared_map)
-
     noise_normalization = aa.util.fit.noise_normalization_from(noise_map=noise_map)
-
     log_likelihood = aa.util.fit.log_likelihood_from(
         chi_squared=chi_squared, noise_normalization=noise_normalization
     )
 
     # chi squared = (1.0/1.0)**2, (0.0), (-1.0/3.0)**2.0, (2.0/4.0)**2.0
-
-    chi_squared = 1.0 + (1.0 / (3.0**2.0)) + 0.25
-    noise_normalization = (
+    expected_chi_squared = 1.0 + (1.0 / (3.0**2.0)) + 0.25
+    expected_noise_normalization = (
         np.log(2 * np.pi * (1.0**2.0))
         + np.log(2 * np.pi * (2.0**2.0))
         + np.log(2 * np.pi * (3.0**2.0))
@@ -328,11 +352,16 @@ def test__log_likelihood_from():
     )
 
     assert log_likelihood == pytest.approx(
-        -0.5 * (chi_squared + noise_normalization), 1e-4
+        -0.5 * (expected_chi_squared + expected_noise_normalization), 1e-4
     )
 
 
-def test__log_likelihood_from__with_mask():
+# ---------------------------------------------------------------------------
+# log_likelihood_from with mask
+# ---------------------------------------------------------------------------
+
+
+def test__log_likelihood_from__with_1d_mask__excludes_masked_pixels_from_chi_squared_and_noise_normalization():
     data = np.array([10.0, 10.0, 10.0, 10.0])
     mask = np.array([True, False, False, True])
     noise_map = np.array([1.0, 2.0, 3.0, 4.0])
@@ -341,34 +370,31 @@ def test__log_likelihood_from__with_mask():
     residual_map = aa.util.fit.residual_map_with_mask_from(
         data=data, mask=mask, model_data=model_data
     )
-
     chi_squared_map = aa.util.fit.chi_squared_map_with_mask_from(
         residual_map=residual_map, mask=mask, noise_map=noise_map
     )
-
     chi_squared = aa.util.fit.chi_squared_with_mask_from(
         mask=mask, chi_squared_map=chi_squared_map
     )
-
     noise_normalization = aa.util.fit.noise_normalization_with_mask_from(
         mask=mask, noise_map=noise_map
     )
-
     log_likelihood = aa.util.fit.log_likelihood_from(
         chi_squared=chi_squared, noise_normalization=noise_normalization
     )
 
-    # chi squared = 0, 0.25, (0.25 and 1.0 are masked)
-
-    chi_squared = 0.0 + (1.0 / 3.0) ** 2.0
-    noise_normalization = np.log(2 * np.pi * (2.0**2.0)) + np.log(
+    # chi squared = 0, (1/3)**2  (pixels at idx 0 and 3 are masked)
+    expected_chi_squared = 0.0 + (1.0 / 3.0) ** 2.0
+    expected_noise_normalization = np.log(2 * np.pi * (2.0**2.0)) + np.log(
         2 * np.pi * (3.0**2.0)
     )
 
     assert log_likelihood == pytest.approx(
-        -0.5 * (chi_squared + noise_normalization), 1e-4
+        -0.5 * (expected_chi_squared + expected_noise_normalization), 1e-4
     )
 
+
+def test__log_likelihood_from__with_2d_mask__excludes_masked_pixels_from_chi_squared_and_noise_normalization():
     data = np.array([[10.0, 10.0], [10.0, 10.0]])
     mask = np.array([[True, False], [False, True]])
     noise_map = np.array([[1.0, 2.0], [3.0, 4.0]])
@@ -377,70 +403,69 @@ def test__log_likelihood_from__with_mask():
     residual_map = aa.util.fit.residual_map_with_mask_from(
         data=data, mask=mask, model_data=model_data
     )
-
     chi_squared_map = aa.util.fit.chi_squared_map_with_mask_from(
         residual_map=residual_map, mask=mask, noise_map=noise_map
     )
-
     chi_squared = aa.util.fit.chi_squared_with_mask_from(
         mask=mask, chi_squared_map=chi_squared_map
     )
-
     noise_normalization = aa.util.fit.noise_normalization_with_mask_from(
         mask=mask, noise_map=noise_map
     )
-
     log_likelihood = aa.util.fit.log_likelihood_from(
         chi_squared=chi_squared, noise_normalization=noise_normalization
     )
 
-    # chi squared = 0, 0.25, (0.25 and 1.0 are masked)
-
-    chi_squared = 0.0 + (1.0 / 3.0) ** 2.0
-    noise_normalization = np.log(2 * np.pi * (2.0**2.0)) + np.log(
+    # chi squared = 0, (1/3)**2  (corner pixels are masked)
+    expected_chi_squared = 0.0 + (1.0 / 3.0) ** 2.0
+    expected_noise_normalization = np.log(2 * np.pi * (2.0**2.0)) + np.log(
         2 * np.pi * (3.0**2.0)
     )
 
     assert log_likelihood == pytest.approx(
-        -0.5 * (chi_squared + noise_normalization), 1e-4
+        -0.5 * (expected_chi_squared + expected_noise_normalization), 1e-4
     )
 
 
-def test__log_likelihood_from__complex_data():
+# ---------------------------------------------------------------------------
+# log_likelihood_from with complex data
+# ---------------------------------------------------------------------------
+
 
+def test__log_likelihood_from__complex_data__chi_squared_sums_real_and_imaginary_components():
     data = np.array([10.0 + 10.0j, 10.0 + 10.0j])
     noise_map = np.array([2.0 + 1.0j, 2.0 + 1.0j])
     model_data = np.array([9.0 + 12.0j, 9.0 + 12.0j])
 
     residual_map = aa.util.fit.residual_map_from(data=data, model_data=model_data)
-
     chi_squared_map = aa.util.fit.chi_squared_map_complex_from(
         residual_map=residual_map, noise_map=noise_map
     )
-
     chi_squared = aa.util.fit.chi_squared_complex_from(chi_squared_map=chi_squared_map)
-
     noise_normalization = aa.util.fit.noise_normalization_complex_from(
         noise_map=noise_map
     )
-
     log_likelihood = aa.util.fit.log_likelihood_from(
         chi_squared=chi_squared, noise_normalization=noise_normalization
     )
 
     # chi squared = 0.25 and 4.0
-
-    chi_squared = 4.25
-    noise_normalization = np.log(2 * np.pi * (2.0**2.0)) + np.log(
+    expected_chi_squared = 4.25
+    expected_noise_normalization = np.log(2 * np.pi * (2.0**2.0)) + np.log(
         2 * np.pi * (1.0**2.0)
     )
 
     assert log_likelihood == pytest.approx(
-        -0.5 * 2.0 * (chi_squared + noise_normalization), 1e-4
+        -0.5 * 2.0 * (expected_chi_squared + expected_noise_normalization), 1e-4
     )
 
 
-def test__log_evidence_from():
+# ---------------------------------------------------------------------------
+# log_likelihood_with_regularization_from / log_evidence_from
+# ---------------------------------------------------------------------------
+
+
+def test__log_likelihood_with_regularization_from__adds_regularization_to_chi_squared_and_noise():
     likelihood_with_regularization_terms = (
         aa.util.fit.log_likelihood_with_regularization_from(
             chi_squared=3.0, regularization_term=6.0, noise_normalization=2.0
@@ -449,6 +474,8 @@ def test__log_evidence_from():
 
     assert likelihood_with_regularization_terms == -0.5 * (3.0 + 6.0 + 2.0)
 
+
+def test__log_evidence_from__includes_curvature_and_regularization_determinant_terms():
     evidences = aa.util.fit.log_evidence_from(
         chi_squared=3.0,
         regularization_term=6.0,
@@ -460,22 +487,28 @@ def test__log_evidence_from():
     assert evidences == -0.5 * (3.0 + 6.0 + 9.0 - 10.0 + 30.0)
 
 
-def test__residual_flux_fraction_map_from():
+# ---------------------------------------------------------------------------
+# residual_flux_fraction_map_from
+# ---------------------------------------------------------------------------
+
+
+def test__residual_flux_fraction_map_from__identical_data_and_model__all_zero():
     data = np.array([10.0, 10.0, 10.0, 10.0])
     model_data = np.array([10.0, 10.0, 10.0, 10.0])
 
     residual_map = aa.util.fit.residual_map_from(data=data, model_data=model_data)
-
     residual_flux_fraction_map = aa.util.fit.residual_flux_fraction_map_from(
         residual_map=residual_map, data=data
     )
 
     assert (residual_flux_fraction_map == np.array([0.0, 0.0, 0.0, 0.0])).all()
 
+
+def test__residual_flux_fraction_map_from__different_data_and_model__correct_fractional_residuals():
+    data = np.array([10.0, 10.0, 10.0, 10.0])
     model_data = np.array([11.0, 10.0, 9.0, 8.0])
 
     residual_map = aa.util.fit.residual_map_from(data=data, model_data=model_data)
-
     residual_flux_fraction_map = aa.util.fit.residual_flux_fraction_map_from(
         residual_map=residual_map, data=data
     )
@@ -483,7 +516,12 @@ def test__residual_flux_fraction_map_from():
     assert (residual_flux_fraction_map == np.array([-0.1, 0.0, 0.1, 0.2])).all()
 
 
-def test__residual_flux_fraction_map_with_mask_from():
+# ---------------------------------------------------------------------------
+# residual_flux_fraction_map_with_mask_from
+# ---------------------------------------------------------------------------
+
+
+def test__residual_flux_fraction_map_with_mask_from__masked_data__zero_at_masked_locations():
     data = np.array([10.0, 10.0, 10.0, 10.0])
     mask = np.array([True, False, False, True])
     model_data = np.array([11.0, 10.0, 9.0, 8.0])
@@ -491,19 +529,21 @@ def test__residual_flux_fraction_map_with_mask_from():
     residual_map = aa.util.fit.residual_map_with_mask_from(
         data=data, mask=mask, model_data=model_data
     )
-
     residual_flux_fraction_map = aa.util.fit.residual_flux_fraction_map_with_mask_from(
         residual_map=residual_map, mask=mask, data=data
     )
 
     assert (residual_flux_fraction_map == np.array([0.0, 0.0, 0.1, 0.0])).all()
 
+
+def test__residual_flux_fraction_map_with_mask_from__different_model__correct_unmasked_fractions():
+    data = np.array([10.0, 10.0, 10.0, 10.0])
+    mask = np.array([True, False, False, True])
     model_data = np.array([11.0, 9.0, 8.0, 8.0])
 
     residual_map = aa.util.fit.residual_map_with_mask_from(
         data=data, mask=mask, model_data=model_data
     )
-
     residual_flux_fraction_map = aa.util.fit.residual_flux_fraction_map_with_mask_from(
         residual_map=residual_map, mask=mask, data=data
     )