cesaraustralia / Dispersal.jl

Tools for simulating organism dispersal
Other
21 stars 3 forks source link

working with Tuple is so much easier! #65

Closed virgile-baudrot closed 3 years ago

virgile-baudrot commented 3 years ago

I made the equivalent function replacing structure (Landevariable) by two Tuples. I feel now that structure is no necessary actually. It was certainly like this since the beginning :-). It's indeed so complex to handle dispersal with structure compared to working with Tuple. I also decomposed a function considering Tuple, which is again more generic.

virgile-baudrot commented 3 years ago

I don't know if you added some functions to handle structure in DynamicsGrids.jl or if you just provide test example. I think it's just a test. Also, Maybe something like this could help explaining how to use Tuples (rather that structure): image

EDIT: thinking about it, it's not clear... I try to make a better scheme

rafaqz commented 3 years ago

What do you mean by Tuple ?

Do you mean a single grid of StaticArray? Or using multiple grids which combined pass in a Tuple to the rule?

You can't actually run a grid of Tuple unless you define zero(::Tuple). DynamicGrids.jl is not allowed to define zero for Tuple as that would define it for the whole language, which is type piracy (which is bad).

There is probably also some performance overhead to using two grids instead of a single grid holding StaticArrays.

But a dagram like that is great, maybe we should even put it in this paper I'm writing.

virgile-baudrot commented 3 years ago

I mean a Tuple of grids. Like you did Tuple{ :pop1,:pop2} where pop1 and pop2 are grids.

virgile-baudrot commented 3 years ago

Arff yes, you are right for performance. Here are performance results between the two versions one with a structure named LandeVariable and the other with a Tuple of grids Tuple{:pfreq,:lc50}. The later seems to be longer EDIT: no difference actually

pfreqinit = [0.1 0.5;
            0.5 0.8]
lc50init = [50.0 50.0; 
            5.0 5.0]
exposure = [50.0 50.0;
            10.0 10.0]

# 1 STRUCTURE WITH STATIC ARRAY
landeinit_str = [LandeVariable(pfreqinit[i,j], lc50init[i,j]) for i in 1:2, j in 1:2]
output_str = ArrayOutput(landeinit_str; tspan=1:3, aux=(exposure=exposure,))
rule_str = SelectionGradient1locusSurvMap(
    exposurekey=:exposure,
    hillcoefficient=2.0,  
    deviation_phenotype=10.0,
    dominance_degree=-1.0
)
# 2 TUPLE OF GRIDS
landeinit = (pfreq = pfreqinit,
            lc50 = lc50init,)
output_tuple = ArrayOutput(landeinit; tspan=1:3, aux=(exposure=exposure,))
grids = Tuple{:pfreq,:lc50}
rule_tuple = SelectionGradientMapTuple{grids,grids}(
    exposurekey=:exposure,
    hillcoefficient=2.0,  
    deviation_phenotype=10.0,
    dominance_degree=-1.0
)

Then:

using BenchmarkTools
julia> @btime sim!(output_tuple, rule_tuple)
  13.200 μs (183 allocations: 11.92 KiB)
julia> @btime sim!(output_tuple, rule_tuple)
  13.399 μs (193 allocations: 12.14 KiB)

And with bigger grids:

gridSize = 1000
pfreqinit = rand(gridSize,gridSize)
lc50init =  rand(gridSize,gridSize).*50
exposure =  rand(gridSize,gridSize).*50
julia> @btime sim!(output_str, rule_str)
  116.758 ms (2000111 allocations: 61.04 MiB)
julia> @btime sim!(output_tuple, rule_tuple)
  118.264 ms (2000195 allocations: 61.05 MiB)
rafaqz commented 3 years ago

Ok that's interesting, I haven't actually checked for ages, we probably need a benchmark suite to track all of these things.

Although, I'm not totally sure those are good timings, it's better to run 100 steps in the simulation than 3 or it has a lot of setup cost. Also have some type instability somewhere, the allocations that suggest that: (2000195 allocations: 61.05 MiB). The simulation shouldn't really allocate like that while it runs - allocations should be in thousands or tens of thousands, not millions.

But don't worry I can fix that later, its only an issue when you try to run sensitivity analysis or optimise the model.

Also, if you're interested in some detail: Tuple in Tuple{:grid,:grid} is not really a tuple of grids, it's just a way we use the type system to specify that you want to receive a NamedTuple in the rule, containing values from the grids you specify. We are just using Tuple as a container for Symbol because you cant dispatch on the contents of (:grid1, grid2) as a type, but you can with Tuple{:grid1,:grid2}. Notice that Tuple{:grid1,:grid2} doesn't describe an actual object that could be constructed!

The grids are actually stored in a NamedTuple inside a SimData object, and we use the grid names you specify to get those grids and rebuild a reduced SimData object with just the grids you want.

Then for every cell we build a NamedTuple of values to pass to the rule for each specific index. If you take it as a single argument you will find it's a NamedTuple, not a Tuple! you can index into it with the grid name like state[:grid1]. See here: https://github.com/cesaraustralia/DynamicGrids.jl/blob/master/src/maprules.jl#L410

But what you are labelling Tuple, we have been calling multiple grids

rafaqz commented 3 years ago

I don't know if you added some functions to handle structure in DynamicsGrids.jl or if you just provide test example. I think it's just a test.

There were a few minor changes to use zero correctly, but otherwise it just worked

virgile-baudrot commented 3 years ago

But what you are labelling Tuple, we have been calling multiple grids

ok, I change naming I did

virgile-baudrot commented 3 years ago

Although, I'm not totally sure those are good timings, it's better to run 100 steps in the simulation than 3 or it has a lot of setup cost. Also have some type instability somewhere, the allocations that suggest that: (2000195 allocations: 61.05 MiB). The simulation shouldn't really allocate like that while it runs - allocations should be in thousands or tens of thousands, not millions

Ok. Indeded, in this sense, I also see the benefit of using Val() as for instance Val(:exposure) which reduce allocations to 193 and 20MiB.

rafaqz commented 3 years ago

Ah ok thats an easy fix. Maybe DynamicGrids should force you to use Val so that cant happen. The rule is accessing it ~1 million times per frame so it needs to be type stable.