Repository: devitocodes/devito
Issue: ConditionalDimension misbehaving a little?

The following code appears to misbehave.

```
from devito import Grid, TimeFunction, ConditionalDimension, Operator, Constant, Eq, configuration

# deci = Constant(name='deci', dtype=np.int32)
deci=8
g = Grid(shape=(50, 50, 50))
time_sub = ConditionalDimension(name='time_sub', parent=g.time_dim, factor=deci)
u = TimeFunction(name='u', grid=g, time_order=2, space_order=4)
usave = TimeFunction(name='usave', grid=g, time_dim=time_sub, time_order=0, space_order=0, save=100)

# Observe that the laplace ends up inside the %8 conditional
actions = [Eq(usave, u), Eq(u.forward, u.laplace+1)]
op = Operator(actions, dle="advanced")
print("---- broken")
print(op.ccode)

# But not in this case
actions = [Eq(u.forward, u.laplace+1), Eq(usave, u)]
op = Operator(actions, dle="advanced")
print("---- working")
print(op.ccode)
```

The "broken" C code here ends up looking like:

```
int Kernel(struct dataobj *restrict u_vec, struct dataobj *restrict usave_vec, const float h_x, const float h_y, const float h_z, const int time_M, const int time_m, struct profiler * timers, const int x0_blk_size, const int x_M, const int x_m, const int y0_blk_size, const int y_M, const int y_m, const int z_M, const int z_m)
{
  for (int time = time_m, t0 = (time)%(3), t1 = (time + 1)%(3); time <= time_M; time += 1, t0 = (time)%(3), t1 = (time + 1)%(3))
  {
    if ((time)%(8) == 0)
    {
      struct timeval start_section0, end_section0;
      gettimeofday(&start_section0, NULL);
      /* Begin section0 */
      bf0(h_x,h_y,h_z,u_vec,usave_vec,t0,t1,time,x_M - (x_M - x_m + 1)%(x0_blk_size),x_m,x0_blk_size,y_M - (y_M - y_m + 1)%(y0_blk_size),y_m,y0_blk_size,z_M,z_m);
      bf0(h_x,h_y,h_z,u_vec,usave_vec,t0,t1,time,x_M - (x_M - x_m + 1)%(x0_blk_size),x_m,x0_blk_size,y_M,y_M - (y_M - y_m + 1)%(y0_blk_size) + 1,(y_M - y_m + 1)%(y0_blk_size),z_M,z_m);
      bf0(h_x,h_y,h_z,u_vec,usave_vec,t0,t1,time,x_M,x_M - (x_M - x_m + 1)%(x0_blk_size) + 1,(x_M - x_m + 1)%(x0_blk_size),y_M - (y_M - y_m + 1)%(y0_blk_size),y_m,y0_blk_size,z_M,z_m);
      bf0(h_x,h_y,h_z,u_vec,usave_vec,t0,t1,time,x_M,x_M - (x_M - x_m + 1)%(x0_blk_size) + 1,(x_M - x_m + 1)%(x0_blk_size),y_M,y_M - (y_M - y_m + 1)%(y0_blk_size) + 1,(y_M - y_m + 1)%(y0_blk_size),z_M,z_m);
      /* End section0 */
      gettimeofday(&end_section0, NULL);
      timers->section0 += (double)(end_section0.tv_sec-start_section0.tv_sec)+(double)(end_section0.tv_usec-start_section0.tv_usec)/1000000;
    }
  }
  return 0;
}

void bf0(const float h_x, const float h_y, const float h_z, struct dataobj *restrict u_vec, struct dataobj *restrict usave_vec, const int t0, const int t1, const int time, const int x0_blk_M, const int x0_blk_m, const int x0_blk_size, const int y0_blk_M, const int y0_blk_m, const int y0_blk_size, const int z_M, const int z_m)
{
  float (*restrict u)[u_vec->size[1]][u_vec->size[2]][u_vec->size[3]] __attribute__ ((aligned (64))) = (float (*)[u_vec->size[1]][u_vec->size[2]][u_vec->size[3]]) u_vec->data;
  float (*restrict usave)[usave_vec->size[1]][usave_vec->size[2]][usave_vec->size[3]] __attribute__ ((aligned (64))) = (float (*)[usave_vec->size[1]][usave_vec->size[2]][usave_vec->size[3]]) usave_vec->data;
  for (int x0_blk = x0_blk_m; x0_blk <= x0_blk_M; x0_blk += x0_blk_size)
  {
    for (int y0_blk = y0_blk_m; y0_blk <= y0_blk_M; y0_blk += y0_blk_size)
    {
      for (int x = x0_blk; x <= x0_blk + x0_blk_size - 1; x += 1)
      {
        for (int y = y0_blk; y <= y0_blk + y0_blk_size - 1; y += 1)
        {
          #pragma omp simd aligned(u,usave:16)
          for (int z = z_m; z <= z_M; z += 1)
          {
            usave[time / 8][x][y][z] = u[t0][x + 4][y + 4][z + 4];
            u[t1][x + 4][y + 4][z + 4] = -2.5F*(u[t0][x + 4][y + 4][z + 4]/((h_z*h_z)) + u[t0][x + 4][y + 4][z + 4]/((h_y*h_y)) + u[t0][x + 4][y + 4][z + 4]/((h_x*h_x))) - 8.33333333e-2F*(u[t0][x + 4][y + 4][z + 2]/((h_z*h_z)) + u[t0][x + 4][y + 4][z + 6]/((h_z*h_z)) + u[t0][x + 4][y + 2][z + 4]/((h_y*h_y)) + u[t0][x + 4][y + 6][z + 4]/((h_y*h_y)) + u[t0][x + 2][y + 4][z + 4]/((h_x*h_x)) + u[t0][x + 6][y + 4][z + 4]/((h_x*h_x))) + 1.33333333F*(u[t0][x + 4][y + 4][z + 3]/((h_z*h_z)) + u[t0][x + 4][y + 4][z + 5]/((h_z*h_z)) + u[t0][x + 4][y + 3][z + 4]/((h_y*h_y)) + u[t0][x + 4][y + 5][z + 4]/((h_y*h_y)) + u[t0][x + 3][y + 4][z + 4]/((h_x*h_x)) + u[t0][x + 5][y + 4][z + 4]/((h_x*h_x))) + 1;
          }
        }
      }
    }
  }
}
```

Note how `u` is only set every 8 iterations, when it should be every time.
The order of these two `Eq`s should not matter, but it appears to.

Thanks!
Tim

Changed files:
tests/test_dimension.py, devito/ir/stree/algorithms.py

diff --git a/devito/ir/stree/algorithms.py b/devito/ir/stree/algorithms.py
index 07a813451..b6195dc0d 100644
--- a/devito/ir/stree/algorithms.py
+++ b/devito/ir/stree/algorithms.py
@@ -66,11 +66,16 @@ def st_schedule(clusters):
         NodeExprs(c.exprs, c.ispace, c.dspace, c.ops, c.traffic, root)
 
         # Add in Conditionals
-        for k, v in mapper.items():
+        drop_guarded = None
+        for k, v in list(mapper.items()):
+            if drop_guarded:
+                mapper.pop(k)
             if k.dim in c.guards:
                 node = NodeConditional(c.guards[k.dim])
                 v.last.parent = node
                 node.parent = v
+                # Drop nested guarded sub-trees
+                drop_guarded = True
 
     return stree
 
diff --git a/tests/test_dimension.py b/tests/test_dimension.py
index fd4ad11af..69342efbe 100644
--- a/tests/test_dimension.py
+++ b/tests/test_dimension.py
@@ -483,6 +483,46 @@ class TestConditionalDimension(object):
         assert np.all([np.allclose(usave.data[i], i*factor)
                       for i in range((nt+factor-1)//factor)])
 
+    def test_basic_shuffles(self):
+        """
+        Like ``test_basic``, but with different equation orderings. Nevertheless,
+        we assert against the same exact values as in ``test_basic``, since we
+        save `u`, not `u.forward`.
+        """
+        nt = 19
+        grid = Grid(shape=(11, 11))
+        time = grid.time_dim
+
+        u = TimeFunction(name='u', grid=grid)
+
+        u2 = TimeFunction(name='u2', grid=grid, save=nt)
+
+        factor = 4
+        time_subsampled = ConditionalDimension('t_sub', parent=time, factor=factor)
+        usave = TimeFunction(name='usave', grid=grid, save=(nt+factor-1)//factor,
+                             time_dim=time_subsampled)
+
+        # Shuffle 1
+        eqns = [Eq(usave, u), Eq(u.forward, u + 1.), Eq(u2.forward, u2 + 1.)]
+        op = Operator(eqns)
+        op.apply(t_M=nt-2)
+        assert np.all(np.allclose(u.data[(nt-1) % 3], nt-1))
+        assert np.all([np.allclose(u2.data[i], i) for i in range(nt)])
+        assert np.all([np.allclose(usave.data[i], i*factor)
+                      for i in range((nt+factor-1)//factor)])
+
+        # Shuffle 2
+        usave.data[:] = 0.
+        u.data[:] = 0.
+        u2.data[:] = 0.
+        eqns = [Eq(u.forward, u + 1.), Eq(usave, u), Eq(u2.forward, u2 + 1.)]
+        op = Operator(eqns)
+        op.apply(t_M=nt-2)
+        assert np.all(np.allclose(u.data[(nt-1) % 3], nt-1))
+        assert np.all([np.allclose(u2.data[i], i) for i in range(nt)])
+        assert np.all([np.allclose(usave.data[i], i*factor)
+                      for i in range((nt+factor-1)//factor)])
+
     def test_spacial_subsampling(self):
         """
         Test conditional dimension for the spatial ones.