Add more tests and fix bugs

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

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)

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

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]

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

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.

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