cesaraustralia / dispersal.jl Goto Github PK

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

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

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

# 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

# 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))

 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, 

@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

@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