tests working

meteor/diffmaps.py

@@ -9,26 +9,11 @@
 
 from .rsmap import Map, assert_is_map
 from .settings import DEFAULT_KPARAMS_TO_SCAN
-from .utils import filter_common_indices
+from .utils import assert_isomorphous, filter_common_indices
 from .validate import ScalarMaximizer, map_negentropy
 
 
-def set_common_crystallographic_metadata(map1: Map, map2: Map, *, output: Map) -> None:
-    if hasattr(map1, "cell"):
-        if hasattr(map2, "cell") and (map1.cell != map2.cell):
-            msg = f"`map1.cell` {map1.cell} != `map2.cell` {map2.cell}"
-            raise AttributeError(msg)
-        output.cell = map1.cell
-
-    if hasattr(map1, "spacegroup"):
-        if hasattr(map2, "spacegroup") and (map1.spacegroup != map2.spacegroup):
-            msg = f"`map1.spacegroup` {map1.spacegroup} != "
-            msg += f"`map2.spacegroup` {map2.spacegroup}"
-            raise AttributeError(msg)
-        output.spacegroup = map1.spacegroup
-
-
-def compute_difference_map(derivative: Map, native: Map) -> Map:
+def compute_difference_map(derivative: Map, native: Map, *, check_isomorphous: bool = True) -> Map:
     """
     Computes amplitude and phase differences between native and derivative structure factor sets.
 
@@ -43,23 +28,30 @@ def compute_difference_map(derivative: Map, native: Map) -> Map:
     ----------
     derivative: Map
         the derivative amplitudes, phases, uncertainties
+
     native: Map
         the native amplitudes, phases, uncertainties
 
+    check_isomorphous: bool
+        perform a check to ensure the two datasets are isomorphous; recommended. Default: True.
+
     Returns
     -------
     diffmap: Map
         map corresponding to the complex difference (derivative - native)
     """
     assert_is_map(derivative, require_uncertainties=False)
     assert_is_map(native, require_uncertainties=False)
+    if check_isomorphous:
+        assert_isomorphous(derivative=derivative, native=native)
 
     derivative, native = filter_common_indices(derivative, native)
 
     delta_complex = derivative.to_structurefactor() - native.to_structurefactor()
     delta = Map.from_structurefactor(delta_complex)
 
-    set_common_crystallographic_metadata(derivative, native, output=delta)
+    delta.cell = native.cell
+    delta.spacegroup = native.spacegroup
 
     if derivative.has_uncertainties and native.has_uncertainties:
         prop_uncertainties = np.sqrt(derivative.uncertainties**2 + native.uncertainties**2)
@@ -76,6 +68,7 @@ def compute_kweights(difference_map: Map, *, k_parameter: float) -> rs.DataSerie
     ----------
     difference_map: Map
         A map of structure factor differences (DeltaF).
+
     k_parameter: float
         A scaling factor applied to the squared `df` values in the weight calculation.
 
@@ -95,7 +88,9 @@ def compute_kweights(difference_map: Map, *, k_parameter: float) -> rs.DataSerie
     return 1.0 / inverse_weights
 
 
-def compute_kweighted_difference_map(derivative: Map, native: Map, *, k_parameter: float) -> Map:
+def compute_kweighted_difference_map(
+    derivative: Map, native: Map, *, k_parameter: float, check_isomorphous: bool = True
+) -> Map:
     """
     Compute k-weighted derivative - native structure factor map.
 
@@ -108,9 +103,13 @@ def compute_kweighted_difference_map(derivative: Map, native: Map, *, k_paramete
     ----------
     derivative: Map
         the derivative amplitudes, phases, uncertainties
+
     native: Map
         the native amplitudes, phases, uncertainties
 
+    check_isomorphous: bool
+        perform a check to ensure the two datasets are isomorphous; recommended. Default: True.
+
     Returns
     -------
     diffmap: Map
@@ -119,8 +118,10 @@ def compute_kweighted_difference_map(derivative: Map, native: Map, *, k_paramete
     # require uncertainties at the beginning
     assert_is_map(derivative, require_uncertainties=True)
     assert_is_map(native, require_uncertainties=True)
+    if check_isomorphous:
+        assert_isomorphous(derivative=derivative, native=native)
 
-    difference_map = compute_difference_map(derivative, native)
+    difference_map = compute_difference_map(derivative, native, check_isomorphous=check_isomorphous)
     weights = compute_kweights(difference_map, k_parameter=k_parameter)
 
     difference_map.amplitudes *= weights
@@ -134,6 +135,7 @@ def max_negentropy_kweighted_difference_map(
     native: Map,
     *,
     k_parameter_values_to_scan: np.ndarray | Sequence[float] = DEFAULT_KPARAMS_TO_SCAN,
+    check_isomorphous: bool = True,
 ) -> rs.DataSet:
     """
     Compute k-weighted differences between native and derivative amplitudes and phases.
@@ -146,11 +148,16 @@ def max_negentropy_kweighted_difference_map(
     ----------
     derivative: Map
         the derivative amplitudes, phases, uncertainties
+
     native: Map
         the native amplitudes, phases, uncertainties
+
     k_parameter_values_to_scan : np.ndarray | Sequence[float]
         The values to scan to optimize the k-weighting parameter, by default is 0.00, 0.01 ... 1.00
 
+    check_isomorphous: bool
+        perform a check to ensure the two datasets are isomorphous; recommended. Default: True.
+
     Returns
     -------
     kweighted_dataset: rs.DataSet
@@ -161,12 +168,15 @@ def max_negentropy_kweighted_difference_map(
     """
     assert_is_map(derivative, require_uncertainties=True)
     assert_is_map(native, require_uncertainties=True)
+    if check_isomorphous:
+        assert_isomorphous(derivative=derivative, native=native)
 
     def negentropy_objective(k_parameter: float) -> float:
         kweighted_map = compute_kweighted_difference_map(
             derivative,
             native,
             k_parameter=k_parameter,
+            check_isomorphous=check_isomorphous,
         )
         return map_negentropy(kweighted_map)
 
@@ -178,6 +188,7 @@ def negentropy_objective(k_parameter: float) -> float:
         derivative,
         native,
         k_parameter=opt_k_parameter,
+        check_isomorphous=check_isomorphous,
     )
 
     return kweighted_dataset, opt_k_parameter

meteor/iterative.py

@@ -14,7 +14,7 @@
     ITERATIVE_TV_MAX_ITERATIONS,
 )
 from .tv import TvDenoiseResult, tv_denoise_difference_map
-from .utils import CellType, SpacegroupType, average_phase_diff_in_degrees
+from .utils import CellType, SpacegroupType, assert_isomorphous, average_phase_diff_in_degrees
 
 log = structlog.get_logger()
 
@@ -214,6 +214,7 @@ def __call__(
         *,
         derivative: Map,
         native: Map,
+        check_isomorphous: bool = True,
     ) -> tuple[Map, pd.DataFrame]:
         """
         Denoise by estimating new, low-TV phases for the `derivative` dataset.
@@ -226,6 +227,9 @@ def __call__(
         native: Map
             the native amplitudes, phases
 
+        check_isomorphous: bool
+            perform a check to ensure the two datasets are isomorphous; recommended. Default: True.
+
         Returns
         -------
         updated_derivative: Map
@@ -236,10 +240,8 @@ def __call__(
             the tv_weight used, the negentropy (after the TV step), and the average phase change in
             degrees.
         """
-        # TODO: do we need this?
-        # initial_derivative, native = filter_common_indices(initial_derivative, native)
-
-        # TODO: check isomorphous
+        if check_isomorphous:
+            assert_isomorphous(derivative=derivative, native=native)
 
         it_tv_complex_derivative, metadata = self._iteratively_denoise_sf_amplitudes(
             native=native.to_structurefactor(),

meteor/tv.py

@@ -126,10 +126,12 @@ def tv_denoise_difference_map(
     difference_map : Map
         The input dataset containing the difference map coefficients (amplitude and phase)
         that will be used to compute the difference map.
+
     full_output : bool, optional
         If `True`, the function returns both the denoised map coefficients and a `TvDenoiseResult`
          object containing the optimal weight and the associated negentropy. If `False`, only
          the denoised map coefficients are returned. Default is `False`.
+
     weights_to_scan : Sequence[float] | None, optional
         A sequence of weight values to explicitly scan for determining the optimal value. If
         `None`, the function uses the golden-section search method to determine the optimal

meteor/utils.py

@@ -19,6 +19,17 @@
 class ShapeMismatchError(Exception): ...
 
 
+class NotIsomorphousError(RuntimeError): ...
+
+
+def assert_isomorphous(*, derivative: rs.DataSet, native: rs.DataSet) -> None:
+    if not native.is_isomorphous(derivative):
+        msg = "`derivative` and `native` datasets are not similar enough; "
+        msg += f"they have cell/spacegroup: {derivative.cell}/{native.cell} and "
+        msg += f"{derivative.spacegroup}/{native.spacegroup} respectively"
+        raise NotIsomorphousError(msg)
+
+
 def filter_common_indices(df1: DataSet, df2: DataSet) -> tuple[DataSet, DataSet]:
     common_indices = df1.index.intersection(df2.index)
     df1_common = df1.loc[common_indices].copy()

test/unit/test_diffmaps.py

@@ -11,9 +11,9 @@
     compute_kweighted_difference_map,
     compute_kweights,
     max_negentropy_kweighted_difference_map,
-    set_common_crystallographic_metadata,
 )
 from meteor.rsmap import Map
+from meteor.utils import NotIsomorphousError
 from meteor.validate import map_negentropy
 
 
@@ -40,80 +40,66 @@ def dummy_native() -> Map:
     return Map(native, index=index).infer_mtz_dtypes()
 
 
-def test_set_common_crystallographic_metadata(dummy_native: Map) -> None:
-    dummy_native.cell = (10.0, 10.0, 10.0, 90.0, 90.0, 90.0)
-    dummy_native.spacegroup = 1
-    map1 = dummy_native
-    map2 = dummy_native.copy()
-    output = dummy_native.copy()
-
-    # ensure things are copied when missing initially
-    del output._cell
-    del output._spacegroup
-    assert not hasattr(output, "cell")
-    assert not hasattr(output, "spacegroup")
-
-    set_common_crystallographic_metadata(map1, map2, output=output)
-
-    assert output.cell == map1.cell
-    assert output.spacegroup == map2.spacegroup
-
-    map2.spacegroup = 19
-    with pytest.raises(AttributeError):
-        set_common_crystallographic_metadata(map1, map2, output=output)
-
-    map2.spacegroup = 1
-    set_common_crystallographic_metadata(map1, map2, output=output)
-
-    dummy_native.cell = (15.0, 15.0, 15.0, 90.0, 90.0, 90.0)
-    with pytest.raises(AttributeError):
-        set_common_crystallographic_metadata(map1, map2, output=output)
-
-
 def test_compute_difference_map_vs_analytical(dummy_derivative: Map, dummy_native: Map) -> None:
     # Manually calculated expected amplitude and phase differences
     expected_amplitudes = np.array([3.0, 5.0])
     expected_phases = np.array([-180.0, 0.0])
     assert isinstance(dummy_native, Map)
     assert isinstance(dummy_derivative, Map)
 
-    result = compute_difference_map(dummy_derivative, dummy_native)
+    result = compute_difference_map(dummy_derivative, dummy_native, check_isomorphous=False)
     assert_almost_equal(result.amplitudes, expected_amplitudes, decimal=4)
     assert_almost_equal(result.phases, expected_phases, decimal=4)
 
 
 @pytest.mark.parametrize(
     "diffmap_fxn",
+    # lambdas to make the call signatures for these functions match `compute_difference_map`
     [
-        compute_difference_map,
-        # lambdas to make the call signatures for these functions match `compute_difference_map`
-        lambda d, n: compute_kweighted_difference_map(d, n, k_parameter=0.5),
-        lambda d, n: max_negentropy_kweighted_difference_map(d, n)[0],
+        lambda d, n, check: compute_difference_map(d, n, check_isomorphous=check),
+        lambda d, n, check: compute_kweighted_difference_map(
+            d, n, k_parameter=0.5, check_isomorphous=check
+        ),
+        lambda d, n, check: max_negentropy_kweighted_difference_map(d, n, check_isomorphous=check)[
+            0
+        ],
     ],
 )
+@pytest.mark.parametrize("check_isomorphous", [True, False])
 def test_cell_spacegroup_propogation(
     diffmap_fxn: Callable,
     dummy_derivative: Map,
     dummy_native: Map,
+    check_isomorphous: bool,
 ) -> None:
+    # these should all cast to gemmi objects
     dummy_derivative.cell = (10.0, 10.0, 10.0, 90.0, 90.0, 90.0)
-    dummy_derivative.spacegroup = 1  # will cast to gemmi.SpaceGroup
-    dummy_native.cell = (10.0, 10.0, 10.0, 90.0, 90.0, 90.0)
+    dummy_derivative.spacegroup = 1
+    dummy_native.cell = (10.01, 10.01, 10.01, 90.01, 90.01, 90.01)
     dummy_native.spacegroup = 1
 
-    result = diffmap_fxn(dummy_derivative, dummy_native)
+    # ensure the native cell is propogated
+    result = diffmap_fxn(dummy_derivative, dummy_native, check_isomorphous)
     assert result.cell == dummy_native.cell
     assert result.spacegroup == dummy_native.spacegroup
 
+    # check we raise or dont with a spacegroup mismatch
     dummy_native.spacegroup = 19
-    with pytest.raises(AttributeError):
-        result = diffmap_fxn(dummy_derivative, dummy_native)
+    if check_isomorphous:
+        with pytest.raises(NotIsomorphousError):
+            _ = diffmap_fxn(dummy_derivative, dummy_native, check_isomorphous)
+    else:
+        _ = diffmap_fxn(dummy_derivative, dummy_native, check_isomorphous)
 
+    # check we raise or dont with a cell mismatch
     dummy_native.spacegroup = 1
-    _ = compute_difference_map(dummy_derivative, dummy_native)
-    dummy_native.cell = (15.0, 15.0, 15.0, 90.0, 90.0, 90.0)
-    with pytest.raises(AttributeError):
-        result = diffmap_fxn(dummy_derivative, dummy_native)
+    _ = diffmap_fxn(dummy_derivative, dummy_native, check_isomorphous)
+    dummy_native.cell = (20.0, 10.0, 10.0, 90.0, 90.0, 90.0)
+    if check_isomorphous:
+        with pytest.raises(NotIsomorphousError):
+            _ = diffmap_fxn(dummy_derivative, dummy_native, check_isomorphous)
+    else:
+        _ = diffmap_fxn(dummy_derivative, dummy_native, check_isomorphous)
 
 
 def test_compute_kweights_vs_analytical() -> None:
@@ -133,7 +119,9 @@ def test_compute_kweighted_difference_map_vs_analytical(
     dummy_derivative: Map,
     dummy_native: Map,
 ) -> None:
-    kwt_diffmap = compute_kweighted_difference_map(dummy_derivative, dummy_native, k_parameter=0.5)
+    kwt_diffmap = compute_kweighted_difference_map(
+        dummy_derivative, dummy_native, k_parameter=0.5, check_isomorphous=False
+    )
     expected_weighted_amplitudes = np.array([1.3247, 1.8280])  # calculated by hand
     expected_weighted_uncertainties = np.array([0.3122, 0.2585])
     assert_almost_equal(kwt_diffmap.amplitudes, expected_weighted_amplitudes, decimal=4)