In [1]:
using Distributions
using ProfileView
using StatsFuns
using Gadfly
using ProgressMeter
using DataFrames
using Colors
using Viridis
using DataStructures


INFO: Recompiling stale cache file /home/diogro/.julia/lib/v0.5/Distributions.ji for module Distributions.
INFO: Recompiling stale cache file /home/diogro/.julia/lib/v0.5/ProfileView.ji for module ProfileView.
WARNING: The call to compilecache failed to create a usable precompiled cache file for module Colors. Got:
WARNING: Module Reexport uuid did not match cache file.
WARNING: eval from module Main to ProfileView:    
Expr(:call, Expr(:., :Base, :include_from_node1)::Any, "/home/diogro/.julia/v0.5/Colors/src/Colors.jl")::Any
  ** incremental compilation may be broken for this module **

INFO: Precompiling module Gadfly.
WARNING: Module Colors with uuid 254316556476747 is missing from the cache.
This may mean module Colors does not support precompilation but is imported by a module that does.
ERROR: LoadError: Declaring __precompile__(false) is not allowed in files that are being precompiled.
 in require(::Symbol) at ./loading.jl:385
 in include_from_node1(::String) at ./loading.jl:488
 in macro expansion; at ./none:2 [inlined]
 in anonymous at ./<missing>:?
 in eval(::Module, ::Any) at ./boot.jl:234
 in process_options(::Base.JLOptions) at ./client.jl:242
 in _start() at ./client.jl:321
while loading /home/diogro/.julia/v0.5/Gadfly/src/Gadfly.jl, in expression starting on line 5
Failed to precompile Gadfly to /home/diogro/.julia/lib/v0.5/Gadfly.ji.

 in compilecache(::String) at ./loading.jl:593
 in require(::Symbol) at ./loading.jl:422
 in include_string(::String, ::String) at ./loading.jl:441

In [2]:
type Individual
    m::Int64
    p::Int64
    mu::Float64
    mu_b::Float64
    mu_sigma::Float64
    y::Array{Float64, 1}
    b::Array{Int64, 2}
    z::Array{Float64, 1}
    Individual(m, p, mu, mu_b, mu_sigma) = new(2*m, p, mu, mu_b, mu_sigma)
end  

type Population
    n::Int64
    base_ind::Individual
    generation::Int64
    fitness::Array{Float64, 1}
    surface::Distributions.FullNormal
    pop::Array{Individual, 1}
    moments::Array{Any, 1}
    Population(n, base_ind, generation) = new(n, base_ind, generation)
end

In [3]:
function Population(n::Int, base_ind::Individual, pop_array::Array{Individual, 1})
    new_pop = Population(n, base_ind, 0)
    new_pop.generation = 0
    new_pop.pop = pop_array
    new_pop.moments = []
    new_pop
end

function RandomInd!(ind::Individual, base_ind::Individual)
    ind.y = rand(Normal(0, base_ind.mu_sigma), ind.m)
    ind.b = isdefined(base_ind, :b) ? copy(base_ind.b) : rand([-1, 0, 1], ind.p, ind.m)
    ind.z = ind.b * ind.y
end

function RandomInd(base_ind::Individual)
    new_ind = Individual(base_ind.m/2, base_ind.p, base_ind.mu, base_ind.mu_b, base_ind.mu_sigma)
    RandomInd!(new_ind, base_ind)
    new_ind
end

function RandomPop!(pop::Population)
    pop.pop = Array(Individual, pop.n)
    pop.generation = 0
    for i = 1:pop.n
        pop.pop[i] = RandomInd(pop.base_ind)
    end
    pop.moments = []
    pop.fitness = zeros(Float64, pop.n)
end


Out[3]:
RandomPop! (generic function with 1 method)

In [24]:
import Base.string
function string(pop::Population)
    return "a Population with $(pop.n) individuals, at generation $(pop.generation)"
end

import Base.print
print(io::IO, pop::Population) = print(io, string(pop))

import Base.show
show(io::IO, pop::Population) = print(io, "This is ", pop)



import Base.getindex
function getindex(pop::Population, i::Integer)
    getindex(pop.pop, i)
end

function getindex(pop::Population, s::UnitRange)
    Population(length(s), pop.base_ind, getindex(pop.pop, s))
end

function append!(pop::Population, ind::Individual)
    pop.pop = [pop.pop; ind]
    pop.n = length(pop.pop)
end

function append!(pop1::Population, pop2::Population)
    # COMPARE BASE IND!
    pop1.pop = [pop1.pop; pop2.pop]
    pop.n = length(pop1.pop)
end

function join(pop1::Population, pop2::Population)
    # TODO: COMPARE BASE IND!
    new_pop = Population(pop1.n + pop2.n, pop1.base_ind)
    new_pop.pop = [pop1.pop; pop2.pop]
    new_pop
end

function copy!(source::Population, sink::Population)
    sink.n          = source.n
    sink.base_ind   = source.base_ind
    sink.generation = source.generation
    sink.fitness    = source.fitness
    sink.surface    = source.surface
    sink.pop        = copy(source.pop)
end

import Base.copy
function copy(source::Population)
    new_pop = Population(source.n, source.base_ind, source.generation)
    copy!(source, new_pop)
    new_pop
end
    
function mutation!(ind::Individual, bin_y, bin_b)
    mutation_y = rand(bin_y)
    if(mutation_y > 0)
        d_uni_y = DiscreteUniform(1, ind.m)
        norm_sigma = Normal(0, ind.mu_sigma)
        for k = range(1, mutation_y)
            i = rand(d_uni_y)
            ind.y[i] = ind.y[i] + rand(norm_sigma)
        end
    end
    mutation_b = rand(bin_b)
    if(mutation_b > 0)
        d_uni_p = DiscreteUniform(1, ind.p)
        d_uni_m = DiscreteUniform(1, ind.m)
        for k = range(1, mutation_b)
            i = rand(d_uni_p)
            j = rand(d_uni_m)
            ind.b[i, j] = ind.b[i, j] == 0 ? rand([-1, 1]) : 0
        end
    end
end

function mutation!(pop::Population)
    bin_y = Binomial(pop.base_ind.m, pop.base_ind.mu)
    bin_b = Binomial(pop.base_ind.m * pop.base_ind.p, pop.base_ind.mu_b)
    for k = 1:pop.n
        mutation!(pop.pop[k], bin_y, bin_b)
    end
end

function fitness!(pop::Population)
    logfit = Float64[logpdf(pop.surface, pop[i].z) for i in 1:pop.n]
    #logfit = Float64[logpdf(pop.surface, pop[i].z) * (1-0.95*sum(abs(pop[i].b))/(size(pop[i].b)[1] * size(pop[i].b)[2])) for i in 1:pop.n]
    pop.fitness = exp(logfit - logsumexp(logfit))
end

function changeSurface!(pop::Population, theta::Array{Float64, 1}, omega::Array{Float64, 2})
    pop.surface = MvNormal(theta, omega)
end

function twoModuleMatrix(vars::Array{Float64, 1}, cor1, cor2)
    n = size(vars)[1]
    n_1 = Int(floor(n/2))
    module_1 = [ones(Float64, n_1); zeros(Float64, n - n_1)]
    module_2 = [zeros(Float64, n_1); ones(Float64, n - n_1)]
    omega = (module_1 * module_1') * cor1 + module_2 * module_2' * cor2
    [omega[i, i] = 1 for i = 1:n]
    omega .* (sqrt(vars) * sqrt(vars)')
end

function cross!(ind_1::Individual, ind_2::Individual, new_ind::Individual,
                d_uni, alele_1::Array{Int64, 1}, alele_2::Array{Int64, 1})
    rand!(d_uni, alele_1)
    rand!(d_uni, alele_2)
    for locus = range(2, 2, convert(Int64, ind_1.m/2))
        new_ind.y[   locus - 1] = ind_1.y[   (locus - 1) + alele_1[convert(Int64, locus/2)]]
        new_ind.y[   locus]     = ind_2.y[   (locus - 1) + alele_2[convert(Int64, locus/2)]]
        for j in range(1, ind_1.p)
            new_ind.b[j, locus - 1] = ind_1.b[j, (locus - 1) + alele_1[convert(Int64, locus/2)]]
            new_ind.b[j, locus]     = ind_2.b[j, (locus - 1) + alele_2[convert(Int64, locus/2)]]
        end
    end
    new_ind.z = new_ind.b * new_ind.y
end

function choose_mates!(pop::Population, mates::Array{Int64, 1})
    matings = rand(Multinomial(pop.n, pop.fitness), 1)
    l = 1
    for k = 1:pop.n
        if(matings[k] > 0)
            for i = 1:matings[k]
                mates[l] = k
                l = l + 1
            end
        end
    end
    round(Int64, shuffle!(mates))
end

function next_generation!(pop::Population, holder_pop::Population, sires, dames, d_uni, alele_1, alele_2; selective = true)
    holder_pop.surface = pop.surface
    if (selective)
        fitness!(pop)
    else
        fill!(pop.fitness, 1./pop.n)
    end
    holder_pop.fitness = pop.fitness
    mutation!(pop)
    choose_mates!(pop, sires)
    choose_mates!(pop, dames)
    for i in 1:pop.n
        cross!(pop[sires[i]], pop[dames[i]], holder_pop.pop[i], d_uni, alele_1, alele_2)
    end
    holder_pop.generation = pop.generation + 1
    copy!(holder_pop, pop)
end 

function moments!(pop::Population)
    pop.moments = [pop.moments; moments(pop)]
end

function moments(pop::Population)
    ys = convert(Array{Float64, 2}, reshape(Float64[ind.y[i] for ind in pop.pop, i in 1:pop.base_ind.m], pop.n, pop.base_ind.m))
    zs = convert(Array{Float64, 2}, reshape(Float64[ind.z[i] for ind in pop.pop, i in 1:pop.base_ind.p], pop.n, pop.base_ind.p))
    mean_b = zeros(Float64, pop.base_ind.p, pop.base_ind.m)
    for i in 1:pop.n
        mean_b = mean_b + pop[i].b
    end
    mean_y = squeeze(mean(ys, 1), 1)
    mean_z = squeeze(mean(zs, 1), 1)
    mean_b = mean_b / pop.n
    
    count_b = countBclass(classifyBj(convert(Array{Int64, 2}, round(Int64, mean_b))))    
    count_b[:gen] = pop.generation
    
    G = cov(zs)
    corrG = cor(zs)
    
    
    
    Dict([("mean_y", mean_y), 
          ("mean_b", mean_b), 
          ("mean_z", mean_z),
          ("zs", zs), 
          ("theta", pop.surface.μ),
          ("gen", pop.generation),
          ("G", G), 
          ("corrG", corrG),
          ("count_b", count_b)])
end


WARNING: Method definition string(Main.Population) in module Main at In[4]:3 overwritten at In[24]:3.
WARNING: Method definition print(IO, Main.Population) in module Main at In[4]:7 overwritten at In[24]:7.
WARNING: Method definition show(IO, Main.Population) in module Main at In[4]:10 overwritten at In[24]:10.
WARNING: Method definition getindex(Main.Population, Integer) in module Main at In[4]:16 overwritten at In[24]:16.
WARNING: Method definition getindex(Main.Population, Base.UnitRange) in module Main at In[4]:20 overwritten at In[24]:20.
WARNING: Method definition append!(Main.Population, Main.Individual) in module Main at In[4]:24 overwritten at In[24]:24.
WARNING: Method definition append!(Main.Population, Main.Population) in module Main at In[4]:30 overwritten at In[24]:30.
WARNING: Method definition join(Main.Population, Main.Population) in module Main at In[4]:36 overwritten at In[24]:36.
WARNING: Method definition copy!(Main.Population, Main.Population) in module Main at In[4]:42 overwritten at In[24]:42.
WARNING: Method definition copy(Main.Population) in module Main at In[4]:52 overwritten at In[24]:52.
WARNING: Method definition mutation!(Main.Individual, Any, Any) in module Main at In[4]:58 overwritten at In[24]:58.
WARNING: Method definition mutation!(Main.Population) in module Main at In[4]:80 overwritten at In[24]:80.
WARNING: Method definition fitness!(Main.Population) in module Main at In[4]:88 overwritten at In[24]:88.
WARNING: Method definition changeSurface!(Main.Population, Array{Float64, 1}, Array{Float64, 2}) in module Main at In[4]:94 overwritten at In[24]:94.
WARNING: Method definition twoModuleMatrix(Array{Float64, 1}, Any, Any) in module Main at In[4]:98 overwritten at In[24]:98.
WARNING: Method definition cross!(Main.Individual, Main.Individual, Main.Individual, Any, Array{Int64, 1}, Array{Int64, 1}) in module Main at In[4]:109 overwritten at In[24]:109.
WARNING: Method definition choose_mates!(Main.Population, Array{Int64, 1}) in module Main at In[4]:123 overwritten at In[24]:123.
WARNING: Method definition next_generation!(Main.Population, Main.Population, Any, Any, Any, Any, Any) in module Main at In[4]:137 overwritten at In[24]:137.
WARNING: Method definition #next_generation!(Array{Any, 1}, Main.#next_generation!, Main.Population, Main.Population, Any, Any, Any, Any, Any) in module Main overwritten.
WARNING: Method definition moments!(Main.Population) in module Main at In[4]:155 overwritten at In[24]:155.
WARNING: Method definition moments(Main.Population) in module Main at In[4]:159 overwritten at In[24]:159.
Out[24]:
moments (generic function with 1 method)

In [25]:
function lowerTri(mat)
    p = size(mat)[1]
    lower = zeros(Float64, Int64((p * p - p)/2))
    k = 1
    for i = 1:p, j = 1:p
        if(i < j) 
            lower[k]= mat[i, j]
            k = k + 1
        end
    end
    lower
end

function AVGcorr(mat; modules = Array[collect(1:Int(floor(size(mat)[1]/2))), 
                                      collect(((Int(floor(size(mat)[1]/2)))+1):size(mat)[1])])
    p = size(mat)[1]
    avg_plus = [mean(lowerTri(mat[mod, mod])) for mod in modules]
    modules_array = Array(Int64, p, p, length(modules))
    for j in 1:length(modules)
        mod_vec = [any(i .== modules[j]) ? 1 : 0 for i in 1:p]
        modules_array[:, :, j] = mod_vec * mod_vec'
    end
    avg_plus, mean(mat[find(sum(modules_array, 3) .== 0)])
end

function plotSurfacePop(pop::Population; gen = length(pop.moments))
    zs = pop.moments[gen]["zs"] 
    theta = pop.moments[gen]["theta"] 
    zs_eig = (eig(pop.surface.Σ.mat)[2][:, end-1:end]' * zs')'

    zs_df = DataFrame(x = zs_eig[:,1], y = zs_eig[:,2], fit = pop.fitness)
    sort!(zs_df,  cols = [:fit])

    s_theta = eig(pop.surface.Σ.mat)[2][:, end-1:end]' * theta
    s_omega = diagm(eig(pop.surface.Σ.mat)[1])[end-1:end,end-1:end]

    limits_x = (reduce(min, [s_theta[1] - 2*sqrt(s_omega[1, 1]); zs_eig[:,1]]),
                reduce(max, [s_theta[1] + 2*sqrt(s_omega[1, 1]); zs_eig[:,1]]))
    limits_y = (reduce(min, [s_theta[2] - 2*sqrt(s_omega[2, 2]); zs_eig[:,2]]),
                reduce(max, [s_theta[2] + 2*sqrt(s_omega[2, 2]); zs_eig[:,2]]))

    plot(layer(z = (x,y) -> pdf(MvNormal(s_theta, s_omega), [x; y]),
               x = linspace(limits_x[1], limits_x[2], 150), y = linspace(limits_y[1], limits_y[2], 150), 
               Geom.contour),
         layer(zs_df, x = "x", y = "y", color = "fit", Geom.point))
end

function plotCorrelations(pop::Population, start = 1)
    n = length(pop.moments)
    df_P = DataFrame(W_1 = [AVGcorr(pop.moments[i]["corrG"])[1][1] for i in start:n],
                     W_2 = [AVGcorr(pop.moments[i]["corrG"])[1][2] for i in start:n],
                     B = [AVGcorr(pop.moments[i]["corrG"])[2] for i in start:n],
                     G = collect(start:n) )
    plot(df_P, layer(y="W_1", x="G", Geom.line, Theme(default_color=colorant"red")),
               layer(y="W_2", x="G", Geom.line, Theme(default_color=colorant"green")),
               layer(y="B", x="G", Geom.line),
    Guide.manual_color_key("Legend", ["Within 1", "Whithin 2", "Between"], ["red", "green", "deepskyblue"]))
end

function plotAVGRatio(pop::Population, start = 1)
    n = length(pop.moments)
    df_P = DataFrame(W = [mean(AVGcorr(pop.moments[i]["corrG"])[1]) / AVGcorr(pop.moments[i]["corrG"])[2] for i in start:n],
                     G = collect(start:n) )
    plot(df_P, layer(y="W", x="G", Geom.line))
end


function plotTraits(pop::Population, start = 1)
    n = size(pop.moments)[1]
    p2 = convert(Int64, pop.base_ind.p/2)
    df_trait = DataFrame(mod1   = [mean(pop.moments[i]["mean_z"][1:p2])     for i in start:n],
                         theta1 = [mean(pop.moments[i]["theta"][1:p2])     for i in start:n],
                         mod2   = [mean(pop.moments[i]["mean_z"][(p2+1):end]) for i in start:n],
                         theta2 = [mean(pop.moments[i]["theta"][(p2+1):end]) for i in start:n],
                         gen    = collect(start:n) )
    plot(df_trait, layer(y="mod1", x="gen", Geom.line, Theme(default_color=colorant"red")),
                   layer(y="theta1", x="gen", Geom.line, Theme(default_color=colorant"darkred")),
                   layer(y="theta2", x="gen", Geom.line, Theme(default_color=colorant"darkblue")),
                   layer(y="mod2", x="gen", Geom.line),
         Guide.manual_color_key("Legend", ["Trait Module 1", "Trait Module 2", "Theta Module 1", "Theta module 2"], 
                                ["red", "deepskyblue", "darkred", "darkblue"]))
end

function plotPleitropy(pop::Population, start = 1)
    df = copy(pop.moments[start]["count_b"])
    for i in (start+1):pop.generation
        df = [df; pop.moments[i]["count_b"]]
    end
    plot(df, x = "gen", y = "count", color="class", Geom.line)
end


WARNING: Method definition lowerTri(Any) in module Main at In[5]:2 overwritten at In[25]:2.
WARNING: Method definition AVGcorr(Any) in module Main at In[5]:16 overwritten at In[25]:16.
WARNING: Method definition #AVGcorr(Array{Any, 1}, Main.#AVGcorr, Any) in module Main overwritten.
WARNING: Method definition plotSurfacePop(Main.Population) in module Main at In[5]:27 overwritten at In[25]:27.
WARNING: Method definition #plotSurfacePop(Array{Any, 1}, Main.#plotSurfacePop, Main.Population) in module Main overwritten.
WARNING: Method definition plotCorrelations(Main.Population) in module Main at In[5]:49 overwritten at In[25]:49.
WARNING: Method definition plotCorrelations(Main.Population, Any) in module Main at In[5]:49 overwritten at In[25]:49.
WARNING: Method definition plotAVGRatio(Main.Population) in module Main at In[5]:61 overwritten at In[25]:61.
WARNING: Method definition plotAVGRatio(Main.Population, Any) in module Main at In[5]:61 overwritten at In[25]:61.
WARNING: Method definition plotTraits(Main.Population) in module Main at In[5]:69 overwritten at In[25]:69.
WARNING: Method definition plotTraits(Main.Population, Any) in module Main at In[5]:69 overwritten at In[25]:69.
WARNING: Method definition plotPleitropy(Main.Population) in module Main at In[5]:85 overwritten at In[25]:85.
WARNING: Method definition plotPleitropy(Main.Population, Any) in module Main at In[5]:85 overwritten at In[25]:85.
Out[25]:
plotPleitropy (generic function with 2 methods)

In [26]:
function classifyBj(b_j::Array{Int64, 1})
    if abs(b_j) == [1, 1, 0, 0] || abs(b_j) == [0, 0, 1, 1] 
        "Modular"
    elseif b_j == [1, 1, -1, -1] || b_j == [-1, -1, 1, 1]
        "Antagonistic"
    elseif b_j == [-1, 1, 1, 1] || b_j == [1, -1, 1, 1] || b_j == [1, 1, -1, 1] || b_j == [1, 1, 1, -1]
        "Antagonistic"
    elseif b_j == [1, -1, -1, -1] || b_j == [-1, 1, -1, -1] || b_j == [-1, -1, 1, -1] || b_j == [-1, -1, -1, 1]
        "Antagonistic"
    elseif b_j == [0, 1, -1, -1] || b_j == [1, 0, -1, -1] || b_j == [1, 1, 0, -1] || b_j == [1, 1, -1, 0] 
        "Antagonistic"
    elseif b_j == [0, -1, 1, 1] || b_j == [-1, 0, 1, 1] || b_j == [-1, -1, 0, 1] || b_j == [-1, -1, 1, 0]
        "Antagonistic"
    elseif b_j == [0, -1, 0, 1] || b_j == [-1, 0, 1, 0] || b_j == [1, 0, -1, 0] || b_j == [0, -1, 0, 1]
        "Antagonistic"
    elseif b_j == [0, -1, 1, 0] || b_j == [-1, 0, 0, 1] || b_j == [1, 0, 0, -1] || b_j == [0, 1, -1, 0]
        "Antagonistic"
    elseif sum(abs(b_j)) == 1
        "Local"
    elseif b_j == [1, 1, 1, 1] || b_j == [-1, -1, -1, -1] 
        "Integrated"
    elseif b_j == [0, 1, 1, 1] || b_j == [1, 0, 1, 1] || b_j == [1, 1, 0, 1] || b_j == [1, 1, 1, 0]
        "Integrated"
    elseif b_j == [0, -1, -1, -1] || b_j == [-1, 0, -1, -1] || b_j == [-1, -1, 0, -1] || b_j == [-1, -1, -1, 0]
        "Integrated"
    elseif b_j == [0, 1, 0, 1] || b_j == [1, 0, 1, 0] || b_j == [0, 1, 1, 0] || b_j == [1, 0, 0, 1]
        "Integrated"
    elseif b_j == [0, -1, 0, -1] || b_j == [-1, 0, -1, 0] || b_j == [0, -1, -1, 0] || b_j == [-1, 0, 0, -1]
        "Integrated"
    elseif b_j == [0, 0, 0, 0]
        "Neutral"
    else
        "Other"
    end
end
function classifyBj(b::Array{Int64, 2})
    m = size(b)[2]
    [classifyBj(b[:,j]) for j in 1:m]
end
function countBclass(b_class)
    counter_b = counter(b_class)
    DataFrame( class = ["Modular", "Antagonistic", "Local", "Integrated", "Neutral", "Other"],
               count = [counter_b[i] for i in ["Modular", "Antagonistic", "Local", "Integrated", "Neutral", "Other"]])
end


WARNING: Method definition classifyBj(Array{Int64, 1}) in module Main at In[6]:2 overwritten at In[26]:2.
WARNING: Method definition classifyBj(Array{Int64, 2}) in module Main at In[6]:37 overwritten at In[26]:37.
WARNING: Method definition countBclass(Any) in module Main at In[6]:41 overwritten at In[26]:41.
Out[26]:
countBclass (generic function with 1 method)

In [27]:
function run_pop(ind::Individual,
                 n_e::Int64,
                 selectionRegimes::Array{String,1},
                 regimeGenerations::Array{Int64, 1};
                 theta::Array{Float64, 1} = zeros(Float64, ind.p), 
                 delta_theta::Array{Float64, 1} = zeros(Float64, ind.p), 
                 omega::Array{Float64, 2} = diagm(ones(Float64, ind.p)), 
                 thin = 100)
    
    pop = Population(n_e, ind, 0)
    RandomPop!(pop)

    changeSurface!(pop, theta, omega)
    fitness!(pop)
    
    run_pop(pop, selectionRegimes, regimeGenerations; 
            theta = theta, delta_theta = delta_theta, omega = omega, thin = thin)
end

function run_pop(pop::Population,
                 selectionRegimes::Array{String,1},
                 regimeGenerations::Array{Int64, 1};
                 theta::Array{Float64, 1} = pop.surface.μ, 
                 delta_theta::Array{Float64, 1} = zeros(Float64, ind.p), 
                 omega::Array{Float64, 2} = pop.surface.Σ,
                 thin = 100)
    
    changeSurface!(pop, theta, omega)
    fitness!(pop)
    
    holder_pop = copy(pop)
    sires = zeros(Int64, pop.n)
    dames = zeros(Int64, pop.n)
    d_uni = DiscreteUniform(0, 1)
    alele_1 = zeros(Int64, Int64(pop.base_ind.m/2))
    alele_2 = zeros(Int64, Int64(pop.base_ind.m/2))
    
    omega_var = omega[1, 1]

    regimes = [i for i in zip(selectionRegimes, regimeGenerations)]
    
    for current_regime in regimes
        @showprogress 2 current_regime[1] for i = 1:current_regime[2]
            if(i % thin == 0) 
                moments!(pop)
            end
            if(startswith(normalize_string(strip(current_regime[1]), casefold = true), "dire"))
                theta += delta_theta
                changeSurface!(pop, theta, omega)
            end
            next_generation!(pop, holder_pop, sires, dames, d_uni, alele_1, alele_2, 
            selective = normalize_string(strip(current_regime[1]), casefold = true) != "drift")
        end
    end
    pop
end


WARNING: Method definition run_pop(Main.Individual, Int64, Array{String, 1}, Array{Int64, 1}) in module Main at In[7]:10 overwritten at In[27]:10.
WARNING: Method definition #run_pop(Array{Any, 1}, Main.#run_pop, Main.Individual, Int64, Array{String, 1}, Array{Int64, 1}) in module Main overwritten.
WARNING: Method definition run_pop(Main.Population, Array{String, 1}, Array{Int64, 1}) in module Main at In[7]:28 overwritten at In[27]:28.
WARNING: Method definition #run_pop(Array{Any, 1}, Main.#run_pop, Main.Population, Array{String, 1}, Array{Int64, 1}) in module Main overwritten.
Out[27]:
run_pop (generic function with 2 methods)

In [28]:
function all_plots(pop::Population, file)
    plot_corr = plotCorrelations(pop, 2000)
    plot_traits = plotTraits(pop, 2000)
    plot_pleio = plotPleitropy(pop, 2000)
    draw(PDF("./figures/$file", 6inch, 6inch), vstack(plot_corr, plot_traits, plot_pleio))
end


WARNING: Method definition all_plots(Main.Population, Any) in module Main at In[17]:2 overwritten at In[28]:2.
Out[28]:
all_plots (generic function with 2 methods)

In [29]:
n_e  = 5000
m    = 100
p    = 4
mu   = 0.0005
mu_b = 0.0001
mu_σ = 0.02

global random_ind = Individual(m, p, mu, mu_b, mu_σ)

delta_theta_speed = 40/1000
delta_theta = delta_theta_speed * [ones(Float64, Int64(random_ind.p/2)); -1*ones(Float64, Int64(random_ind.p/2))]
delta_theta_1 = delta_theta_speed * [ones(Float64, Int64(random_ind.p/2)); zeros(Float64, Int64(random_ind.p/2))]
delta_theta_2 = delta_theta_speed * [zeros(Float64, Int64(random_ind.p/2)); -1 * ones(Float64, Int64(random_ind.p/2))]

omega_var = 2
omega = twoModuleMatrix(omega_var*ones(Float64, random_ind.p), 0.8, 0.8)


Out[29]:
4×4 Array{Float64,2}:
 2.0  1.6  0.0  0.0
 1.6  2.0  0.0  0.0
 0.0  0.0  2.0  1.6
 0.0  0.0  1.6  2.0

In [33]:
corridor_pop = run_pop(random_ind, n_e, ["Drift", "Stab", "Direct", "Stab"], [2000, 2000, 2000, 1000]; 
                     delta_theta = delta_theta_1, omega = omega, thin = 1)
#divergent_pop = run_pop(random_ind, n_e, ["Drift", "Stab", "Direct", "Stab"], [2000, 2000, 2000, 1000]; 
#                     delta_theta = delta_theta, omega = omega, thin = 1)
#random_pop = run_pop(random_pop, ["Stab"], [1000]; 
#                     delta_theta = delta_theta_2, omega = omega, thin = 1)
#random_pop = run_pop(random_pop, ["Direct"], [1000]; 
#                     delta_theta = delta_theta_1, omega = omega, thin = 1)
#random_pop = run_pop(random_pop, ["Direct"], [1000]; 
#                     delta_theta = delta_theta_2, omega = omega, thin = 1)


Drift100% Time: 0:05:34
Stab100% Time: 0:05:48
Direct100% Time: 0:05:56
Stab100% Time: 0:02:57
Out[33]:
This is a Population with 5000 individuals, at generation 7000

In [32]:
all_plots(corridor_pop, "corridor.pdf")
all_plots(divergent_pop, "divergent.pdf")

In [9]:
random_pop.moments[end]["corrG"]


UndefVarError: random_pop not defined

In [22]:
plot_corr = plotCorrelations(divergent_pop, 2000)


Out[22]:
G -4.0×10³ -3.0×10³ -2.0×10³ -1.0×10³ 0 1.0×10³ 2.0×10³ 3.0×10³ 4.0×10³ 5.0×10³ 6.0×10³ 7.0×10³ 8.0×10³ 9.0×10³ 1.0×10⁴ 1.1×10⁴ 1.2×10⁴ 1.3×10⁴ -3.00×10³ -2.80×10³ -2.60×10³ -2.40×10³ -2.20×10³ -2.00×10³ -1.80×10³ -1.60×10³ -1.40×10³ -1.20×10³ -1.00×10³ -8.00×10² -6.00×10² -4.00×10² -2.00×10² 0 2.00×10² 4.00×10² 6.00×10² 8.00×10² 1.00×10³ 1.20×10³ 1.40×10³ 1.60×10³ 1.80×10³ 2.00×10³ 2.20×10³ 2.40×10³ 2.60×10³ 2.80×10³ 3.00×10³ 3.20×10³ 3.40×10³ 3.60×10³ 3.80×10³ 4.00×10³ 4.20×10³ 4.40×10³ 4.60×10³ 4.80×10³ 5.00×10³ 5.20×10³ 5.40×10³ 5.60×10³ 5.80×10³ 6.00×10³ 6.20×10³ 6.40×10³ 6.60×10³ 6.80×10³ 7.00×10³ 7.20×10³ 7.40×10³ 7.60×10³ 7.80×10³ 8.00×10³ 8.20×10³ 8.40×10³ 8.60×10³ 8.80×10³ 9.00×10³ 9.20×10³ 9.40×10³ 9.60×10³ 9.80×10³ 1.00×10⁴ 1.02×10⁴ 1.04×10⁴ 1.06×10⁴ 1.08×10⁴ 1.10×10⁴ 1.12×10⁴ 1.14×10⁴ 1.16×10⁴ 1.18×10⁴ 1.20×10⁴ -5.0×10³ 0 5.0×10³ 1.0×10⁴ 1.5×10⁴ -3.00×10³ -2.50×10³ -2.00×10³ -1.50×10³ -1.00×10³ -5.00×10² 0 5.00×10² 1.00×10³ 1.50×10³ 2.00×10³ 2.50×10³ 3.00×10³ 3.50×10³ 4.00×10³ 4.50×10³ 5.00×10³ 5.50×10³ 6.00×10³ 6.50×10³ 7.00×10³ 7.50×10³ 8.00×10³ 8.50×10³ 9.00×10³ 9.50×10³ 1.00×10⁴ 1.05×10⁴ 1.10×10⁴ 1.15×10⁴ 1.20×10⁴ Within 1 Whithin 2 Between Legend -2.5 -2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 2.5 3.0 -2.00 -1.95 -1.90 -1.85 -1.80 -1.75 -1.70 -1.65 -1.60 -1.55 -1.50 -1.45 -1.40 -1.35 -1.30 -1.25 -1.20 -1.15 -1.10 -1.05 -1.00 -0.95 -0.90 -0.85 -0.80 -0.75 -0.70 -0.65 -0.60 -0.55 -0.50 -0.45 -0.40 -0.35 -0.30 -0.25 -0.20 -0.15 -0.10 -0.05 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.45 0.50 0.55 0.60 0.65 0.70 0.75 0.80 0.85 0.90 0.95 1.00 1.05 1.10 1.15 1.20 1.25 1.30 1.35 1.40 1.45 1.50 1.55 1.60 1.65 1.70 1.75 1.80 1.85 1.90 1.95 2.00 2.05 2.10 2.15 2.20 2.25 2.30 2.35 2.40 2.45 2.50 -2 0 2 4 -2.0 -1.9 -1.8 -1.7 -1.6 -1.5 -1.4 -1.3 -1.2 -1.1 -1.0 -0.9 -0.8 -0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2.0 2.1 2.2 2.3 2.4 2.5 W_1

In [216]:
plot_traits = plotTraits(random_pop, 2000)


Out[216]:
gen -2×10³ -1×10³ 0 1×10³ 2×10³ 3×10³ 4×10³ 5×10³ 6×10³ 7×10³ 8×10³ 9×10³ -1000 -900 -800 -700 -600 -500 -400 -300 -200 -100 0 100 200 300 400 500 600 700 800 900 1000 1100 1200 1300 1400 1500 1600 1700 1800 1900 2000 2100 2200 2300 2400 2500 2600 2700 2800 2900 3000 3100 3200 3300 3400 3500 3600 3700 3800 3900 4000 4100 4200 4300 4400 4500 4600 4700 4800 4900 5000 5100 5200 5300 5400 5500 5600 5700 5800 5900 6000 6100 6200 6300 6400 6500 6600 6700 6800 6900 7000 7100 7200 7300 7400 7500 7600 7700 7800 7900 8000 -5×10³ 0 5×10³ 1×10⁴ -1000 -800 -600 -400 -200 0 200 400 600 800 1000 1200 1400 1600 1800 2000 2200 2400 2600 2800 3000 3200 3400 3600 3800 4000 4200 4400 4600 4800 5000 5200 5400 5600 5800 6000 6200 6400 6600 6800 7000 7200 7400 7600 7800 8000 Trait Module 1 Trait Module 2 Theta Module 1 Theta module 2 Legend -250 -200 -150 -100 -50 0 50 100 150 200 250 300 -200 -195 -190 -185 -180 -175 -170 -165 -160 -155 -150 -145 -140 -135 -130 -125 -120 -115 -110 -105 -100 -95 -90 -85 -80 -75 -70 -65 -60 -55 -50 -45 -40 -35 -30 -25 -20 -15 -10 -5 0 5 10 15 20 25 30 35 40 45 50 55 60 65 70 75 80 85 90 95 100 105 110 115 120 125 130 135 140 145 150 155 160 165 170 175 180 185 190 195 200 205 210 215 220 225 230 235 240 245 250 -200 0 200 400 -200 -190 -180 -170 -160 -150 -140 -130 -120 -110 -100 -90 -80 -70 -60 -50 -40 -30 -20 -10 0 10 20 30 40 50 60 70 80 90 100 110 120 130 140 150 160 170 180 190 200 210 220 230 240 250 mod1

In [23]:
plot_pleio = plotPleitropy(divergent_pop, 2000)


Out[23]:
gen -1.0×10⁴ -8.0×10³ -6.0×10³ -4.0×10³ -2.0×10³ 0 2.0×10³ 4.0×10³ 6.0×10³ 8.0×10³ 1.0×10⁴ 1.2×10⁴ 1.4×10⁴ 1.6×10⁴ 1.8×10⁴ -8.00×10³ -7.50×10³ -7.00×10³ -6.50×10³ -6.00×10³ -5.50×10³ -5.00×10³ -4.50×10³ -4.00×10³ -3.50×10³ -3.00×10³ -2.50×10³ -2.00×10³ -1.50×10³ -1.00×10³ -5.00×10² 0 5.00×10² 1.00×10³ 1.50×10³ 2.00×10³ 2.50×10³ 3.00×10³ 3.50×10³ 4.00×10³ 4.50×10³ 5.00×10³ 5.50×10³ 6.00×10³ 6.50×10³ 7.00×10³ 7.50×10³ 8.00×10³ 8.50×10³ 9.00×10³ 9.50×10³ 1.00×10⁴ 1.05×10⁴ 1.10×10⁴ 1.15×10⁴ 1.20×10⁴ 1.25×10⁴ 1.30×10⁴ 1.35×10⁴ 1.40×10⁴ 1.45×10⁴ 1.50×10⁴ 1.55×10⁴ 1.60×10⁴ -1×10⁴ 0 1×10⁴ 2×10⁴ -8.00×10³ -7.50×10³ -7.00×10³ -6.50×10³ -6.00×10³ -5.50×10³ -5.00×10³ -4.50×10³ -4.00×10³ -3.50×10³ -3.00×10³ -2.50×10³ -2.00×10³ -1.50×10³ -1.00×10³ -5.00×10² 0 5.00×10² 1.00×10³ 1.50×10³ 2.00×10³ 2.50×10³ 3.00×10³ 3.50×10³ 4.00×10³ 4.50×10³ 5.00×10³ 5.50×10³ 6.00×10³ 6.50×10³ 7.00×10³ 7.50×10³ 8.00×10³ 8.50×10³ 9.00×10³ 9.50×10³ 1.00×10⁴ 1.05×10⁴ 1.10×10⁴ 1.15×10⁴ 1.20×10⁴ 1.25×10⁴ 1.30×10⁴ 1.35×10⁴ 1.40×10⁴ 1.45×10⁴ 1.50×10⁴ 1.55×10⁴ 1.60×10⁴ Modular Antagonistic Local Integrated Neutral Other class