support multiple instances per feature description

src/MAP_elites.jl

@@ -208,7 +208,7 @@ function population_from_archive(P,χ)
    return χ[findall(!iszero,P)]
 end
 
-function plot_MAP_archive(p)
+function plot_MAP_archive(p::Array{Float64, 2})
     X = []
     Y = []
     Z = []
@@ -223,7 +223,7 @@ function plot_MAP_archive(p)
     Z_ = Z[non_zero_inds]
     X_ = X[non_zero_inds]
     Y_ = Y[non_zero_inds]
-    fig = scatter(X_,Y_,log10.(Z_),legend=false)
+    fig = scatter(X_,Y_,log10.(Z_),legend=false,color=:green)
     xlabel!("Integer-order count")
     ylabel!("Fractional-order count")
     zlabel!("Log fitness")
@@ -247,7 +247,7 @@ function annotated_archive(p,x) #find a way to avoid the overlay plotting, maybe
         return a
 end
 
-function plot_MAP_archive_next(p)
+function plot_MAP_archive_next(p::Array{Float64, 2})
     X = []
     Y = []
     Z = []
@@ -262,15 +262,15 @@ function plot_MAP_archive_next(p)
     Z_ = Z[non_zero_inds]
     X_ = X[non_zero_inds]
     Y_ = Y[non_zero_inds]
-    fig = scatter(X_,Y_,log10.(Z_),legend=false)
+    fig = scatter(X_,Y_,log10.(Z_),legend=false,color=:green)
     xlabel!("Integer-order count")
     ylabel!("Fractional-order count")
     title!("Feature-performance plot")
     return fig
 end
 
 
-function plot_MAP_archive(p, cam = (40, 30))
+function plot_MAP_archive(p::Array{Float64, 2}, cam = (40, 30))
     X = []
     Y = []
     Z = []
@@ -284,7 +284,7 @@ function plot_MAP_archive(p, cam = (40, 30))
     Z_ = Z[findall(!iszero,Z)]
     X_ = X[findall(!iszero,Z)]
     Y_ = Y[findall(!iszero,Z)]
-    fig = scatter(X_,Y_,log10.(Z_),legend=false,camera = cam)  
+    fig = scatter(X_,Y_,log10.(Z_),legend=false,camera = cam,color=:green)  
     xlabel!("Integer-order count")
     ylabel!("Fractional-order count")
     zlabel!("log10 fitness")
@@ -440,7 +440,8 @@ function Archive_nyquist(measurements,frequencies,P,χ,cutoff=1e-3)
     archive_nyquist_plot = scatter(real(measurements),-imag(measurements),label="Measured",legend = :outertopright)
 # get population from archive 
     circuit_population = filter_fitness(population_from_archive(P,χ),cutoff)
-    circuit_population_sorted = sort_by_complexity(circuit_population)
+    circuit_population_unique = make_unique(circuit_population)
+    circuit_population_sorted = sort_by_complexity(circuit_population_unique)
 # plot nyquist for each circuit in archive
     for circ in circuit_population_sorted
         simulations = simulateimpedance_noiseless(circ,frequencies)
@@ -460,3 +461,144 @@ function filter_fitness(circuit_population,cutoff)
     filter(x->x.fitness<cutoff,circuit_population)
 end
 
+### Version with multiple entries (n) at each point in the feature space.
+
+function circuit_map_elites_R_m(measurements,frequencies;population_size= 30,iterations=10,feat_desc=fractional_descriptor,convergence_threshold = 5e-5,bounds = nothing,terminals="RCLP",head=8,initial_population=nothing, pe = 3)
+    parameter_bounds = Dict('R'=>[0,1.0e9],'C'=>[0,10],'L'=>[0,5],'P'=>[[0,0],[1.0e9,1]],'W'=>[0,1.0e9],'+'=>[0,0],'-'=>[0,0])
+    if typeof(bounds) == Dict{Char, Vector}
+        for key in keys(bounds)
+            parameter_bounds[key] = bounds[key]
+        end
+    end
+    if isnothing(initial_population)  
+        population = initializepopulation_R(parameter_bounds,population_size,head,terminals) 
+    elseif typeof(initial_population) == Vector{Circuit}
+        population = initial_population
+    elseif typeof(initial_population) in [Vector{String},Vector{Tuple{String, Vector{Float64}}}]
+        population = population_from_string(initial_population, head, population_size, terminals)
+    end
+    simplifypopulation!(population,terminals)
+    @suppress evaluate_population_fitness!(population,measurements,frequencies,parameter_bounds)
+    init_descriptions = [feat_desc(circ) for circ in population]
+    init_fitnesses = [circ.fitness for circ in population]
+    min_fit, iter = 1,0;
+    P = fill(0.0,(12,12,pe)); #Change dimensions
+    χ = Array{Circuit}(undef,(12,12,pe));
+    for (desc,fit,circ) in zip(init_descriptions,init_fitnesses,population) 
+        archive_update!(desc,fit,circ,P,χ)
+    end
+    while (iter ≤ iterations) && (min_fit ≥ convergence_threshold)
+        iter += 1
+        for e in 1:population_size
+            x′ = simplifycircuit(circuit_offspring_R(random_selection(χ,P),random_selection(χ,P),terminals)) #Random Selection also needs to be updated
+            b′ = feat_desc(x′) #als input.
+                @suppress evaluate_circuit_fitness!(x′,measurements,frequencies,parameter_bounds)
+                p′ = x′.fitness
+                archive_update!(b′,p′,x′,P,χ)
+        end
+        min_fit = minimum(P[P .> 0])
+    end
+    return P,χ,iter
+end
+
+function archive_update!(desc,fit,circ,P,χ)
+    # Check if archive contains at least one zero entry at feature descriptor location.
+    firstzero_index = findfirst(x -> x == 0, P[desc[1],desc[2],:])
+    if !isnothing(firstzero_index)
+        P[desc[1],desc[2],firstzero_index] = fit
+        χ[desc[1],desc[2],firstzero_index] = circ
+    elseif maximum(P[desc[1],desc[2],:]) > fit
+        P[desc[1],desc[2],findmax(P[desc[1],desc[2],:])[2]] = fit
+        χ[desc[1],desc[2],findmax(P[desc[1],desc[2],:])[2]] = circ
+    end
+end
+
+
+function circuit_map_elites_continue_R_m(iterations,P,χ,measurements,frequencies,terminals="RCLP",feat_desc=fractional_descriptor,bounds=nothing) #continue evolving starting from a given archive
+    parameter_bounds = Dict('R'=>[0,1.0e9],'C'=>[0,10],'L'=>[0,5],'P'=>[[0,0],[1.0e9,1]],'W'=>[0,1.0e9],'+'=>[0,0],'-'=>[0,0])
+    if typeof(bounds) == Dict{Char, Vector}
+        for key in keys(bounds)
+            parameter_bounds[key] = bounds[key]
+        end
+    end
+    iter = 0
+    population_size = length(population_from_archive(P,χ))
+    while (iter ≤ iterations) 
+        iter += 1
+        for e in 1:population_size
+            x′ = simplifycircuit(circuit_offspring_R(random_selection(χ,P),random_selection(χ,P),terminals))
+            b′ = feat_desc(x′) #als input.
+                @suppress evaluate_circuit_fitness!(x′,measurements,frequencies,parameter_bounds)
+                p′ = x′.fitness
+            archive_update!(b′,p′,x′,P,χ)
+        end
+    end
+    return P,χ
+end
+
+
+
+function plot_MAP_archive(p)
+
+    map_arch = plot_MAP_archive(p[:,:,1])
+    for i in 2:size(p)[3]
+        map_arch = plot_MAP_archive!(p[:,:,i])
+    end
+    return map_arch
+end
+
+function Archive_nyquist_m(measurements,frequencies,P,χ,cutoff=1e-3)
+    # plot ground truth
+        archive_nyquist_plot = scatter(real(measurements),-imag(measurements),label="Measured",legend = :outertopright)
+    # get population from archive 
+        circuit_population = filter_fitness(population_from_archive(P,χ),cutoff)
+        circuit_population_sorted = sort_by_complexity(circuit_population)
+    # plot nyquist for each circuit in archive
+        for circ in circuit_population_sorted
+            simulations = simulateimpedance_noiseless(circ,frequencies)
+            plot!(real(simulations),-imag(simulations),label=readablecircuit(circ))
+        end
+        xlabel!("Real{Z}/Ω")
+        ylabel!("-Imag{Z}/Ω")
+        return archive_nyquist_plot
+    end
+
+
+    function population_from_archive_m(P,χ)
+        return χ[findall(!iszero,P)]
+     end
+
+
+    function plot_MAP_archive!(p::Array{Float64, 2}, cam = (40, 30))
+        X = []
+        Y = []
+        Z = []
+        for x in 1:size(p)[1]
+            for y in 1:size(p)[2]
+                push!(X,x)
+                push!(Y,y)
+                push!(Z,p[x,y])
+            end
+        end
+        Z_ = Z[findall(!iszero,Z)]
+        X_ = X[findall(!iszero,Z)]
+        Y_ = Y[findall(!iszero,Z)]
+        fig = scatter!(X_,Y_,log10.(Z_),legend=false,camera = cam,color=:green)  
+        xlabel!("Integer-order count")
+        ylabel!("Fractional-order count")
+        zlabel!("log10 fitness")
+        title!("Feature-performance plot")
+        return fig
+    end
+
+    function make_unique(circuitpopulation)
+        newpop = []
+        newpop_strings = []
+        for circ in circuitpopulation
+            if !(readablecircuit(circ) in newpop_strings)
+                push!(newpop_strings,readablecircuit(circ))
+                push!(newpop,circ)
+            end
+        end
+        return newpop
+    end