cesaraustralia / dispersal.jl Goto Github PK
View Code? Open in Web Editor NEWTools for simulating organism dispersal
License: Other
Tools for simulating organism dispersal
License: Other
This issue is used to trigger TagBot; feel free to unsubscribe.
If you haven't already, you should update your TagBot.yml
to include issue comment triggers.
Please see this post on Discourse for instructions and more details.
If you'd like for me to do this for you, comment TagBot fix
on this issue.
I'll open a PR within a few hours, please be patient!
We need to be more explicit about the origins of methodologies and include math representations in the docs.
Following discussion on Population type in #70 (comment)_
I start reviewing the possibilities, given quite big "Comment". There are 2 parts: one about Genetics, a second about Phenotypes
For the Genetics, we have to define Locus
, a Zygote_unit
and a Zygote
(an individual) and Deme
(Population).
I use Locus
definition like https://en.wikipedia.org/wiki/Locus_(genetics). The argument ploidy
which is additional to the definition is to facilitate the handling of polysomy (like trisomy) but don't know if it's relevant for know.
struct Locus # internal
chromosome::Int32 # to specify the number of the chromosome
arm::Bool # 0 and 1 are short and long arm (i.e., p-arm and q-arm)
position::Float32 # for distance and sub-band
ploidy::Int32 # ploidy for this chormosome number
end
# no operation are possible on Locus
chromosome(l::Locus) = l.chromosome
arm(l::Locus) = l.arm
position(l::Locus) = l.position
ploidy(l::Locus) = l.ploidy
L1a = Locus(1,0,1,2) # locus on diploid
chromosome(L1a)
arm(L1a)
position(L1a)
ploidy(L1a)
L1b = Locus(1,0,1,2) # it can have many locus on chromosome 1
L1a == L1b # TRUE
L2 = Locus(2,0,1,2)
L1a == L2 # FALSE
The Zygote_unit
is used to define the zygosity for a single Locus
struct Zygote_unit
locus::Locus
zygosity::Array{Char, 1} # Char should be enough since we have Locus!
# check Ploidy
function Zygote_unit(l,z)
if length(z) != ploidy(l)
throw(ArgumentError("zygosity does not match locus ploidy"))
end
new(l,z)
end
end
locus(zu::Zygote_unit) = zu.locus
zygosity(zu::Zygote_unit) = zu.zygosity
ploidy(zu::Zygote_unit) = ploidy(locus(zu))
chromosome(zu::Zygote_unit) = chromosome(locus(zu))
Lz = Locus(1,0,1,2)
Zygote_unit(Lz,['R','S']) # ok
Zygote_unit(Lz,['R']) # error
Zygote_unit(Lz,['R', 'S', 'R']) # error
Zz = Zygote_unit(Lz,['R','S']) # ok
locus(Zz)
zygosity(Zz)
ploidy(Zz)
chromosome(Zz)
import Base: push!
function push!(zu1::Zygote_unit, zu2::Zygote_unit)
init = [zu1]
if chromosome(zu1) == chromosome(zu2)
if ploidy(zu1) == ploidy(zu2)
push!(init,zu2)
else
throw(ArgumentError("ploidies does not match"))
end
else
push!(init,zu2)
end
end
L1 = Locus(1,0,1,2)
L2 = Locus(2,0,1,2)
L3 = Locus(1,0,1,3)
Z1 = Zygote_unit(L1,['R','R'])
Z2 = Zygote_unit(L2,['R','S'])
push!(Z1,Z2) # should be validated
Z3 = Zygote_unit(L3,['R','S', 'A'])
push!(Z1,Z3) # should be error
Then, a Zygote
is a pool of Zygote_unit
use to define an individual
abstract type AbstractIndividual end
# Individual can be add in a population, multiply (reproduction?), remove from a population, ...
struct Zygote <: AbstractIndividual
genotype::Array{Zygote_unit, 1}
# internal constructor
# # TODO: (1) respect ploidy of chromosome
# # TODO: (2) respect unicity of Locus
# function Zygote(g)
# ...
# end
end
# What is it to make a dispersal of an individual?
# Base.:+(z::Zygote, x::Number) = ... # no sense
# Base.:+(x::Number, z::Zygote) = ... # no sense
# Base.:+(z1::Zygote, z2::Zygote) = ... # could be concat, but concat in Julia is with '*'
# Base.:-(z::Zygote, x::Number) = ... # like +
# Base.:-(x::Number, z::Zygote) = ...
# Base.:-(z1::Zygote, z2::Zygote) = ...
# Base.:*(z::Deme, x::Number) = ... # like reproduction ?
# Base.:*(x::Number, z::Deme) = ... # no sense
# Base.:*(x::Deme, z::Deme) = ... # no sense
# Base.:/(z::Deme, x::Number) = ... # no sense... in-dividual
# Base.:/(x::Number, z::Deme) = ... # no sense... in-dividual
# Base.:/(x::Deme, z::Deme) = ... # no sense... in-dividual
Z1RR = Zygote_unit(L1,['R','R'])
Z1RS = Zygote_unit(L1,['R','S'])
Z1SS = Zygote_unit(L1,['S','S'])
Z2RR = Zygote_unit(L2,['R','R'])
Z2RS = Zygote_unit(L2,['R','S'])
Z2SS = Zygote_unit(L2,['S','S'])
# From this 3 zygote_unit we can have 9 genotypes, each defining a Zygote
Zygote_1RR_2RR = Zygote([Z1RR,Z2RR])
Zygote_1RR_2RS = Zygote([Z1RR,Z2RS])
Zygote_1RR_2SS = Zygote([Z1RR,Z2SS])
Zygote_1RS_2RR = Zygote([Z1RS,Z2RR])
Zygote_1RS_2RS = Zygote([Z1RS,Z2RS])
Zygote_1RS_2SS = Zygote([Z1RS,Z2SS])
Zygote_1SS_2RR = Zygote([Z1SS,Z2RR])
Zygote_1SS_2RS = Zygote([Z1SS,Z2RS])
Zygote_1SS_2SS = Zygote([Z1SS,Z2SS])
# indexin(chromosome(zu),chromosome(z))[1] != nothing
# function push!(z::Zygote,zu::Zygote_unit)
# index = indexin(chromosome(zu),chromosome(z))[1]
# if index != nothing
# if ploidy(zu) == ploidy(z[index])
# push!(zu1,zu2)
# else
# throw(ArgumentError("ploidy does not match"))
# end
# else
# push!(zu1,zu2)
# end
# end
Finally, come the definition for Population. A Deme
is a pool of Zygote
.
# We should be able to define all operation (+,-,*,/, zero and dot) on AbstractPopulation object in order to be used with DynamicGrids.
abstract type AbstractPopulation end
struct Deme <: AbstractPopulation
genotypes::Array{Zygote, 1}
end
Deme([Zygote([Z1RR]), Zygote([Z1RS]), Zygote([Z1SS]), Zygote([Z1SS])])
import Base
# Base.:+(d::Deme, x::Number) = ... # no sense
# Base.:+(x::Number, d::Deme) = ... # no sense
# Base.:+(d1::Deme, d2::Deme) = ... # could be concat, but concat in Julia is with '*'
# Base.:-(d::Deme, x::Number) = ...
# Base.:-(x::Number, d::Deme) = ...
# Base.:-(d1::Deme, d2::Deme) = ...
# Base.:*(d::Deme, x::Number) = ... # Concatenation ?
# Base.:*(x::Number, d::Deme) = ... # a repeat concatenation, would be more relevant with '^'
# Base.:*(x::Deme, d::Deme) = ... # a repeat concatenation, would be more relevant with '^'
# Base.zero(::Type{<:Deme{T}}) where {T} = Deme(zero(T)) # should look like this
Since some difficulties on the Deme
object, I found usefull to define a StatisticsPopulation
. The StatisticsPopulation
is then super easy to use for phenotypes.
# Maybe every thing make sense for a StatisticsPopulation rather than an AbstractPopulation
abstract type StatisticsPopulation <: AbstractPopulation end
struct StatsDeme <: StatisticsPopulation
size::Float64 # size of the population. Could be a Int64 but hard to use frequency/statistics on it
genotypes::Array{Zygote, 1}
frequency::Array{Float64,1}
# TODO: 1. elements in `genotype` must be unique!
# TODO: 2. length of `genotype` must = length frequency => and used to defined size
# TODO: 3. frequency is not neccessary sum to 1, but scaling as to be done!
function StatsDeme(s,g,f)
new(s,g,f./sum(f))
end
end
size(sd::StatsDeme) = sd.size
genotypes(sd::StatsDeme) = sd.genotypes
frequency(sd::StatsDeme) = sd.frequency
SD1 = StatsDeme(
3,
[Zygote([Z1RR]), Zygote([Z1RS]), Zygote([Z1SS])], # the possible genotypes
[1, 2, 7] # the frequency of each
)
size(SD1)
frequency(SD1)
genotypes(SD1)
import Base
Base.:*(a::Number, sd::StatsDeme) = StatsDeme(a*sd.size, sd.genotypes, sd.frequency)
Base.:*(sd::StatsDeme, a::Number) = StatsDeme(a*sd.size, sd.genotypes, sd.frequency)
size(SD1*4)
size(SD1*0.2)
# Base.:+(sd::StatsDeme, a::Number) = ... # no sense ?
# Base.:+(a::Number, sd::StatsDeme) = ... # no sense ?
# Base.:+(sd1::StatsDeme, sd2::StatsDeme) = ... # Concat... but in Julia it's '*'
# Base.:-(sd::StatsDeme, a::Number) = ... # no sense ?
# Base.:-(a::Number, sd::StatsDeme) = ... # no sense ?
# Base.:-(sd1::StatsDeme, sd2::StatsDeme) = ... # Concat... but in Julia it's '*'
# Base.zero(::Type{<:StatsDeme{T1,T2,T3}}) where {T1,T2,T3} = StatsDeme(zero(T1), zero(T2), zero(T3))
These are legacy from when there where two methods for each. Now it seems pretty verbose.
I think LogisticGrowth
is better than ExactLogisticGrowth
, as for ExactExponentialGrowth
and especially ExactExponentialGrowthMap
is really getting too long. @jamesmaino any thoughts? I think you named these.
Hey there,
Dispersal.jl sounds like a cool name, but depending on your background it may mean totally different things. E.g. I come from electron transport and there dispersal really has no connection to what it represents here.
If you want to consider a title more relevant to what the package does, it might help you gain more visibility for this work, which seems to be very good work by looking at the plots. "OrganismDispersal" is something I just thought on the fly, but you are more qualified to choose :P
As we move towards getting Dispersal.jl registered, we will need to clarify what the plan is for the new genetic methods.
I think it would be great if at least some of them were, but I'm not sure what the intent is here.
There are a bunch of things I would like to simplify at some stage, but I also don't want to get in the way of virgile getting work done - I can refactor things once the methods are all finalised.
@jamesmaino probably good if you had a look at these
https://github.com/cesaraustralia/Dispersal.jl/blob/master/src/rules/kernel/common.jl#L108-L191
I feel like I could solve this one, but I failed...
no method matching precalc_auxtimeindex(::GeoArray...
Closest candidates are:
precalc_auxtimeindex(::AbstractArray{var"#s18",3} where var"#s18", ::Rule, ::DynamicGrids.SimData) at C:\Users\virgi\.julia\packages\Dispersal\592tM\src\utils.jl:18
precalc_auxtimeindex(::AbstractArray{var"#s16",3} where var"#s16", ::Rule, ::DynamicGrids.SimData, ::Any) at C:\Users\virgi\.julia\packages\Dispersal\592tM\src\utils.jl:18
precalc_auxtimeindex(::AbstractArray{T,2} where T,
Well, I did simulation by recovering reshaping Arrays includes in GeoArray.data
... but i'm sure you have it
The downsampling code we use to simplify the human dispersal problem is too much of a utility to be in this package. We should use aggregate
from GeoData instead, or move it to DimensionalData.jl (as I've been meaning to for a while) so we don't have to depend on GeoData.jl here - we already depend on DimensionalData.
Adding a weight
parameter/scalar to DispersalKernel
would allow multiplying a size-structured population by a weights vector prior to dispersal to simulate size specific dispersal patterns. As it's a parameter, it could also come from a Grid
or Aux
to use heterogeneous dispersal based on local conditions, for mixed or sized structure populations. Or whatever.
But I think it would be useful. The default value would be nothing
, which wouldn't have any overhead.
It's not ready to include in a registered package - too specific to a particular study.
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
Dear Rafael,
I wanted to test Dispersal, but I had the following issue when precompiling the package. I'm new using Julia, so maybe it's not linked to your package.
Best
julia> using Dispersal
[ Info: Precompiling Dispersal [8797f018-3445-11e9-3969-c13eb94a0661]
WARNING: Method definition similar(DimensionalData.AbstractDimensionalArray{T, N, D} where D<:Tuple where N where T, Type{T}, Tuple{Int64, Vararg{Int64, N} where N}) where {T} in module DimensionalData at /home/vbaudrot/.julia/packages/DimensionalData/RX4n2/src/array.jl:47 overwritten at /home/vbaudrot/.julia/packages/DimensionalData/RX4n2/src/array.jl:61.
** incremental compilation may be fatally broken for this module **
WARNING: Method definition similar(DimensionalData.AbstractDimensionalArray{T, N, D} where D<:Tuple where N where T, Type{T}, Tuple{Union{Integer, Base.OneTo{T} where T<:Integer}, Vararg{Union{Integer, Base.OneTo{T} where T<:Integer}, N}}) where {T, N} in module DimensionalData at /home/vbaudrot/.julia/packages/DimensionalData/RX4n2/src/array.jl:49 overwritten at /home/vbaudrot/.julia/packages/DimensionalData/RX4n2/src/array.jl:63.
** incremental compilation may be fatally broken for this module **
WARNING: Method definition similar(DimensionalData.AbstractDimensionalArray{T, N, D} where D<:Tuple where N where T, Type{T}, Tuple{Int64, Vararg{Int64, N} where N}) where {T} in module DimensionalData at /home/vbaudrot/.julia/packages/DimensionalData/RX4n2/src/array.jl:61 overwritten at /home/vbaudrot/.julia/packages/DimensionalData/RX4n2/src/array.jl:47.
** incremental compilation may be fatally broken for this module **
WARNING: Method definition similar(DimensionalData.AbstractDimensionalArray{T, N, D} where D<:Tuple where N where T, Type{T}, Tuple{Int64, Vararg{Int64, N} where N}) where {T} in module DimensionalData at /home/vbaudrot/.julia/packages/DimensionalData/RX4n2/src/array.jl:47 overwritten at /home/vbaudrot/.julia/packages/DimensionalData/RX4n2/src/array.jl:61.
** incremental compilation may be fatally broken for this module **
WARNING: Method definition similar(DimensionalData.AbstractDimensionalArray{T, N, D} where D<:Tuple where N where T, Type{T}, Tuple{Union{Integer, Base.OneTo{T} where T<:Integer}, Vararg{Union{Integer, Base.OneTo{T} where T<:Integer}, N}}) where {T, N} in module DimensionalData at /home/vbaudrot/.julia/packages/DimensionalData/RX4n2/src/array.jl:63 overwritten at /home/vbaudrot/.julia/packages/DimensionalData/RX4n2/src/array.jl:49.
** incremental compilation may be fatally broken for this module **
WARNING: Method definition similar(DimensionalData.AbstractDimensionalArray{T, N, D} where D<:Tuple where N where T, Type{T}, Tuple{Union{Integer, Base.OneTo{T} where T<:Integer}, Vararg{Union{Integer, Base.OneTo{T} where T<:Integer}, N}}) where {T, N} in module DimensionalData at /home/vbaudrot/.julia/packages/DimensionalData/RX4n2/src/array.jl:49 overwritten at /home/vbaudrot/.julia/packages/DimensionalData/RX4n2/src/array.jl:63.
** incremental compilation may be fatally broken for this module **
WARNING: could not import DynamicGrids.applyinteraction into Dispersal
WARNING: could not import DynamicGrids.applyinteraction! into Dispersal
WARNING: could not import DynamicGrids.neighborhood into Dispersal
WARNING: could not import DynamicGrids.setneighbor! into Dispersal
WARNING: could not import DynamicGrids.mapreduceneighbors into Dispersal
WARNING: could not import DynamicGrids.currenttimestep into Dispersal
WARNING: could not import DynamicGrids.starttime into Dispersal
WARNING: could not import DynamicGrids.stoptime into Dispersal
WARNING: could not import DynamicGrids.tspan into Dispersal
ERROR: LoadError: LoadError: UndefVarError: starttime not defined
Hi Raf,
this is an implementation of a discrete growth with 3 subpopulations (here sharing carying capacity) which work pretty well. If you have any idea to implement something generic for N populations:
@mix @Timestep struct InstrinsicGrowthRate3{GR1, GR2, GR3}
intrinsicrate1::GR1 | 0.1 | true | (0.0, 10.0) | "Intrinsic rate of growth per timestep"
intrinsicrate2::GR2 | 0.1 | true | (0.0, 10.0) | "Intrinsic rate of growth per timestep"
intrinsicrate3::GR3 | 0.1 | true | (0.0, 10.0) | "Intrinsic rate of growth per timestep"
end
@InstrinsicGrowthRate3 @CarryCap struct DiscreteGrowth3{R,W} <: GrowthRule{R,W} end
@inline applyrule(data, rule::DiscreteGrowth3, (population1, population2,population3), args...) = begin
sumPop = population1+population2+population3
( rule.intrinsicrate1 * population1 * (1-sumPop/rule.carrycap),
rule.intrinsicrate2 * population2 * (1-sumPop/rule.carrycap),
rule.intrinsicrate3 * population3 * (1-sumPop/rule.carrycap) )
end
A generic could potentially look like this:
@Timestep struct InstrinsicGrowthRateN{GR}
intrinsicrate::Vector{GR}
numberSupPop::Int
numberSupPop = length(intrinsicrate)
end
# I feel like popState would be like the "hood" is a dispersal context?
@inline applyrule(data, rule::DiscreteGrowth, popState, args...) = begin
sumPop = sum(popState)
...
end
EDIT: I edit typos
You can pass in aux data and a human dispersal rule with the wrong size arrays, and it does weird things and throws indecipherable errors. This bit me a number of times switching model scales between australia and sydney regions.
We should always check aux array size matches gridsize(simdata)
in precalc
wherever we use them.
It's really a gravity model - and could be useful for non-human processes with similar properties.
DateTime of layer is not always matching the DateTime it shoudl be in the context of the simulation.
Write more tests of this.
Layers should be able to be specified with DateTime based Time dimension (using DimensionalData.jl) so that Month() etc can be used, even exact date-time for the whole simulation.
This will be expensive, so we should implement precalc
method for all rules to generate arbitrary precalc data outside of the main loop. It will be functional, not stateful - data objects will be rebuilt for all rules to contain the precalc they want, whatever it is, and passed in during the inner loop.
The upside is it will make numbers/unitful methods faster as well - calculating the indices to interpolate isn't so cheap anyway so they may as well be precalculated.
A declarative, efficient, and flexible JavaScript library for building user interfaces.
๐ Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
An Open Source Machine Learning Framework for Everyone
The Web framework for perfectionists with deadlines.
A PHP framework for web artisans
Bring data to life with SVG, Canvas and HTML. ๐๐๐
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
Some thing interesting about web. New door for the world.
A server is a program made to process requests and deliver data to clients.
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
Some thing interesting about visualization, use data art
Some thing interesting about game, make everyone happy.
We are working to build community through open source technology. NB: members must have two-factor auth.
Open source projects and samples from Microsoft.
Google โค๏ธ Open Source for everyone.
Alibaba Open Source for everyone
Data-Driven Documents codes.
China tencent open source team.