function run_evolution(n_gen::Int64,evo_config::EMAPConfig,cell_config::CellularConfig,
        exp_config::ExperimentalConfig,grn_config::FGRNConfig,save_mid_run::Bool)

    n_species = 16
    n_families = 4
    
    fittest = Int(floor((evo_config.adult_niche_size*evo_config.elitism_percentage)))
    
    map_float = [[zeros(evo_config.n_total_genes,8)] for i in 1:n_species]
    map_int = [[vec(vcat(evo_config.genome_start,rand(2:4,evo_config.n_total_genes - evo_config.genome_preserve,1)))] for i in 1:n_species]
    map_fitnesses = [[0.] for i in 1:n_species]
    map_families = zeros(Int,n_species)
    map_lineage = [[[0,0,0]] for i in 1:n_species]
    map_discovered = [zeros(Int64,n_species)] 
    
    fitness_records_cells = [[0.] for i in 1:n_species]
    fitness_records_families = [[0.] for i in 1:n_families]
    
    for i in 1:n_species
        map_families[i] = sum(digits(i - 1, base=2,pad = 4) |> reverse)
    end
    
    valid_species = findall(x->x in evo_config.families,map_families)

    pop_fingerprints = zeros(Int,evo_config.initial_pop_size)
    pop_families = zeros(Int,evo_config.initial_pop_size) 
    pop_fitnesses = zeros(Float64,evo_config.initial_pop_size)

    pop_float = [zeros(Float64,evo_config.n_total_genes,8) for i in 1:evo_config.initial_pop_size]
    pop_int = [zeros(Int64,evo_config.n_total_genes) for i in 1:evo_config.initial_pop_size]

    for i in 1:evo_config.initial_pop_size
        pop_float[i] = rand(evo_config.n_total_genes,8)
        pop_int[i] = vec(vcat(evo_config.genome_start,rand(2:4,evo_config.n_total_genes - evo_config.genome_preserve,1)))
    end

    @sync for i in 1:evo_config.initial_pop_size

        @spawn pop_fitnesses[i],pop_fingerprints[i] = evaluate_individual_fitness(pop_float[i],
            pop_int[i],evo_config,cell_config,exp_config,grn_config)
    end

    for cell in valid_species
        candidates = findall(x-> x==cell, pop_fingerprints)

        if size(candidates,1) > evo_config.adult_niche_size

            best = partialsortperm(pop_fitnesses[candidates], 1:evo_config.adult_niche_size, rev=true)

            map_fitnesses[cell] = pop_fitnesses[candidates[best]]
            map_float[cell] = pop_float[candidates[best]]
            map_int[cell] = pop_int[candidates[best]]
            map_lineage[cell] = [[17,17,0] for i in 1:evo_config.adult_niche_size]

        elseif size(candidates,1) != 0

            map_fitnesses[cell] = pop_fitnesses[candidates]
            map_float[cell] = pop_float[candidates]
            map_int[cell] = pop_int[candidates]
            map_lineage[cell] = [[17,17,0] for i in 1:size(candidates,1)]

        else

            map_fitnesses[cell] = []
            map_float[cell] = []
            map_int[cell] = []
            map_lineage[cell] = []
        end

        if size(map_fitnesses[cell],1) > 0
            push!(fitness_records_cells[cell], mean(map_fitnesses[cell]))
        else
            push!(fitness_records_cells[cell], 0.) 
        end

    end

    for i in 1:n_families
        fitness_records_families[i] = mean(fitness_records_cells[map_families .== i])
    end


    @printf "\nFinished initial evaluation"
    
    for g in 1:n_gen
        

        if (save_mid_run) && (g % 10 == 0)
            save("EMAPf_mid_run_backup_map_float.jld", "data", map_float)
            save("EMAPf_mid_run_backup_map_int.jld", "data", map_int)
            save("EMAPf_mid_run_backup_map_fit.jld", "data", map_fitnesses) 
            save("EMAPf_mid_run_backup_map_fam.jld", "data", map_families) 
            save("EMAPf_mid_run_backup_map_lin.jld", "data", map_lineage) 
            save("EMAPf_mid_run_backup_map_disc.jld", "data", map_discovered)
            save("EMAPf_mid_run_backup_fitness_records_families.jld", "data", fitness_records_families)
            save("EMAPf_mid_run_backup_fitness_records_cells.jld", "data", fitness_records_cells) 
        end

        @printf "\nFitness generation %i : 1-family %g , 2-family %g , 3-family %g  , 4-family %g" g fitness_records_families[1][end] fitness_records_families[2][end] fitness_records_families[3][end] fitness_records_families[4][end]
            

        child_map_float = []
        child_map_int = []
        child_map_lineage = []
        
        niche_operating = findall(x-> size(x,1) > 0, map_float)
            
        for cell in valid_species
            
            if size(map_float[cell],1) > 0
                
                candidates = findall(x-> (map_families[x] == map_families[cell]) && (x in niche_operating) && (x != cell) , [i for i in 1:16])
            
                for s in 1:evo_config.child_map_size

                    if (size(candidates,1) > 1) && rand() < evo_config.inter_species_crossover # 20% chance inter-species crossover - elitist
                        
                        choice = sample(candidates,1)[1]

                        best_niche = argmax(map_fitnesses[cell])
                        best_choice = argmax(map_fitnesses[choice])

                        ci, cf  = create_child(map_int[cell][best_niche], map_float[cell][best_niche],map_int[choice][best_choice],map_float[choice][best_choice])

                        push!(child_map_float,cf)
                        push!(child_map_int,ci)
                        push!(child_map_lineage,[cell,choice,g])
                        
                    else
                            
                        if size(map_float[cell],1) == 1 # mutate if not enough for cross-over
                            
                            ci, cf = create_child(map_int[cell][1],map_float[cell][1],evo_config)
                            
                            push!(child_map_float,cf)
                            push!(child_map_int,ci)
                            push!(child_map_lineage,[cell,cell,g])
                            
                        elseif rand() < evo_config.normal_crossover # otherwise 80% chance of crossover
                        
                            if size(map_float[cell],1) <= fittest

                                parents = sample([i for i in 1:size(map_float[cell],1)], 2, replace = false)

                                ci, cf  = create_child(map_int[cell][parents[1]], map_float[cell][parents[1]],map_int[cell][parents[2]],map_float[cell][parents[2]])

                                push!(child_map_float,cf)
                                push!(child_map_int,ci)
                                push!(child_map_lineage,[cell,cell,g])

                            else
                                
                                if rand() < evo_config.non_elite_crossover # 20% chance of non-elitist crossover
                                    parents = sample([i for i in 1:size(map_float[cell],1)], 2, replace = false)

                                    ci, cf  = create_child(map_int[cell][parents[1]], map_float[cell][parents[1]],map_int[cell][parents[2]],map_float[cell][parents[2]])

                                    push!(child_map_float,cf)
                                    push!(child_map_int,ci)
                                    push!(child_map_lineage,[cell,cell,g])
                                else # 80% chance of elitist crossover
                                    best = partialsortperm(map_fitnesses[cell], 1:fittest, rev=true)
                                    parents = sample(best, 2, replace = false)

                                    ci, cf  = create_child(map_int[cell][parents[1]], map_float[cell][parents[1]],map_int[cell][parents[2]],map_float[cell][parents[2]])

                                    push!(child_map_float,cf)
                                    push!(child_map_int,ci)
                                    push!(child_map_lineage,[cell,cell,g])
                                end

                            end
                            
                        else # 20% chance of mutation - not elitist
                            choice = rand(1:size(map_float[cell],1)) 
                            
                            ci, cf = create_child(map_int[cell][choice],map_float[cell][choice],evo_config)
                            
                            push!(child_map_float,cf)
                            push!(child_map_int,ci)
                            push!(child_map_lineage,[cell,cell,g])
                        end
                            
                    end
                end
            end
        end
                    
        child_map_fitnesses = zeros(Float64,size(child_map_float,1))
        child_map_fingerprints = zeros(Int64,size(child_map_float,1))
            
        @sync for child_id in 1:size(child_map_float,1)
            
            @spawn child_map_fitnesses[child_id],child_map_fingerprints[child_id] = evaluate_individual_fitness(child_map_float[child_id],
                child_map_int[child_id],evo_config,cell_config,exp_config,grn_config)
            
            
        end
        @printf "\n"
        for cell in valid_species
                
            candidates = findall(x-> x==cell, child_map_fingerprints)
                         
            if size(map_float[cell],1) >= evo_config.adult_niche_size
                
                @printf "%i : C | " cell
                
                for child_id in candidates
                    weakest = argmin(map_fitnesses[cell])    
                    if child_map_fitnesses[child_id] > map_fitnesses[cell][weakest]
                        map_fitnesses[cell][weakest] = child_map_fitnesses[child_id]
                        map_float[cell][weakest] = child_map_float[child_id]
                        map_int[cell][weakest] = child_map_int[child_id]
                        map_lineage[cell][weakest] = child_map_lineage[child_id]
                    end
                end
                            
            elseif size(candidates,1) != 0
                
                @printf "%i : F | " cell

                for child_id in candidates
                    if (size(map_float[cell],1) < evo_config.adult_niche_size)
                        push!(map_fitnesses[cell],child_map_fitnesses[child_id])
                        push!(map_float[cell],child_map_float[child_id])
                        push!(map_int[cell],child_map_int[child_id])
                        push!(map_lineage[cell],child_map_lineage[child_id])
                    end
                end
                        
            else
                nothing           
            end
                
            if size(map_fitnesses[cell],1) > 0
                push!(fitness_records_cells[cell], mean(map_fitnesses[cell]))
            else
                push!(fitness_records_cells[cell], 0.) 
            end
        end
            
        
        for i in 1:n_families
            fitness_records_families[i] = mean(fitness_records_cells[map_families .== i])
        end
        
        discovered = zeros(Int,n_species)

        discovered[niche_operating] .= 1
        
        push!(map_discovered,discovered)
        
    end
    
    results = Dict{String, Any}("EMAP_genomes_float" => map_float,
        "EMAP_genomes_int" => map_int,
        "EMAP_genomes_species_fitnesses" => map_fitnesses,
        "EMAP_genomes_family_id_map" => map_families,
        "EMAP_genomes_lineage" => map_lineage,
        "EMAP_genomes_discovered" => map_discovered,
        "EMAP_genomes_family_fitness_records" => fitness_records_families,
        "EMAP_genomes_species_fitness_records" => fitness_records_cells
    )

    return results
end

# An in-place version of MAP-Flags - this allows evolution to continue on existing populations 

function run_evolution!(n_gen::Int64,fitness_records_cells::Vector{Vector{Float64}},
        fitness_records_families::Vector{Vector{Float64}},map_float::Vector{Vector{Matrix{Float64}}},map_int::Vector{Vector{Vector{Int64}}},
        map_fitnesses::Vector{Vector{Float64}},map_families::Vector{Int64},map_lineage::Vector{Vector{Vector{Int64}}},
        map_discovered::Vector{Vector{Int64}},evo_config::EMAPConfig,cell_config::CellularConfig,
        exp_config::ExperimentalConfig,grn_config::FGRNConfig,save_mid_run::Bool)

    n_species = 16
    n_families = 4
    
    fittest = Int(floor((evo_config.adult_niche_size*evo_config.elitism_percentage)))
    
    valid_species = findall(x->x in evo_config.families,map_families)
    
    for g in 1:n_gen
        
        if (save_mid_run) && (g % 10 == 0)
            save("EMAPf_mid_run_backup_map_float.jld", "data", map_float)
            save("EMAPf_mid_run_backup_map_int.jld", "data", map_int)
            save("EMAPf_mid_run_backup_map_fit.jld", "data", map_fitnesses) 
            save("EMAPf_mid_run_backup_map_fam.jld", "data", map_families) 
            save("EMAPf_mid_run_backup_map_lin.jld", "data", map_lineage) 
            save("EMAPf_mid_run_backup_map_disc.jld", "data", map_discovered)
            save("EMAPf_mid_run_backup_fitness_records_families.jld", "data", fitness_records_families)
            save("EMAPf_mid_run_backup_fitness_records_cells.jld", "data", fitness_records_cells) 
        end

        @printf "\nFitness generation %i : 1-family %g , 2-family %g , 3-family %g  , 4-family %g" g fitness_records_families[1][end] fitness_records_families[2][end] fitness_records_families[3][end] fitness_records_families[4][end]
            

        child_map_float = []
        child_map_int = []
        child_map_lineage = []
        
        niche_operating = findall(x-> size(x,1) > 0, map_float)
            
        for cell in valid_species
            
            if size(map_float[cell],1) > 0
                
                candidates = findall(x-> (map_families[x] == map_families[cell]) && (x in niche_operating) && (x != cell) , [i for i in 1:16])
            
                for s in 1:evo_config.child_map_size

                    if (size(candidates,1) > 1) && rand() < evo_config.inter_species_crossover # 20% chance inter-species crossover - elitist
                        
                        choice = sample(candidates,1)[1]

                        best_niche = argmax(map_fitnesses[cell])
                        best_choice = argmax(map_fitnesses[choice])

                        ci, cf  = create_child(map_int[cell][best_niche], map_float[cell][best_niche],map_int[choice][best_choice],map_float[choice][best_choice])

                        push!(child_map_float,cf)
                        push!(child_map_int,ci)
                        push!(child_map_lineage,[cell,choice,g])
                        
                    else
                            
                        if size(map_float[cell],1) == 1 # mutate if not enough for cross-over
                            
                            ci, cf = create_child(map_int[cell][1],map_float[cell][1],evo_config)
                            
                            push!(child_map_float,cf)
                            push!(child_map_int,ci)
                            push!(child_map_lineage,[cell,cell,g])
                            
                        elseif rand() < evo_config.normal_crossover # otherwise 80% chance of crossover
                        
                            if size(map_float[cell],1) <= fittest

                                parents = sample([i for i in 1:size(map_float[cell],1)], 2, replace = false)

                                ci, cf  = create_child(map_int[cell][parents[1]], map_float[cell][parents[1]],map_int[cell][parents[2]],map_float[cell][parents[2]])

                                push!(child_map_float,cf)
                                push!(child_map_int,ci)
                                push!(child_map_lineage,[cell,cell,g])

                            else
                                
                                if rand() < evo_config.non_elite_crossover # 20% chance of non-elitist crossover
                                    parents = sample([i for i in 1:size(map_float[cell],1)], 2, replace = false)

                                    ci, cf  = create_child(map_int[cell][parents[1]], map_float[cell][parents[1]],map_int[cell][parents[2]],map_float[cell][parents[2]])

                                    push!(child_map_float,cf)
                                    push!(child_map_int,ci)
                                    push!(child_map_lineage,[cell,cell,g])
                                else # 80% chance of elitist crossover
                                    best = partialsortperm(map_fitnesses[cell], 1:fittest, rev=true)
                                    parents = sample(best, 2, replace = false)

                                    ci, cf  = create_child(map_int[cell][parents[1]], map_float[cell][parents[1]],map_int[cell][parents[2]],map_float[cell][parents[2]])

                                    push!(child_map_float,cf)
                                    push!(child_map_int,ci)
                                    push!(child_map_lineage,[cell,cell,g])
                                end

                            end
                            
                        else # 20% chance of mutation - not elitist
                            choice = rand(1:size(map_float[cell],1)) 
                            
                            ci, cf = create_child(map_int[cell][choice],map_float[cell][choice],evo_config)
                            
                            push!(child_map_float,cf)
                            push!(child_map_int,ci)
                            push!(child_map_lineage,[cell,cell,g])
                        end
                            
                    end
                end
            end
        end
                    
        child_map_fitnesses = zeros(Float64,size(child_map_float,1))
        child_map_fingerprints = zeros(Int64,size(child_map_float,1))
            
        @sync for child_id in 1:size(child_map_float,1)
            
            @spawn child_map_fitnesses[child_id],child_map_fingerprints[child_id] = evaluate_individual_fitness(child_map_float[child_id],
                child_map_int[child_id],evo_config,cell_config,exp_config,grn_config)
            
            
        end
        @printf "\n"
        for cell in valid_species
                
            candidates = findall(x-> x==cell, child_map_fingerprints)
                         
            if size(map_float[cell],1) >= evo_config.adult_niche_size
                
                @printf "%i : C | " cell
                
                for child_id in candidates
                    weakest = argmin(map_fitnesses[cell])    
                    if child_map_fitnesses[child_id] > map_fitnesses[cell][weakest]
                        map_fitnesses[cell][weakest] = child_map_fitnesses[child_id]
                        map_float[cell][weakest] = child_map_float[child_id]
                        map_int[cell][weakest] = child_map_int[child_id]
                        map_lineage[cell][weakest] = child_map_lineage[child_id]
                    end
                end
                            
            elseif size(candidates,1) != 0
                
                @printf "%i : F | " cell

                for child_id in candidates
                    if (size(map_float[cell],1) < evo_config.adult_niche_size)
                        push!(map_fitnesses[cell],child_map_fitnesses[child_id])
                        push!(map_float[cell],child_map_float[child_id])
                        push!(map_int[cell],child_map_int[child_id])
                        push!(map_lineage[cell],child_map_lineage[child_id])
                    end
                end
                        
            else
                nothing           
            end
                
            if size(map_fitnesses[cell],1) > 0
                push!(fitness_records_cells[cell], mean(map_fitnesses[cell]))
            else
                push!(fitness_records_cells[cell], 0.) 
            end
        end
            
        
        for i in 1:n_families
            fitness_records_families[i] = mean(fitness_records_cells[map_families .== i])
        end
        
        discovered = zeros(Int,n_species)

        discovered[niche_operating] .= 1
        
        push!(map_discovered,discovered)
        
    end

end

function run_evolution(n_gen::Int,evo_config::POOLConfig,cell_config::CellularConfig,
    exp_config::ExperimentalConfig,grn_config::FGRNConfig)

    @printf "\rInitialise"
        
    shared_pool_float::Vector{Matrix{Float64}} = []
    shared_pool_int::Vector{Vector{Int64}} = []
    shared_pool_fitnesses::Vector{Float64} = []

    for i in 1:evo_config.n_pools
        push!(shared_pool_float,rand(evo_config.n_total_genes,8))
        push!(shared_pool_int,vec(vcat(evo_config.genome_start,rand(2:4,evo_config.n_total_genes - evo_config.genome_preserve,1))))
        push!(shared_pool_fitnesses,0)
    end

    fitness_records::Vector{Vector{Float64}} = []
    best_float::Vector{Matrix{Float64}} = []
    best_int::Vector{Vector{Int64}} = []
    best_fitnesses::Vector{Float64} = []

    for i in 1:evo_config.n_pools
        push!(fitness_records,zeros(Float64,n_gen))
        push!(best_float,rand(evo_config.n_total_genes,8))
        push!(best_int,vec(vcat(evo_config.genome_start,rand(2:4,evo_config.n_total_genes - evo_config.genome_preserve,1))))
        push!(best_fitnesses,0)
    end
    
    @printf "\rCreate Pools"

    Threads.@threads for i in 1:evo_config.n_pools
        
        adult_pop_float = [rand(evo_config.n_total_genes,8) for i in 1:evo_config.adult_pop_size]
        adult_pop_int = [vec(vcat(evo_config.genome_start,rand(2:4,evo_config.n_total_genes - evo_config.genome_preserve,1))) for i in 1:evo_config.adult_pop_size]
        adult_fitnesses = zeros(evo_config.adult_pop_size); 
        
        child_pop_float = [zeros(evo_config.n_total_genes,8) for i in 1:evo_config.child_pop_size]
        child_pop_int = [zeros(Int,evo_config.n_total_genes) for i in 1:evo_config.child_pop_size]
        child_fitnesses = zeros(evo_config.child_pop_size)

        adult_ages = zeros(evo_config.adult_pop_size)
        adult_fitnesses = zeros(evo_config.adult_pop_size)

        topk = Int64(floor(evo_config.adult_pop_size*evo_config.elitism_percentage))
        topm = Int64(floor(evo_config.adult_pop_size*evo_config.elitism_migration_percentage))

        for i in 1:evo_config.adult_pop_size
            
            adult_fitnesses[i] = evaluate_individual_fitness(adult_pop_float[i],adult_pop_int[i],evo_config,cell_config,
                                                                    exp_config,grn_config)

            adult_ages[i] = 1
        end

        av_fitness = zeros(n_gen)

        for g in 1:n_gen
            
            if Threads.threadid() == 1

                @printf "\rGeneration %i" g
            end
            
            av_fitness[g] = mean(adult_fitnesses)

            @. adult_ages = adult_ages + 1
            
            if rand() < evo_config.emmigration_rate
                
                migrate = partialsortperm(adult_fitnesses, 1:topm, rev=true)
                traveller = sample(migrate,1)[1]
                
                shared_pool_float[Threads.threadid()] = adult_pop_float[traveller]
                shared_pool_int[Threads.threadid()] = adult_pop_int[traveller]
                shared_pool_fitnesses[Threads.threadid()] = adult_fitnesses[traveller]
                
                adult_pop_float[traveller] = rand(evo_config.n_total_genes,8)
                
                adult_pop_int[traveller] = vec(vcat(evo_config.genome_start,rand(2:4,evo_config.n_total_genes - evo_config.genome_preserve,1)))
                
                adult_fitnesses[traveller] = evaluate_individual_fitness(adult_pop_float[traveller],
                    adult_pop_int[traveller],evo_config,cell_config,exp_config,grn_config)
                
                adult_ages[traveller] = 1
            end
                        
            if rand() < evo_config.immigration_rate
                
                foreigner_id = [i for i in 1:evo_config.n_pools if i != Threads.threadid()][rand(1:evo_config.n_pools-1)]

                adult_pop_float[argmin(adult_fitnesses)] = shared_pool_float[foreigner_id]
                adult_pop_int[argmin(adult_fitnesses)] = shared_pool_int[foreigner_id]
                adult_fitnesses[argmin(adult_fitnesses)] = shared_pool_fitnesses[foreigner_id]
                adult_ages[argmin(adult_fitnesses)] = 1
            end

            breed = partialsortperm(adult_fitnesses, 1:topk, rev=true)

            for s in 1:evo_config.child_pop_size
                
                if rand() < evo_config.crossover

                    parents = sample(breed, 2, replace = false)  

                    child_pop_int[s], child_pop_float[s] = create_child(adult_pop_int[parents[1]],
                        adult_pop_float[parents[1]],
                        adult_pop_int[parents[2]],
                        adult_pop_float[parents[2]])
                else
                    parents = sample(breed, 1, replace = false)  

                    child_pop_int[s], child_pop_float[s] = create_child(adult_pop_int[parents[1]],
                        adult_pop_float[parents[1]],evo_config)
                end
            end

            for s in 1:evo_config.child_pop_size

                child_fitnesses[s] = evaluate_individual_fitness(child_pop_float[s],child_pop_int[s],evo_config,cell_config,
                                                                    exp_config,grn_config)
            end

            for s in 1:evo_config.child_pop_size

                if child_fitnesses[s] > adult_fitnesses[argmin(adult_fitnesses)]

                    adult_pop_float[argmin(adult_fitnesses)] = child_pop_float[s]
                    adult_pop_int[argmin(adult_fitnesses)] = child_pop_int[s]
                    adult_fitnesses[argmin(adult_fitnesses)] = child_fitnesses[s]
                    adult_ages[argmin(adult_fitnesses)] = 1

                end
            end
        end

        fitness_records[Threads.threadid()] = av_fitness
        best_float[Threads.threadid()]  = adult_pop_float[argmax(adult_fitnesses)]
        best_int[Threads.threadid()]  = adult_pop_int[argmax(adult_fitnesses)]
        best_fitnesses[Threads.threadid()]  = adult_fitnesses[argmax(adult_fitnesses)]
    end

    results = Dict{String, Any}("POOL_fitness_record" => fitness_records,
                                "POOL_genome_float" => best_float,
                                "POOL_genome_int" => best_int,
                                "POOL_genome_fitnesses"  => best_fitnesses)

    return results
end