feat: add --check-increase option

news/sort-squeezed-x.rst

@@ -0,0 +1,23 @@
+**Added:**
+
+* Add ``--check-increase`` option for squeeze morph.
+
+**Changed:**
+
+* <news item>
+
+**Deprecated:**
+
+* <news item>
+
+**Removed:**
+
+* <news item>
+
+**Fixed:**
+
+* <news item>
+
+**Security:**
+
+* <news item>

news/strictly-increasing-squeeze.rst

@@ -0,0 +1,23 @@
+**Added:**
+
+* Raise ``ValueError`` if ``x_squeezed`` is not strictly increasing.
+
+**Changed:**
+
+* <news item>
+
+**Deprecated:**
+
+* <news item>
+
+**Removed:**
+
+* <news item>
+
+**Fixed:**
+
+* <news item>
+
+**Security:**
+
+* <news item>

src/diffpy/morph/morph_io.py

@@ -408,7 +408,7 @@ def tabulate_results(multiple_morph_results):
     return tabulated_results
 
 
-def handle_warnings(squeeze_morph):
+def handle_extrapolation_warnings(squeeze_morph):
     if squeeze_morph is not None:
         extrapolation_info = squeeze_morph.extrapolation_info
         is_extrap_low = extrapolation_info["is_extrap_low"]
@@ -443,3 +443,22 @@ def handle_warnings(squeeze_morph):
                 wmsg,
                 UserWarning,
             )
+
+
+def handle_check_increase_warning(squeeze_morph):
+    if squeeze_morph is not None:
+        if squeeze_morph.strictly_increasing:
+            wmsg = None
+        else:
+            wmsg = (
+                "Warning: The squeeze morph has interpolated your morphed "
+                "function from a non-monotonically increasing grid. "
+                "This can result in strange behavior in the non-uniqe "
+                "grid regions. To disable this setting, "
+                "please enable --check-increasing."
+            )
+        if wmsg:
+            warnings.warn(
+                wmsg,
+                UserWarning,
+            )

src/diffpy/morph/morphapp.py

@@ -209,6 +209,16 @@ def custom_error(self, msg):
             "See online documentation for more information."
         ),
     )
+    group.add_option(
+        "--check-increase",
+        action="store_true",
+        dest="check_increase",
+        help=(
+            "Disable squeeze morph to interpolat morphed function "
+            "from a non-monotonically increasing grid."
+        ),
+    )
+
     group.add_option(
         "--smear",
         type="float",
@@ -573,7 +583,7 @@ def single_morph(
                 except ValueError:
                     parser.error(f"{coeff} could not be converted to float.")
         squeeze_poly_deg = len(squeeze_dict_in.keys())
-        squeeze_morph = morphs.MorphSqueeze()
+        squeeze_morph = morphs.MorphSqueeze(check_increase=opts.check_increase)
         chain.append(squeeze_morph)
         config["squeeze"] = squeeze_dict_in
         # config["extrap_index_low"] = None
@@ -707,9 +717,10 @@ def single_morph(
         chain(x_morph, y_morph, x_target, y_target)
 
     # THROW ANY WARNINGS HERE
-    io.handle_warnings(squeeze_morph)
-    io.handle_warnings(shift_morph)
-    io.handle_warnings(stretch_morph)
+    io.handle_extrapolation_warnings(squeeze_morph)
+    io.handle_check_increase_warning(squeeze_morph)
+    io.handle_extrapolation_warnings(shift_morph)
+    io.handle_extrapolation_warnings(stretch_morph)
 
     # Get Rw for the morph range
     rw = tools.getRw(chain)

src/diffpy/morph/morphpy.py

@@ -51,6 +51,7 @@ def __get_morph_opts__(parser, scale, stretch, smear, plot, **kwargs):
         "reverse",
         "diff",
         "get-diff",
+        "check-increase",
     ]
     opts_to_ignore = ["multiple-morphs", "multiple-targets"]
     for opt in opts_storing_values:

src/diffpy/morph/morphs/morphsqueeze.py

@@ -1,6 +1,7 @@
 """Class MorphSqueeze -- Apply a polynomial to squeeze the morph
 function."""
 
+import numpy
 from numpy.polynomial import Polynomial
 from scipy.interpolate import CubicSpline
 
@@ -68,8 +69,58 @@ class MorphSqueeze(Morph):
     squeeze_cutoff_low = None
     squeeze_cutoff_high = None
 
-    def __init__(self, config=None):
+    def __init__(self, config=None, check_increase=False):
         super().__init__(config)
+        self.check_increase = check_increase
+
+    def _ensure_strictly_increase(self, x, x_sorted):
+        if list(x) != list(x_sorted):
+            if self.check_increase:
+                raise ValueError(
+                    "Error: The polynomial applied by the squeeze morph has "
+                    "resulted in a grid that is no longer strictly increasing"
+                    ", likely due to a convergence issue. A strictly "
+                    "increasing grid is required for diffpy.morph to compute "
+                    "the morphed function through cubic spline interpolation. "
+                    "Here are some suggested methods to resolve this:\n"
+                    "(1) Please decrease the order of your polynomial and "
+                    "try again.\n"
+                    "(2) If you are using initial guesses of all 0, please "
+                    "ensure your objective function only requires a small "
+                    "polynomial squeeze to match your reference. (In other "
+                    "words, there is good agreement between the two functions"
+                    ".)\n"
+                    "(3) If you expect a large polynomial squeeze to be needed"
+                    ", please ensure your initial parameters for the "
+                    "polynomial morph result in good agreement between your "
+                    "reference and objective functions. One way to obtain "
+                    "such parameters is to first apply a --hshift and "
+                    "--stretch morph. Then, use the hshift parameter for a0 "
+                    "and stretch parameter for a1."
+                )
+            else:
+                self.strictly_increasing = False
+        else:
+            self.strictly_increasing = True
+
+    def _sort_squeeze(self, x, y):
+        """Sort x,y according to the value of x."""
+        xy = list(zip(x, y))
+        xy_sorted = sorted(xy, key=lambda pair: pair[0])
+        x_sorted, y_sorted = list(zip(*xy_sorted))
+        return x_sorted, y_sorted
+
+    def _handle_duplicates(self, x, y):
+        """Remove duplicated x and use the mean value of y corresponded
+        to the duplicated x."""
+        unq_x, unq_inv = numpy.unique(x, return_inverse=True)
+        if len(unq_x) == len(x):
+            return x, y
+        else:
+            y_avg = numpy.zeros_like(unq_x)
+            for i in range(len(unq_x)):
+                y_avg[i] = numpy.array(y)[unq_inv == i].mean()
+            return unq_x, y_avg
 
     def morph(self, x_morph, y_morph, x_target, y_target):
         """Apply a polynomial to squeeze the morph function.
@@ -78,13 +129,19 @@ def morph(self, x_morph, y_morph, x_target, y_target):
         data.
         """
         Morph.morph(self, x_morph, y_morph, x_target, y_target)
-
         coeffs = [self.squeeze[f"a{i}"] for i in range(len(self.squeeze))]
         squeeze_polynomial = Polynomial(coeffs)
         x_squeezed = self.x_morph_in + squeeze_polynomial(self.x_morph_in)
-        self.y_morph_out = CubicSpline(x_squeezed, self.y_morph_in)(
+        x_squeezed_sorted, y_morph_sorted = self._sort_squeeze(
+            x_squeezed, self.y_morph_in
+        )
+        self._ensure_strictly_increase(x_squeezed, x_squeezed_sorted)
+        x_squeezed_sorted, y_morph_sorted = self._handle_duplicates(
+            x_squeezed_sorted, y_morph_sorted
+        )
+        self.y_morph_out = CubicSpline(x_squeezed_sorted, y_morph_sorted)(
             self.x_morph_in
         )
-        self.set_extrapolation_info(x_squeezed, self.x_morph_in)
+        self.set_extrapolation_info(x_squeezed_sorted, self.x_morph_in)
 
         return self.xyallout

tests/test_morphsqueeze.py

@@ -170,3 +170,175 @@ def test_morphsqueeze_extrapolate(user_filesystem, squeeze_coeffs, wmsg_gen):
     )
     with pytest.warns(UserWarning, match=expected_wmsg):
         single_morph(parser, opts, pargs, stdout_flag=False)
+
+
+@pytest.mark.parametrize(
+    "squeeze_coeffs, x_morph",
+    [
+        ({"a0": 0.01, "a1": 0.01, "a2": -0.1}, np.linspace(0, 10, 101)),
+    ],
+)
+def test_non_strictly_increasing_squeeze(squeeze_coeffs, x_morph):
+    x_target = x_morph
+    y_target = np.sin(x_target)
+    coeffs = [squeeze_coeffs[f"a{i}"] for i in range(len(squeeze_coeffs))]
+    squeeze_polynomial = Polynomial(coeffs)
+    x_squeezed = x_morph + squeeze_polynomial(x_morph)
+    # non-strictly-increasing
+    assert not np.all(np.sign(np.diff(x_squeezed)) > 0)
+    y_morph = np.sin(x_squeezed)
+    # all zero initial guess
+    morph_results = morphpy.morph_arrays(
+        np.array([x_morph, y_morph]).T,
+        np.array([x_target, y_target]).T,
+        squeeze=[0, 0, 0],
+        apply=True,
+    )
+    _, y_morph_actual = morph_results[1].T  # noqa: F841
+    y_morph_expected = np.sin(x_morph)  # noqa: F841
+    # squeeze morph extrapolates.
+    #   Need to extract extrap_index from morph_results to examine
+    #   the convergence.
+    # assert np.allclose(y_morph_actual, y_morph_expected, atol=1e-3)
+    # Raise warning when called without --check-increase
+    with pytest.warns() as w:
+        morph_results = morphpy.morph_arrays(
+            np.array([x_morph, y_morph]).T,
+            np.array([x_target, y_target]).T,
+            squeeze=[0.01, 0.01, -0.1],
+            apply=True,
+        )
+    assert w[0].category is UserWarning
+    actual_wmsg = " ".join([str(w[i].message) for i in range(len(w))])
+    expected_wmsg = (
+        "Warning: The squeeze morph has interpolated your morphed "
+        "function from a non-monotonically increasing grid. "
+    )
+    assert expected_wmsg in actual_wmsg
+    _, y_morph_actual = morph_results[1].T  # noqa: F841
+    y_morph_expected = np.sin(x_morph)  # noqa: F841
+    # squeeze morph extrapolates.
+    #   Need to extract extrap_index from morph_results to examine
+    #   the convergence.
+    # assert np.allclose(y_morph_actual, y_morph_expected, atol=1e-3)
+    # System exits when called with --check-increase
+    with pytest.raises(SystemExit) as excinfo:
+        morphpy.morph_arrays(
+            np.array([x_morph, y_morph]).T,
+            np.array([x_target, y_target]).T,
+            squeeze=[0.01, 0.009, -0.1],
+            check_increase=True,
+        )
+    actual_emsg = str(excinfo.value)
+    expected_emsg = "2"
+    assert expected_emsg == actual_emsg
+
+
+@pytest.mark.parametrize(
+    "squeeze_coeffs, x_morph",
+    [
+        ({"a0": -1, "a1": -1, "a2": 2}, np.linspace(-1, 1, 101)),
+        (
+            {"a0": -1, "a1": -1, "a2": 0, "a3": 0, "a4": 2},
+            np.linspace(-1, 1, 101),
+        ),
+    ],
+)
+def test_sort_squeeze_bad(user_filesystem, squeeze_coeffs, x_morph):
+    # call in .py without --check-increase
+    x_target = x_morph
+    y_target = np.sin(x_target)
+    coeffs = [squeeze_coeffs[f"a{i}"] for i in range(len(squeeze_coeffs))]
+    squeeze_polynomial = Polynomial(coeffs)
+    x_squeezed = x_morph + squeeze_polynomial(x_morph)
+    y_morph = np.sin(x_squeezed)
+    morph = MorphSqueeze()
+    morph.squeeze = squeeze_coeffs
+    with pytest.warns() as w:
+        morphpy.morph_arrays(
+            np.array([x_morph, y_morph]).T,
+            np.array([x_target, y_target]).T,
+            squeeze=coeffs,
+            apply=True,
+        )
+    assert len(w) == 1
+    assert w[0].category is UserWarning
+    actual_wmsg = str(w[0].message)
+    expected_wmsg = (
+        "Warning: The squeeze morph has interpolated your morphed "
+        "function from a non-monotonically increasing grid. "
+    )
+    assert expected_wmsg in actual_wmsg
+
+    # call in .py with --check-increase
+    with pytest.raises(ValueError) as excinfo:
+        morphpy.morph_arrays(
+            np.array([x_morph, y_morph]).T,
+            np.array([x_target, y_target]).T,
+            squeeze=coeffs,
+            check_increase=True,
+            apply=True,
+        )
+    actual_emsg = str(excinfo.value)
+    expected_emsg = (
+        "Error: The polynomial applied by the squeeze morph has "
+        "resulted in a grid that is no longer strictly increasing"
+    )
+    assert expected_emsg in actual_emsg
+
+    # call in CLI without --check-increase
+    morph_file, target_file = create_morph_data_file(
+        user_filesystem / "cwd_dir", x_morph, y_morph, x_target, y_target
+    )
+    parser = create_option_parser()
+    (opts, pargs) = parser.parse_args(
+        [
+            "--squeeze",
+            ",".join(map(str, coeffs)),
+            f"{morph_file.as_posix()}",
+            f"{target_file.as_posix()}",
+            "--apply",
+            "-n",
+        ]
+    )
+    with pytest.warns(UserWarning) as w:
+        single_morph(parser, opts, pargs, stdout_flag=False)
+    assert len(w) == 1
+    actual_wmsg = str(w[0].message)
+    assert expected_wmsg in actual_wmsg
+
+    # call in CLI with --check-increase
+    parser = create_option_parser()
+    (opts, pargs) = parser.parse_args(
+        [
+            "--squeeze",
+            ",".join(map(str, coeffs)),
+            f"{morph_file.as_posix()}",
+            f"{target_file.as_posix()}",
+            "--apply",
+            "-n",
+            "--check-increase",
+        ]
+    )
+    with pytest.raises(ValueError) as excinfo:
+        single_morph(parser, opts, pargs, stdout_flag=False)
+    actual_emsg = str(excinfo.value)
+    assert expected_emsg in actual_emsg
+
+
+def test_handle_duplicates():
+    unq_x = np.linspace(0, 11, 10)
+    iter = 10
+    morph = MorphSqueeze()
+    for i in range(iter):
+        actual_x = np.random.choice(unq_x, size=20)
+        actual_y = np.sin(actual_x)
+        actual_handled_x, actual_handled_y = morph._handle_duplicates(
+            actual_x, actual_y
+        )
+        expected_handled_x = np.unique(actual_x)
+        expected_handled_y = np.array(
+            [actual_y[actual_x == x].mean() for x in expected_handled_x]
+        )
+        assert np.allclose(actual_handled_x, expected_handled_x)
+        assert np.allclose(actual_handled_y, expected_handled_y)