Fix bugs for creating networks with more than 3 chemical species

Project.toml

@@ -1,7 +1,7 @@
 name = "ReactionNetworkEvolution"
 uuid = "61220d1e-4cdb-4fc3-9e6c-02bd7f830b85"
 authors = ["really-lilly <[email protected]>"]
-version = "0.1.4"
+version = "0.1.5"
 
 [deps]
 ArgParse = "c7e460c6-2fb9-53a9-8c5b-16f535851c63"

README.md

@@ -204,7 +204,7 @@ Any settings that are not specified in the JSON file will be set to the default
     <td>If generating networks to fit data, path to .csv file with data to fit - if no path is given, ReactionNetworkEvolution.jl will attempt to generate oscillators</td>
   </tr>
   <tr>
-    <td>track_metadata</td>
+    <td>track_fitness</td>
     <td>true</td>
     <td>Track information about evolution such as top fitness, number of species, etc for each generation</td>
   </tr>
@@ -216,7 +216,7 @@ Any settings that are not specified in the JSON file will be set to the default
   <tr>
     <td>process_output_oscillators</td>
     <td>true</td>
-    <td>If true, utomatically test if evolved networks are oscillators or not</td>
+    <td>If true, automatically test if evolved networks are oscillators or not</td>
   </tr>
 
 </table>

run_evolution.jl

@@ -3,6 +3,23 @@ Pkg.activate(".")
 using ReactionNetworkEvolution
 using ArgParse
 
+
+
+function parse_initial_concentrations(initial_concentrations::String)
+    # Strip whitespace
+    initial_concentrations = filter(x -> !isspace(x), initial_concentrations)
+    # Strip brackets
+    initial_concentrations = initial_concentrations[2:end-1]
+    # Separate by commas 
+    initial_concentrations = split(initial_concentrations, ",")
+    new_array = Vector{Float64}()
+    for num in initial_concentrations
+        push!(new_array, parse(Float64, num))
+    end
+    return new_array
+end
+
+
 s = ArgParseSettings()
 
 @add_arg_table s begin
@@ -12,6 +29,12 @@ s = ArgParseSettings()
     "--ntrials"
         arg_type = Int
         default = -1
+    "--nspecies"
+        arg_type = Int
+        default = 3
+    "--initial_concentrations"
+        arg_type = String
+        default = "[1., 5., 9.]"
     "--population_size"
         arg_type = Int
         default = -1
@@ -31,8 +54,12 @@ end
 
 parsed_args = parse_args(ARGS, s)
 
-@time ReactionNetworkEvolution.run_evolution(ngenerations=parsed_args["ngenerations"],
+parsed_args["initial_concentrations"] = parse_initial_concentrations(parsed_args["initial_concentrations"])
+
+ReactionNetworkEvolution.run_evolution(ngenerations=parsed_args["ngenerations"],
               ntrials=parsed_args["ntrials"],
+              nspecies=parsed_args["nspecies"],
+              initial_concentrations = parsed_args["initial_concentrations"],
               population_size=parsed_args["population_size"],
               pathtosettings=parsed_args["pathtosettings"],
               outputpath=parsed_args["outputpath"],

settings.json

@@ -5,7 +5,8 @@
 "percent_rateconstant_change": 20,			
 "p_new_rateconstant": 0.15,				
 "population_size": 100,						
-"ngenerations": 800,					
+"ngenerations": 800,
+"nspecies": 3,					
 "nreactions": 5,					
 "max_offspring_portion": 0.1,				
 "writeout_threshold": 0.05,				
@@ -25,9 +26,10 @@
 "chemical_species_names": ["S0", "S1", "S2"],		
 "initial_concentrations": [1.0, 5.0, 9.0],			
 "enable_speciation": true,				
-"track_metadata": true,					
+"track_fitness": true,					
 "average_fitness": false,				
 "same_fitness_crossover": false,			
 "same_fitness_percent_range": 5,		
-"process_output_oscillators": true,			
+"process_output_oscillators": false,
+"objective_data_path": "DEFAULT",			
 "note": ""}						

src/ReactionNetworkEvolution.jl

@@ -7,25 +7,11 @@ import Dates: now, format
 include("evo_utils.jl")
 include("process_output.jl")
 
-function evolve_networks(batchnum::Int64, parentdir::String, settings::Settings, objfunct::ObjectiveFunction)
+function evolve_networks(parentdir::String, settings::Settings, objfunct::ObjectiveFunction)
 
-    starttime = format(now(), "YYYYmmdd_HHMMSS")
-    starttime = "$starttime" * randstring(3)  # for race conditions
-    starttime = joinpath(parentdir, starttime)
-    mkdir(starttime)
-    if settings.track_metadata
+    if settings.track_fitness
         tracker = Dict{String, Any}(
-            "batch_num" => batchnum,
             "top_individual_fitness" => Vector{Float64}(),
-            "num_unique_networks" => Vector{Int64}(),
-            "num_individuals" => Vector{Int64}(),
-            "num_best_network" => Vector{Int64}(),
-            "total_num_species" => Vector{Int64}(),
-            "avg_species_distances" => Vector{Float64}(),
-            "min_species_distances" => Vector{Float64}(),
-            "max_species_distances" => Vector{Float64}(),
-            "delta" => Vector{Float64}(),
-            "avg_intraspecies_distance" => Vector{Float64}()
         )
     end
 
@@ -40,32 +26,16 @@ function evolve_networks(batchnum::Int64, parentdir::String, settings::Settings,
         species_by_IDs, DELTA = speciate(species_by_IDs, population, DELTA, settings)
     end
 
-    if settings.track_metadata
-        num_unique = length(Set(population))
-        push!(tracker["num_unique_networks"], num_unique)
-        push!(tracker["num_individuals"], length(population))
-        push!(tracker["total_num_species"], length(keys(species_by_IDs)))
-        avg, mn, mx = get_diversity_stats(species_by_IDs, settings)
-        push!(tracker["avg_species_distances"], avg)
-        push!(tracker["min_species_distances"], mn)
-        push!(tracker["max_species_distances"], mx)
-        push!(tracker["delta"], DELTA)
-        push!(tracker["avg_intraspecies_distance"], get_intraspecies_distances(species_by_IDs, settings))
-    end
 
-
     for i in 1:settings.ngenerations
 
-        # TODO: save time by also returning the top fitness and network?
         species_by_IDs, total_fitness = evaluate_population_fitness(objfunct, species_by_IDs, settings)
-        species_by_IDs, total_offspring = calculate_num_offspring(species_by_IDs, total_fitness, settings, writeoutdir=joinpath(starttime, "stalled_models"))
+        species_by_IDs, total_offspring = calculate_num_offspring(species_by_IDs, total_fitness, settings, writeoutdir=joinpath(parentdir, "stalled_models"))
 
         bestnetwork, maxfitness = gettopmodel(species_by_IDs)
 
-        if settings.track_metadata
-            push!(tracker["top_individual_fitness"],maxfitness)
-            bestnetwork_count = count_best_networks(bestnetwork, population)
-            push!(tracker["num_best_network"], bestnetwork_count)
+        if settings.track_fitness
+            push!(tracker["top_individual_fitness"], maxfitness)
         end
 
         if maxfitness > settings.writeout_threshold 
@@ -74,18 +44,6 @@ function evolve_networks(batchnum::Int64, parentdir::String, settings::Settings,
 
         population = reproduce_networks(species_by_IDs, settings, ng, objfunct, total_offspring)
 
-        if settings.track_metadata
-            num_unique = length(Set(population))
-            push!(tracker["num_unique_networks"], num_unique)
-            push!(tracker["num_individuals"], length(population))
-            push!(tracker["total_num_species"], length(keys(species_by_IDs)))
-            avg, mn, mx = get_diversity_stats(species_by_IDs, settings)
-            push!(tracker["avg_species_distances"], avg)
-            push!(tracker["min_species_distances"], mn)
-            push!(tracker["max_species_distances"], mx)
-            push!(tracker["avg_intraspecies_distance"], get_intraspecies_distances(species_by_IDs, settings))
-        end
-
         if settings.enable_speciation
             species_by_IDs, DELTA = speciate(species_by_IDs, population, DELTA, settings)
         else
@@ -96,37 +54,26 @@ function evolve_networks(batchnum::Int64, parentdir::String, settings::Settings,
     species_by_IDs, total_fitness = evaluate_population_fitness(objfunct, species_by_IDs, settings)
     bestnetwork, maxfitness = gettopmodel(species_by_IDs)
 
-    if settings.track_metadata
+    if settings.track_fitness
         push!(tracker["top_individual_fitness"],maxfitness)
-        bestnetwork_count = count_best_networks(bestnetwork, population)
-        push!(tracker["num_best_network"], bestnetwork_count)
-        push!(tracker["delta"], DELTA)
-        push!(tracker["avg_intraspecies_distance"], get_intraspecies_distances(species_by_IDs, settings))
     end
 
-    writeoutnetwork(bestnetwork, "bestmodel_$(bestnetwork.ID)", directory=joinpath(starttime,"final_models"))
-
-    for ID in keys(species_by_IDs)
-        species = species_by_IDs[ID]
-        topnetwork = species.topnetwork
-        if topnetwork.ID != bestnetwork.ID
-            writeoutnetwork(topnetwork, "$(species.ID)", directory=joinpath(starttime, "final_models"))
-        end
-    end
+    writeoutnetwork(bestnetwork, "$(bestnetwork.ID).ant", parentdir)
 
-    if settings.track_metadata
+    if settings.track_fitness
         stringtracker = json(tracker)
-        open(joinpath(starttime, "datatracker.json"), "w") do f
+        open(joinpath(parentdir, "$(bestnetwork.ID)_fitness.json"), "w") do f
             write(f, stringtracker)
         end
     end
 
-    writeout_settings(settings, joinpath(starttime, "settings.json"))
 end
 
 function run_evolution(;
     ngenerations::Int64=-1,
     ntrials::Int64=-1,
+    nspecies::Int64=3,
+    initial_concentrations::Vector{Float64}=[1., 5., 9.],
     population_size::Int64=-1,
     pathtosettings::String="",
     outputpath::String="",
@@ -164,27 +111,30 @@ function run_evolution(;
         pathtosettings = "DEFAULT"
     end
 
-    settings, objectivefunction = read_usersettings(pathtosettings, ngenerations=ngenerations, population_size=population_size, seed=seed, note=note)
+    settings, objectivefunction = read_usersettings(pathtosettings, ngenerations=ngenerations, nspecies=nspecies, initial_concentrations=initial_concentrations, population_size=population_size, seed=seed, note=note)
 
     println("Beginning $ntrials trials with $(settings.ngenerations) generations each")
 
-    for i in 1:ntrials
+    @time for i in 1:ntrials
         if settings.verbose && (i%10 == 0 || i == ntrials)
             println(". trial $i")
         elseif settings.verbose
             print(".")
         end
-        evolve_networks(i, path, settings, objectivefunction)
+        evolve_networks(path, settings, objectivefunction)
     end
 
+    writeout_settings(settings, joinpath(outputpath, "settings.json"))
+
     if settings.process_output_oscillators
         println("Sorting oscillators...")
-        process_oscillators(outputpath)
+        process_oscillators(path, track_fitness=settings.track_fitness)
     else
         println("Output written to $path")
     end
 
     println("Done")
 end
 
-end #module
+end #module
+

src/evo_utils.jl

@@ -750,7 +750,7 @@ function general_crossover(network1::ReactionNetwork, network2::ReactionNetwork)
 
 end
 
-function writeoutnetwork(network::ReactionNetwork, filename::String; directory::String="stalled_models")
+function writeoutnetwork(network::ReactionNetwork, filename::String, directory::String)
     astr = convert_to_antimony(network)
     astr *= "\n#fitness: $(network.fitness)"