dsl: Patch issue #1578 (Check legality of symbolic derivatives)

[Leitevmd]

Fix #1578

[codecov]

Codecov Report

Merging #1598 (4a10c37) into master (cfc7a26) will decrease coverage by 8.60%.
The diff coverage is 88.88%.

@@            Coverage Diff             @@
##           master    #1598      +/-   ##
==========================================
- Coverage   86.68%   78.07%   -8.61%     
==========================================
  Files         218      201      -17     
  Lines       33389    31872    -1517     
  Branches     4342     4147     -195     
==========================================
- Hits        28943    24884    -4059     
- Misses       3959     6454    +2495     
- Partials      487      534      +47     

Impacted Files	Coverage Δ
tests/test_derivatives.py	98.74% <86.66%> (-1.26%)	⬇️
devito/types/dimension.py	91.65% <90.00%> (-0.22%)	⬇️
devito/operator/operator.py	88.53% <100.00%> (ø)
devito/types/dense.py	83.18% <100.00%> (-9.88%)	⬇️
tests/test_mpi.py	0.93% <0.00%> (-98.13%)	⬇️
devito/mpi/routines.py	28.25% <0.00%> (-67.35%)	⬇️
tests/test_data.py	40.69% <0.00%> (-56.93%)	⬇️
devito/data/utils.py	36.24% <0.00%> (-56.34%)	⬇️
devito/data/data.py	49.01% <0.00%> (-45.43%)	⬇️
tests/test_gpu_openacc.py	56.69% <0.00%> (-40.95%)	⬇️
... and 82 more

---

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update cfc7a26...4a10c37. Read the comment docs.

[FabioLuporini]

I'm not sure this is the sort of check we're after... and I'll need @mloubout 's input . So, here's my first comments:

• is this patch really going to catch for example a situation like u.dx.dx.dx.dx... ? and, if so, cna we add it to the test case
• I'm not sure looking at deriv_order and space_order is enough. For example, what if we use left or right derivatives (eg u.dxl.dxl.dxl...), which tend to create unbalanced stencils?
• this patch is running an a-priori check . What if we insted run an a-postieriori check? ie, we evaluate, and before returning we simply check that all of the Functions (e.g., u(t, x+3, y)) accesses fall within the halo (e.g., for the aforementioned u(t, x+3, y), we would have to guarantee that 3 < halo_x_left + halo_x_right. Do you see what I mean?

Note: there's a non-negligible probability I'm talking nonsense above. Let's discuss this

[FabioLuporini]

Mumble mumble...

In fact, here's another idea: what if we just ran the "stay within halo" check at indexify() time? that is, here: https://github.com/devitocodes/devito/blob/master/devito/types/basic.py#L992
essentially when the mathematical object u(t, x, y) gets transformed into a lower level object u[t, x, y]. This way we would also intercept misues of the indexed notation: https://github.com/devitocodes/devito/wiki/FAQ#what-is-the-indexed-notation

thoughts?

[mloubout]

I think I like the second one better indexify yes. It would cat indexify/subs/...... cases where an out of halo index may have been introduced out side of FD as well. I don't think testing it in derivative is the way to go indeed.

[Leitevmd]

is this patch really going to catch for example a situation like u.dx.dx.dx.dx... ? and, if so, can we add it to the test case

It doesn't! I wasn't aware of that. MFE updated!

I'm not sure looking at deriv_order and space_order is enough. For example, what if we use left or right derivatives (eg u.dxl.dxl.dxl...), which tend to create unbalanced stencils?

You're right, space_order is not enough, we need to check _size_nodomain sizes for each dimension side.

I was trying to check indices on indexify, and things worked well in common cases. The problem is when we have functions whose indices have been replaced (e.g. f.subs(x, ...)). It can be a new dimension, sub-dimension, symbol, constant, etc. Some of them we just don't know what values they will get, right? In such cases, should we just ignore them? And when I know its limits, e.g. subdimension offsets, should I take this into account (N < subdim_offsets + halos + paddings)?

[FabioLuporini]

@Leitevmd you're right

Some of them we just don't know what values they will get, right

So, I think there are two relevant use cases, and we should focus on them only:

• Dimension
• SubDimension

Other types of Dimensions may be used, but they're not so relevant for the sort of checks we are envisaging here.

I would also observe that:

• SubDimension is just a special Dimension. With a SubDimension I think we should only consider the worst case scenario, i.e. when the offset (what we call "thickness", and what I think you called subdim_offsets?) are 0, that is when the SubDimension coincides (at runtime) with the whole Dimension.

Further:

• I don't think we should take padding into account, because that's not a value useful to the compiler. You could have an Operator op which gets run with different Functions, ie op(u=u1, ...) op(u=u2, ...) op(u=u3, ...) ... then u1, u2, u3 must have the same halo as u (the runtime otherwise will raise an exception), but the padding can be anything. Padding is just a performance feature.

So... this made me realise that we could more easily implement this check within the compiler directly, since there's already a piece of code that checks for OOB accesses, but that was mainly conceived for the SteppingDimensions (in essence, I need to know for various reasons what the max offset is in t, t-1, t+1, t+2...). The check occurs here: https://github.com/devitocodes/devito/blob/master/devito/ir/support/utils.py#L60 , and the return value of this function is used here: https://github.com/devitocodes/devito/blob/master/devito/ir/equations/equation.py#L140. And so I think a possible (much simpler) solution for this feature consists of checking that all of the Dimensions in oobs are is_NonlinearDerived (ie, SteppingDimension or ModuloDimension or ConditionalDimension), while if there's at least one "true" Dimension (Dimension or SubDimension) then we error out?

[Leitevmd]

@FabioLuporini, I think I need some help on this issue. Please check if I’m going in the right direction here. Note that, to avoid several errors, the proposed oob check has the following behavior:

• Arrays ignored
• Subdimensions with thickness != 0 ignored
• Only space dimensions considered

Even with these considerations, there are some errors still occurring, due to the use of Indexeds with out of bound indices on these tests.

[FabioLuporini]

@Leitevmd

first of all, thanks for the update!

Now, your questions:

• Arrays ignored

yes

• Subdimensions with thickness != 0 ignored

I'm going to quote an excerpt from my previous comment:

SubDimension is just a special Dimension. With a SubDimension I think we should only consider the worst case scenario, i.e. when the offset (what we call "thickness", and what I think you called subdim_offsets?) are 0, that is when the SubDimension coincides (at runtime) with the whole Dimension.

do u see what I mean here? I think SubDimensions should be handled

• Only space dimensions considered

OK, this is a good question. Which makes me wonder the following. When we have an OOB, what is the content of the oobs dictionary here ? Shouldn't it have entries != 0 along the (space) dimension which would trigger an OOB?

I don't know why this isn't happening anymore, but now that I think about it, I had designed the equation lowering such that if there were an OOB along one of the dimensions, then the iteration space would have shrunk accordingly (through the non-zero extremes in oobs)

Do you see what I mean here?

While, what we're trying to do with this PR -- and I'm happy with it -- is more like: scratch this shrinking thing non sense (which btw doesn't seem to be working anymore?) and let's raise an exception if your halo isn't enough. But -- as your question points out -- only do that if it's a SpaceDimension (or, similarly, don't do that if it's a TimeDimension)

due to the use of Indexeds with out of bound indices on these tests.

this is OK for now I think, though I'm not sure yet why we're not catching the Indexeds in the current implementation?

[FabioLuporini]

ah, yeah @Leitevmd , you are checking the SubDimension parent, so it seems OK?

[Leitevmd]

Please check the current patch and tests.

Summary

• Dimensions with OOB indices are handled by shrinking the iteration space. No segfaults then.

  • This doesn’t cover SubDimensions

• The only problem of segfault I've seen so far is when using SubDimensions with OOB indices.
  We have two ways to solve it:

  1. Raising an exception for SubDimensions with OOB index
     a) Only for thickness == 0 (CURRENT PATCH)
     b) For all (only if indices exceed thickness)
  2. Increase thickness based on stencil indices

TL;DR

I had designed the equation lowering such that if there were an OOB along one of the dimensions, then the iteration space would have shrunk accordingly (through the non-zero extremes in oobs)

By the way, Isn’t this shrinking mechanism too restrictive? It doesn't consider the halo. For example if we have a function f with space_order=1 we can access f[x-1]+f[x+1] without shrinking x_m and X_M, but if I access f[x-1]+f[x+2], instead of decreasing x_M by 1, it’s decreased by 3 (I doesn't consider halo shift of 2).

I don't know why this isn't happening anymore

It's happening! I couldn’t reproduce segfaults with the example u.biharmonic(1/m) from #1578, nor with those dt.dt.dt.dt, since the iteration space is shrunk to prevent OOB accesses. The only way I could reproduce a segfault was with SubDimensions OOB indices (see why below).

I'm not sure yet why we're not catching the Indexeds in the current implementation?

Yes, we were catching them. We were raising exceptions for these indexed OOB indices. That’s why many tests with Indexeds were failing in my patch. I think these tests have OOB indices relying on the fact that they are handled by shrinking the iteration space.

When we have an OOB, what is the content of the oobs dictionary here ?

oobs is a set containing time and space dimensions with oob indices detected from expression stencils. Here is an example: For a function f with space_order=1, accessing f[t,x+2,y+1] give us a set oobs={t, time, x} because x+2 is out of bound. By the way, why is t always included in oobs?

you are checking the SubDimension parent, so it seems OK?

That’s what the original detect_oobs does. It just checks the parent dimension. If we use a SubDimension as index, the stencil indices will be checked against its parent dimension halo (f._size_nodomain[p].left + f._size_halo[p].right). I think this is OK (although too restrictive since it doesn’t consider thickness when checking OOB indices). I think the problem actually is that, if a SubDimension has an OOB index, only the parent dimension will be included in oobs. That doesn’t help us to handle the OOB when using SubDimensions. They are just ignored, thus incurring in segfaults.

SubDimension is just a special Dimension. With a SubDimension I think we should only consider the worst case scenario, i.e. when the offset (what we call "thickness", and what I think you called subdim_offsets?) are 0, that is when the SubDimension coincides (at runtime) with the whole Dimension.

I think SubDimensions should be handled

By these quotes I understand you proposed that SubDimensions should have the stencil indices checked only if thickness is equal to 0. We would just ignore OOB when SubDimension thickness is not 0. This would still not prevent OOB segfaults with SubDimensions that have thickness != 0.

[FabioLuporini]

aaaahhhhh I see, so the issue is with SubDimensions only... and this is in line with what I have seen recently in other codes making use of SubDimensions. That makes sense.

OK, first of all, one orthogonal thing that I just spotted:

f._size_nodomain[p].left

I think this should rather be the _size_halo[p].left . The _size_nodomain also takes the padding into account, but that's a runtime argument (ie, the padding of the Functions used to construct an Operator doesn't (shouldn't) really matter, because a runtime Function (when u do op.apply(f=f1, ...) could have a different padding , thus there could be a mismatch and faulty indices might be generated. Do u see what I mean?

I understand you proposed that SubDimensions should have the stencil indices checked only if thickness is equal to 0

not quite. I propose to check that as if the thickness were equal to 0 (that's the "worst case scenario" -- the SubDimension would span the same exact iteration space as its parent).

I think the problem actually is that, if a SubDimension has an OOB index, only the parent dimension will be included in oobs. That doesn’t help us to handle the OOB when using SubDimensions. They are just ignored, thus incurring in segfaults.

is this where we currently (master) fail?

only the parent dimension will be included in oobs.

would it help to add to found the whole hierarchy of dimensions, ie here we rather do found.udpate(p._defines) ?

EDIT: never mind, we already do it here

so I'm not entirely sure why is the SubDimension check failing in master yet?

And finally, one thought... for which I'll need both yours and @mloubout 's feedback: I now start thinking that it might just be better to raise an exception in case of OOB rather than "silently"shrinking the iteration space. So we detect the OOB along the non-stepping dimensions, and if there's at least one we raise exception with a suitable message. What do you guys think? So this is clearly orthogonal wrt my non-understanding of why SubDimension OOBs aren't currently detected ☝️

[Leitevmd]

I now start thinking that it might just be better to raise an exception in case of OOB rather than "silently"shrinking the iteration space.

If we're going to raise exceptions for OOB accesses instead of shrinking iteration limits, we’ll need to change several tests, as they rely on this behavior. What if we keep shrinking the iteration space, but warning the user about it? Or if we give an option to the user for this feature?

so I'm not entirely sure why is the SubDimension check failing in master yet?

Master fails if we use something like that Eq(so0[t+1,xi,y], so0[t,xi-2,y] + so0[t,xi+2,y] ), when [xi-2] exceeds the so0 function halo (0 in this example) plus xi subdimension thickness. That's because we don’t take any action if a SubDimention is used in index here, like we do for Dimensions here.

If we want to make something similar to SubDimensions as we do for Dimension today, we could increase ltkn, rtkn based on interval.lower and interval.upper here.

I propose to check that as if the thickness were equal to 0

Now I get it! Not “only if”, but “as if”. The problem with this is we’re going to be too restrictive and some tests which rely on thickness would fail: test_bcs, test_flow_detection_interior and test_subdimmiddle_parallel see here.

[FabioLuporini]

If we're going to raise exceptions for OOB accesses instead of shrinking iteration limits, we’ll need to change several tests, as they rely on this behavior. What if we keep shrinking the iteration space, but warning the user about it? Or if we give an option to the user for this feature?

How many tests would we have to change? Perhaps one or two of these are tests that are actually testing this behaviour... so we could either drop them or change them st they check for the exception to be raised?

Shrinking the iteration space feels like hiding a way a bug in the user code. I think raising a warning is a good start, but not 100% sure it's the best way to go (instead of raising an exception). Any thoughts @mloubout ?

Now I get it! Not “only if”, but “as if”. The problem with this is we’re going to be too restrictive and some tests which rely on thickness would fail: test_bcs, test_flow_detection_interior and test_subdimmiddle_parallel

Ah! yes I see now... and I also remember why those tests have space order 0 (basically to get halo==0, which is something you can get away with, thus saving on memory, if you are using SubDimensions/SubDomains and not using MPI). Good spot.

If we want to make something similar to SubDimensions as we do for Dimension today, we could increase ltkn, rtkn based on interval.lower and interval.upper here.

We can't unfortunately, as that could lead to all sort of bugs. For instance, the thickness is often used to define the size of the domain boundary, thus altering its size would lead to applying boundary conditions over a different region than originally planned.

So, proposal for the detection of a potential OOB: for a SubDimension -- could we rather run the same check we run for Dimension, but rather than checking only against the halo, we check against halo+thickness (ie, thus treating the thickness as a compile-time constant, the very same way we do for the halo)?

Note that in this PR we're talking about two different things:

• fixing the detection of OOB for SubDimensions
• how to use the OOB information -- shrinking iteration space vs warning vs raising an exception

[FabioLuporini]

@Leitevmd

So, proposal for the detection of a potential OOB: for a SubDimension -- could we rather run the same check we run for Dimension, but rather than checking only against the halo, we check against halo+thickness (ie, thus treating the thickness as a compile-time constant, the very same way we do for the halo)?

I guess this won't work, and in any case I don't like it, because we can't formally consider the thickness a compile time parameter.

Maybe you're right, we could rather run this checks at runtime, ie upon op.apply(...). So, when we pass the dspace here to check for the in-boundedness of the array accesses, could we improve SubDimension._arg_check to run these checks? The dspace should (can u confirm this?) contain enough information about "how large" the stencil is, ie how much we might go OOB, plus we know the runtime value of the thickness... is that what you had in mind?

[Leitevmd]

@FabioLuporini please check this again. Now I'm checking thickness in order to put the subdimension into oobs. I noticed thickness is used to express the outside offset in middle subdimensions, but it refers to the inside width in left/right subdimensions. Because of that, on the former, the right limit for OOB is just halo + thickness, while on the latter I’m using halo + (shape - thickness) here. By this I’m assuming a function with shape (2,20,20) and a left x SubDimention with thickness 4 could have a stencil which accesses [t,x+16,y] and that won't cause any illegal memory access. I think that makes sense, and you?

Maybe you're right, we could rather run this checks at runtime, ie upon op.apply(...). So, when we pass the dspace here to check for the in-boundedness of the array accesses, could we improve SubDimension._arg_check to run these checks? The dspace should (can u confirm this?) contain enough information about "how large" the stencil is, ie how much we might go OOB, plus we know the runtime value of the thickness... is that what you had in mind?

I think detect_oobs is the right place to check this since we have all the info we need to catch the OOB. Even if we want to take action in arg_check, we need to detect OOB first in onder fill to dspace with the full interval of the stencil here.

Note that in this PR we're talking about two different things:

• fixing the detection of OOB for SubDimensions
• how to use the OOB information -- shrinking iteration space vs warning vs raising an exception

At this point, we are cathing OOBs for SubDimensions into oobs and then raising an exception if any. We are still shrinking iteration space for OOB with Dimensions. I’m going to check all tests we need to fix in order to change that to raising exceptions.