From ebfd3edb2d710a8a13b8c5ae53e49361f6b09b5f Mon Sep 17 00:00:00 2001 From: Sparks29032 Date: Sat, 14 Dec 2024 20:02:10 -0500 Subject: [PATCH 1/2] nsinterp --- doc/source/examples/resampleexample.rst | 10 ++++++- news/nsinterp.rst | 23 ++++++++++++++ src/diffpy/utils/resampler.py | 40 +++++++++++++++++++++++++ tests/test_resample.py | 22 +++++++++++++- 4 files changed, 93 insertions(+), 2 deletions(-) create mode 100644 news/nsinterp.rst diff --git a/doc/source/examples/resampleexample.rst b/doc/source/examples/resampleexample.rst index f7ba5e18..2757b791 100644 --- a/doc/source/examples/resampleexample.rst +++ b/doc/source/examples/resampleexample.rst @@ -58,7 +58,7 @@ given enough datapoints. nickel_resample = wsinterp(grid, nickel_grid, nickel_func) target_resample = wsinterp(grid, target_grid, target_func) - We can now plot the difference to see that these two functions are quite similar.: + We can now plot the difference to see that these two functions are quite similar.:: plt.plot(grid, target_resample) plt.plot(grid, nickel_resample) @@ -78,3 +78,11 @@ given enough datapoints. In the case of our dataset, our band-limit is ``qmax=25.0`` and our function spans :math:`r \in (0.0, 60.0)`. Thus, our original grid requires :math:`25.0 * 60.0 / \pi < 478`. Since our grid has :math:`601` datapoints, our reconstruction was perfect as shown from the comparison between ``Nickel.gr`` and ``NiTarget.gr``. + + This computation is implemented in the function ``nsinterp``.:: + + from diffpy.utils.resampler import nsinterp + qmin = 0 + qmax = 25 + nickel_resample = (nickel_grid, nickel_func, qmin, qmax) + diff --git a/news/nsinterp.rst b/news/nsinterp.rst new file mode 100644 index 00000000..9b716b87 --- /dev/null +++ b/news/nsinterp.rst @@ -0,0 +1,23 @@ +**Added:** + +* Function nsinterp for automatic interpolation onto the Nyquist-Shannon grid. + +**Changed:** + +* + +**Deprecated:** + +* + +**Removed:** + +* + +**Fixed:** + +* + +**Security:** + +* diff --git a/src/diffpy/utils/resampler.py b/src/diffpy/utils/resampler.py index 5f6062d2..76108404 100644 --- a/src/diffpy/utils/resampler.py +++ b/src/diffpy/utils/resampler.py @@ -77,6 +77,46 @@ def wsinterp(x, xp, fp, left=None, right=None): return fp_at_x +def nsinterp(xp, fp, qmin=0, qmax=25, left=None, right=None): + """One-dimensional Whittaker-Shannon interpolation onto the Nyquist-Shannon grid. + + Takes a band-limited function fp and original grid xp and resamples fp on the NS grid. + Uses the minimum number of points N required by the Nyquist sampling theorem. + N = (qmax-qmin)(rmax-rmin)/pi, where rmin and rmax are the ends of the real-space ranges. + fp must be finite, and the user inputs qmin and qmax of the frequency-domain. + + Parameters + ---------- + xp: ndarray + The array of known x values. + fp: ndarray + The array of y values associated with xp. + qmin: float + The lower band limit in the frequency domain. + qmax: float + The upper band limit in the frequency domain. + + Returns + ------- + x: ndarray + The Nyquist-Shannon grid computed for the given qmin and qmax. + fp_at_x: ndarray + The interpolated values at points x. Returns a single float if x is a scalar, + otherwise returns a numpy.ndarray. + """ + # Ensure numpy array + xp = np.array(xp) + rmin = np.min(xp) + rmax = np.max(xp) + + nspoints = int(np.round((qmax-qmin)*(rmax-rmin)/np.pi)) + + x = np.linspace(rmin, rmax, nspoints) + fp_at_x = wsinterp(x, xp, fp) + + return x, fp_at_x + + def resample(r, s, dr): """Resample a PDF on a new grid. diff --git a/tests/test_resample.py b/tests/test_resample.py index 784932bb..e43af070 100644 --- a/tests/test_resample.py +++ b/tests/test_resample.py @@ -3,7 +3,7 @@ import numpy as np import pytest -from diffpy.utils.resampler import wsinterp +from diffpy.utils.resampler import wsinterp, nsinterp def test_wsinterp(): @@ -30,3 +30,23 @@ def test_wsinterp(): 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)) + + +def test_nsinterp(): + # Create a cosine function cos(2x) for x \in [0, 3pi] + xp = np.linspace(0, 3*np.pi, 100) + fp = np.cos(2 * xp) + + # Want to resample onto the grid {0, pi, 2pi, 3pi} + x = np.array([0, np.pi, 2*np.pi, 3*np.pi]) + + # Get wsinterp result + ws_f = wsinterp(x, xp, fp) + + # Use nsinterp with qmin-qmax=4/3 + qmin = np.random.uniform() + qmax = qmin + 4/3 + ns_x, ns_f = nsinterp(xp, fp, qmin, qmax) + + assert np.allclose(x, ns_x) + assert np.allclose(ws_f, ns_f) From b1213b1f557de11f9aa20b6ad20974fbe50fb598 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Sun, 15 Dec 2024 01:07:07 +0000 Subject: [PATCH 2/2] [pre-commit.ci] auto fixes from pre-commit hooks --- doc/source/examples/resampleexample.rst | 1 - src/diffpy/utils/resampler.py | 2 +- tests/test_resample.py | 8 ++++---- 3 files changed, 5 insertions(+), 6 deletions(-) diff --git a/doc/source/examples/resampleexample.rst b/doc/source/examples/resampleexample.rst index 2757b791..3166e42b 100644 --- a/doc/source/examples/resampleexample.rst +++ b/doc/source/examples/resampleexample.rst @@ -85,4 +85,3 @@ given enough datapoints. qmin = 0 qmax = 25 nickel_resample = (nickel_grid, nickel_func, qmin, qmax) - diff --git a/src/diffpy/utils/resampler.py b/src/diffpy/utils/resampler.py index 76108404..890f29fb 100644 --- a/src/diffpy/utils/resampler.py +++ b/src/diffpy/utils/resampler.py @@ -109,7 +109,7 @@ def nsinterp(xp, fp, qmin=0, qmax=25, left=None, right=None): rmin = np.min(xp) rmax = np.max(xp) - nspoints = int(np.round((qmax-qmin)*(rmax-rmin)/np.pi)) + nspoints = int(np.round((qmax - qmin) * (rmax - rmin) / np.pi)) x = np.linspace(rmin, rmax, nspoints) fp_at_x = wsinterp(x, xp, fp) diff --git a/tests/test_resample.py b/tests/test_resample.py index e43af070..6e6294ed 100644 --- a/tests/test_resample.py +++ b/tests/test_resample.py @@ -3,7 +3,7 @@ import numpy as np import pytest -from diffpy.utils.resampler import wsinterp, nsinterp +from diffpy.utils.resampler import nsinterp, wsinterp def test_wsinterp(): @@ -34,18 +34,18 @@ def test_wsinterp(): def test_nsinterp(): # Create a cosine function cos(2x) for x \in [0, 3pi] - xp = np.linspace(0, 3*np.pi, 100) + xp = np.linspace(0, 3 * np.pi, 100) fp = np.cos(2 * xp) # Want to resample onto the grid {0, pi, 2pi, 3pi} - x = np.array([0, np.pi, 2*np.pi, 3*np.pi]) + x = np.array([0, np.pi, 2 * np.pi, 3 * np.pi]) # Get wsinterp result ws_f = wsinterp(x, xp, fp) # Use nsinterp with qmin-qmax=4/3 qmin = np.random.uniform() - qmax = qmin + 4/3 + qmax = qmin + 4 / 3 ns_x, ns_f = nsinterp(xp, fp, qmin, qmax) assert np.allclose(x, ns_x)