TV weights equal or less than zero

meteor/rsmap.py

@@ -12,6 +12,7 @@
 from .utils import (
     canonicalize_amplitudes,
     complex_array_to_rs_dataseries,
+    numpy_array_to_map,
 )
 
 
@@ -249,11 +250,14 @@ def set_uncertainties(self, values: rs.DataSeries, column_name: str = "SIGF") ->
         else:
             # otherwise, create a new column
             self._uncertainty_column = column_name
-            position = len(self.columns)
-            if position != 2:  # noqa: PLR2004, should be 2: just amplitudes & phases
+            number_of_columns = len(self.columns)
+            number_of_columns_with_just_amplitudes_and_phases = 2
+            if number_of_columns != number_of_columns_with_just_amplitudes_and_phases:
                 msg = "Misconfigured columns"
                 raise RuntimeError(msg)
-            super().insert(position, self._uncertainty_column, values, allow_duplicates=False)
+            super().insert(
+                number_of_columns, self._uncertainty_column, values, allow_duplicates=False
+            )
 
     @property
     def complex_amplitudes(self) -> np.ndarray:
@@ -308,6 +312,48 @@ def from_gemmi(
             uncertainty_column=uncertainty_column,
         )
 
+    @classmethod
+    def from_3d_numpy_map(
+        cls, map_grid: np.ndarray, *, spacegroup: Any, cell: Any, high_resolution_limit: float
+    ) -> Map:
+        """
+        Create a `Map` from a 3d grid of voxel values stored in a numpy array.
+
+        Parameters
+        ----------
+        map_grid: np.ndarray
+            The array, laid out in Gemmi format
+        spacegroup: Any
+            Specifies which spacegroup, can be an int, gemmi.SpaceGroup, ...
+        cell
+            Specifies cell, can be a tuple, gemmi.Cell, ...
+        high_resolution_limit: float
+            The resolution of the map, irregardless of the sampling; we need this to infer the map
+            sampling
+
+        Returns
+        -------
+        map: Map
+            The map coefficients
+
+        See Also
+        --------
+        For information about Gemmi data layout: https://gemmi.readthedocs.io/en/latest/grid.html
+        """
+        number_of_dimensions_in_universe = 3
+        if len(map_grid.shape) != number_of_dimensions_in_universe:
+            msg = "`map_grid` should be a 3D array representing a realspace map"
+            raise ValueError(msg)
+        ccp4 = numpy_array_to_map(
+            map_grid,
+            spacegroup=spacegroup,
+            cell=cell,
+        )
+        return cls.from_ccp4_map(
+            ccp4_map=ccp4,
+            high_resolution_limit=high_resolution_limit,
+        )
+
     def to_ccp4_map(self, *, map_sampling: int) -> gemmi.Ccp4Map:
         map_coefficients_gemmi_format = self.to_gemmi()
         ccp4_map = gemmi.Ccp4Map()

meteor/tv.py

@@ -13,9 +13,6 @@
     TV_MAX_NUM_ITER,
     TV_STOP_TOLERANCE,
 )
-from .utils import (
-    numpy_array_to_map,
-)
 from .validate import ScalarMaximizer, negentropy
 
 
@@ -29,6 +26,8 @@ class TvDenoiseResult:
 
 def _tv_denoise_array(*, map_as_array: np.ndarray, weight: float) -> np.ndarray:
     """Closure convienence function to generate more readable code."""
+    if weight <= 0.0:
+        return map_as_array
     return denoise_tv_chambolle(
         map_as_array,
         weight=weight,
@@ -131,17 +130,17 @@ def negentropy_objective(tv_lambda: float):
         map_as_array=realspace_map_array,
         weight=maximizer.argument_optimum,
     )
-    final_realspace_map_as_ccp4 = numpy_array_to_map(
+    final_map = Map.from_3d_numpy_map(
         final_realspace_map_as_array,
         spacegroup=difference_map.spacegroup,
         cell=difference_map.cell,
-    )
-
-    final_map = Map.from_ccp4_map(
-        ccp4_map=final_realspace_map_as_ccp4,
         high_resolution_limit=difference_map.resolution_limits[1],
     )
 
+    # propogate uncertainties
+    if difference_map.has_uncertainties:
+        final_map.set_uncertainties(difference_map.uncertainties)
+
     # sometimes `from_ccp4_map` adds reflections -- systematic absences or
     # reflections just beyond the resolution limt; remove those
     extra_indices = final_map.index.difference(difference_map.index)

test/unit/test_rsmap.py

@@ -274,6 +274,28 @@ def test_from_ccp4_map(ccp4_map: gemmi.Ccp4Map) -> None:
     assert len(rsmap) > 0
 
 
+def test_from_numpy(noise_free_map: Map) -> None:
+    _, resolution = noise_free_map.resolution_limits
+    array = np.array(noise_free_map.to_ccp4_map(map_sampling=3).grid)
+    new_map = Map.from_3d_numpy_map(
+        array,
+        spacegroup=noise_free_map.spacegroup,
+        cell=noise_free_map.cell,
+        high_resolution_limit=resolution,
+    )
+
+    pd.testing.assert_series_equal(new_map.amplitudes, noise_free_map.amplitudes, atol=1e-4)
+    assert_phases_allclose(new_map.phases, noise_free_map.phases)
+
+    with pytest.raises(ValueError, match="`map_grid`"):
+        Map.from_3d_numpy_map(
+            array[0, :, :],  # only has 2 dimensions
+            spacegroup=noise_free_map.spacegroup,
+            cell=noise_free_map.cell,
+            high_resolution_limit=resolution,
+        )
+
+
 @pytest.mark.parametrize("map_sampling", [1, 2, 2.25, 3, 5])
 def test_ccp4_map_round_trip(
     map_sampling: int,

test/unit/test_tv.py

@@ -3,8 +3,8 @@
 from typing import Sequence
 
 import numpy as np
+import pandas as pd
 import pytest
-import reciprocalspaceship as rs
 
 from meteor import tv
 from meteor.rsmap import Map
@@ -28,14 +28,14 @@ def rms_between_coefficients(map1: Map, map2: Map) -> float:
     "lambda_values_to_scan",
     [
         None,
-        [0.01],
+        [-1.0, 0.0, 1.0],
     ],
 )
 @pytest.mark.parametrize("full_output", [False, True])
 def test_tv_denoise_map_smoke(
     lambda_values_to_scan: None | Sequence[float],
     full_output: bool,
-    random_difference_map: rs.DataSet,
+    random_difference_map: Map,
 ) -> None:
     output = tv.tv_denoise_difference_map(
         random_difference_map,
@@ -44,10 +44,22 @@ def test_tv_denoise_map_smoke(
     )  # type: ignore[call-overload]
     if full_output:
         assert len(output) == 2
-        assert isinstance(output[0], rs.DataSet)
+        assert isinstance(output[0], Map)
         assert isinstance(output[1], tv.TvDenoiseResult)
     else:
-        assert isinstance(output, rs.DataSet)
+        assert isinstance(output, Map)
+
+
+def test_tv_denoise_zero_weight(random_difference_map: Map) -> None:
+    weight = 0.0
+    output = tv.tv_denoise_difference_map(
+        random_difference_map,
+        lambda_values_to_scan=[weight],
+        full_output=False,
+    )
+    random_difference_map.canonicalize_amplitudes()
+    output.canonicalize_amplitudes()
+    pd.testing.assert_frame_equal(random_difference_map, output, rtol=1e-3)
 
 
 @pytest.mark.parametrize("lambda_values_to_scan", [None, DEFAULT_LAMBDA_VALUES_TO_SCAN])