dsl: Implementation of new interpolation methods #1932

guaacoelho · 2022-06-08T16:19:40Z

Hi Guys @FabioLuporini @mloubout @navjotk @georgebisbas ,

I've been working with @nogueirapeterson to include new interpolation techniques on Devito. We implement the cubic and sinc interpolation methods. Thus, we create two new classes: CubicInterpolator and SincInterpolator, which contains, respectively, the cubic and sinc interpolation methods.

The cubic interpolator was implemented using the cubic convolution technique (10.1109/TASSP.1981.1163711) and the sinc interpolator was implemented using Kaiser windowed sinc functions (https://doi.org/10.1190/1.1451454).

We hope you appreciate this contribution and help us to improve this implementation with your great review!

…cInterpolator class): - Implementation of _interpolation_indices method - Creation of relative_position_map, in sparse.py - Call to interpolation_indices method inside interpolation's callback

…ic interpolator

…ng 2D and 3D cubic interpolation

…zation with the sinc interpolator

…ng 2D and 3D sinc interpolation - Add the lookuptable.h file to the operator's default import list

mloubout

Great! Thanks a lot. I'll have a thorough walk-through.

So SincInterpolation is not testable, right. It would be good to make it testable in some way. Also, there is a hardcoded path to lookuptable.h that will break the rest of the code.

mloubout

Ok made a first pass. I think it's quite nice and provides some pretty cool features. I left some comments and I'm happy to discuss if not quite clear.

mloubout · 2022-06-08T18:02:57Z

devito/types/sparse.py

+        mathematically equivalent expanded form `coord/spacing -
+        origin/spacing`. This particular form is problematic when a sparse
+        point is in close proximity of the grid origin, since due to a larger
+        machine precision error it may cause a +-1 error in the computation of


That +-1 used to be an error and we made changes to the API to make sure it would not happen anymore. Can you point us to a small MFE with the breaking change.

Actually this was my mistake, I took the same comments from position_map function and forgot to update after the rebase. So there is no error in here

mloubout · 2022-06-08T18:04:12Z

devito/types/sparse.py

+            return tuple(product(range(-1, 3), repeat=self.grid.dim))
+        elif self.sinc:
+            r = self.interpolator.r
+            return tuple(product(range(-r + 1, r+1), repeat=self.grid.dim))


Isn't this case true for all r? I.e for cubic you have r=2 which is range(-1, 3) and for the linear case r=1 whci his range(0, 2)

You're right! range(-r + 1, r+1) works for all r.

I'll change this

So in the end we won't have this conditional, I hope?

I updated the PR with this fix.

mloubout · 2022-06-08T18:04:28Z

tests/test_interpolation_cubic.py

+                    Function, TimeFunction)
+
+
+def unit_box(name='a', shape=(11, 11), grid=None):


Redefinition of existing fixture.
same below

mloubout · 2022-06-08T18:05:39Z

tests/test_interpolation_cubic.py

+    ((101, 101), [(3.4, 95.), (3., 92)]),
+    ((101, 101, 101), [(3.4, 95.), (3., 92), (10., 10.)])
+])
+def test_interpolate(shape, coords, npoints=20):


Can this all be gathered under a common test class?

Also I don't think it needs a new file, test_interpolation.py exists for this

mloubout · 2022-06-08T18:06:06Z

tests/test_interpolation_cubic.py

+    # print(p.coordinates.data)
+    op(a=a)
+
+    # Precisão 3D nao bate em 10-6 mas sim em 10-5


There is leftover Portugues throughout

mloubout · 2022-06-08T18:08:20Z

devito/operations/interpolators.py

+
+        idx_subs = []
+        for i, idx in enumerate(index_matrix):
+            # Introduce ConditionalDimension so that we don't go OOB


This should only be necessary if r > variables.space_order

So if r <= variables.space_order there is no need to use ConditionalDimension in this case?

I adjust this part of the code and updated here.

You'll see that I added a conditionalif hasattr(x, 'space_order') else False. I needed to do this because in some cases variables does not have the space_order attribute

mloubout · 2022-06-08T18:12:16Z

devito/operations/interpolators.py

+            lhs = self.sfunction.subs(self_subs)
+            last = [Inc(lhs, rhs)] if increment else [Eq(lhs, rhs)]
+
+            # Creates the symbolic equation that calls the populate function


What's this?

To implement the sinc interpolation I used a look up table that holds the values of the coefficients. Before it could be accessed, it was necessary to initialize this table with the values. This symbolic Function populate() creates a call to the function responsible for initialize the look up table.

This function is defined inside look up table.h, that I apparently forgot to add here. I'm attaching the lookuptable.h in this comment, this way you can see how it works

lookuptable.zip

mloubout · 2022-06-08T18:14:06Z

devito/operations/interpolators.py

+            dim_pos = self.sfunction.grid.dimensions
+
+            # Checks if the data is 3D or 2D
+            if len(dim_pos) == 3:


this can be setup at __init_ with self._eqs = self.bicubic if self.sfunction.grid.dim = 2 else ...

...also lets try to setup as much at init as possible?

mloubout · 2022-06-08T18:15:05Z

devito/operations/interpolators.py

+                                                                   idx_subs,
+                                                                   dim_pos=dim_pos)]
+            else:
+                eqs = [np.sum(v) for v in self._bicubic_equations(_expr,


why doesn't ..._equations directly return the sum?
Also why np.sum? Isn't it still partially symbolic?

Not returning directly the sum, the return of ..._equations can be used for both injection and interpolation. In interpolation it sum the values and in injection just assign the values for its respective grid point.

About the use of numpy: Yeah, it is still a symbolic. I wanted to generate a symbolic equation that sum all the individual values, that's why i used np.sum. Is there some problem in using it?

Ok makes sense for sum. But should not use numpy as it is meant for numerical values. Should use Devito.Add(*v)

Does this actually work? Numpy summing over sympy symbols?

Yes, it works summing over sympy symbols. But I will take a look at Devito.Add(*v) to use it

I updated the code and it is using Add(*v) now

mloubout · 2022-06-08T18:16:40Z

devito/operations/interpolators.py

+                                                          field_offset=field_offset)
+
+            rhs = Symbol(name='sum', dtype=self.sfunction.dtype)
+            summands = [Eq(rhs, 0., implicit_dims=self.sfunction.dimensions)]


everything till here is the same as in LinearInterpolation right? and after the bicubic_equation too. So may be nice to gather those under a PolynomialInterpolation with order=1 or order=3 that specialize on the equations only since it's the only thing that changes

I added a PolynomialInterpolator class that works as a parent class. It's not exactly what you just said, but I think it solves the code reuse problem

FabioLuporini · 2022-06-09T08:23:56Z

devito/operations/interpolators.py

+class CubicInterpolator(GenericInterpolator):
+    """
+    Concrete implementation of GenericInterpolator implementing a cubic interpolation
+    scheme.


FabioLuporini · 2022-06-09T08:24:13Z

devito/operations/interpolators.py

+    scheme.
+    Parameters
+    ----------
+    sfunction: The SparseFunction that this Interpolator operates on.


this is not legal numpydoc format

Indeed this is not legal, but I guess it is our "fault" that we have not numpy-doc-ized this file.... i.e. look here:
https://github.com/devitocodes/devito/pull/1932/files/c068956a657999c10d1d977347f0f52f40132f1f#diff-2c4e29d657bce68c3dc96bb61bed3d1038b983ab72845761bf2ead45e30e222aR25

FabioLuporini · 2022-06-09T08:25:07Z

devito/operations/interpolators.py

+
+        return eqs
+
+    def _interpolation_indices(self, variables, offset=0, field_offset=0):


isn't there anything we could factorize in a common superclass?

isn't there anything we could factorize in a common superclass?

Yes, I think it has. I'm working on it

I updated the code by adding the superclass you suggested

FabioLuporini · 2022-06-09T08:26:01Z

devito/operations/interpolators.py

+
+    def __init__(self, sfunction):
+        self.sfunction = sfunction
+        self.r = 4


This should probably be a class property

Is it a problem to let the user choose?

Is it a problem to let the user choose?

The Bessel function was implemented inside the .h file considering r = 4. So it would be a problem to allow the user to choose the value of r

FabioLuporini

Thank you for this contribution. Looking great! I just left a bunch of comments

FabioLuporini · 2022-06-09T08:27:25Z

tests/test_interpolation_cubic.py

@@ -0,0 +1,269 @@
+import numpy as np


IMHO this file should be squashed with test_interpolation, and as much code and existing tests as possible should be reused, through pytest.parametrize for example

EdCaunt · 2022-06-09T13:04:54Z

tests/test_interpolation_cubic.py

+    ((101, 101, 101), [(3.4, 95.), (3., 92), (10., 10.)])
+])
+def test_interpolate(shape, coords, npoints=20):
+    """Test generic point interpolation testing the x-coordinate of an


Opening triple quotes want to be on line above for consistency with the rest of the codebase, i.e.

""" Test generic... """

EdCaunt · 2022-06-09T13:05:45Z

devito/operations/interpolators.py

+        ----------
+        expr : expr-like
+            Input expression to interpolate.
+        position: List


Should be list (all lowercase)

EdCaunt · 2022-06-09T13:09:50Z

devito/operations/interpolators.py

+               for ii in range(4)]
+        return eqs
+
+    def _tricubic_equations(self, expr, idx_subs, dim_pos=None):


Could this (and _bicubic_equations) be combined into some _ncubic_equations for N-dimensional interpolation?

Yes, it could. I changed the code and made this change you suggested

EdCaunt · 2022-06-09T13:13:32Z

devito/operations/interpolators.py

+        idx_subs: dict
+            Structure responsible for mapping the order of substitution of dimensions
+            in expr.
+        idx_2d: Integer


Should be int I think?

EdCaunt · 2022-06-09T13:15:51Z

devito/operations/interpolators.py

+
+        return eqs
+
+    def _sinc_equations3D(self, expr, idx_subs, dim_pos=None):


Again, can these be consolidated into _sinc_equationsND()? Seems like a decent amount of lines are the same and would be nice to have a generalization to N-D

EdCaunt · 2022-06-09T13:19:30Z

tests/test_interpolation_cubic.py

@@ -0,0 +1,269 @@
+import numpy as np


Would it be worth having a test to check the convergence of the cubic interpolator? Not sure what convergence would be expected for the sinc one

Does this seem like a decent idea @mloubout ?

It isn't as I pointed out above

Maybe I'm just being stupid and missing it, but where?

#1932 (comment)

Sorry completely missed what you were talking about only saw the numpy.

Convergence may be tricky, it's usually easier to just do an analytical check where you set the field to polynomial values and check the interpolation is exact.

navjotk

Left some minor comments.

Also have two higher-level questions:

Is SincInterpolator the best name for what this is? (Sinc is only a part of the interpolation function)
How does this approach compare to using the machinery provided in PrecomputedInterpolator, i.e. the Sinc class could inherit from that, precompute the coefficients (which I think you're doing in some way anyway), and hand them off to PrecomputedInterpolator for the rest?

navjotk · 2022-06-10T11:08:16Z

devito/operations/interpolators.py

+            option = 1, ponto imediatamente vizinho ao ponto de intepolação.
+            Onde 0 < |arg| < 1
+            option = 2, pontos extremos da vizinhança.
+            Onde 1 < |arg| < 2


Am I seeing an old version? We probably want our documentation in English :)

navjotk · 2022-06-10T11:10:32Z

devito/operations/interpolators.py

+            dim_pos = self.sfunction.grid.dimensions
+
+            # Checks if the data is 3D or 2D
+            if len(dim_pos) == 3:


...also lets try to setup as much at init as possible?

navjotk · 2022-06-10T11:11:23Z

devito/operations/interpolators.py

+                                                                   idx_subs,
+                                                                   dim_pos=dim_pos)]
+            else:
+                eqs = [np.sum(v) for v in self._bicubic_equations(_expr,


Does this actually work? Numpy summing over sympy symbols?

navjotk · 2022-06-10T11:12:24Z

devito/operations/interpolators.py

+                eqs = [v for v in self._tricubic_equations(_expr,
+                                                           idx_subs,
+                                                           dim_pos=dim_pos)]
+            else:


Can we do this with inheritance? Or in a way that generalises beyond 2/3 dimensions?

navjotk · 2022-06-10T11:12:57Z

devito/operations/interpolators.py

+
+    def __init__(self, sfunction):
+        self.sfunction = sfunction
+        self.r = 4


Is it a problem to let the user choose?

navjotk · 2022-06-10T11:16:26Z

devito/operator/operator.py

@@ -133,7 +134,8 @@ class Operator(Callable):
    """

    _default_headers = [('_POSIX_C_SOURCE', '200809L')]
-    _default_includes = ['stdlib.h', 'math.h', 'sys/time.h']
+    _default_includes = ['stdlib.h', 'math.h', 'sys/time.h',
+                         f'{os.path.expanduser("~")}/lookuptable.h']


I need help understanding what this header file is doing. Also, is this a static header file that should become a standard part of Devito? Alternately, can this be generated using the JIT infra?

To implement the sinc interpolation I used a look up table that holds the values of the coefficients. So this header contains the look up table and the methods responsible for initialize and get the values from it.

Access to these functions is defined from the definition of symbolic functions during the generation of symbolic equations.

I updated the PR adding the lookuptable.h file, this way you can take a look if you want to.

About using JIT to generate the code, I don't know how it could be using JIT but I will take a look.

And yes, the static header file should become a standard part of Devito.

navjotk · 2022-06-10T11:19:07Z

devito/types/sparse.py

+            return tuple(product(range(-1, 3), repeat=self.grid.dim))
+        elif self.sinc:
+            r = self.interpolator.r
+            return tuple(product(range(-r + 1, r+1), repeat=self.grid.dim))


So in the end we won't have this conditional, I hope?

…ay the Interpolator initialization can use it

…o it

…r, so it was added a conditional checking this condition

mloubout · 2022-06-14T13:05:02Z

The .h file shouldn't really be necessary. Bessel function are defined in math.h so should be able to make a symbolic expression directly calling it rather than the lookup table

navjotk · 2022-06-14T15:11:34Z

Also, I want to understand if this is significantly different from using PrecomputedInterpolator (i.e. this class could extend and use the functionality there)

guaacoelho · 2022-06-14T19:45:26Z

The .h file shouldn't really be necessary. Bessel function are defined in math.h so should be able to make a symbolic expression directly calling it rather than the lookup table

@mloubout Before the decision of create this '.h' file it was made a search for some library that could already has the implementation of the Bessel function that we needed (the zero-order modified Bessel function of the first kind).

I found some implementations of the bessel functions, but none of the zero-order modified Bessel function of the first kind.

guaacoelho · 2022-06-14T20:08:41Z

Also, I want to understand if this is significantly different from using PrecomputedInterpolator (i.e. this class could extend and use the functionality there)

@navjotk At the same way that linear interpolation, both cubic and sinc interpolation also can be done by using PrecomputedInterpolator, if the right coefficients are used. So I can say that the main difference is that this new classes removes the responsibility of obtaining the interpolation coefficients from the user.

navjotk · 2022-06-23T12:43:09Z

@guaacoelho But LinearInterpolator does not do any pre computation. If this class is doing pre computation of coefficients, it could store them in the properties provided by PrecomputedInterpolator and use the rest of the infrastructure there. i.e. This class precomputes the coefficients and hands them off to PrecomputedInterpolator, not the user.

…e between the interpolators classes

georgebisbas

Thanks for this contribution!

georgebisbas · 2022-07-18T07:50:24Z

devito/operations/interpolators.py

@@ -1,5 +1,6 @@
 from abc import ABC, abstractmethod
-
+# Polly interpolator é a que contem as duas interpolações


Translate and avoid comment in this area

georgebisbas · 2022-07-18T07:54:40Z

devito/operations/interpolators.py


    def interpolate(self, expr, offset=0, increment=False, self_subs={}):
        """
        Generate equations interpolating an arbitrary expression into ``self``.
-


add line back

georgebisbas · 2022-07-18T07:54:58Z

devito/operations/interpolators.py

+            variables = list(retrieve_function_carriers(_expr))
+
+            # Need to get origin of the field in case it is staggered
+            # TODO: handle each variable staggereing spearately


georgebisbas · 2022-07-18T07:57:53Z

devito/operator/operator.py

@@ -133,7 +134,8 @@ class Operator(Callable):
    """

    _default_headers = [('_POSIX_C_SOURCE', '200809L')]
-    _default_includes = ['stdlib.h', 'math.h', 'sys/time.h']
+    _default_includes = ['stdlib.h', 'math.h', 'sys/time.h',
+                         f'{os.path.expanduser("~/devito")}/lookuptable.h']


Ι think this .h is not something we encourage

mloubout · 2023-08-03T13:57:02Z

FYI sparse interpolation is beeing revamp in #2128 which should help updae this PR.

mloubout · 2024-04-30T19:13:59Z

sinc interpolation added in #2342 , closing.

Feel free to open a new PR for cubic with the new API.

guaacoelho added 12 commits June 8, 2022 10:49

dsl: Add CubicInterpolator class with some of your main methods

0bf8d43

dsl: Partial implementation of interpolations's callback method (Cubi…

94f6bc1

…cInterpolator class): - Implementation of _interpolation_indices method - Creation of relative_position_map, in sparse.py - Call to interpolation_indices method inside interpolation's callback

dsl: Add the possibility to initialize a sparse function with the cub…

feab04a

…ic interpolator

dsl: Add methods that create the equations responsible for implementi…

7ff95ed

…ng 2D and 3D cubic interpolation

dsl: Add implementation of injection method

a9fdea1

dsl: Add SincInterpolator class with some of your main methods

f4b5ea5

dsl: Alterations in SparseFunction creation, which allow the initiali…

ca59ffc

…zation with the sinc interpolator

dsl: Add methods that create the equations responsible for implementi…

35099a9

…ng 2D and 3D sinc interpolation - Add the lookuptable.h file to the operator's default import list

dsl: Add implementation of injection method (SincInterpolator)

4b089b2

misc: Fix flake8 identation problems

964820b

tests: Add the cubic interpolation test script

872cbec

dsl: Changes the _radius attribute to allow MPI execution

c068956

guaacoelho changed the title ~~Interpolation methods~~ dsl: Implementation of new interpolation methods Jun 8, 2022

mloubout added the API api (symbolics, types, ...) label Jun 8, 2022

mloubout reviewed Jun 8, 2022

View reviewed changes

guaacoelho added 2 commits June 8, 2022 16:02

misc: Removes some comments that not supposed to be in the code

d6684e8

dsl: Changes in the code to make it more readable and clean

72736cc

FabioLuporini reviewed Jun 9, 2022

View reviewed changes

FabioLuporini requested changes Jun 9, 2022

View reviewed changes

EdCaunt reviewed Jun 9, 2022

View reviewed changes

guaacoelho added 2 commits June 9, 2022 15:41

misc: Updating comments

e0d9308

misc: Updating comments

3dd2817

navjotk reviewed Jun 10, 2022

View reviewed changes

guaacoelho added 3 commits June 10, 2022 18:05

dsl: Defines self._radius before the interpolator creation. In this w…

db65094

…ay the Interpolator initialization can use it

dsl: Adds the lookuptable header inside Devito and updates the path t…

60c9c89

…o it

dsl: Replacing np.sum with Sympy.Add

468b403

guaacoelho added 2 commits June 13, 2022 15:03

misc: Fix flake8 identation problems

f584986

dsl: Need conditional dimension only if radius > variables.space_orde…

ad98770

…r, so it was added a conditional checking this condition

dsl: Makes r a class property

42325fc

guaacoelho added 2 commits June 28, 2022 16:23

dsl: Add a parent class (PolynomialInterpolator) providing inheritanc…

7630861

…e between the interpolators classes

misc: Small fixes

08439d5

georgebisbas reviewed Jul 18, 2022

View reviewed changes

misc: Removes some comments that not supposed to be in the code

98835de

mloubout closed this Apr 30, 2024

		Function, TimeFunction)


		def unit_box(name='a', shape=(11, 11), grid=None):


		return eqs

		def _interpolation_indices(self, variables, offset=0, field_offset=0):


		return eqs

		def _sinc_equations3D(self, expr, idx_subs, dim_pos=None):

		@@ -1,5 +1,6 @@
		from abc import ABC, abstractmethod

		# Polly interpolator é a que contem as duas interpolações

dsl: Implementation of new interpolation methods #1932

dsl: Implementation of new interpolation methods #1932

Conversation

guaacoelho commented Jun 8, 2022 • edited Loading

mloubout left a comment

Choose a reason for hiding this comment

mloubout left a 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

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

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

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

navjotk left a 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 Jun 14, 2022

navjotk commented Jun 14, 2022

guaacoelho commented Jun 14, 2022

guaacoelho commented Jun 14, 2022 • edited Loading

navjotk commented Jun 23, 2022

georgebisbas left a 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 Aug 3, 2023

guaacoelho commented Jun 8, 2022 •

edited

Loading

guaacoelho commented Jun 14, 2022 •

edited

Loading