diff --git a/news/morphsqueeze.rst b/news/morphsqueeze.rst new file mode 100644 index 00000000..750fbcf1 --- /dev/null +++ b/news/morphsqueeze.rst @@ -0,0 +1,23 @@ +**Added:** + +* Polynomial squeeze of x-axis of morphed data + +**Changed:** + +* + +**Deprecated:** + +* + +**Removed:** + +* + +**Fixed:** + +* + +**Security:** + +* diff --git a/src/diffpy/morph/morphs/morphsqueeze.py b/src/diffpy/morph/morphs/morphsqueeze.py new file mode 100644 index 00000000..e63adc8b --- /dev/null +++ b/src/diffpy/morph/morphs/morphsqueeze.py @@ -0,0 +1,74 @@ +import numpy as np +from numpy.polynomial import Polynomial +from scipy.interpolate import CubicSpline + +from diffpy.morph.morphs.morph import LABEL_GR, LABEL_RA, Morph + + +class MorphSqueeze(Morph): + """Apply a polynomial to squeeze the morph function. The morphed + data is returned on the same grid as the unmorphed data.""" + + # Define input output types + summary = "Squeeze morph by polynomial shift" + xinlabel = LABEL_RA + yinlabel = LABEL_GR + xoutlabel = LABEL_RA + youtlabel = LABEL_GR + parnames = ["squeeze"] + # extrap_index_low: last index before interpolation region + # extrap_index_high: first index after interpolation region + extrap_index_low = None + extrap_index_high = None + + def morph(self, x_morph, y_morph, x_target, y_target): + """Squeeze the morph function. + + This applies a polynomial to squeeze the morph non-linearly. + + Configuration Variables + ----------------------- + squeeze : list + The polynomial coefficients [a0, a1, ..., an] for the squeeze + function where the polynomial would be of the form + a0 + a1*x + a2*x^2 and so on. The order of the polynomial is + determined by the length of the list. + + Returns + ------- + A tuple (x_morph_out, y_morph_out, x_target_out, y_target_out) + where the target values remain the same and the morph data is + shifted according to the squeeze. The morphed data is returned on + the same grid as the unmorphed data. + + Example + ------- + Import the squeeze morph function: + >>> from diffpy.morph.morphs.morphsqueeze import MorphSqueeze + Provide initial guess for squeezing coefficients: + >>> squeeze_coeff = [0.1, -0.01, 0.005] + Run the squeeze morph given input morph array (x_morph, y_morph) + and target array (x_target, y_target): + >>> morph = MorphSqueeze() + >>> morph.squeeze = squeeze_coeff + >>> x_morph_out, y_morph_out, x_target_out, y_target_out = morph( + ... x_morph, y_morph, x_target, y_target) + To access parameters from the morph instance: + >>> x_morph_in = morph.x_morph_in + >>> y_morph_in = morph.y_morph_in + >>> x_target_in = morph.x_target_in + >>> y_target_in = morph.y_target_in + >>> squeeze_coeff_out = morph.squeeze + """ + Morph.morph(self, x_morph, y_morph, x_target, y_target) + + squeeze_polynomial = Polynomial(self.squeeze) + x_squeezed = self.x_morph_in + squeeze_polynomial(self.x_morph_in) + self.y_morph_out = CubicSpline(x_squeezed, self.y_morph_in)( + self.x_morph_in + ) + low_extrap = np.where(self.x_morph_in < x_squeezed[0])[0] + high_extrap = np.where(self.x_morph_in > x_squeezed[-1])[0] + self.extrap_index_low = low_extrap[-1] if low_extrap.size else None + self.extrap_index_high = high_extrap[0] if high_extrap.size else None + return self.xyallout diff --git a/tests/test_morphsqueeze.py b/tests/test_morphsqueeze.py new file mode 100644 index 00000000..0de38e81 --- /dev/null +++ b/tests/test_morphsqueeze.py @@ -0,0 +1,86 @@ +import numpy as np +import pytest +from numpy.polynomial import Polynomial + +from diffpy.morph.morphs.morphsqueeze import MorphSqueeze + +squeeze_coeffs_list = [ + # The order of coefficients is [a0, a1, a2, ..., an] + # Negative cubic squeeze coefficients + [-0.01, -0.0005, -0.0005, -1e-6], + # Positive cubic squeeze coefficients + [0.2, 0.01, 0.001, 0.0001], + # Positive and negative cubic squeeze coefficients + [0.2, -0.01, 0.002, -0.0001], + # Quadratic squeeze coefficients + [-0.2, 0.005, -0.0004], + # Linear squeeze coefficients + [0.1, 0.3], + # 4th order squeeze coefficients + [0.2, -0.01, 0.001, -0.001, 0.0001], + # Zeros and non-zeros, the full polynomial is applied + [0, 0.03, 0, -0.0001], + # Testing zeros, expect no squeezing + [0, 0, 0, 0, 0, 0], +] +morph_target_grids = [ + # UCs from issue 181: https://github.com/diffpy/diffpy.morph/issues/181 + # UC2: Same range and same grid density + (np.linspace(0, 10, 101), np.linspace(0, 10, 101)), + # UC4: Target range wider than morph, same grid density + (np.linspace(0, 10, 101), np.linspace(-2, 20, 221)), + # UC6: Target range wider than morph, target grid density finer than morph + (np.linspace(0, 10, 101), np.linspace(-2, 20, 421)), + # UC8: Target range wider than morph, morph grid density finer than target + (np.linspace(0, 10, 401), np.linspace(-2, 20, 200)), + # UC10: Morph range starts and ends earlier than target, same grid density + (np.linspace(-2, 10, 121), np.linspace(0, 20, 201)), + # UC12: Morph range wider than target, same grid density + (np.linspace(-2, 20, 201), np.linspace(0, 10, 101)), +] + + +@pytest.mark.parametrize("x_morph, x_target", morph_target_grids) +@pytest.mark.parametrize("squeeze_coeffs", squeeze_coeffs_list) +def test_morphsqueeze(x_morph, x_target, squeeze_coeffs): + y_target = np.sin(x_target) + squeeze_polynomial = Polynomial(squeeze_coeffs) + x_squeezed = x_morph + squeeze_polynomial(x_morph) + y_morph = np.sin(x_squeezed) + low_extrap = np.where(x_morph < x_squeezed[0])[0] + high_extrap = np.where(x_morph > x_squeezed[-1])[0] + extrap_index_low_expected = low_extrap[-1] if low_extrap.size else None + extrap_index_high_expected = high_extrap[0] if high_extrap.size else None + x_morph_expected = x_morph + y_morph_expected = np.sin(x_morph) + morph = MorphSqueeze() + morph.squeeze = squeeze_coeffs + x_morph_actual, y_morph_actual, x_target_actual, y_target_actual = morph( + x_morph, y_morph, x_target, y_target + ) + extrap_index_low = morph.extrap_index_low + extrap_index_high = morph.extrap_index_high + if extrap_index_low is None: + extrap_index_low = 0 + elif extrap_index_high is None: + extrap_index_high = -1 + assert np.allclose( + y_morph_actual[extrap_index_low + 1 : extrap_index_high], + y_morph_expected[extrap_index_low + 1 : extrap_index_high], + atol=1e-6, + ) + assert np.allclose( + y_morph_actual[:extrap_index_low], + y_morph_expected[:extrap_index_low], + atol=1e-3, + ) + assert np.allclose( + y_morph_actual[extrap_index_high:], + y_morph_expected[extrap_index_high:], + atol=1e-3, + ) + assert morph.extrap_index_low == extrap_index_low_expected + assert morph.extrap_index_high == extrap_index_high_expected + assert np.allclose(x_morph_actual, x_morph_expected) + assert np.allclose(x_target_actual, x_target) + assert np.allclose(y_target_actual, y_target)