Commit 867a8c8f authored by Ben Farr's avatar Ben Farr
Browse files

Add option for different out-of-bounds values in 1D interp

parent 2287620b
......@@ -13,7 +13,7 @@ def interp_from_txt(infile, flow=None, fhigh=None, asd_file=False):
freqs, psd = np.loadtxt(infile, unpack=True)
if asd_file:
psd = np.square(psd)
return BoundedInterp1D(freqs, psd, low=flow, high=fhigh)
return BoundedInterp1D(freqs, psd, low=flow, high=fhigh, below_domain_val=np.inf, above_domain_val=np.inf)
def interp_bayesline_from_txt(infile, ignore_lines=False, flow=None, fhigh=None):
"""Read a bayesline PSD from file, with the option of ignoring lines"""
......@@ -22,4 +22,4 @@ def interp_bayesline_from_txt(infile, ignore_lines=False, flow=None, fhigh=None)
psd = spline if ignore_lines else spline + line
psd /= 2.
return BoundedInterp1D(freqs, psd, low=flow, high=fhigh)
return BoundedInterp1D(freqs, psd, low=flow, high=fhigh, below_domain_val=np.inf, above_domain_val=np.inf)
......@@ -12,7 +12,7 @@ from scipy.interpolate import interp1d
class BoundedInterp1D(object):
"""A bounded 1-D interpolant that returns the specified value outside the
input domain"""
def __init__(self, x, y, low=None, high=None, outside_domain_val=np.inf, **kwargs):
def __init__(self, x, y, low=None, high=None, below_domain_val=np.nan, above_domain_val=np.nan, **kwargs):
self._dtype = y.dtype
self._complex = np.iscomplexobj(y)
......@@ -26,7 +26,8 @@ class BoundedInterp1D(object):
else:
self._high = x.max()
self._outside_domain_val = outside_domain_val
self._below_domain_val = below_domain_val
self._above_domain_val = above_domain_val
self._interp_re = interp1d(x, np.real(y), **kwargs)
if self._complex:
......@@ -43,16 +44,24 @@ class BoundedInterp1D(object):
return self._high
@property
def outside_domain_val(self):
"""Value returned when outside of the domain"""
return self._outside_domain_val
def below_domain_val(self):
"""Value returned when below the lower domain boundary"""
return self._below_domain_val
@property
def above_domain_val(self):
"""Value returned when above the upper domain boundary"""
return self._above_domain_val
def __call__(self, pts):
pts = np.atleast_1d(pts)
result = np.empty_like(pts, dtype=self._dtype)
in_bounds = (pts > self.low) & (pts < self.high)
result[~in_bounds] = self.outside_domain_val
below = pts < self.low
above = pts > self.high
in_bounds = (~below) & (~above)
result[below] = self.below_domain_val
result[above] = self.above_domain_val
result[in_bounds] = self._interp_re(pts[in_bounds])
if self._complex:
result[in_bounds] += 1j*self._interp_im(pts[in_bounds])
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment