ConditionalDimension tutorial¶This tutorial explains how to create equations that only get executed under given conditions.
from devito import Dimension, Function, Grid
import numpy as np
# We define a 10x10 grid, dimensions are x, y
shape = (10, 10)
grid = Grid(shape = shape)
x, y = grid.dimensions
# Define function 𝑓. We will initialize f's data with ones on its diagonal.
f = Function(name='f', grid=grid)
f.data[:] = np.eye(10)
f.data
Data([[1., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
[0., 1., 0., 0., 0., 0., 0., 0., 0., 0.],
[0., 0., 1., 0., 0., 0., 0., 0., 0., 0.],
[0., 0., 0., 1., 0., 0., 0., 0., 0., 0.],
[0., 0., 0., 0., 1., 0., 0., 0., 0., 0.],
[0., 0., 0., 0., 0., 1., 0., 0., 0., 0.],
[0., 0., 0., 0., 0., 0., 1., 0., 0., 0.],
[0., 0., 0., 0., 0., 0., 0., 1., 0., 0.],
[0., 0., 0., 0., 0., 0., 0., 0., 1., 0.],
[0., 0., 0., 0., 0., 0., 0., 0., 0., 1.]], dtype=float32)
We begin by constructing an Operator that increases the values of f, a Function defined on a two-dimensional Grid, by one.
#NBVAL_IGNORE_OUTPUT
from devito import Eq, Operator
op0 = Operator(Eq(f, f + 1))
op0.apply()
f.data
Operator `Kernel` ran in 0.01 s
Data([[2., 1., 1., 1., 1., 1., 1., 1., 1., 1.],
[1., 2., 1., 1., 1., 1., 1., 1., 1., 1.],
[1., 1., 2., 1., 1., 1., 1., 1., 1., 1.],
[1., 1., 1., 2., 1., 1., 1., 1., 1., 1.],
[1., 1., 1., 1., 2., 1., 1., 1., 1., 1.],
[1., 1., 1., 1., 1., 2., 1., 1., 1., 1.],
[1., 1., 1., 1., 1., 1., 2., 1., 1., 1.],
[1., 1., 1., 1., 1., 1., 1., 2., 1., 1.],
[1., 1., 1., 1., 1., 1., 1., 1., 2., 1.],
[1., 1., 1., 1., 1., 1., 1., 1., 1., 2.]], dtype=float32)
Every value has been updated by one. You should be able to see twos on the diagonal and ones everywhere else.
#print(op0.ccode) # Print the generated code
In order to construct a ConditionalDimension, one needs to get familiar with relationals. The Devito API for relationals currently supports the following classes of relation, which inherit from their SymPy counterparts:
Relationals are used to define a condition and are passed as an argument to ConditionalDimension at construction time.
A ConditionalDimension defines a non-convex iteration sub-space derived from a parent Dimension, implemented by the compiler generating conditional "if-then" code within the parent Dimension's iteration space.
A ConditionalDimension is used in two typical cases:
To constrain the execution of loop iterations to make sure certain conditions are honoured
To subsample a Function
from devito import ConditionalDimension
print(ConditionalDimension.__doc__)
Symbol defining a non-convex iteration sub-space derived from a ``parent``
Dimension, implemented by the compiler generating conditional "if-then" code
within the parent Dimension's iteration space.
Parameters
----------
name : str
Name of the dimension.
parent : Dimension, optional
The parent Dimension.
factor : int, optional
The number of iterations between two executions of the if-branch. If None
(default), ``condition`` must be provided.
condition : expr-like, optional
An arbitrary SymPy expression, typically involving the ``parent``
Dimension. When it evaluates to True, the if-branch is executed. If None
(default), ``factor`` must be provided.
indirect : bool, optional
If True, use `self`, rather than the parent Dimension, to
index into arrays. A typical use case is when arrays are accessed
indirectly via the ``condition`` expression. Defaults to False.
Examples
--------
Among the other things, ConditionalDimensions are indicated to implement
Function subsampling. In the following example, an Operator evaluates the
Function ``g`` and saves its content into ``f`` every ``factor=4`` iterations.
>>> from devito import Dimension, ConditionalDimension, Function, Eq, Operator
>>> size, factor = 16, 4
>>> i = Dimension(name='i')
>>> ci = ConditionalDimension(name='ci', parent=i, factor=factor)
>>> g = Function(name='g', shape=(size,), dimensions=(i,))
>>> f = Function(name='f', shape=(int(size/factor),), dimensions=(ci,))
>>> op = Operator([Eq(g, 1), Eq(f, g)])
The Operator generates the following for-loop (pseudocode)
.. code-block:: C
for (int i = i_m; i <= i_M; i += 1) {
g[i] = 1;
if (i%4 == 0) {
f[i / 4] = g[i];
}
}
Another typical use case is when one needs to constrain the execution of
loop iterations so that certain conditions are honoured. The following
artificial example uses ConditionalDimension to guard against out-of-bounds
accesses in indirectly accessed arrays.
>>> from sympy import And
>>> ci = ConditionalDimension(name='ci', parent=i,
... condition=And(g[i] > 0, g[i] < 4, evaluate=False))
>>> f = Function(name='f', shape=(int(size/factor),), dimensions=(ci,))
>>> op = Operator(Eq(f[g[i]], f[g[i]] + 1))
The Operator generates the following for-loop (pseudocode)
.. code-block:: C
for (int i = i_m; i <= i_M; i += 1) {
if (g[i] > 0 && g[i] < 4) {
f[g[i]] = f[g[i]] + 1;
}
}
Now it's time to show a more descriptive example. We define a conditional that increments by one all values of a Function that are larger than zero. We name our ConditionalDimension ci.
In this example we want to update again by one all the elements in f that are larger than zero. Before updating the elements we reinitialize our data to the eye function.
#NBVAL_IGNORE_OUTPUT
from devito import Gt
f.data[:] = np.eye(10)
# Define the condition f(x, y) > 0
condition = Gt(f, 0)
# Construct a ConditionalDimension
ci = ConditionalDimension(name='ci', parent=y, condition=condition)
op1 = Operator(Eq(f, f + 1))
op1.apply()
f.data
# print(op1.ccode) # Uncomment to view code
Operator `Kernel` ran in 0.01 s
Data([[2., 1., 1., 1., 1., 1., 1., 1., 1., 1.],
[1., 2., 1., 1., 1., 1., 1., 1., 1., 1.],
[1., 1., 2., 1., 1., 1., 1., 1., 1., 1.],
[1., 1., 1., 2., 1., 1., 1., 1., 1., 1.],
[1., 1., 1., 1., 2., 1., 1., 1., 1., 1.],
[1., 1., 1., 1., 1., 2., 1., 1., 1., 1.],
[1., 1., 1., 1., 1., 1., 2., 1., 1., 1.],
[1., 1., 1., 1., 1., 1., 1., 2., 1., 1.],
[1., 1., 1., 1., 1., 1., 1., 1., 2., 1.],
[1., 1., 1., 1., 1., 1., 1., 1., 1., 2.]], dtype=float32)
We've constructed ci, but it's still unused, so op1 is still identical to op0. We need to pass ci as an "implicit Dimension" to the equation that we want to be conditionally executed as shown in the next cell.
#NBVAL_IGNORE_OUTPUT
f.data[:] = np.eye(10)
op2 = Operator(Eq(f, f + 1, implicit_dims=ci))
print(op2.body.body[-1])
op2.apply()
assert (np.count_nonzero(f.data - np.diag(np.diagonal(f.data)))==0)
assert (np.count_nonzero(f.data) == 10)
f.data
Operator `Kernel` ran in 0.01 s
/* Begin section0 */
START_TIMER(section0)
for (int x = x_m; x <= x_M; x += 1)
{
#pragma omp simd aligned(f:32)
for (int y = y_m; y <= y_M; y += 1)
{
if (f[x + 1][y + 1] > 0)
{
f[x + 1][y + 1] = f[x + 1][y + 1] + 1;
}
}
}
STOP_TIMER(section0,timers)
/* End section0 */
Data([[2., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
[0., 2., 0., 0., 0., 0., 0., 0., 0., 0.],
[0., 0., 2., 0., 0., 0., 0., 0., 0., 0.],
[0., 0., 0., 2., 0., 0., 0., 0., 0., 0.],
[0., 0., 0., 0., 2., 0., 0., 0., 0., 0.],
[0., 0., 0., 0., 0., 2., 0., 0., 0., 0.],
[0., 0., 0., 0., 0., 0., 2., 0., 0., 0.],
[0., 0., 0., 0., 0., 0., 0., 2., 0., 0.],
[0., 0., 0., 0., 0., 0., 0., 0., 2., 0.],
[0., 0., 0., 0., 0., 0., 0., 0., 0., 2.]], dtype=float32)
The generated code is as expected and only the elements that were greater than zero were incremented by one.
Let's now create a new Function g, initialized with ones along its diagonal.
We want to add f(x,y) to g(x,y) only if (i) g != 0 and (ii) y < 5 (i.e., over the first five columns). To join the two conditions we can use sympy.And.
#NBVAL_IGNORE_OUTPUT
from sympy import And
from devito import Ne, Lt
f.data[:] = np.eye(10)
g = Function(name='g', grid=grid)
g.data[:] = np.eye(10)
ci = ConditionalDimension(name='ci', parent=y, condition=And(Ne(g, 0), Lt(y, 5)))
op3 = Operator(Eq(f, f + g, implicit_dims=ci))
op3.apply()
print(op3.body.body[-1])
assert (np.count_nonzero(f.data - np.diag(np.diagonal(f.data)))==0)
assert (np.count_nonzero(f.data) == 10)
assert np.all(f.data[np.nonzero(f.data[:5,:5])] == 2)
assert np.all(f.data[5:,5:] == np.eye(5))
f.data
Operator `Kernel` ran in 0.01 s
/* Begin section0 */
START_TIMER(section0)
for (int x = x_m; x <= x_M; x += 1)
{
#pragma omp simd aligned(f,g:32)
for (int y = y_m; y <= y_M; y += 1)
{
if (y < 5 && g[x + 1][y + 1] != 0)
{
f[x + 1][y + 1] = f[x + 1][y + 1] + g[x + 1][y + 1];
}
}
}
STOP_TIMER(section0,timers)
/* End section0 */
Data([[2., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
[0., 2., 0., 0., 0., 0., 0., 0., 0., 0.],
[0., 0., 2., 0., 0., 0., 0., 0., 0., 0.],
[0., 0., 0., 2., 0., 0., 0., 0., 0., 0.],
[0., 0., 0., 0., 2., 0., 0., 0., 0., 0.],
[0., 0., 0., 0., 0., 1., 0., 0., 0., 0.],
[0., 0., 0., 0., 0., 0., 1., 0., 0., 0.],
[0., 0., 0., 0., 0., 0., 0., 1., 0., 0.],
[0., 0., 0., 0., 0., 0., 0., 0., 1., 0.],
[0., 0., 0., 0., 0., 0., 0., 0., 0., 1.]], dtype=float32)
You can see that f has been updated only for the first five columns with the f+g expression.
A ConditionalDimension can be also used at Function construction time. Let's use ci from the previous cell to explicitly construct a Function h. You will notice that in this case there is no need to pass implicit_dims to the Eq as the ConditionalDimension ci is now a dimension in h. h still has the size of the full domain, not the size of the domain that satisfies the condition.
#NBVAL_IGNORE_OUTPUT
h = Function(name='h', shape=grid.shape, dimensions=(x, ci))
op4 = Operator(Eq(h, h + g))
op4.apply()
print(op4.body.body[-1])
assert (np.count_nonzero(h.data) == 5)
h.data
Operator `Kernel` ran in 0.01 s
/* Begin section0 */
START_TIMER(section0)
for (int x = x_m; x <= x_M; x += 1)
{
for (int y = y_m; y <= y_M; y += 1)
{
if (y < 5 && g[x + 1][y + 1] != 0)
{
h[x + 1][y + 1] = g[x + 1][y + 1] + h[x + 1][y + 1];
}
}
}
STOP_TIMER(section0,timers)
/* End section0 */
Data([[1., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
[0., 1., 0., 0., 0., 0., 0., 0., 0., 0.],
[0., 0., 1., 0., 0., 0., 0., 0., 0., 0.],
[0., 0., 0., 1., 0., 0., 0., 0., 0., 0.],
[0., 0., 0., 0., 1., 0., 0., 0., 0., 0.],
[0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
[0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
[0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
[0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
[0., 0., 0., 0., 0., 0., 0., 0., 0., 0.]], dtype=float32)
ConditionalDimensions are additionally indicated to implement Function subsampling. In the following example, an Operator subsamples Function g and saves its content into f every factor=4 iterations.
#NBVAL_IGNORE_OUTPUT
size, factor = 16, 4
i = Dimension(name='i')
ci = ConditionalDimension(name='ci', parent=i, factor=factor)
g = Function(name='g', shape=(size,), dimensions=(i,))
# Intialize g
g.data[:,]= list(range(size))
f = Function(name='f', shape=(int(size/factor),), dimensions=(ci,))
op5 = Operator([Eq(f, g)])
print(op5.body.body[-1])
op5.apply()
assert np.all(f.data[:] == g.data[0:-1:4])
print("\n Data in g \n", g.data)
print("\n Subsampled data in f \n", f.data)
Operator `Kernel` ran in 0.01 s
/* Begin section0 */
START_TIMER(section0)
for (int i = i_m; i <= i_M; i += 1)
{
if ((i)%(4) == 0)
{
f[i / 4] = g[i];
}
}
STOP_TIMER(section0,timers)
/* End section0 */
Data in g
[ 0. 1. 2. 3. 4. 5. 6. 7. 8. 9. 10. 11. 12. 13. 14. 15.]
Subsampled data in f
[ 0. 4. 8. 12.]
Finally, we can notice that f holds as expected the subsampled values of g.
ConditionalDimension can be used as an alternative to save snaps to disk as illustrated in 08_snapshotting. The example subsamples the TimeDimension, capturing snaps of the our wavefield at equally spaced snaps.
Here we present an example of using ConditionalDimension by counting the nonzero elements of an array.
from devito.types import Scalar
from devito import Inc
grid = Grid(shape=shape)
# Define function 𝑓. We will initialize f's data with ones on its diagonal.
f = Function(name='f', shape=shape, dimensions=grid.dimensions, space_order=0)
f.data[:] = np.eye(shape[0])
f.data[:]
Data([[1., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
[0., 1., 0., 0., 0., 0., 0., 0., 0., 0.],
[0., 0., 1., 0., 0., 0., 0., 0., 0., 0.],
[0., 0., 0., 1., 0., 0., 0., 0., 0., 0.],
[0., 0., 0., 0., 1., 0., 0., 0., 0., 0.],
[0., 0., 0., 0., 0., 1., 0., 0., 0., 0.],
[0., 0., 0., 0., 0., 0., 1., 0., 0., 0.],
[0., 0., 0., 0., 0., 0., 0., 1., 0., 0.],
[0., 0., 0., 0., 0., 0., 0., 0., 1., 0.],
[0., 0., 0., 0., 0., 0., 0., 0., 0., 1.]], dtype=float32)
#NBVAL_IGNORE_OUTPUT
# ConditionalDimension with Ne(f, 0) condition
ci = ConditionalDimension(name='ci', parent=f.dimensions[-1],
condition=Ne(f, 0))
# g is our counter. g needs to be a Function to write in
g = Function(name='g', shape=(1,), dimensions=(ci,), dtype=np.int32)
# Extra Scalar used to avoid reduction in gcc-5
eq0 = Inc(g[0], 1, implicit_dims=(f.dimensions + (ci,)))
op0 = Operator([eq0])
op0.apply()
assert (np.count_nonzero(f.data) == g.data[0])
print(op0.body.body[-1])
g.data[0]
Operator `Kernel` ran in 0.01 s
/* Begin section0 */
START_TIMER(section0)
for (int x = x_m; x <= x_M; x += 1)
{
for (int y = y_m; y <= y_M; y += 1)
{
if (f[x][y] != 0)
{
int r0 = 1;
g[1] += r0;
}
}
}
STOP_TIMER(section0,timers)
/* End section0 */
10
Here we present another example of using ConditionalDimension. An Operator saves the nonzero positions of an array. We know that we have g.data[0] nonzero elements from Example A.
#NBVAL_IGNORE_OUTPUT
count = g.data[0] # Number of nonzeros
# Dimension used only to nest different size of Functions under the same dim
id_dim = Dimension(name='id_dim')
# Conditional for nonzero element
ci = ConditionalDimension(name='ci', parent=f.dimensions[-1],
condition=Ne(f, 0))
g = Function(name='g', shape=(count, len(f.dimensions)),
dimensions=(ci, id_dim), dtype=np.int32, space_order=0)
eqs = []
# Extra Scalar used to avoid reduction in gcc-5
k = Scalar(name='k', dtype=np.int32)
eqi = Eq(k, -1)
eqs.append(eqi)
eqii = Inc(k, 1, implicit_dims=(f.dimensions + (ci,)))
eqs.append(eqii)
for n, i in enumerate(f.dimensions):
eqs.append(Eq(g[k, n], f.dimensions[n], implicit_dims=(f.dimensions + (ci,))))
# TODO: Must be openmp=False for now due to issue #2061
op0 = Operator(eqs, openmp=False)
op0.apply()
print(op0.body.body[-1])
g.data[:]
ids_np = np.nonzero(f.data[:])
for i in range(len(shape)-1):
assert all(g.data[:, i] == ids_np[i])
Operator `Kernel` ran in 0.01 s
int k = -1;
/* Begin section0 */
START_TIMER(section0)
for (int x = x_m; x <= x_M; x += 1)
{
for (int y = y_m; y <= y_M; y += 1)
{
if (f[x][y] != 0)
{
k += 1;
g[k][0] = x;
g[k][1] = y;
}
}
}
STOP_TIMER(section0,timers)
/* End section0 */