From 3dc1e0f8e6c42387cb7ea93a95e75a92fbc409ab Mon Sep 17 00:00:00 2001 From: "Ali R. Vahdati" Date: Wed, 31 Aug 2022 17:50:31 +0200 Subject: [PATCH] Add more tests and fix bugs --- src/EvoDynamics.jl | 5 ++ src/api.jl | 2 - src/interactions.jl | 20 ++++++- src/reproduction.jl | 6 +- src/simulation.jl | 4 +- test/params.jl | 4 +- test/simulation_tests.jl | 116 +++++++++++++++++++++++++++++++++++++-- 7 files changed, 140 insertions(+), 17 deletions(-) diff --git a/src/EvoDynamics.jl b/src/EvoDynamics.jl index c657289..4bc22a6 100644 --- a/src/EvoDynamics.jl +++ b/src/EvoDynamics.jl @@ -9,6 +9,11 @@ using StaticArrays using YAML using RandomNumbers +export runmodel, model_initiation, load_parameters +# export abiotic_survive!, burn_energy!, consume_food!, interact!, remove_agent!, reproduce!, migrate!, target_species_ids, eat!, abiotic_distance, abiotic_fitness, interaction_power, phenotypic_distance, get_abiotic_phenotype, get_biotic_phenotype, update_fitness!, vertex2coord, coord2vertex, mutate!, adjust_fitness!, create_gamete, crossing_overs, return_opt_phenotype +# export check_site, pick_site, get_migration_trait +# export step! + include("simulation.jl") include("interactions.jl") include("migration.jl") diff --git a/src/api.jl b/src/api.jl index c3a5cb3..cfbce7f 100644 --- a/src/api.jl +++ b/src/api.jl @@ -1,5 +1,3 @@ -export runmodel - "Returns a tuple whose entries are the mean fitness of each species." function mean_fitness_per_species(model::ABM) mean_fitness = Array{Float64}(undef, model.nspecies) diff --git a/src/interactions.jl b/src/interactions.jl index 9ecc00e..cd5f728 100644 --- a/src/interactions.jl +++ b/src/interactions.jl @@ -32,7 +32,7 @@ Distance is a probability and is calculated by the number of standard deviations Variance can stand for selection coefficient, with a larger variance, less sever the difference. """ -function phenotypic_distance(phenotype1, phenotype2, interaction_coeff, variance::Real) +function phenotypic_distance_different_species(phenotype1, phenotype2, interaction_coeff, variance::Real) distance = 0.0 counter = 0.0 for ph1 in phenotype1 @@ -50,7 +50,23 @@ function phenotypic_distance(phenotype1, phenotype2, interaction_coeff, variance return distance end -phenotypic_distance(ag1::Ind, ag2::Ind, model) = phenotypic_distance(ag1.biotic_phenotype, ag2.biotic_phenotype, model.interactions[ag1.species, ag2.species], model.biotic_variances[ag1.species]) +function phenotypic_distance_same_species(phenotype1, phenotype2, interaction_coeff, variance::Real) + distance = 0.0 + counter = 0.0 + for (ph1, ph2) in zip(phenotype1, phenotype2) + d = abs(0.5 - cdf(Normal(ph1, variance), ph2)) # 0.5 is the value when ph1 and ph2 are identical + distance += d + counter += 1.0 + end + distance += distance # Double the distance so that it is in range 0 to 1. + distance /= counter # average distance + if interaction_coeff < 0 + distance = 1.0 - distance + end + return distance +end + +phenotypic_distance(ag1::Ind, ag2::Ind, model) = ag1.species == ag2.species ? phenotypic_distance_same_species(ag1.biotic_phenotype, ag2.biotic_phenotype, model.interactions[ag1.species, ag2.species], model.biotic_variances[ag1.species]) : phenotypic_distance_different_species(ag1.biotic_phenotype, ag2.biotic_phenotype, model.interactions[ag1.species, ag2.species], model.biotic_variances[ag1.species]) """ Abiotic distance to the optimal values diff --git a/src/reproduction.jl b/src/reproduction.jl index ced914b..578e1d1 100644 --- a/src/reproduction.jl +++ b/src/reproduction.jl @@ -81,7 +81,7 @@ function create_one_offspring(ag1::Ind, ag2::Ind, model::ABM) W = adjust_fitness(W, ag1.species, model) initial_energy = model.initial_energy[species] - offspring = add_agent!(ag1.pos, model, ag1.species, bph, abph, episMat, pleioMat, q, 1, sex, interaction_history, initial_energy, W, true, 0, -1) + offspring = add_agent!(ag1.pos, model, ag1.species, deepcopy(bph), deepcopy(abph), deepcopy(episMat), deepcopy(pleioMat), deepcopy(q), 1, sex, deepcopy(interaction_history), initial_energy, W, true, 0, -1) return offspring end @@ -93,7 +93,7 @@ function reproduce!(ag1::Ind, ag2::Ind, model::ABM) if model.mating_schemes[ag1.species] == 0.0 reproduction_success = 1.0 else - d = phenotypic_distance(ag1.biotic_phenotype, ag2.biotic_phenotype, model.mating_schemes[ag1.species], model.biotic_variances[ag1.species]) + d = phenotypic_distance(ag1, ag2, model) reproduction_success = interaction_power(ag1, ag2, d, model) end growth_rate = model.growthrates[ag1.species] @@ -117,7 +117,7 @@ function reproduce!(agent::Ind, model::ABM) interaction_history = deepcopy(agent.interaction_history) interaction_history[1:end] .= 0 for c in 1:nchildren - offspring = add_agent!(agent.pos, model, agent.species, agent.biotic_phenotype, agent.abiotic_phenotype, agent.epistasisMat, agent.pleiotropyMat, agent.q, 1, agent.sex, interaction_history, model.initial_energy[agent.species], agent.W, true, 0, -1) + offspring = add_agent!(agent.pos, model, agent.species, deepcopy(agent.biotic_phenotype), deepcopy(agent.abiotic_phenotype), deepcopy(agent.epistasisMat), deepcopy(agent.pleiotropyMat), deepcopy(agent.q), 1, agent.sex, interaction_history, model.initial_energy[agent.species], agent.W, true, 0, -1) mutate!(offspring, model) end elseif model.ploidy[agent.species] == 2 && agent.interaction_history[agent.species] == model.step[1] && agent.time_met_other_sex == model.step[1] diff --git a/src/simulation.jl b/src/simulation.jl index ca7f75a..d1ec617 100644 --- a/src/simulation.jl +++ b/src/simulation.jl @@ -105,7 +105,7 @@ function model_initiation(dd) end interaction_history = MVector{model.nspecies,Int}(fill(-1, model.nspecies)) initial_energy = model.initial_energy[sp] - add_agent!(model.nodes[pos], model, sp, biotic_ph, abiotic_ph, epistasisMat[sp], pleiotropyMat[sp], expressionArrays[sp], 1, sex, interaction_history, initial_energy, W, true, 0, -1) + add_agent!(model.nodes[pos], model, sp, deepcopy(biotic_ph), deepcopy(abiotic_ph), deepcopy(epistasisMat[sp]), deepcopy(pleiotropyMat[sp]), deepcopy(expressionArrays[sp]), 1, sex, deepcopy(interaction_history), initial_energy, W, true, 0, -1) end end end @@ -302,7 +302,7 @@ function mutate!(agent::Ind, model::ABM) end # mutate pleiotropy matrix if rand(model.rng, model.mutProbs[agent.species][2]) - randnumbers = rand(model.rng, model.rng, model.mutMagnitudes[agent.species][2], size(agent.pleiotropyMat)) + randnumbers = rand(model.rng, model.mutMagnitudes[agent.species][2], size(agent.pleiotropyMat)) agent.pleiotropyMat[randnumbers] .= .!agent.pleiotropyMat[randnumbers] mutated = true end diff --git a/test/params.jl b/test/params.jl index 5dfde6e..9d5f3e6 100644 --- a/test/params.jl +++ b/test/params.jl @@ -32,8 +32,8 @@ species1 = Dict( :selection_coefficient => 0.02, :abiotic_variance => 1.0, # variance of a normal distribution used in determining the phenotypic distance of agents to the optimal environmental phenotypes. The larger the variance, the less important is the distance. :biotic_variance => 10.0, # same as above but for determining the biotic phenotypic distance between two individuals (used in any kind of interaction). Larger values mean that all pairs are equally likely to interact, regardless of their phenotype difference. - :mutation_probabilities => [0.9, 0.0, 0.0], # for gene expression array, pleiotropy matrix and epistasis matrix, respectively - :mutation_magnitudes => [0.05, 0.0, 0.01], # same as above + :mutation_probabilities => [0.99, 0.99, 0.99], # for gene expression array, pleiotropy matrix and epistasis matrix, respectively + :mutation_magnitudes => [0.01, 0.1, 0.1], # same as above :N => [1000, 0, 0, 0], # number of individuals per site at time 0 :environmental_noise => 0.0, # variance of a normal distribution with mean 0 # each row is the optimal phenotypes for each site for all abiotic traits. There are as many element as number of sites times number of abiotic traits. The first N elements are for the first abiotic trait, where N is the number of sites, and so on. diff --git a/test/simulation_tests.jl b/test/simulation_tests.jl index 1df1654..55cde5d 100644 --- a/test/simulation_tests.jl +++ b/test/simulation_tests.jl @@ -1,6 +1,33 @@ +@testset "Basic functions" begin + dd = EvoDynamics.load_parameters(param_file) + model = EvoDynamics.model_initiation(dd) + EvoDynamics.remove_agent!(model[1], model) + @test !haskey(model.agents, 1) == true + + EvoDynamics.consume_food!(model[2], model) + @test model.resources == [1999 1980; 1830 1900] + + @test EvoDynamics.coord2vertex((2, 1), model) == 2 + @test EvoDynamics.coord2vertex((1, 2), model) == 3 + @test EvoDynamics.coord2vertex((2, 2), model) == 4 + @test EvoDynamics.vertex2coord(2, model) == (2, 1) + @test EvoDynamics.vertex2coord(3, model) == (1, 2) + @test EvoDynamics.vertex2coord(4, model) == (2, 2) + + @test model[2].W == model[3].W + model[2].q[1] -= 0.05 + EvoDynamics.update_fitness!(model[2], model) + @test model[2].W != model[3].W +end + @testset "Mutation" begin dd = EvoDynamics.load_parameters(param_file) model = EvoDynamics.model_initiation(dd) + + @test EvoDynamics.nagents(model) == 1100 + @test size(model.space) == (2, 2) + @test length(EvoDynamics.ids_in_position((1, 1), model)) == EvoDynamics.nagents(model) + ag1 = model.agents[1] ag2 = model.agents[2] @@ -18,14 +45,91 @@ EvoDynamics.update_fitness!(ag1, model) @test ag1.W != Worg + # no leakage. using deepcopy for creating new agents + @test model[3].q == model[4].q + @test model[3].pleiotropyMat == model[4].pleiotropyMat + @test model[3].epistasisMat == model[4].epistasisMat + EvoDynamics.mutate!(model[3], model) + @test model[3].q != model[4].q + @test model[3].pleiotropyMat != model[4].pleiotropyMat + @test model[3].epistasisMat != model[4].epistasisMat + @test model[5].q == model[4].q + @test model[5].pleiotropyMat == model[4].pleiotropyMat + @test model[5].epistasisMat == model[4].epistasisMat end -# #TODO -# @testset "Migration" begin +@testset "Seeding" begin + dd = EvoDynamics.load_parameters(param_file) + dd[:seed] = 1 + model = model_initiation(dd) + EvoDynamics.step!(model, EvoDynamics.agent_step!, 1) + + dd2 = EvoDynamics.load_parameters(param_file) + dd2[:seed] = 2 + model2 = EvoDynamics.model_initiation(dd2) + EvoDynamics.step!(model2, EvoDynamics.agent_step!, 1) -# end + dd3 = EvoDynamics.load_parameters(param_file) + dd3[:seed] = 1 + model3 = EvoDynamics.model_initiation(dd3) + EvoDynamics.step!(model3, EvoDynamics.agent_step!, 1) -# #TODO -# @testset "Reproduction" begin + @test collect(keys(model.agents)) != collect(keys(model2.agents)) + @test collect(keys(model.agents)) == collect(keys(model3.agents)) +end -# end +@testset "Interactions" begin + dd = EvoDynamics.load_parameters(param_file) + model = EvoDynamics.model_initiation(dd) + + agent = model[1] + bph = EvoDynamics.get_biotic_phenotype(agent, model) + aph = EvoDynamics.get_abiotic_phenotype(agent, model) + + @test length(bph) == length(dd[:species][1][:biotic_phenotypes]) + @test length(aph) == length(dd[:species][1][:abiotic_phenotypes]) + + # identical phenotypes in the same species results in the same distance: + agent2 = model[1001] + d = EvoDynamics.phenotypic_distance(agent, agent2, model) + dsame = EvoDynamics.phenotypic_distance(agent, agent, model) + dsame2 = EvoDynamics.phenotypic_distance(agent2, agent2, model) + @test dsame == dsame2 == 0 + + # identical phenotypes in two different species results in non-zero distance: + agent.biotic_phenotype .= agent2.biotic_phenotype + d2 = EvoDynamics.phenotypic_distance(agent, agent2, model) + @test d2 > 0.0 + @test EvoDynamics.phenotypic_distance(agent, agent2, model) == EvoDynamics.phenotypic_distance(agent2, agent, model) + + optphen = EvoDynamics.return_opt_phenotype(agent2.species, agent.pos, model) + @test EvoDynamics.abiotic_distance(agent2.abiotic_phenotype, optphen, model.abiotic_variances[agent2.species]) > 0 + agent2.abiotic_phenotype .= optphen + @test EvoDynamics.abiotic_distance(agent2.abiotic_phenotype, optphen, model.abiotic_variances[agent2.species]) == 0 + @test EvoDynamics.abiotic_fitness(agent2, model) == 1.0 + + EvoDynamics.interact!(agent, agent2, model) + @test agent.isalive == false + + @test length(EvoDynamics.target_species_ids(agent, model)) == 2 + @test length(EvoDynamics.target_species_ids(agent2, model)) == 2 + + # reproduction + nagentsbefore = EvoDynamics.nagents(model) + EvoDynamics.reproduce!(agent2, model) + EvoDynamics.reproduce!(agent2, model) + EvoDynamics.reproduce!(agent2, model) + nagentsafter = EvoDynamics.nagents(model) + @test nagentsbefore < nagentsafter + + agent3 = model[10] + agent4 = model[11] + agent3.sex = true + agent4.sex = false + + nagentsbefore = EvoDynamics.nagents(model) + EvoDynamics.reproduce!(agent3, agent4, model) + EvoDynamics.reproduce!(agent3, agent4, model) + nagentsafter = EvoDynamics.nagents(model) + @test nagentsbefore < nagentsafter +end