Skip to content

Improve wsinterp docstrings and test function #221

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 6 commits into from
Dec 13, 2024
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion README.rst
Original file line number Diff line number Diff line change
Expand Up @@ -43,7 +43,7 @@ General purpose shared utilities for the diffpy libraries.
The diffpy.utils package provides functions for extracting array data from
variously formatted text files, an interpolation function based on the
Whittaker-Shannon formula that can be used to resample a PDF or other profile
function over a new grid, `Diffraction_objects` for conveniently manipulating
function over a new grid, `DiffractionObject` for conveniently manipulating
diffraction data, and some wx GUI utilities used by the PDFgui
program.

Expand Down
2 changes: 0 additions & 2 deletions src/diffpy/utils/__init__.py
Original file line number Diff line number Diff line change
Expand Up @@ -20,5 +20,3 @@

# silence the pyflakes syntax checker
assert __version__ or True

# End of file
2 changes: 0 additions & 2 deletions src/diffpy/utils/parsers/__init__.py
Original file line number Diff line number Diff line change
Expand Up @@ -15,5 +15,3 @@

"""Various utilities related to data parsing and manipulation.
"""

# End of file
3 changes: 0 additions & 3 deletions src/diffpy/utils/parsers/loaddata.py
Original file line number Diff line number Diff line change
Expand Up @@ -344,6 +344,3 @@ def isfloat(s):
except ValueError:
pass
return False


# End of file
61 changes: 29 additions & 32 deletions src/diffpy/utils/resampler.py
Original file line number Diff line number Diff line change
Expand Up @@ -15,25 +15,26 @@

"""Various utilities related to data parsing and manipulation."""

import numpy
import numpy as np


# NOTE - this should be faster than resample below and conforms more closely to
# numpy.interp. I'm keeping resample for legacy reasons.
def wsinterp(x, xp, fp, left=None, right=None):
"""One-dimensional Whittaker-Shannon interpolation.

This uses the Whittaker-Shannon interpolation formula to interpolate the value of fp (array),
which is defined over xp (array), at x (array or float).
Reconstruct a continuous signal from discrete data points by utilizing sinc functions
as interpolation kernels. This function interpolates the values of fp (array),
which are defined over xp (array), at new points x (array or float).
The implementation is based on E. T. Whittaker's 1915 paper
(https://doi.org/10.1017/S0370164600017806).

Parameters
----------
x: ndarray
Desired range for interpolation.
The x values at which interpolation is computed.
xp: ndarray
Defined range for fp.
The array of known x values.
fp: ndarray
Function to be interpolated.
The array of y values associated with xp.
left: float
If given, set fp for x < xp[0] to left. Otherwise, if left is None (default) or not given,
set fp for x < xp[0] to fp evaluated at xp[-1].
Expand All @@ -43,24 +44,23 @@ def wsinterp(x, xp, fp, left=None, right=None):

Returns
-------
float:
If input x is a scalar (not an array), return the interpolated value at x.
ndarray:
If input x is an array, return the interpolated array with dimensions of x.
ndarray or float
The interpolated values at points x. Returns a single float if x is a scalar,
otherwise returns a numpy.ndarray.
"""
scalar = numpy.isscalar(x)
scalar = np.isscalar(x)
if scalar:
x = numpy.array(x)
x = np.array(x)
x.resize(1)
# shape = (nxp, nx), nxp copies of x data span axis 1
u = numpy.resize(x, (len(xp), len(x)))
u = np.resize(x, (len(xp), len(x)))
# Must take transpose of u for proper broadcasting with xp.
# shape = (nx, nxp), v(xp) data spans axis 1
v = (xp - u.T) / (xp[1] - xp[0])
# shape = (nx, nxp), m(v) data spans axis 1
m = fp * numpy.sinc(v)
m = fp * np.sinc(v)
# Sum over m(v) (axis 1)
fp_at_x = numpy.sum(m, axis=1)
fp_at_x = np.sum(m, axis=1)

# Enforce left and right
if left is None:
Expand Down Expand Up @@ -100,36 +100,33 @@ def resample(r, s, dr):
dr0 = r[1] - r[0] # Constant timestep

if dr0 < dr:
rnew = numpy.arange(r[0], r[-1] + 0.5 * dr, dr)
snew = numpy.interp(rnew, r, s)
rnew = np.arange(r[0], r[-1] + 0.5 * dr, dr)
snew = np.interp(rnew, r, s)
return rnew, snew

elif dr0 > dr:
# Tried to pad the end of s to dampen, but nothing works.
# m = (s[-1] - s[-2]) / dr0
# b = (s[-2] * r[-1] - s[-1] * r[-2]) / dr0
# rpad = r[-1] + numpy.arange(1, len(s))*dr0
# rpad = r[-1] + np.arange(1, len(s))*dr0
# spad = rpad * m + b
# spad = numpy.concatenate([s,spad])
# rnew = numpy.arange(0, rpad[-1], dr)
# snew = numpy.zeros_like(rnew)
# spad = np.concatenate([s,spad])
# rnew = np.arange(0, rpad[-1], dr)
# snew = np.zeros_like(rnew)
# Accommodate for the fact that r[0] might not be 0
# u = (rnew-r[0]) / dr0
# for n in range(len(spad)):
# snew += spad[n] * numpy.sinc(u - n)
# snew += spad[n] * np.sinc(u - n)

# sel = numpy.logical_and(rnew >= r[0], rnew <= r[-1])
# sel = np.logical_and(rnew >= r[0], rnew <= r[-1])

rnew = numpy.arange(0, r[-1], dr)
snew = numpy.zeros_like(rnew)
rnew = np.arange(0, r[-1], dr)
snew = np.zeros_like(rnew)
u = (rnew - r[0]) / dr0
for n in range(len(s)):
snew += s[n] * numpy.sinc(u - n)
sel = numpy.logical_and(rnew >= r[0], rnew <= r[-1])
snew += s[n] * np.sinc(u - n)
sel = np.logical_and(rnew >= r[0], rnew <= r[-1])
return rnew[sel], snew[sel]

# If we got here, then no resampling is required
return r.copy(), s.copy()


# End of file
2 changes: 0 additions & 2 deletions src/diffpy/utils/version.py
Original file line number Diff line number Diff line change
Expand Up @@ -22,5 +22,3 @@
from importlib.metadata import version

__version__ = version("diffpy.utils")

# End of file
2 changes: 0 additions & 2 deletions src/diffpy/utils/wx/__init__.py
Original file line number Diff line number Diff line change
Expand Up @@ -15,5 +15,3 @@

"""Utilities related wx Python GUIs.
"""

# End of file
3 changes: 0 additions & 3 deletions src/diffpy/utils/wx/gridutils.py
Original file line number Diff line number Diff line change
Expand Up @@ -163,6 +163,3 @@ def _indicesToBlocks(indices):
i0 = i
rv = [tuple(ij) for ij in rngs]
return rv


# End of file
28 changes: 15 additions & 13 deletions tests/test_resample.py
Original file line number Diff line number Diff line change
@@ -1,30 +1,32 @@
import random

import numpy as np
import pytest

from diffpy.utils.resampler import wsinterp


def test_wsinterp():
import random

# Check known points are unchanged by interpolation
# FIXME: if another SW interp function exists, run comparisons for interpolated points

# Sampling rate
ssr = 44100**-1 # Standard sampling rate for human-hearable frequencies
t = ssr

# Creating a symmetric set of sample points around zero.
n = 5
xp = np.array([i * t for i in range(-n, n + 1, 1)])
x = np.array([i * t for i in range(-n - 1, n + 2, 1)])
xp = np.array([i * ssr for i in range(-n, n + 1, 1)])
x = np.array([i * ssr for i in range(-n - 1, n + 2, 1)])
assert len(xp) == 11 and len(x) == 13

# Interpolate a few random datasets
# Generate a new set of fp values across 10 trial runs
trials = 10
for trial in range(trials):
fp = np.array([random.random() * ssr for i in range(-n, n + 1, 1)])

for _ in range(trials):
# Create random function values (fp) at the points defined in xp above
fp = np.array([random.random() * ssr for _ in range(-n, n + 1, 1)])
fp_at_x = wsinterp(x, xp, fp)

# Check that the known points are unchanged by interpolation
assert np.allclose(fp_at_x[1:-1], fp)
for i in range(len(x)):
assert fp_at_x[i] == pytest.approx(wsinterp(x[i], xp, fp))


if __name__ == "__main__":
test_wsinterp()
Loading