Performance calculations carry over from previous runs (#32)

* Performance evaluations carry over from previous. Warn user when sim is throttled. Small test added

* v0.1.1, update changelog

* Improve coverage

Author: xoth42

CHANGELOG.md

@@ -1,5 +1,9 @@
 # News
 
+## v0.1.1 - 2025-07-29
+- Performance evaluations carry over from previous runs. 
+- Warning when simulation is throttled (fidelity = 1)
+
 ## v0.1.0 - 2025-06-13
 
 - Simulation capabilities for commond gates.

Project.toml

@@ -1,7 +1,7 @@
 name = "QEPOptimize"
 uuid = "9f67e06f-6394-43be-a72d-a2000081c01d"
 authors = ["Stefan Krastanov <[email protected]>", "QEPOptimize contributors"]
-version = "0.1.0"
+version = "0.1.1"
 
 [deps]
 BPGates = "2c1a90cd-bec4-4415-ae35-245016909a8f"

example/pluto.jl

@@ -1,5 +1,5 @@
 ### A Pluto.jl notebook ###
-# v0.20.10
+# v0.20.13
 
 using Markdown
 using InteractiveUtils
@@ -64,6 +64,8 @@ md"""
 
 * Number of Simulations: $(@bind num_simulations PlutoUI.Slider(100:100:5000, default=1000, show_value=true))
 
+* Max performance calculations per circuit: $(@bind max_perf_calcs PlutoUI.Slider(1:1:50, default=10, show_value=true))
+
 * Population Size: $(@bind pop_size PlutoUI.Slider(10:10:100, default=20, show_value=true))
 
 * Initial Operations: $(@bind start_ops PlutoUI.Slider(5:20, default=10, show_value=true))
@@ -73,7 +75,7 @@ md"""
 
 * Number of Evolution Steps: $(@bind evolution_steps PlutoUI.Slider(50:150, default=80, show_value=true))
 
-* New Mutants: $(@bind new_mutants PlutoUI.Slider(5:5:30, default=10, show_value=true))
+* New Mutants: $(@bind new_mutants PlutoUI.Slider(5:1:30, default=10, show_value=true))
 
 * Drop Probability: $(@bind p_drop PlutoUI.Slider(0.0:0.05:0.5, default=0.1, show_value=true))
 
@@ -92,6 +94,7 @@ begin
 	    code_distance=1, # Not needed to change for now
 	    pop_size=pop_size,
 	    noises=[NetworkFidelity(network_fidelity), PauliNoise(paulix, pauliy, pauliz)],
+		max_performance_calcs=max_perf_calcs
 	)
 
 	init_config = (;

src/QEPOptimize.jl

@@ -1,5 +1,7 @@
 module QEPOptimize
 
+import Base:+,==
+
 using BPGates
 using BPGates: mctrajectory!, continue_stat, PauliNoise # TODO these should be exported by default
 

src/evolve.jl

@@ -19,6 +19,7 @@ function step!(
     p_drop=0.1,
     p_mutate=0.1,
     p_gain=0.1,
+    max_performance_calcs=10
 )
     # Mark existing individuals as survivors
     # Survivors ensure that some individuals are carried over unchanged, maintaining good solutions
@@ -36,7 +37,8 @@ function step!(
         num_simulations,
         purified_pairs,
         number_registers, # TODO (low priority) this should be by-default derived from `indiv`
-        noises=[NetworkFidelity(0.9)] # TODO configurable noise
+        noises=[NetworkFidelity(0.9)], # TODO configurable noise
+        max_performance_calcs
     )
     cull!(population, pop_size)
 end
@@ -50,8 +52,19 @@ function multiple_steps_with_history!(
     fitness_history[1, :] = [i.fitness for i in population.individuals]
     transition_counts = []
 
+    throttling_warned = 0
+
     for i in 1:steps
         step!(population; step_config...)
+
+        # Check to make sure that the optimizer is not in the 'fitness = 1.0' failure mode
+        if population.individuals[1].fitness == 1.0 && throttling_warned < THROTTLE_WARNINGS
+            throttling_warned += 1
+            @warn "Simulation is throttled: Increase simulation count, or decrease new mutants to fix. Top circuit has fitness = 1.0"
+            # If fitness is 1, then all of the simulations done to evaluate a circuit show no errors. This implies the circuit is good, but stops the optimizer from performing well. Increasing simulation count can fix this issue
+            # but, decreasing new mutants will also fix the issue. This is because fewer new circuits need to be simulated, causing more simulations on the current circuits to get a better non-1.0 fidelity. 
+        end
+
         fitness_history[i+1,:] = [i.fitness for i in population.individuals]
         push!(transition_counts, counter([i.history for i in population.individuals]))
         step_callback()
@@ -138,21 +151,26 @@ function simulate_and_sort!(
     purified_pairs::Int=1,
     number_registers::Int=2, # TODO (low priority) this should be by-default derived from `indiv`
     code_distance::Int=1,
-    noises=[NetworkFidelity(0.9)]
+    noises=[NetworkFidelity(0.9)],
+    max_performance_calcs::Int=10
 )
     # calculate and update each individual's performance
     function update!(indiv)
-        calculate_performance!(indiv;
-            num_simulations,
-            purified_pairs,
-            number_registers,
-            code_distance,
-            noises)
+        # Restrict performance calculation if this indiv has already reached the max calcs. 
+        # However, if the calculated fidelity is undetermined (1.0), then it needs more calculations to try and get a non-one value (ie: 0.9999). 
+        # Otherwise, if circuits have a fidelity of 1.0, it is difficult to distinguish between ones that are better (f:0.9999) or worse (f:0.9997).  
+        if indiv.performance.num_calcs < max_performance_calcs || indiv.fitness == 1.0
+            calculate_performance!(indiv;
+                num_simulations,
+                purified_pairs,
+                number_registers,
+                code_distance,
+                noises)
+        end
         indiv.fitness = indiv.performance.purified_pairs_fidelity # TODO make it possible to select the type of fitness to evaluate
     end
 
     # Parallel processing for performance calculations. Max threads will be set by the threads specified when running julia. ex) julia -t 16
-    max_threads::Int = Threads.nthreads()
 
     tmap(update!,population.individuals)
 
@@ -177,7 +195,8 @@ function initialize_pop!(
     num_simulations::Int=100,
     purified_pairs::Int=1,
     code_distance::Int=1,
-    noises=[NetworkFidelity(0.9)]
+    noises=[NetworkFidelity(0.9)],
+    max_performance_calcs=10
 )
     valid_pairs=1:number_registers # TODO (low priority) decouple valid_pairs from number_registers
 
@@ -197,7 +216,8 @@ function initialize_pop!(
         purified_pairs,
         number_registers,
         code_distance,
-        noises
+        noises,
+        max_performance_calcs
     )
     cull!(population,pop_size)
 end

src/performance_eval.jl

@@ -14,6 +14,13 @@ function calculate_performance!(
 
     count_success = 0
     counts_marginals = zeros(Int,purified_pairs) # an array to find F₁, F₂, …, Fₖ (tracks how often each purified bell pair is in the desired state)
+    
+    # Edge case: purified pairs/registers have changed, and the circuit has previously had a performance calculation. In this case, the previous performance must be discarded, as its error probabilities are now invalid.
+    if indiv.performance.num_calcs > 0 && length(indiv.performance.error_probabilities) != purified_pairs+1 
+        # signal the previous performance to be discarded 
+        indiv.performance = Performance() # implies it is a blank performance
+    end
+
     counts_nb_errors = zeros(Int,purified_pairs+1) # an array to find P₀, P₁, …, Pₖ (tracks how often a given total number of errors happens) Careful with indexing it as it includes a P₀!
 
     noisy_purification_circuit = indiv.ops
@@ -39,8 +46,8 @@ function calculate_performance!(
     end
 
     if count_success == 0
-        indiv.performance =  Performance(counts_nb_errors, 0,0, 0, 0) # TODO this is probably going to break the optimization runs and lead to picking low performing individuals for certain cost functions -- do this better
-        return indiv.performance
+        # TODO this is probably going to break the optimization runs and lead to picking low performing individuals for certain cost functions -- do this better
+        return update_performance!(indiv,Performance(counts_nb_errors,0,0, 0, 0,1))
     end
 
     p_success = count_success    / num_simulations # proportion of successful simulations
@@ -50,7 +57,30 @@ function calculate_performance!(
     correctable_errors = div(code_distance  - 1, 2) # maximum number of correctable errors based on code distance after teleportation
     indiv_logical_qubit_fidelity = sum(err_probs[1:min(end, correctable_errors+1)]) # calculates the logical qubit fidelity by summing the probabilities of correctable errors
 
-    indiv.performance =  Performance(err_probs, err_probs[1], indiv_logical_qubit_fidelity, mean(marginals), p_success)
+    return update_performance!(indiv,Performance(
+        err_probs, 
+        err_probs[1], 
+        indiv_logical_qubit_fidelity, 
+        mean(marginals), 
+        p_success,
+        1)
+    )
+end
+
+"""
+    update_performance!(indiv::Individual,new::Performance)
+
+Helper function to deal with circuit init, and performance averaging. The indiv given will have their performance averaged with the new performance, if their current performance is defined (nonzero). Otherwise, the new performance will overrite their current. 
+"""
+function update_performance!(indiv::Individual,new::Performance)
+    @assert new.num_calcs >= 1 # This should always be the case at this point. The new performance should not be empty.
 
+    if indiv.performance.num_calcs == 0
+        # old perf has not been calculated yet/not usable, so use the new one
+        indiv.performance = new
+    else 
+        # old perf is calculated, so average the two
+        indiv.performance += new
+    end
     return indiv.performance
 end

src/types.jl

@@ -1,6 +1,8 @@
 # TODO (low priority) this would be a great place to use an Enum or, even better, an algebraic data type (ADT)
 const HISTORIES = [:manual, :survivor, :random, :child, :drop, :gain, :swap, :opsmutate]
 
+const THROTTLE_WARNINGS = 10 # max amount of warnings per multiple_steps_with_history! call
+
 "A convenient structure to store various purification performance metrics."
 struct Performance
     "a list of probabilities as such `[probability for no errors, probability for one Bell pair to be erroneous, probability for two Bell pairs to be erroneous, ..., probability for all Bell pairs to be erroneous]`"
@@ -15,12 +17,62 @@ struct Performance
     average_marginal_fidelity::Float64
     "the proportion of runs of a given protocol that do not have detected errors (i.e. Alice and Bob do not measure any errors)"
     success_probability::Float64
+    "the number of times that performance has been calculated and averaged for this circuit"
+    num_calcs::Int64
 end
 
-Performance() = Performance(Float64[], 0.0, 0.0, 0.0, 0.0)
+Performance() = Performance(Float64[], 0.0, 0.0, 0.0, 0.0,0)
+
+Base.copy(p::Performance) = Performance(copy(p.error_probabilities), p.purified_pairs_fidelity, p.logical_qubit_fidelity, p.average_marginal_fidelity, p.success_probability,p.num_calcs)
+
+"""
+    +(a::Performance,b::Performance)
+
+Add two performance points, and average them.
+"""
+function +(a::Performance,b::Performance) 
+    # start with using a's
+    new_error_probabilities = a.error_probabilities
+
+    # if error probs are empty on a, use b. if b is also empty, dosen't matter, both are empty, use a.
+    if isempty(a.error_probabilities)
+        new_error_probabilities =  b.error_probabilities
+    elseif !isempty(b.error_probabilities)
+        # make sure error probs are same shape
+        if length(a.error_probabilities) != length(b.error_probabilities)
+            @warn "Error probabilities mismatch, defaulting to new performance"
+            # differing error_prob lengths should have been dealt with. This warning is only reached if there is a bug in the circuit performance calculator
+            return b
+        else 
+            # both error probs are nonempty, and same length. Sum and divide them 
+            for i in 1:length(a.error_probabilities)
+                new_error_probabilities[i] = (a.error_probabilities[i] + b.error_probabilities[i]) / 2
+            end
+        end
+    end
+
+    # Now average the rest of the attributes (floats)
+    new_purified_pairs_fidelity = (a.purified_pairs_fidelity + b.purified_pairs_fidelity) / 2
+    new_logical_qubit_fidelity = (a.logical_qubit_fidelity + b.logical_qubit_fidelity) / 2
+    new_average_marginal_fidelity = (a.average_marginal_fidelity + b.average_marginal_fidelity) / 2
+    new_success_probability = (a.success_probability + b.success_probability) / 2
+
+    # and mark that new calcs were added
+    new_num_calcs = a.num_calcs + b.num_calcs
+    return Performance(new_error_probabilities,new_purified_pairs_fidelity ,new_logical_qubit_fidelity,new_average_marginal_fidelity,new_success_probability,new_num_calcs)
+end
 
-Base.copy(p::Performance) = Performance(copy(p.error_probabilities), p.purified_pairs_fidelity, p.logical_qubit_fidelity, p.average_marginal_fidelity, p.success_probability)
+"""
+    ==(a::Performance, b::Performance)
 
+Compare equality of all internal values of a performance stuct. Used for testing.
+"""
+==(a::Performance, b::Performance) = (a.error_probabilities == b.error_probabilities &&
+                                    a.purified_pairs_fidelity == b.purified_pairs_fidelity &&
+                                    a.logical_qubit_fidelity == b.logical_qubit_fidelity &&
+                                    a.average_marginal_fidelity == b.average_marginal_fidelity &&
+                                    a.success_probability == b.success_probability && 
+                                    a.num_calcs == b.num_calcs)
 
 "An individual (circuit) in the population we are evolving"
 mutable struct Individual

test/test_types.jl

@@ -0,0 +1,73 @@
+using TestItems
+
+@testitem "performance averaging" begin
+    using QEPOptimize
+    using QEPOptimize:Performance
+
+    low = Performance([0,0,1],0,0,0,0,1)
+    high = Performance([1,0,0],1,1,1,1,1)
+    mid = low + high
+    @test mid == Performance([0.5,0,0.5],0.5,0.5,0.5,0.5,2)
+
+    # works with +=
+    low = Performance([0,0,1],0,0,0,0,1)
+    high = Performance([1,0,0],1,1,1,1,1)
+
+    low += high
+    @test low == mid
+end
+
+@testitem "Performance averages deal with low sim count/circuit failures" begin 
+    using QEPOptimize
+    using QEPOptimize: initialize_pop!, step!, NetworkFidelity 
+    using BPGates: PauliNoise 
+
+    config = (;
+        num_simulations=5, # needs to be large enough to resolve circuit noise but you can make smaller too # TODO anneal this and save old results
+        number_registers=4, # do 3 for something faster
+        purified_pairs=1,
+        code_distance=1,
+        pop_size = 20,
+        noises=[NetworkFidelity(0.5), PauliNoise(0.01/3, 0.01/3, 0.01/3)],
+    )
+
+    init_config = (;
+        start_ops = 10,
+        start_pop_size = 1000,
+        config...
+    )
+
+    step_config = (;
+        max_ops = 15, # do 5 for something faster
+        new_mutants = 10,
+        p_drop = 0.1,
+        p_mutate = 0.1,
+        p_gain = 0.1,
+        config...
+    )
+
+    pop = Population()
+
+    initialize_pop!(pop; init_config...)
+
+    _, fitness_history, transition_counts_matrix, transition_counts_keys = multiple_steps_with_history!(pop, 80; step_config...)
+
+end
+
+@testitem "performance averages deal with bad error_probs" begin
+    using QEPOptimize
+    using QEPOptimize:Performance
+
+    # defaults to the second performance if error probs are different
+    normal = Performance([0,0,1],0,0,0,0,1)
+    new_err_probs = Performance([1,0,0,0],1,1,1,1,1)
+    result = normal + new_err_probs
+    @test result == Performance([1,0,0,0],1,1,1,1,1)
+
+    # defaults to the second err_probs if no error probs, still averages
+    normal = Performance([],0,0,0,0,1)
+    new_err_probs = Performance([1,0,0,0],1,1,1,1,1)
+    result = normal + new_err_probs
+    @test result == Performance([1,0,0,0],0.5,0.5,.5,.5,2)
+    
+end