Skip to content

Commit

Permalink
add: use ws interp for nyq grid
Browse files Browse the repository at this point in the history
  • Loading branch information
Tieqiong committed Jan 24, 2025
1 parent 792e22d commit 6436942
Show file tree
Hide file tree
Showing 2 changed files with 75 additions and 17 deletions.
23 changes: 23 additions & 0 deletions news/ws.rst
Original file line number Diff line number Diff line change
@@ -0,0 +1,23 @@
**Added:**

* Added Whittaker-Shannon interpolation option for grid_interpolation.

**Changed:**

* Use WS interpolation for Nyquist grid.

**Deprecated:**

* <news item>

**Removed:**

* <news item>

**Fixed:**

* <news item>

**Security:**

* <news item>
69 changes: 52 additions & 17 deletions src/diffpy/pdfgui/control/fitdataset.py
Original file line number Diff line number Diff line change
Expand Up @@ -22,6 +22,7 @@
from diffpy.pdfgui.control.controlerrors import ControlStatusError
from diffpy.pdfgui.control.parameter import Parameter
from diffpy.pdfgui.control.pdfdataset import PDFDataSet
from diffpy.utils.resampler import wsinterp


class FitDataSet(PDFDataSet):
Expand Down Expand Up @@ -610,13 +611,14 @@ def _updateRcalcSampling(self):
if frmax - (rcalcfirst + nrcalc * frstep) > frstep * 1e-8:
nrcalc += 1
newrcalc = rcalcfirst + frstep * numpy.arange(nrcalc + 1)
tp = self.getFitSamplingType()
# Gcalc:
if len(self._Gcalc) > 0:
newGcalc = grid_interpolation(self._rcalc, self._Gcalc, newrcalc)
newGcalc = grid_interpolation(self._rcalc, self._Gcalc, newrcalc, tp)
self._Gcalc = list(newGcalc)
# dGcalc
if len(self._dGcalc) > 0:
newdGcalc = grid_interpolation(self._rcalc, self._dGcalc, newrcalc)
newdGcalc = grid_interpolation(self._rcalc, self._dGcalc, newrcalc, tp)
self._dGcalc = list(newdGcalc)
# invalidate Gtrunc and dGtrunc
self._Gtrunc = []
Expand Down Expand Up @@ -709,7 +711,8 @@ def _set_dGcalc(self, value):
def _get_Gtrunc(self):
self._updateRcalcSampling()
if not self._Gtrunc:
newGtrunc = grid_interpolation(self.robs, self.Gobs, self.rcalc)
tp = self.getFitSamplingType()
newGtrunc = grid_interpolation(self.robs, self.Gobs, self.rcalc, tp)
self._Gtrunc = list(newGtrunc)
return self._Gtrunc

Expand All @@ -724,13 +727,15 @@ def _set_Gtrunc(self, value):
def _get_dGtrunc(self):
self._updateRcalcSampling()
if not self._dGtrunc:
tp = self.getFitSamplingType()
# use sum to avoid index error for empty arrays
newdGtrunc = grid_interpolation(
self.robs,
self.dGobs,
self.rcalc,
youtleft=sum(self.dGobs[:1]),
youtright=sum(self.dGobs[-1:]),
tp,
left=sum(self.dGobs[:1]),
right=sum(self.dGobs[-1:]),
)
self._dGtrunc = list(newdGtrunc)
return self._dGtrunc
Expand Down Expand Up @@ -775,19 +780,8 @@ def _set_crw(self, value):
##############################################################################
# helper functions
##############################################################################
def _linear_interpolation(x0, y0, x1, youtleft, youtright):


def grid_interpolation(x0, y0, x1, youtleft=0.0, youtright=0.0):
"""Linear interpolation of x0, y0 values to a new grid x1.
x0 -- original x-grid, must be equally spaced
y0 -- original y values
x1 -- new x-grid, it can have any spacing
youtleft -- value for interpolated y1 for x1 below the x0 range
youtright -- value for interpolated y1 for x1 above the x0 range
Return numpy.array of interpolated y1 values.
"""
x0 = numpy.asarray(x0, copy=None, dtype=float)
y0 = numpy.asarray(y0, copy=None, dtype=float)
n0 = len(x0)
Expand Down Expand Up @@ -819,6 +813,47 @@ def grid_interpolation(x0, y0, x1, youtleft=0.0, youtright=0.0):
y1[m1] = w0lo * y0[ilo0] + w0hi * y0[ihi0]
return y1

def grid_interpolation(x0, y0, x1, tp, left=None, right=None):
"""
Interpolate values from one grid onto another using either linear
or Whittaker–Shannon interpolation.
Parameters
----------
x0 : array_like
Original x-grid, must be equally spaced.
y0 : array_like
Original values defined on x0.
x1 : array_like
New x-grid upon which to interpolate.
tp : {'data', 'Nyquist', 'custom'}, optional
Corresponding fit sampling type. Use Whittaker–Shannon interpolation
for Nyquist resampling and linear interpolation otherwise.
If not provided, linear interpolation is used.
left : float, optional
Value for interpolated y1 for x1 below the x0 range.
Default: if tp='Nyquist' then y1[0] is used. Otherwise 0.0 is used.
right : float, optional
Value for interpolated y1 for x1 above the x0 range.
Default: if tp='Nyquist' then y1[-1] is used. Otherwise 0.0 is used.
Returns
-------
numpy.ndarray
Array of interpolated values on the new grid x1.
Notes
-----
When tp='Nyquist', the function calls :func:`wsinterp` to perform Whittaker–Shannon interpolation.
Otherwise it uses the internal :func:`_linear_interpolation` routine.
"""
if tp == 'Nyquist':
return wsinterp(x1, x0, y0, left, right)
else:
left = 0.0 if left is None else left
right = 0.0 if right is None else right
return _linear_interpolation(x0, y0, x1, left, right)


# simple test code
if __name__ == "__main__":
Expand Down

0 comments on commit 6436942

Please sign in to comment.