Skip to content

Squeeze morph: Adding UCs tests #182

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 9 commits into from
Apr 19, 2025
23 changes: 23 additions & 0 deletions news/morphsqueeze.rst
Original file line number Diff line number Diff line change
@@ -0,0 +1,23 @@
**Added:**

* Polynomial squeeze of x-axis of morphed data

**Changed:**

* <news item>

**Deprecated:**

* <news item>

**Removed:**

* <news item>

**Fixed:**

* <news item>

**Security:**

* <news item>
73 changes: 73 additions & 0 deletions src/diffpy/morph/morphs/morphsqueeze.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,73 @@
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"]

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,
min_index, max_index) where the target values remain the same and
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

we don't need to return min and max index. These can be just updated in place. Also, these names are not so intuitive to me. maybe something like extrap_index_low and extrap_index_high? Let's define these and set them to None when MorphSqueeze is instantiated.

the morph data is shifted according to the squeeze. The morphed
data is returned on the same grid as the unmorphed data.
The min_index and max_index are the last index before the
interpolated region and the first index after the interpolated
region, respectively. If there is no extrapolation it returns None.

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
)
left_extrap = np.where(self.x_morph_in < x_squeezed[0])[0]
right_extrap = np.where(self.x_morph_in > x_squeezed[-1])[0]
min_index = left_extrap[-1] if left_extrap.size else None
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

self.extrap_index_low (or left) and don't return it.

Describe this variable in the constructor where it is first defined.

max_index = right_extrap[0] if right_extrap.size else None
return self.xyallout + (min_index, max_index)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

apart from it being a bad idea in general, this breaks backwards compatibility.

78 changes: 78 additions & 0 deletions tests/test_morphsqueeze.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,78 @@
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.2, -0.01, -0.001, -0.0001],
# 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.007],
# Linear squeeze coefficients
[0.1, 0.3],
# 4th order squeeze coefficients
[0.2, -0.01, 0.001, -0.001, 0.0004],
# 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
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is now fantastic! I love these tests. Future us will love you.

(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, 221), 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)
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,
low_extrap_idx,
high_extrap_idx,
) = morph(x_morph, y_morph, x_target, y_target)
if low_extrap_idx is None and high_extrap_idx is None:
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

this seems too complicated to me, and also does it test the extrapolation? How about something like

if low_extrap_idx is None: 
    low_extrap_idx = 0
elif high_extrap_idx is None:
    high_extrap_idx = -1
assert np.allclose(y_morph_actual[low_extrap_index:high_extrap_idx], y_morph_expected[low_extrap_index:high_extrap_idx], atol=1e-6)
assert np.allclose(y_morph_actual[:low_extrap_index-1], y_morph_expected[:low_extrap_index-1], atol=1e-5)
assert np.allclose(y_morph_actual[high_extrap_index+1:], y_morph_expected[high_extrap_index+1:], atol=1e-5)
assert morph.low_extrap_idx == expected_low_extrap_idx
assert morph.high_extrap_idx == expected_high_extrap_idx

assert np.allclose(y_morph_actual, y_morph_expected, atol=1e-6)
else:
interp_start = low_extrap_idx + 1 if low_extrap_idx is not None else 0
interp_end = (
high_extrap_idx
if high_extrap_idx is not None
else len(y_morph_actual)
)
assert np.allclose(
y_morph_actual[interp_start:interp_end],
y_morph_expected[interp_start:interp_end],
atol=1e-6,
)
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)
Loading