Fixing boundary checks on masked grids

[nicolasmerino41]

Hey there again :)

I opened the PR to discuss a "small?" implementation of masked boundary checks in OutwardsDispersal. I haven't had much success yet but we're on it :) I've tried to modify OutwardsDispersal such that it can hold an optional mask:

struct OutwardsDispersal{R,W,S<:Stencils.AbstractKernelStencil} <: SetNeighborhoodRule{R,W}
    stencil::S
    mask::Union{Nothing, Array{Bool}}  # Add an optional mask parameter
end

And then modified its applyrule! so it discards cells outside of the mask.

Perhaps I instead should create a brand new OutwardsDispersalWithMask to keep them separate but since I don't know what other parts of the "family" packages are related with these rules, this feels more tedious to build.

Just for the sake of understanding the package behaviour:

• How come InwardsDispersal works fine while not being defined anywhere in the package nor DG/Stencils/etc? Like this is all it's written:

  
  struct InwardsDispersal{R,W,S<:Stencils.AbstractKernelStencil} <: NeighborhoodRule{R,W}
  stencil::S
  end
  function InwardsDispersal{R,W}(; kw...) where {R,W}
  InwardsDispersal{R,W}(DispersalKernel(; kw...))
  end

@inline applyrule(data, rule::InwardsDispersal, N, I) = kernelproduct(rule)

I know what `kernelproduct() ` is used for, but where did you define its method for an InwardsDispersal rule? 

Edit: Btw, just to know. Are tests passing for you? Because I could not pass them even before making any modification.

[rafaqz]

Ok instead of attaching the mask to the rule, we can just get it from data with mask(data).

We can have the rule field as nothing and some other value? maybe a singleton _Mask()? Then in the constructor for the rule we switch true/false to _Mask()/nothing.

This is better than the field being true or false, because the compiler cant remove the check for true/false, but it can for two different types.

[rafaqz]

kernelproduct is in Stencils.jl now (try @which kernelproduct).

Thats available to Dispersal.jl, so it just works.

(we need kernelproduct because Base.dot is recursive, which is awful)

And maybe you wonder how the constructor for InwardsDispersal works?

There are default constructors, see here: https://github.com/cesaraustralia/DynamicGrids.jl/blob/main/src/rules.jl#L96-L114

[nicolasmerino41]

But how will mask(data) work when it's not missingval's what's in the masked area? And I don't think using the with= keyword would work either cause it'd need to be specified before.

[nicolasmerino41]

Well, I'll leave it for today. I don't even know what I'm doing at this point jajaja. I created what I think it's a singleton struct _Mask end so it can be specified when using OutwardsDispersal, otherwise mask_flag is nothing by default.

Then I added the mask(data) approach within applyrule! but still don't understand why it should work, which it doesn't so far :)

@inline function applyrule!(data, rule::OutwardsDispersal{R,W}, N, I) where {R,W}
    N == zero(N) && return nothing

    # Check if the current cell is masked, skip if it is
    mask_data = if rule.mask_flag === NoMask nothing else mask(data) end
    if !isnothing(mask_data) && !mask_data[I...]
        return nothing
    end

    sum = zero(N)
    for (offset, k) in zip(offsets(rule), kernel(rule))
        target = I .+ offset
        (target_mod, inbounds) = inbounds(data, target)
        if inbounds && (isnothing(mask_data) || mask_data[target_mod...])
            @inbounds propagules = N * k  
            @inbounds add!(data[W], propagules, target_mod...)  
            sum += propagules
        end
    end
    @inbounds sub!(data[W], sum, I...)
    return nothing
end

[rafaqz]

I mean DynamicGrids.mask! That's what you will use inside DynamicGrids.jl, as theres no Rasters.jl dep.

Its unfortunate the name is the same, its from before Rasters existed. It just gets the mask object that you pass in to the output.

[nicolasmerino41]

That makes much more sense!

On Thu, 13 Jun 2024, 19:11 Rafael Schouten, @.***> wrote:

I mean DynamicGrids.mask! That's what you will use inside DynamicGrids.jl, as theres no Rasters.jl dep.

Its unfortunate the name is the same, its from before Rasters existed. It just gets the mask object that you pass in to the output.

— Reply to this email directly, view it on GitHub https://github.com/cesaraustralia/Dispersal.jl/pull/114#issuecomment-2166343318, or unsubscribe https://github.com/notifications/unsubscribe-auth/ALXYX3YAAYR6JS5EDNVTXI3ZHHHE3AVCNFSM6AAAAABJIEWBZWVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDCNRWGM2DGMZRHA . You are receiving this because you authored the thread.Message ID: @.***>

[nicolasmerino41]

Good morning. It's all good now. It works perfectly. I took the freedom to write a test for it but you might find it trivial so feel free to tell me to erase it.

It works perfect within the package but I can't use it in another environment. I guess I'm missing some definition? Any idea on that?

Also, I included Stencils in deps, I don't know if that's correct but I was struggling to make use of it without being there, let me know :)

Edit: Also: Tests don't pass in general but that's fine, right? Like they didn't do before anyway, right?

Edit2: Nevermind about the definition, I just needed to specify the package like:

outdisp_with_mask = OutwardsDispersal(
    formulation=ExponentialKernel(λ=0.0125),
    distancemethod=AreaToArea(30),
    mask_flag=Dispersal.Mask()
)

I never understand why sometime you have to specify the package and sometimes not.

[rafaqz]

Amazing!

And trivial tests are actually great. The main idea is to clearly demonstrate that obvious things do work exactly as we intend. Especially that your problem is solved - the populations are conserved over many iterations in the simulation.

Its fine to have Stencils in the deps, to be more explicit. But you can also do using DynamicGrids.Stencils.

Looks like I will have to work on the tests here too to get them passing, I only did DynamicGrids.jl so far. These rules should work on any version of DynamicGrids.jl too.

The last thing for you now (if you havent already) is to make sure things are still fast, and not too much slower with the mask check. Long simulations @time is fine, but BenchmarkTools.jl @btime is more accurate for short runs. You might want to look at a ProfileView.jl @profview profile as well, to make sure there are no type instability coming from the rule. Paste timings and /or profile images here if you do any.

[nicolasmerino41]

True! I forgot to benchmark it. I'll get to it. I was confused with DynamicGrids.Stencils since I thought that was referring to when stencils was part of DG and since I was using the new DG#dev, I assumed stencils feature would come from Stencils.jl but then it didn't work without also including DynamicGrids.Stencils, so quite confusing to me jeje

[nicolasmerino41]

Hey there,

I just ran the benchmarks, sorry for the delay :) So for the default package (before the changes):

# Define a mask
mask_data = [i == 1 || i == 100 || j == 1 || j == 100 ? false : true for i in 1:100, j in 1:100]

# Create a grid with empty borders matching the mask
init = map(x -> x ? 100.0 : 0.0, mask_data)

outdisp = OutwardsDispersal(
    formulation=Dispersal.ExponentialKernel(λ=0.0125),
    distancemethod=AreaToArea(30)
)
# Run the simulation with a mask
output_with_mask = ArrayOutput(init; tspan=1:1000, mask=mask_data)
@btime sim!(output_with_mask, Ruleset(outdisp; boundary=Reflect()))

# Run the simulation without a mask
output_without_mask = ArrayOutput(init; tspan=1:1000)
@btime sim!(output_without_mask, rule_without_mask)

This gave 51ms for the masked simulation and 43ms for the one without.

Now, I ran the exact same code with the new PR implementation, the same two given mask_flag = NoMask is default and an additional one with Mask(), i.e:

# Define a mask
mask_data = [i == 1 || i == 100 || j == 1 || j == 100 ? false : true for i in 1:100, j in 1:100]

# Create a grid with empty borders matching the mask
init = map(x -> x ? 100.0 : 0.0, mask_data)

# Create OutwardsDispersal without a mask, NoMask is default
outdisp_without_mask = OutwardsDispersal(
    formulation=ExponentialKernel(λ=0.0125),
    distancemethod=AreaToArea(30)
)

# Run the simulation with a mask
output_with_mask = ArrayOutput(init; tspan=1:1000, mask=mask_data)
@btime sim!(output_with_mask, Ruleset(outdisp_without_mask; boundary=Reflect()))

# Run the simulation without a mask
output_without_mask = ArrayOutput(init; tspan=1:1000)
@btime sim!(output_without_mask, Ruleset(outdisp_without_mask; boundary=Reflect()))

## TESTING NEW IMPLEMENTATION
# Create OutwardsDispersal with a mask
outdisp_with_mask = OutwardsDispersal(
    formulation=Dispersal.ExponentialKernel(λ=0.0125),
    distancemethod=AreaToArea(30),
    mask_flag=Dispersal.Mask()
)

# Run the simulation with a mask
output_with_mask_and_masked_rule = ArrayOutput(init; tspan=1:1000, mask=mask_data)
@btime sim!(output_with_mask_and_masked_rule, Ruleset(outdisp_without_mask; boundary=Reflect()))

And I got 59ms, 54ms and 127s, respectively. So it looks like using the new feature does slow down the process a lil bit. Fortunately, the implementation only caused a tiny slowing when you choose NOT to use the new feature. We can notify the perfomance drop caused by Mask() so the user is aware. Here I attach a ProfileView screenshot, I've tried to interprete it but it links to many scripts from who knows where, perhaps the amount of red lines gives you a hint.

Have good weekend :) Nico

Ps: For perspective, the first two trials using InwardsDispersal are 33ms and 15ms, for masked and unmasked grid, respectively.

[rafaqz]

Ah interesting. Yes there should be no cost of not using the mask.

That profile is likely the first run, it looks like compilation? Maybe run it a few times ;)

[nicolasmerino41]

Yes, you're right. Here's after repeating: The only long red line refers to: boot.jl, eval: line 385

[rafaqz]

The profile looks fine. Its really the top lines that matter, those are the ones taking time, the lower ones are the functions that lead to them.

[nicolasmerino41]

Hey, so I applied the changes you suggested. Clearly slow it down even more. I've used a different PC but here are the times now: With mask but OutwardsDispersal NoMask(): 197ms Without mask and OutwardsDispersal NoMask(): 193ms With mask and OutwardsDispersal Mask(): 260ms Without mask and InwardsDispersal: 20ms With mask and InwardsDispersal: 50ms

[nicolasmerino41]

Ok, for perspective here's a full benchmarking for the original package (O), my script (M) and your suggestions(S): With mask but OutwardsDispersal NoMask(): O:75ms, M:190ms, S:197ms Without mask and OutwardsDispersal NoMask(): O:66ms, M:186, S:193ms With mask and OutwardsDispersal Mask(): O:NoFeature, M:260ms, S:260ms Without mask and InwardsDispersal: 20ms With mask and InwardsDispersal: 50ms I guess if you'd like to go back to original performance we'd have to take the route of creating a completely different function, like OutwardsDispersalWithMask? so you can keep the performance of the original OutwardsDispersal through not checking any mask at all

[rafaqz]

I don't think we need a different function, as long as we let the compiler see the differences it will just delete them completely.

That line I've highlighted needs to go inside a conditional that depends on the type of mask_flag. Then it wont run. I can rewrite this with a PR to this PR if thats easier for you.

This should work and have the original performance when mask_data is nothing

  sum = zero(N)
    for (offset, k) in zip(offsets(rule), kernel(rule))
        target = I .+ offset
        inbounds = if isnothing(mask_data)    
            true
        else        
            (target_mod, inbounds) = inbounds(data, target)
            mask_data[target_mod...])
        end
        if inbounds
            propagules = N * k  
            @inbounds add!(data[W], propagules, target_mod...)  
            sum += propagules
        end
    end

[nicolasmerino41]

Ok, kind of fixed? Times now are (N), compared to original (O): With mask but OutwardsDispersal NoMask(): N: 76ms, O:76ms Without mask and OutwardsDispersal NoMask(): N:66ms, O:66ms With mask and OutwardsDispersal Mask(): N:278ms, O:NoFeature

However, now I'm getting losses on real boundaries, which didn't use to happen. For example: A non_masked grid with non-mask OutwardsDispersal causes losses. A masked_grid with masked OutwardsDispersal but with "some" mask borders equal to real edges, also cause losses on those edges (e.g. an L shape mask). If no mask edge touches a real edge, then it works fine.
I guess we're missing some bounds-check with this new code?

@inline function applyrule!(data, rule::OutwardsDispersal{R,W}, N, I) where {R,W}
    N == zero(N) && return nothing

    # Check if the current cell is masked, skip if it is
    mask_data = if rule.mask_flag === NoMask() nothing else DynamicGrids.mask(data) end
    if !isnothing(mask_data) && !mask_data[I...]
        return nothing
    end

    sum = zero(N)
    for (offset, k) in zip(offsets(rule), kernel(rule))
        target = I .+ offset
        inbounds = if isnothing(mask_data)    
            true
        else        
            (target_mod, inbounds) = DynamicGrids.inbounds(data, target)
            mask_data[target_mod...]
        end
        if inbounds
            @inbounds propagules = N * k  
            @inbounds add!(data[W], propagules, target...)  
            sum += propagules
        end
    end
    @inbounds sub!(data[W], sum, I...)
    return nothing
end

[rafaqz]

Nice timings!

A non_masked grid with non-mask OutwardsDispersal causes losses

The original algorithm should be identical to what it was. It looks identical to me?

I assumed there was always some losses at the edges for non-masked grids with Remove boundary conditions?

(And bounds checks are pretty expensive, as you can see in your algorithm timings. We could add another singleton like Mask to opt into them?)

[nicolasmerino41]

I don't know, I'm quite confused now, because the tests I wrote actually checked that when there was no mask and NoMask() OutwardsDispersal there should only be floating error loss (which is below 1.0) and I was passing them (when using Wrap or Reflect of course). But now I don't pass them, neither I pass them using the original Dispersal#dev. InwardsDispersal only has floating error loss when not using a mask. I don't think there should be that large loss when not using a mask (and using Wrap or Reflect), I find it quite inconvenient because if you really need to shut down losses you'll have to create a buffer mask. It doesn't personally affect because my masks barely touch the grid edges, but seems suboptimal to me :)

We could add another singleton like Mask to opt into them? I don't really follow you. Like, if you don't specify mask_flag in OutwardsDispersal it will not use the feature, not bound-check and perform as fast as usual, no?

[rafaqz]

Yes having to use a fake mask sounds annoying. Personally I always have a mask if I use Remove, otherwise I would use Wrap (shouldn't Wrap be loss-free? Don't the writes get wrapped around? maybe thats another bug lol). But others may not.

I don't really follow you. Like, if you don't specify mask_flag in OutwardsDispersal it will not use the feature, not bound-check and perform as fast as usual, no?

Ahh I think I get you, yes we can just use the same Mask singleton when there is no mask array just to mean no out of bounds losses. Maybe we should change the name to Lossless or something to cover both of these cases?

(We will also need to add a docstring to describe what these things do and why you would choose to use them or not)

[nicolasmerino41]

shouldn't Wrap be loss-free? Exactly! I think it should but it's not now unless there's a full mask around borders. I feel that bug didn't use to be there.

we can just use the same Mask singleton when there is no mask array just to mean no out of bounds losses That's a good idea, but it doesn't work because if there's no mask(data), the function won't access the bound-checking part of the code and lose individuals anyway. This for example gives a false:

output_without_mask_and_masked_rule = ArrayOutput(init; tspan=1:1000)

r = sim!(output_without_mask_and_masked_rule, Ruleset(outdisp_with_mask; boundary=Reflect()))
sum(r[1]) ≈ sum(r[1000]) #FALSE

Because it's basically the same problem, if a mask doesn't cover all grid borders, it'll lose individuals.

We will also need to add a docstring to describe what these things do and why you would choose to use them or not. Consider it done :)

[rafaqz]

Yes I meant we would have to adjust the code so that the Mask flag will also trigger checking for bounds if there is no Mask object. And for that reason we would need to change the name too so its not connected to mask.

For Wrap and Reflect to work at the boundaries we need to use the target_mod indices that have been fixed to be wrapped or reflected, while for Remove we just don't do the check at all. So we would need another check, probably on the boundary condition being Remove, or anything else.

Or, we just always do the inbounds check for all cases and take the performance hit to keep the algorithm simple... so only write if we are inbounds and if we do then write to the modified indices.

Maybe that's just simpler and better?

(you can probably see why this was simple but a bit buggy... its pretty hard to get all of these things right and have good performance and simple code all at the same time)

[nicolasmerino41]

I mean, that's your call. If you're ok with a decrease in performance I'm fine with that. And InwardsDispersal is always there in need of efficiency.

And yes, I see the tradeoff :)

[rafaqz]

Ok lets just make it correct for now :)

I can go back in an retrieve some Remove boundaries performance later if need be.

[nicolasmerino41]

Ok, so I brought the code back to always bound-check. We're back to the not-so-good timings we discussed the other day but the tests all pass.

I also added a more detailed description in Outwards.jl so the user clearly knows what mask_flag is for. If there is any other place or any other info you'd like to specifiy, let me know :)

[nicolasmerino41]

Cool, all changes done. I guess by AvoidMasked you meant that masked area is out of bounds so you avoid it, but also sounds like you avoid using the mask? Which is the opposite... jeje So, if it's ok with you, I used CheckMaskEdges and IgnoreMaskEdges instead.

[rafaqz]

Yes it's hard to not get a double negative there!

But those names sound good to me.