api: numpify sinc precompute and add test

devito/operations/interpolators.py

@@ -2,21 +2,18 @@
 from functools import wraps
 
 import sympy
+import numpy as np
 from cached_property import cached_property
 
 from devito.finite_differences.differentiable import Mul
-<<<<<<< HEAD
 from devito.finite_differences.elementary import floor
-=======
-from devito.finite_differences.elementary import floor, sqrt
->>>>>>> 5815a494d (api: support sinc interpolation)
 from devito.symbolics import retrieve_function_carriers, retrieve_functions, INT
 from devito.tools import as_tuple, flatten, filter_ordered
 from devito.types import (ConditionalDimension, Eq, Inc, Evaluable, Symbol,
                           CustomDimension, SubFunction)
 from devito.types.utils import DimensionTuple
 
-__all__ = ['LinearInterpolator', 'PrecomputedInterpolator']
+__all__ = ['LinearInterpolator', 'PrecomputedInterpolator', 'SincInterpolator']
 
 
 def check_radius(func):
@@ -129,6 +126,8 @@ class GenericInterpolator(ABC):
     Abstract base class defining the interface for an interpolator.
     """
 
+    _name = "generic"
+
     @abstractmethod
     def inject(self, *args, **kwargs):
         pass
@@ -137,6 +136,10 @@ def inject(self, *args, **kwargs):
     def interpolate(self, *args, **kwargs):
         pass
 
+    @property
+    def name(self):
+        return self._name
+
     def _arg_defaults(self, **args):
         return {}
 
@@ -149,6 +152,8 @@ class WeightedInterpolator(GenericInterpolator):
     and multiplied at a given point: `w[x, y] = wx[x] * wy[y]`
     """
 
+    _name = 'weighted'
+
     def __init__(self, sfunction):
         self.sfunction = sfunction
 
@@ -377,6 +382,9 @@ class LinearInterpolator(WeightedInterpolator):
     ----------
     sfunction: The SparseFunction that this Interpolator operates on.
     """
+
+    _name = 'linear'
+
     @property
     def _weights(self):
         c = [(1 - p) * (1 - r) + p * r
@@ -410,6 +418,8 @@ class PrecomputedInterpolator(WeightedInterpolator):
     sfunction: The SparseFunction that this Interpolator operates on.
     """
 
+    _name = 'precomp'
+
     def _positions(self, implicit_dims):
         if self.sfunction.gridpoints is None:
             return super()._positions(implicit_dims)
@@ -441,8 +451,10 @@ class SincInterpolator(PrecomputedInterpolator):
 
     """
 
+    _name = 'sinc'
+
     # Table 1
-    _b_table = {1: 0.0, 2: 1.84, 3: 3.04,
+    _b_table = {2: 2.94, 3: 4.53,
                 4: 4.14, 5: 5.26, 6: 6.40,
                 7: 7.51, 8: 8.56, 9: 9.56, 10: 10.64}
 
@@ -473,13 +485,14 @@ def _arg_defaults(self, coords=None, sfunc=None):
             raise ValueError("No coordinates or sparse function provided")
         # Coords to indices
         coords = (coords - np.array(sfunc.grid.origin)) / np.array(sfunc.grid.spacing)
-        coords = np.floor(coords) - coords + 1 - self.r
+        coords = coords - np.floor(coords)
+
         # Precompute sinc
         for j, r in enumerate(self._rdim):
             data = np.zeros((coords.shape[0], 2*self.r), dtype=sfunc.dtype)
             for ri in range(2*self.r):
-                rpos = ri + coords[:, j]
-                num = np.i0(b*sqrt(1 - (rpos/self.r)**2))
+                rpos = ri - self.r + 1 - coords[:, j]
+                num = np.i0(b*np.sqrt(1 - (rpos/self.r)**2))
                 data[:, ri] = num / b0 * np.sinc(rpos)
             args[self.interpolation_coeffs[r].name] = data
 

devito/types/sparse.py

@@ -31,7 +31,7 @@
 
 
 _interpolators = {'linear': LinearInterpolator, 'sinc': SincInterpolator}
-_default_radius = {'linear': 1, 'sinc': 2}
+_default_radius = {'linear': 1, 'sinc': 4}
 
 
 class AbstractSparseFunction(DiscreteFunction):
@@ -808,14 +808,17 @@ class SparseFunction(AbstractSparseFunction):
 
     _sub_functions = ('coordinates',)
 
-    __rkwargs__ = AbstractSparseFunction.__rkwargs__ + ('coordinates', 'interpolator')
+    __rkwargs__ = AbstractSparseFunction.__rkwargs__ + ('coordinates', 'interpolation')
 
     def __init_finalize__(self, *args, **kwargs):
         super().__init_finalize__(*args, **kwargs)
 
-        interp = kwargs.get('interpolator', 'linear')
+        interp = kwargs.get('interpolation', 'linear')
+        self.interpolation = interp
         self.interpolator = _interpolators[interp](self)
         self._radius = kwargs.get('r', _default_radius[interp])
+        if interp == 'sinc' and self._radius < 2:
+            raise ValueError("The 'sinc' interpolator requires a radius of at least 2")
 
         # Set up sparse point coordinates
         coordinates = kwargs.get('coordinates', kwargs.get('coordinates_data'))
@@ -838,7 +841,7 @@ def _arg_defaults(self, alias=None):
         defaults = super()._arg_defaults(alias=alias)
 
         key = alias or self
-        coords = 
+        coords = defaults.get(key.coordinates.name, key.coordinates.data)
         defaults.update(key.interpolator._arg_defaults(coords=coords,
                                                        sfunc=key))
         return defaults

examples/seismic/acoustic/acoustic_example.py

@@ -4,7 +4,7 @@
 from devito.logger import info
 from devito import Constant, Function, smooth, norm
 from examples.seismic.acoustic import AcousticWaveSolver
-from examples.seismic import demo_model, setup_geometry, seismic_args, Receiver
+from examples.seismic import demo_model, setup_geometry, seismic_args
 
 
 def acoustic_setup(shape=(50, 50, 50), spacing=(15.0, 15.0, 15.0),
@@ -15,7 +15,7 @@ def acoustic_setup(shape=(50, 50, 50), spacing=(15.0, 15.0, 15.0),
                        fs=fs, **kwargs)
 
     # Source and receiver geometries
-    geometry = setup_geometry(model, tn)
+    geometry = setup_geometry(model, tn, **kwargs)
 
     # Create solver object to provide relevant operators
     solver = AcousticWaveSolver(model, geometry, kernel=kernel,
@@ -65,16 +65,19 @@ def run(shape=(50, 50, 50), spacing=(20.0, 20.0, 20.0), tn=1000.0,
 
 @pytest.mark.parametrize('shape', [(101,), (51, 51), (16, 16, 16)])
 @pytest.mark.parametrize('k', ['OT2', 'OT4'])
-def test_isoacoustic_stability(shape, k):
+@pytest.mark.parametrize('interp', ['linear', 'sinc'])
+def test_isoacoustic_stability(shape, k, interp):
     spacing = tuple([20]*len(shape))
-    _, _, _, [rec, _] = run(shape=shape, spacing=spacing, tn=20000.0, nbl=0, kernel=k)
+    _, _, _, [rec, _] = run(shape=shape, spacing=spacing, tn=20000.0, nbl=0, kernel=k,
+                            interpolation=interp)
     assert np.isfinite(norm(rec))
 
 
 @pytest.mark.parametrize('fs, normrec, dtype', [(True, 369.955, np.float32),
                                                 (False, 459.1678, np.float64)])
-def test_isoacoustic(fs, normrec, dtype):
-    _, _, _, [rec, _] = run(fs=fs, dtype=dtype)
+@pytest.mark.parametrize('interp', ['linear', 'sinc'])
+def test_isoacoustic(fs, normrec, dtype, interp):
+    _, _, _, [rec, _] = run(fs=fs, dtype=dtype, interpolation=interp)
     assert np.isclose(norm(rec), normrec, rtol=1e-3, atol=0)
 
 

examples/seismic/acoustic/operators.py

@@ -1,6 +1,5 @@
 from devito import Eq, Operator, Function, TimeFunction, Inc, solve, sign
 from devito.symbolics import retrieve_functions, INT
-from examples.seismic import PointSource, Receiver
 
 
 def freesurface(model, eq):
@@ -119,11 +118,8 @@ def ForwardOperator(model, geometry, space_order=4,
     u = TimeFunction(name='u', grid=model.grid,
                      save=geometry.nt if save else None,
                      time_order=2, space_order=space_order)
-    src = PointSource(name='src', grid=geometry.grid, time_range=geometry.time_axis,
-                      npoint=geometry.nsrc)
-
-    rec = Receiver(name='rec', grid=geometry.grid, time_range=geometry.time_axis,
-                   npoint=geometry.nrec)
+    src = geometry.src
+    rec = geometry.rec
 
     s = model.grid.stepping_dim.spacing
     eqn = iso_stencil(u, model, kernel)
@@ -160,10 +156,8 @@ def AdjointOperator(model, geometry, space_order=4,
 
     v = TimeFunction(name='v', grid=model.grid, save=None,
                      time_order=2, space_order=space_order)
-    srca = PointSource(name='srca', grid=model.grid, time_range=geometry.time_axis,
-                       npoint=geometry.nsrc)
-    rec = Receiver(name='rec', grid=model.grid, time_range=geometry.time_axis,
-                   npoint=geometry.nrec)
+    srca = geometry.new_src(name='srca', src_type=None)
+    rec = geometry.rec
 
     s = model.grid.stepping_dim.spacing
     eqn = iso_stencil(v, model, kernel, forward=False)
@@ -206,8 +200,7 @@ def GradientOperator(model, geometry, space_order=4, save=True,
                      else None, time_order=2, space_order=space_order)
     v = TimeFunction(name='v', grid=model.grid, save=None,
                      time_order=2, space_order=space_order)
-    rec = Receiver(name='rec', grid=model.grid, time_range=geometry.time_axis,
-                   npoint=geometry.nrec)
+    rec = geometry.rec
 
     s = model.grid.stepping_dim.spacing
     eqn = iso_stencil(v, model, kernel, forward=False)
@@ -244,11 +237,8 @@ def BornOperator(model, geometry, space_order=4,
     m = model.m
 
     # Create source and receiver symbols
-    src = Receiver(name='src', grid=model.grid, time_range=geometry.time_axis,
-                   npoint=geometry.nsrc)
-
-    rec = Receiver(name='rec', grid=model.grid, time_range=geometry.time_axis,
-                   npoint=geometry.nrec)
+    src = geometry.src
+    rec = geometry.rec
 
     # Create wavefields and a dm field
     u = TimeFunction(name="u", grid=model.grid, save=None,

examples/seismic/source.py

@@ -117,8 +117,6 @@ def __args_setup__(cls, *args, **kwargs):
                 raise TypeError("Need either `npoint` or `coordinates`")
             kwargs['npoint'] = coordinates.shape[0]
 
-        kwargs.setdefault('r', 1)
-        kwargs.setdefault('interpolator', 'linear')
         return args, kwargs
 
     def __init_finalize__(self, *args, **kwargs):

examples/seismic/utils.py

@@ -3,13 +3,14 @@
 
 from devito import error, configuration, warning
 from devito.tools import Pickable
+from devito.types.sparse import _default_radius
 
 from .source import *
 
 __all__ = ['AcquisitionGeometry', 'setup_geometry', 'seismic_args']
 
 
-def setup_geometry(model, tn, f0=0.010):
+def setup_geometry(model, tn, f0=0.010, interpolation='linear', **kwargs):
     # Source and receiver geometries
     src_coordinates = np.empty((1, model.dim))
     src_coordinates[0, :] = np.array(model.domain_size) * .5
@@ -18,8 +19,10 @@ def setup_geometry(model, tn, f0=0.010):
 
     rec_coordinates = setup_rec_coords(model)
 
+    r = kwargs.get('r', _default_radius[interpolation])
     geometry = AcquisitionGeometry(model, rec_coordinates, src_coordinates,
-                                   t0=0.0, tn=tn, src_type='Ricker', f0=f0)
+                                   t0=0.0, tn=tn, src_type='Ricker', f0=f0,
+                                   interpolation=interpolation, r=r)
 
     return geometry
 
@@ -58,7 +61,7 @@ class AcquisitionGeometry(Pickable):
     """
 
     __rargs__ = ('grid', 'rec_positions', 'src_positions', 't0', 'tn')
-    __rkwargs__ = ('f0', 'src_type')
+    __rkwargs__ = ('f0', 'src_type', 'interpolation', 'r')
 
     def __init__(self, model, rec_positions, src_positions, t0, tn, **kwargs):
         """
@@ -85,6 +88,8 @@ def __init__(self, model, rec_positions, src_positions, t0, tn, **kwargs):
         self._dt = model.critical_dt
         self._t0 = t0
         self._tn = tn
+        self.interpolation = kwargs.get('interpolation', 'linear')
+        self._r = kwargs.get('r', _default_radius[self.interpolation])
 
     def resample(self, dt):
         self._dt = dt
@@ -140,14 +145,23 @@ def nsrc(self):
     def dtype(self):
         return self.grid.dtype
 
+    @property
+    def r(self):
+        return self._r
+
+    @property
+    def interpolation(self):
+        return self._interpolation
+
     @property
     def rec(self):
         return self.new_rec()
 
     def new_rec(self, name='rec'):
         return Receiver(name=name, grid=self.grid,
                         time_range=self.time_axis, npoint=self.nrec,
-                        coordinates=self.rec_positions)
+                        coordinates=self.rec_positions,
+                        interpolation=self.interpolation, r=self._r)
 
     @property
     def adj_src(self):
@@ -157,7 +171,8 @@ def adj_src(self):
         adj_src = sources[self.src_type](name='rec', grid=self.grid, f0=self.f0,
                                          time_range=self.time_axis, npoint=self.nrec,
                                          coordinates=self.rec_positions,
-                                         t0=self._t0w, a=self._a)
+                                         t0=self._t0w, a=self._a,
+                                         interpolation=self.interpolation, r=self._r)
         # Revert time axis to have a proper shot record and not compute on zeros
         for i in range(self.nrec):
             adj_src.data[:, i] = adj_src.wavelet[::-1]
@@ -172,12 +187,14 @@ def new_src(self, name='src', src_type='self'):
             warning("No source type defined, returning uninitiallized (zero) source")
             return PointSource(name=name, grid=self.grid,
                                time_range=self.time_axis, npoint=self.nsrc,
-                               coordinates=self.src_positions)
+                               coordinates=self.src_positions,
+                               interpolation=self.interpolation, r=self._r)
         else:
             return sources[self.src_type](name=name, grid=self.grid, f0=self.f0,
                                           time_range=self.time_axis, npoint=self.nsrc,
                                           coordinates=self.src_positions,
-                                          t0=self._t0w, a=self._a)
+                                          t0=self._t0w, a=self._a,
+                                          interpolation=self.interpolation, r=self._r)
 
 
 sources = {'Wavelet': WaveletSource, 'Ricker': RickerSource, 'Gabor': GaborSource}