A rule depending on 2 layers and chaining rules using arguments

[virgile-baudrot]

Hi Raf,

the 2 questions are related to previous on DynamicsGrids, but here is to define it using applyrule. I defined this rule "SurvLogLogisticMap", which is almost a copy/paste of growth function your already defined:

@mix @columns struct LC50{LC}
    # Field           | Default  | Flat  | Bounds      | Description
    LC50::LC          | 10.0     | true  | (0.0, 1e9)  | "Lethal concentration for 50% of individuals."
end

@mix @columns struct HillCoefficient{HC}
    # Field           | Default  | Flat  | Bounds      | Description
    hillcoefficient::HC         | 0.1      | true  | (0.0, 1e9) | "Hill coefficient, or shape of a log logistic function"
end

@mix @columns struct Exposure{XP}
    # Field           | Default  | Flat  | Bounds      | Description
    exposure::XP      | 100.0      | true  | (0.0, 1e9) | "Exposure: environmental concentration to which the population is exposed to"
end

abstract type SurvRule{R,W} <: CellRule{R,W} end
abstract type SurvMapRule{R,W} <: SurvRule{R,W} end

DynamicGrids.precalcrules(rule::SurvMapRule, data) = begin
    precalclayer(layer(rule, data), rule, data)
end

@Layers @LC50 @HillCoefficient struct SurvLogLogisticMap{R,W} <: SurvMapRule{R,W} end

@inline applyrule(data, rule::SurvLogLogisticMap, population, index, args...) = begin
    population > zero(population) || return zero(population)
    exposure = layer(rule, data, index)
    @fastmath population * ( 1/ (1 + (exposure / rule.LC50)^rule.hillcoefficient) )

The Map version depend on a single layer "exposure". How would you define it if it was depending on 2 layers, "exposure" and "LC50"?

I generate an output of LC50 using a selection gradient functions (see on the fork: https://github.com/virgile-baudrot/Dispersal.jl/blob/master/src/rules/selectionGradient.jl):

@inline applyrule(data, rule::SelectionGradientSurvMap, LC50, index, args...) = begin
    LC50 > zero(LC50) || return zero(LC50)
    exposure = layer(rule, data, index)
    @fastmath LC50 + rule.additiveGeneticVariance * rule.hillcoefficient*(exposure/LC50)^rule.hillcoefficient / ( LC50 *((exposure/LC50)^rule.hillcoefficient +1)^2)
end

I saw you did it with growthmap but did not get how would you pipe the rules?

Best, Virgile

[virgile-baudrot]

I found how to deal with it which is actually very very simple :-).

I'm going to add simple example in the documentation to build layer series with cat(A,B,C,...; dims = 3), like this one:

popSizeInit = [ 1.0 4.0 7.0;
                2.0 5.0 8.0;
                3.0 6.0 9.0]

intrinsicRate = cat([ 1.0 1.0 1.0;
                    1.0 1.0 1.0;
                    1.0 1.0 1.0],
                    [ 2.0 2.0 2.0;
                    2.0 2.0 2.0;
                    2.0 2.0 2.0],
                    [ 1.0 1.0 1.0;
                    1.0 1.0 1.0;
                    1.0 1.0 1.0]; dims=3)

popSizeGrids = ArrayOutput(popSizeInit; tspan=1:6, aux=(intrinsicRate=intrinsicRate,));
growthRule = Ruleset(DiscreteGrowthMap(layerkey=:intrinsicRate));
sim!(popSizeGrids, growthRule);
@test popSizeGrids[1] == intrinsicRate[:,:,1] .* popSizeInit
@test popSizeGrids[2] == intrinsicRate[:,:,2] .* popSizeGrids[1]
@test popSizeGrids[3] == intrinsicRate[:,:,3] .* popSizeGrids[2]
@test popSizeGrids[4] == intrinsicRate[:,:,1] .* popSizeGrids[3]
@test popSizeGrids[5] == intrinsicRate[:,:,2] .* popSizeGrids[4]
@test popSizeGrids[6] == intrinsicRate[:,:,3] .* popSizeGrids[5]

[rafaqz]

You fork is looking good by the way, great to see you could work with those examples we had. I have a plan to get rid of a lot of those macros too!

With the example, you could useDimArray from https://github.com/rafaqz/DimensionalData.jl/ instead of regular arrays so we can use real DateTime in the index. It's already a dependency of Dispersal.jl. Otherwise you have to use Int for tspan - with DimArray it will find the right dates to cycle over if you use e.g. tspan=DateTime(2001, 1, 1):Month(1):DateTime(2004, 1, 1) in the output.

But maybe that's too complicated!

[rafaqz]

And sorry, for some reason I missed the original issue - sometimes I get a lot of github emails (like hundreds) so don't hessitate to ask again or email me directly.

[virgile-baudrot]

Don't worry, it's perfect. And thank you for this additional way of doing times, which would have real interest when applied to a more real situation.

I actually look deeper into layer object and I expanded to handle several layers in the same time with this functionality:

layer(l::AbstractArray{T,4}, I, timeindex, layerindex) where T =
    l[I..., timeindex, layerindex]  

precalclayer(::AbstractArray{<:Any,4}, rule::Rule, data) =
    @set rule.timeindex = precalc_timeindex(layer(rule, data), rule, data)

DynamicGrids.applyrule(data, rule::LayerCopy, state, index, layerindex, args...) =
    layer(rule, data, index, layerindex)

It's quite cool, still have to be sure it's working in all situation. Basically, it works like this:

@Layers @Timestep struct ExactLogisticGrowthMap2{R,W} <: GrowthMapRule{R,W} end

@inline applyrule(data, rule::ExactLogisticGrowthMap2, population, timeindex, args...) = begin
    population > zero(population) || return zero(population)
    intrinsicrate = layer(rule, data, timeindex, 1)
    carrycap = layer(rule, data, timeindex, 2)
    if intrinsicrate > zero(intrinsicrate)
        @fastmath (population * carrycap) / (population + (carrycap - population) *
                                             exp(-intrinsicrate * rule.nsteps))
    else
        @fastmath population * exp(intrinsicrate * rule.nsteps)
    end
end

popSizeInit = [ 1.0 4.0 7.0;
                    2.0 5.0 8.0;
                    3.0 6.0 9.0]

    intrinsicRate = cat([ 1.0 1.0 1.0;
                        1.0 1.0 1.0;
                        1.0 1.0 1.0],
                        [ 2.0 2.0 2.0;
                        2.0 2.0 2.0;
                        2.0 2.0 2.0],
                        [ 1.0 1.0 1.0;
                        1.0 1.0 1.0;
                        1.0 1.0 1.0]; dims=3)

    carryingCapacity = cat([ 10.0 10.0 10.0;
                            10.0 10.0 10.0;
                            10.0 10.0 10.0],
                            [ 10.0 10.0 10.0;
                            10.0 10.0 10.0;
                            10.0 10.0 10.0],
                            [ 10.0 10.0 10.0;
                            10.0 10.0 10.0;
                            10.0 10.0 10.0]; dims=3)

    popParameter = cat(intrinsicRate, carryingCapacity; dims = 4)

    popSizeGrids = ArrayOutput(popSizeInit; tspan=1:6, aux=(popParameter=popParameter,));
    growthRule = Ruleset(ExactLogisticGrowthMap2(layerkey=:popParameter));
    sim!(popSizeGrids, growthRule);

[rafaqz]

Great, we have talked about carrycap being environmentally driven for ages. It made things just a bit too complicated in the last paper but really it makes much more sense ecologically.

Edit: actually looking at that you shouldn't need any of that for accessing timeindex. It runs automatically for Rule <: GrowthMapRule! So You shouldn't need the LayerCopy things, you should be able to do:

@inline applyrule(data, rule::ExactLogisticGrowthMap2, population, index, args...) = begin
    population > zero(population) || return zero(population)
    intrinsicrate = layer(rule, data, rule.timeindex, 1)
    carrycap = layer(rule, data, rule.timeindex, 2)
    if intrinsicrate > zero(intrinsicrate)
        @fastmath (population * carrycap) / (population + (carrycap - population) *
                                             exp(-intrinsicrate * rule.nsteps))
    else
        @fastmath population * exp(intrinsicrate * rule.nsteps)
    end
end

By the way if you put "```julia" for your code block it will highlight, I'm slow to understand the code without that

[rafaqz]

It may be better if you write out exactly what you are trying to do with the layers, I can probably work out a simpler way for you.

I am also realising I need to document things a lot better.

[virgile-baudrot]

Edit: actually looking at that you shouldn't need any of that for accessing timeindex.

Indeed in this example it's not necessary. But the idea was to have several layers that change in time as parameter.

[rafaqz]

Can you lay out what its for and what the use case is? I don't understand it well enough from the code.

What is that modelling ecologically?

[rafaqz]

I dont understand this bit:

popParameter = cat(intrinsicRate, carryingCapacity; dims = 4)

You can use as many separater aux layers as you want! Just give them different names. Is there a reason they are combined?

It's seems like that isn't very clear from how I've written this.

[virgile-baudrot]

Ok, but how do you get which layer you need in this function

@inline applyrule(data, rule::ExactLogisticGrowthMap2, population, index, args...) = begin
    population > zero(population) || return zero(population)
    intrinsicrate = layer(rule, data, rule.timeindex, 1)
    carrycap = layer(rule, data, rule.timeindex, 2)
    if intrinsicrate > zero(intrinsicrate)
        @fastmath (population * carrycap) / (population + (carrycap - population) *
                                             exp(-intrinsicrate * rule.nsteps))
    else
        @fastmath population * exp(intrinsicrate * rule.nsteps)
    end
end

[rafaqz]

Ok, I see that this is confusing. Layers was written long before aux(data) layers existed. I'll work out a solution for you, this should be much simpler.

[virgile-baudrot]

The rough idea is:

1. a stack of layers defining variation of exposure profile of a chemical
2. computing the evolutionnary response in phenotype of the dose-response trait (assuming infinite population size)
3. using the 2 layers exposure and phenotype response to get population dyanmics

I'm going to write out equations and simulation, I send it to you asap.

[rafaqz]

Just to be clear on the language we are using, it will help these conversations a lot:

"layers" are time-slices of static Arrays attached to Outputs with aux=somedata

"grids" are the dynamic Arrays in the simulation, attached to Outputs with init=(a=inita, b=initb).

Is your 1. static aux layers and 3. dynamic grids?

This is also confusing me. Here you have timeindex:

@inline applyrule(data, rule::ExactLogisticGrowthMap2, population, timeindex, args...) = begin

This argument is the index - its the spatial index in the grid plane for the current cell, coming from DynamicGrids.jl, it doesn't relate to time.

[rafaqz]

The way to do what you want is something like (untested):

# Define the struct from scratch without all the macros, just Base@kwdef for default values
# instead of `layerkey` we have two keys `ratekey` and `carrycapkey` so there are two
# datasets you can retrieve. 
@Base.defkw struct ExactLogisticGrowthMap2{R,W,TI,N} <: GrowthMapRule{R,W,RK,CKTI} 
    ratekey::RK             = Val{:intrisicrate}
    carrycapkey::CK     = Val{:carrycap}
    timeindex::TI           = 1
    nsteps::N                = 1
end

@inline applyrule(data, rule::ExactLogisticGrowthMap2, population, index, args...) = begin
    population > zero(population) || return zero(population)

    # I'll think of a cleaner way to do this part
    intrinsicrate = aux(data)[unwrap(rule.ratekey)][index..., timeindex(rule)]
    carrycap = aux(data)[unwrap(rule.carrycapkey)][index..., timeindex(rule)]

    if intrinsicrate > zero(intrinsicrate)
        @fastmath (population * carrycap) / (population + (carrycap - population) *
                                             exp(-intrinsicrate * rule.nsteps))
    else
        @fastmath population * exp(intrinsicrate * rule.nsteps)
    end
end

Edit: made some updates to the above code.

Sorry the code in growth.jl is confusing. It should be much easier to add fields to structs.

I've just written this in the last week to get rid of all the macros so this code is much easier to write. https://github.com/rafaqz/ModelParameters.jl

[rafaqz]

I'll write a PR this weeks to clean all of this up in Dispersal.jl.

[rafaqz]

But don't worry, keep doing what you are doing for now and get something working. We will clean up later :)

[virgile-baudrot]

Yes, this way to implement is a lot clearer.

I just looked at ModelParameters.jl, it's seems seriously an amazing tool. I did a lot of Bayesian inference and was curious about finding a way to make some inference with Dispersal.jl... it's like you grant my wishes :-D