Merge branch 'master' of https://github.com/MaximeVH/EquivalentCircui…

…ts.jl

Project.toml

@@ -7,8 +7,10 @@ version = "0.3.1"
 BlackBoxOptim = "a134a8b2-14d6-55f6-9291-3336d3ab0209"
 Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"
 DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"
+Distributed = "8ba89e20-285c-5b6f-9357-94700520ee1b"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 GeneralizedGenerated = "6b9d7cbe-bcb9-11e9-073f-15a7a543e2eb"
+Hwloc = "0e44f5e4-bd66-52a0-8798-143a42290a1d"
 Optim = "429524aa-4258-5aef-a3af-852621145aeb"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 

examples/example.jl

@@ -1,23 +1,30 @@
-using EquivalentCircuits, CSV
+using CSV
+using DataFrames
+using EquivalentCircuits
 
 """
-Circuitidentification is done with the circuitevolution function. This function can take several keyword arguments,
-allowing users to tune the circuitidentification if necessary. These can be found in the function documentation.
-The measurements and frequencies should be provided either as a CSV file, or as arrays.
-
+Circuit identification is done with the `circuit_evolution` function. This function can take
+several keyword arguments, allowing users to tune the circuit identification if necessary.
+These can be found in the function documentation. The measurements and frequencies should
+be provided either as a CSV file, or as arrays.
 """
-#Measurements input as CSV file with as columns: real impedances, imaginary impedances and frequencies.
-circuitevolution("example_measurements.csv")
+# Example 1: Input is CSV file with as columns: Re(Z), Im(Z), frequencies measurements
+circuit = circuit_evolution("example_measurements.csv")
 
-#Example of input as arrays.
-data = readdlm(raw"C:\Users\Dell\Documents\EquivalentCircuits.jl\example_measurements.csv",',')
-measurements = data[:,1] .+ data[:,2]im
-frequencies = data[:,3]
-circuitevolution(measurements,frequencies)
+# Example 2: Input is an array with as columns: Re(Z), Im(Z), frequencies measurements
+data = CSV.read("example_measurements.csv", DataFrame)
+Z = data[:, 1] .+ data[:, 2]im
+freq = data[:, 3]
+@time circuit = circuit_evolution(Z, freq)
 
 """
-When the equivalent circuit to be used is known, its parameters can be fit using the parameteroptimisation function.
-
+When the equivalent circuit to be used is known, its parameters can be fit using the
+`parameteroptimisation` function.
 """
+p = parameteroptimisation("[[C1,R2]-R3,P4]", Z, freq)
 
-parameteroptimisation("[[C1,R2]-R3,P4]", measurements, frequencies)
+"""
+If you want to generate multiple candidate circuits, you can use the `circuit_evolution_batch`
+function. This function is parallelised on all available cores.
+"""
+circuits = circuit_evolution_batch(Z, freq; generations=30, population_size=100, iters=6);

src/CircuitEvolution.jl

@@ -1,6 +1,8 @@
+using Hwloc
+
 """
-EquivalentCircuit(circuitstring::String, Parameters::NamedTuple)   
-   
+EquivalentCircuit(circuitstring::String, Parameters::NamedTuple)
+
 A structure used as the output of the evolutionary algorithm. It's fields are 'circuitstring', which is the string representation of a circuit (e.g. "R1-[C2,R3]-P4")
 and 'Parameters', a NamedTuple with the circuit's components and their corresponding parameter values.
 
@@ -30,13 +32,13 @@ function parametertuple(circuit,parameters)
 end
 
 # function initializecircuit(head=8,terminals="RCLP")
-#     karva = generatekarva(head,terminals) 
+#     karva = generatekarva(head,terminals)
 #     parameters = karva_parameters(karva)
 #     return Circuit(karva,parameters,nothing)
 # end
 
 function initializecircuit(bounds,head=8,terminals="RCLP")
-    karva = generatekarva(head,terminals) 
+    karva = generatekarva(head,terminals)
     parameters = karva_parameters_(karva,bounds)
     return Circuit(karva,parameters,nothing)
 end
@@ -80,8 +82,8 @@ function circuitfitness(circuit,measurements,frequencies,bounds)
     objective = objectivefunction_noweight(circfunc,measurements,frequencies)
     params = max.(min.(params,upper),lower.+10^-15)
     optparams,fitness = optimizeparameters(objective,params,lower,upper)
-    optparams = max.(min.(optparams,upper),lower.+10^-15) #workaround of optim.jl's bug 
-    return deflatten_parameters(optparams,tree,param_inds), fitness, param_inds 
+    optparams = max.(min.(optparams,upper),lower.+10^-15) #workaround of optim.jl's bug
+    return deflatten_parameters(optparams,tree,param_inds), fitness, param_inds
 end
 
 function circuitfitness(circuit,measurements,frequencies)
@@ -90,7 +92,7 @@ function circuitfitness(circuit,measurements,frequencies)
     objective = objectivefunction_noweight(circfunc,measurements,frequencies)
     optparams,fitness = optimizeparameters(objective,params,upper)
     optparams = max.(min.(optparams,upper),0) #workaround of optim.jl's bug
-    return deflatten_parameters(optparams,tree,param_inds), fitness, param_inds 
+    return deflatten_parameters(optparams,tree,param_inds), fitness, param_inds
 end
 
 function generate_offspring(population,terminals="RCLP")
@@ -121,19 +123,19 @@ end
 function evaluate_fitness!(population,measurements,frequencies,bounds)
     params = []
     param_inds = []
-    for circuit in population 
+    for circuit in population
         params,circuit.fitness,param_inds = @suppress circuitfitness(circuit,measurements,frequencies,bounds)
         circuit.parameters[param_inds] = params
-    end 
+    end
 end
 
 function evaluate_fitness!(population,measurements,frequencies)
     params = []
     param_inds = []
-    for circuit in population 
+    for circuit in population
         params,circuit.fitness,param_inds = circuitfitness(circuit,measurements,frequencies)
         circuit.parameters[param_inds] = params
-    end 
+    end
 end
 
 function population_from_string(circuitstring_list, head=8, size=30, terminals="RCLP")
@@ -156,7 +158,7 @@ end
 """
 circuit_evolution(measurements::Array{Complex{Float64},1},frequencies::Array{Float64,1}; <keyword arguments>)
 
-Identify an equivalent electrical circuit that fits a given set of electrochemical impedance spectroscopy measurements. 
+Identify an equivalent electrical circuit that fits a given set of electrochemical impedance spectroscopy measurements.
 
 The inputs are a circuit (e.g. "R1-[C2,R3]-P4") an array of complex-valued impedance measurements and their corresponding
 frequencies. The output is an EquivalentCircuit object containing a field circuitstring, which is a string denoting the identified circuit
@@ -169,7 +171,7 @@ and a field Parameters, which is a NamedTuple of the circuit's components with t
 - `head::Integer=8`: a hyperparameter than controls the maximum considered complexity of the circuits.
 - `cutoff::Float64=0.8`: a hyperparameter that controls the circuit complexity by removing redundant components.
 - `initial_population::Array{Circuit,1}=nothing`:the option to provide an initial population of circuits
-(obtained by using the loadpopulation function) with which the algorithm starts. 
+(obtained by using the loadpopulation function) with which the algorithm starts.
  Alternatively, users can provide a custom list of circuits which can either be a list of one or more circuit strings or a list of tuples
  where each tuple has the circuit string as first value and the parameters as second value.
  - `convergence_threshold::Float64=5e-4`: Convergence threshold for circuit identification. Increase this to reduce strictness, e.g., in case of noisy measurements.
@@ -181,7 +183,7 @@ julia> using EquivalentCircuits, Random
 
 julia> Random.seed!(25);
 
-julia> measurements = [5919.9 - 15.7, 5918.1 - 67.5im, 5887.1 - 285.7im, 5428.9 - 997.1im, 3871.8 - 978.9im, 
+julia> measurements = [5919.9 - 15.7, 5918.1 - 67.5im, 5887.1 - 285.7im, 5428.9 - 997.1im, 3871.8 - 978.9im,
 3442.9 - 315.4im, 3405.5 - 242.5im, 3249.6 - 742.0im, 1779.4 - 1698.9im,  208.2 - 777.6im, 65.9 - 392.5im];
 
 julia> frequencies = [0.10, 0.43, 1.83, 7.85, 33.60, 143.84, 615.85,  2636.65, 11288.38, 48329.30, 100000.00];
@@ -199,18 +201,18 @@ function circuit_evolution(measurements,frequencies;generations::Real=10,populat
     end
     # Either initialize a new population, or work with a provided initial population.
     # @assert typeof(initial_population) in [nothing,Vector{Circuit},Vector{String},Vector{Tuple{String, Vector{Float64}}}] "The initial population must be a list."
-    if isnothing(initial_population)  
+    if isnothing(initial_population)
         population = initializepopulation(parameter_bounds,population_size,head,terminals) #initializevariedpopulation(population_size,head)
     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
     # Theoretical simplification of the initial population.
-        simplifypopulation!(population,terminals) 
+        simplifypopulation!(population,terminals)
         evaluate_fitness!(population,measurements,frequencies,parameter_bounds)
         sort!(population)
-        generation = 0 
+        generation = 0
         # Keep track of the fittest individual circuit.
         elite = minimum(population)
         min_fitness = elite.fitness
@@ -236,16 +238,16 @@ function circuit_evolution(measurements,frequencies;generations::Real=10,populat
         @info "Algorithm did not converge"
     else
         best_circuit = readablecircuit(population[1])
-        return EquivalentCircuit(best_circuit,parameteroptimisation(best_circuit,measurements,frequencies)) #readablecircuit.(population[1:top_n]) 
+        return EquivalentCircuit(best_circuit,parameteroptimisation(best_circuit,measurements,frequencies)) #readablecircuit.(population[1:top_n])
     end
 end
 
 """
     circuit_evolution(filepath::String; <keyword arguments>)
 
- Identify an equivalent electrical circuit that fits a given set of electrochemical impedance spectroscopy measurements. 
+ Identify an equivalent electrical circuit that fits a given set of electrochemical impedance spectroscopy measurements.
 
- The inputs are a circuit (e.g. "R1-[C2,R3]-P4") and a filepath to a CSV file containing the three following columns: 
+ The inputs are a circuit (e.g. "R1-[C2,R3]-P4") and a filepath to a CSV file containing the three following columns:
  the real part of the impedance, the imaginary part of the impedance, and the frequencies corresponding to the measurements.
  The output is NamedTuple of the circuit's components with their corresponding parameter values.
 
@@ -257,7 +259,7 @@ end
 - `head::Integer=8`: a hyperparameter than controls the maximum considered complexity of the circuits.
 - `cutoff::Float64=0.8`: a hyperparameter that controls the circuit complexity by removing redundant components.
 - `initial_population::Array{Circuit,1}=nothing`:the option to provide an initial population of circuits
-(obtained by using the loadpopulation function) with which the algorithm starts. 
+(obtained by using the loadpopulation function) with which the algorithm starts.
  Alternatively, users can provide a custom list of circuits which can either be a list of one or more circuit strings or a list of tuples
  where each tuple has the circuit string as first value and the parameters as second value.
  - `convergence_threshold::Float64=5e-4`: Convergence threshold for circuit identification. Increase this to reduce strictness, e.g., in case of noisy measurements.
@@ -280,18 +282,18 @@ function circuit_evolution(filepath::String;generations::Real=10,population_size
     end
     # Either initialize a new population, or work with a provided initial population.
     # @assert typeof(initial_population) in [nothing,Vector{Circuit},Vector{String},Vector{Tuple{String, Vector{Float64}}}] "The initial population must be a list."
-    if isnothing(initial_population)  
+    if isnothing(initial_population)
         population = initializepopulation(parameter_bounds,population_size,head,terminals) #initializevariedpopulation(population_size,head)
     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
     # Theoretical simplification of the initial population.
-        simplifypopulation!(population,terminals) 
+        simplifypopulation!(population,terminals)
         evaluate_fitness!(population,measurements,frequencies,parameter_bounds)
         sort!(population)
-        generation = 0 
+        generation = 0
         # Keep track of the fittest individual circuit.
         elite = minimum(population)
         min_fitness = elite.fitness
@@ -317,6 +319,83 @@ function circuit_evolution(filepath::String;generations::Real=10,population_size
         @info "Algorithm did not converge"
     else
         best_circuit = readablecircuit(population[1])
-        return EquivalentCircuit(best_circuit,parameteroptimisation(best_circuit,measurements,frequencies)) #readablecircuit.(population[1:top_n]) 
+        return EquivalentCircuit(best_circuit,parameteroptimisation(best_circuit,measurements,frequencies)) #readablecircuit.(population[1:top_n])
     end
 end
+
+
+"""
+    circuit_evolution_batch(args...; kwargs...)
+
+Similar to the `circuit_evolution` function, but to generate multiple candidate circuits.
+
+# Arguments
+
+- `measurements::Array{Complex{Float64},1}`: Array of complex-valued impedance measurements.
+- `frequencies::Array{Float64,1}`: Frequencies corresponding to the measurements.
+
+# Keywords
+
+- `generations::Integer=10`: Maximum number of iterations of the evolutionary algorithm.
+- `population_size::Integer=30`: Number of individuals in the population during each iteration.
+- `terminals::String="RCLP"`: Allowed components to be included in the circuit.
+- `head::Integer=8`: Controls the maximum considered complexity of the circuits.
+- `cutoff::Float64=0.8`: Controls the circuit complexity by removing redundant components.
+- `initial_population::Array{Circuit,1}=nothing`: The option to provide an initial
+    population of circuits (obtained by using the loadpopulation function) with which the
+    algorithm starts. Alternatively, users can provide a custom list of circuits which can
+    either be a list of one or more circuit strings or a list of tuples where each tuple
+    has the circuit string as first value and the parameters as second value.
+- `convergence_threshold::Float64=5e-4`: Convergence threshold for circuit identification.
+    Increase to reduce strictness, e.g., in case of noisy measurements.
+- `bounds`::Dict{Char, Vector}: Upper/lower bounds for the circuit component parameter values.
+- `numprocs::Integer=num_physical_cores()`: Number of processes to use for parallelisation.
+- `iters`: Number of candidate circuits to generate.
+
+# Returns
+
+- `results::Vector{EquivalentCircuit}`: A vector of EquivalentCircuit objects.
+
+"""
+function circuit_evolution_batch(
+    measurements,
+    frequencies;
+    generations::Real=10,
+    population_size=30,
+    terminals="RCLP",
+    head=8,
+    cutoff=0.8,
+    initial_population=nothing,
+    convergence_threshold=5e-4,
+    bounds=nothing,
+    numprocs=num_physical_cores(),
+    iters,
+)
+    # Set up workers
+    _workers = addprocs(min(numprocs, iters))
+    import_module_on_workers(_workers)
+
+    # Run the circuit evolution in parallel on all workers
+    @info "Starting circuit evolution on $(length(_workers)) workers"
+    results = pmap(
+        _ -> circuit_evolution(
+            measurements,
+            frequencies;
+            generations=generations,
+            population_size=population_size,
+            terminals=terminals,
+            head=head,
+            cutoff=cutoff,
+            initial_population=initial_population,
+            convergence_threshold=convergence_threshold,
+            bounds=bounds,
+        ),
+        WorkerPool(_workers),
+        1:iters,
+    )
+
+    # Remove workers
+    rmprocs(_workers)
+
+    return results
+end

src/EquivalentCircuits.jl

@@ -1,6 +1,6 @@
 module EquivalentCircuits
 
-    export circuit_evolution, circuit_search
+    export circuit_evolution, circuit_search, circuit_evolution_batch
     export parameteroptimisation, circuitfunction
     export loadpopulation, get_parameter_upper_bound
     export simulateimpedance_noiseless
@@ -9,7 +9,7 @@ module EquivalentCircuits
     using BlackBoxOptim
     import Base: isless, length
 
-    include("Circuits.jl") 
+    include("Circuits.jl")
     include("CircuitFunction.jl")
     include("EvolutionOperators.jl")
     include("ObjectiveFunction.jl")

src/Miscellaneous.jl

@@ -1,3 +1,5 @@
+using Distributed
+
 # circuitstring conversion to and from impedance.py
 
 function impy_to_ec(circuitstring::String)
@@ -49,10 +51,10 @@ function ReadGamry(fn_DTA::AbstractString)
     else
         # --- original: Pt          Time	Freq	 Zreal	   Zimag	 Zsig	      Zmod	   Zphz	    Idc	    Vdc	    IERange
         s_header = ["", "Point",    "Time", "Frequ", "Z_real", "Z_imag", "Amplitude", "Z_mod", "Z_phz", "I_dc", "V_dc", "IE_range"]
-        filecontent = CSV.File(fn_data; 
+        filecontent = CSV.File(fn_data;
             header = s_header,
-            skipto = idx_header_line, 
-            decimal = char_dicimal_delim, 
+            skipto = idx_header_line,
+            decimal = char_dicimal_delim,
             delim = "\t");
     end
     # --- check column missmatch:
@@ -105,4 +107,12 @@ macro suppress(block)
             end
         end
     end
-end
+end
+
+
+# https://github.com/MilesCranmer/SymbolicRegression.jl/blob/master/src/Configure.jl
+function import_module_on_workers(procs)
+    @everywhere procs begin
+        Base.MainInclude.eval(:(using EquivalentCircuits))
+    end
+end