diff --git a/news/sort-squeezed-x.rst b/news/sort-squeezed-x.rst new file mode 100644 index 0000000..9fe0f2a --- /dev/null +++ b/news/sort-squeezed-x.rst @@ -0,0 +1,23 @@ +**Added:** + +* Add ``--check-increase`` option for squeeze morph. + +**Changed:** + +* + +**Deprecated:** + +* + +**Removed:** + +* + +**Fixed:** + +* + +**Security:** + +* diff --git a/news/strictly-increasing-squeeze.rst b/news/strictly-increasing-squeeze.rst new file mode 100644 index 0000000..f14fabf --- /dev/null +++ b/news/strictly-increasing-squeeze.rst @@ -0,0 +1,23 @@ +**Added:** + +* Raise ``ValueError`` if ``x_squeezed`` is not strictly increasing. + +**Changed:** + +* + +**Deprecated:** + +* + +**Removed:** + +* + +**Fixed:** + +* + +**Security:** + +* diff --git a/src/diffpy/morph/morph_io.py b/src/diffpy/morph/morph_io.py index 84d5c15..fbe06ec 100644 --- a/src/diffpy/morph/morph_io.py +++ b/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, + ) diff --git a/src/diffpy/morph/morphapp.py b/src/diffpy/morph/morphapp.py index 327ed2d..65a41f3 100755 --- a/src/diffpy/morph/morphapp.py +++ b/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) diff --git a/src/diffpy/morph/morphpy.py b/src/diffpy/morph/morphpy.py index bac6eee..3889f82 100644 --- a/src/diffpy/morph/morphpy.py +++ b/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: diff --git a/src/diffpy/morph/morphs/morphsqueeze.py b/src/diffpy/morph/morphs/morphsqueeze.py index 3d0250d..2bfbf90 100644 --- a/src/diffpy/morph/morphs/morphsqueeze.py +++ b/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 diff --git a/tests/test_morphsqueeze.py b/tests/test_morphsqueeze.py index 07b9937..e35c89e 100644 --- a/tests/test_morphsqueeze.py +++ b/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)