api: Add Hicks (sinc) interpolation api #2342

mloubout · 2024-04-01T16:57:46Z

Add support for kaiser-windowed sinc interpolation

Fixes #595

codecov · 2024-04-01T17:26:24Z

Codecov Report

Attention: Patch coverage is 85.71429% with 28 lines in your changes are missing coverage. Please review.

Project coverage is 79.41%. Comparing base (311cd4e) to head (353af4f).

Files	Patch %	Lines
conftest.py	27.58%	19 Missing and 2 partials ⚠️
devito/operations/interpolators.py	90.74%	4 Missing and 1 partial ⚠️
examples/seismic/acoustic/acoustic_example.py	87.50%	1 Missing ⚠️
examples/seismic/utils.py	93.33%	1 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##           master    #2342      +/-   ##
==========================================
- Coverage   86.66%   79.41%   -7.26%     
==========================================
  Files         232      232              
  Lines       43445    43578     +133     
  Branches     8063     8072       +9     
==========================================
- Hits        37653    34606    -3047     
- Misses       5087     8215    +3128     
- Partials      705      757      +52

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

EdCaunt · 2024-04-02T15:56:59Z

devito/operations/interpolators.py

+
+    def _arg_defaults(self, coords=None, sfunc=None):
+        args = {}
+        b = self._b_table[self.r]


Does this need a catch for r > 10?

It will just throw a key-error but can add a check (I doubt anyone will run it this high that's a 20x20x20 cube in 3d super expensive)

agree to add a check. Users can be surprising...

FabioLuporini

Can we also add an ABox test in PRO?

FabioLuporini · 2024-04-02T07:06:21Z

.github/workflows/pytest-core-mpi.yml

@@ -48,6 +48,11 @@ jobs:
        python3 scripts/clear_devito_cache.py
        python3 -m pytest --cov --cov-config=.coveragerc --cov-report=xml -m parallel tests/

+    - name: Test examples


for uniformity, "Test examples with MPI" ?

FabioLuporini · 2024-04-02T07:08:54Z

devito/operations/interpolators.py

+    def name(self):
+        return self._name
+
+    def _arg_defaults(self, **args):


if we want this method here it would be cleaner to have this class inhereting from ArgProvider.

Perhaps we move it into a subclass?

No, the interpolator isn't an arg provide, its sparse function is and is calling it only.

FabioLuporini · 2024-04-03T07:21:49Z

devito/operations/interpolators.py

+        # Position only replacement, not radius dependent.
+        # E.g src.inject(vp(x)*src) needs to use vp[posx] at all points
+        # not vp[posx + rx]
+        if pos_only is not None:


if you turn pos_only into pos_only=() (instead of None) you don't need this if I think?

Convention throughout is None for kwargs default (for tuple/dict/...)

no, the convention is to use immutable objects only, and tuple in that sense is fine ; dict/list/set aren't instead obv

for example, I'm using =() in ir/iet/nodes.py

alternatively, you may use as_tuple

FabioLuporini · 2024-04-03T07:22:22Z

devito/operations/interpolators.py

+
+    _name = 'sinc'
+
+    # Table 1


What's Table 2? u mean Table 1 in the paper?

Yes Table 1 in the paper

FabioLuporini · 2024-04-03T07:23:10Z

devito/operations/interpolators.py

+
+    def _arg_defaults(self, coords=None, sfunc=None):
+        args = {}
+        b = self._b_table[self.r]


agree to add a check. Users can be surprising...

FabioLuporini · 2024-04-03T07:23:55Z

devito/operations/interpolators.py

+        coords = coords / np.array(sfunc.grid.spacing)
+        coords = coords - np.floor(coords)
+
+        # Precompute sinc


I think this wants a NOTE: <comment about performance> ? if you measured it's slowish? or is it super quick since the loops are very short?

I haven't seen any drastic slow compute since it's numpified so not sure want to force user to be like "hey here is a warning should I nitpick here"

FabioLuporini · 2024-04-03T07:24:19Z

devito/types/sparse.py

+        if interp == 'sinc' and self._radius < 2:
+            raise ValueError("The 'sinc' interpolator requires a radius of at least 2")
+
+        # Set space ordert to `r` for safety


but not sure why this comment ended up here

FabioLuporini · 2024-04-03T07:25:48Z

examples/seismic/acoustic/acoustic_example.py

@@ -60,21 +60,26 @@ def run(shape=(50, 50, 50), spacing=(20.0, 20.0, 20.0), tn=1000.0,
    solver.jacobian(dm, autotune=autotune)
    info("Applying Gradient")
    solver.jacobian_adjoint(rec, u, autotune=autotune, checkpointing=checkpointing)
-    return summary.gflopss, summary.oi, summary.timings, [rec, u.data]
+    return summary.gflopss, summary.oi, summary.timings, [rec, u.data._local]


this will pick the (dirty) padding as well I think? why not just data or data_with_halo?

for MPI debug the _local is a lot more useful

FabioLuporini · 2024-04-03T07:27:14Z

examples/seismic/acoustic/operators.py

    z = zfs.parent

-    # Functions present in the stencil
-    funcs = retrieve_functions(rhs)
+    # Retrieve vertical derivatives


what's all these new additions?

it does the same, except it only evaluates z derivatives instead of everything since the rest isn't needed which is much cheaper especially for things like tti.

FabioLuporini · 2024-04-03T07:27:21Z

examples/seismic/acoustic/operators.py

-    # Functions present in the stencil
-    funcs = retrieve_functions(rhs)
+    # Retrieve vertical derivatives
+    Dz = {d for d in retrieve_derivatives(rhs) if z in d.dims}


FabioLuporini · 2024-04-03T14:47:04Z

devito/operations/interpolators.py

+        # Position only replacement, not radius dependent.
+        # E.g src.inject(vp(x)*src) needs to use vp[posx] at all points
+        # not vp[posx + rx]
+        if pos_only is not None:


no, the convention is to use immutable objects only, and tuple in that sense is fine ; dict/list/set aren't instead obv

for example, I'm using =() in ir/iet/nodes.py

alternatively, you may use as_tuple

FabioLuporini · 2024-04-03T14:47:17Z

.github/workflows/pytest-core-mpi.yml

+      - name: Test examples with MPI
+        run: |
+          docker run --rm -e DEVITO_SAFE_MATH=1 -e DEVITO_MPI=1 -e OMP_NUM_THREADS=1 --name examplerun devito_img mpiexec -n 2 pytest examples/seismic/acoustic


comment why safe-math needed

FabioLuporini · 2024-04-03T14:48:12Z

devito/operations/interpolators.py

        # List of indirection indices for all adjacent grid points
-        idx_subs, temps = self._interp_idx(variables, implicit_dims=implicit_dims)
+        idx_subs, temps = self._interp_idx(list(fields), implicit_dims=implicit_dims,


why does it need to be casted to a list? it was a tuple already, so it feels pretty odd

mloubout added API api (symbolics, types, ...) feature-request labels Apr 1, 2024

mloubout force-pushed the sinc-interp-final branch 4 times, most recently from 77bd1b2 to 554aece Compare April 1, 2024 17:19

mloubout force-pushed the sinc-interp-final branch 3 times, most recently from f708213 to 0749d57 Compare April 1, 2024 17:57

mloubout marked this pull request as draft April 1, 2024 19:39

mloubout marked this pull request as ready for review April 1, 2024 19:39

mloubout force-pushed the sinc-interp-final branch 4 times, most recently from 2010409 to ef27943 Compare April 2, 2024 16:42

EdCaunt reviewed Apr 2, 2024

View reviewed changes

mloubout force-pushed the sinc-interp-final branch 5 times, most recently from 82cced0 to 8172aab Compare April 3, 2024 03:26

FabioLuporini reviewed Apr 3, 2024

View reviewed changes

mloubout added 7 commits April 3, 2024 08:40

api: support sinc interpolation

cf4b65e

api: numpify sinc precompute and add test

83bbef3

mpi: revamp conftest

4aadbf8

mpi: fix sinc setup with mpi and add back missing mpi example tests

c4ac44e

api: fix sparse position index for negative location

4db1120

examples: improved fs only check z derivs

99448f1

misc: review cleanup

a83fac7

mloubout force-pushed the sinc-interp-final branch 2 times, most recently from 25721fe to 2705f26 Compare April 3, 2024 13:26

FabioLuporini reviewed Apr 3, 2024

View reviewed changes

api: switch to recommended scipy i0 when available

353af4f

mloubout force-pushed the sinc-interp-final branch from 2705f26 to 353af4f Compare April 3, 2024 14:54

mloubout merged commit 364997c into master Apr 3, 2024
31 checks passed

mloubout deleted the sinc-interp-final branch April 3, 2024 17:08

mloubout mentioned this pull request Apr 30, 2024

dsl: Implementation of new interpolation methods #1932

Closed

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

api: Add Hicks (sinc) interpolation api #2342

api: Add Hicks (sinc) interpolation api #2342

mloubout commented Apr 1, 2024 •

edited

codecov bot commented Apr 1, 2024 •

edited

EdCaunt Apr 2, 2024

mloubout Apr 2, 2024

FabioLuporini Apr 3, 2024

FabioLuporini left a comment

FabioLuporini Apr 2, 2024

FabioLuporini Apr 2, 2024

mloubout Apr 3, 2024 •

edited

FabioLuporini Apr 3, 2024

mloubout Apr 3, 2024

FabioLuporini Apr 3, 2024

FabioLuporini Apr 3, 2024

mloubout Apr 3, 2024

FabioLuporini Apr 3, 2024

FabioLuporini Apr 3, 2024

mloubout Apr 3, 2024

FabioLuporini Apr 3, 2024

FabioLuporini Apr 3, 2024

FabioLuporini Apr 3, 2024

mloubout Apr 3, 2024

FabioLuporini Apr 3, 2024

mloubout Apr 3, 2024

FabioLuporini Apr 3, 2024

FabioLuporini Apr 3, 2024

FabioLuporini Apr 3, 2024

FabioLuporini Apr 3, 2024

api: Add Hicks (sinc) interpolation api #2342

api: Add Hicks (sinc) interpolation api #2342

Conversation

mloubout commented Apr 1, 2024 • edited

codecov bot commented Apr 1, 2024 • edited

Codecov Report

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

FabioLuporini left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

mloubout Apr 3, 2024 • edited

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

mloubout commented Apr 1, 2024 •

edited

codecov bot commented Apr 1, 2024 •

edited

mloubout Apr 3, 2024 •

edited