Merge pull request #29 from Kinetica-jl/adjustable-iterative

Optional undirected levels at start of iterative exploration

Project.toml

@@ -1,8 +1,7 @@
 name = "Kinetica"
 uuid = "3847ce11-53ec-444b-aa85-3b6606472139"
 authors = ["joegilkes <[email protected]>"]
-version = "0.5.10"
-
+version = "0.5.11"
 
 [deps]
 BSON = "fbb218c0-5317-5bc6-957e-2ee96dd4b1f0"

src/exploration/explore_utils.jl

@@ -295,11 +295,15 @@ end
 """
     identify_next_seeds(sol, sd::SpeciesData, seed_conc<:AbstractFloat[, elim_small_na::Int=0, 
                         ignore::Vector{String}=[], saveto::Union{String, Nothing}=nothing])
+    identify_next_seeds(sol, sd::SpeciesData[, elim_small_na::Int=0, ignore::Vector{String}=[],
+                        saveto::Union{String, Nothing}=nothing])
 
 Selects seed species for the next level of network exploration.
 
 Identifies species in `sol` with a maximum concentration above
 `seed_conc` and returns an array of the SMILES of these species.
+If `seed_conc` is not provided, assumes all species in `sd` should
+become seeds.
 
 Species meeting selection criteria can be manually ignored
 by including their SMILES in the `ignore` argument. Similarly,
@@ -347,5 +351,41 @@ function identify_next_seeds(sol, sd::SpeciesData, seed_conc::AbstractFloat;
         end
     end
 
+    return next_seeds
+end
+
+function identify_next_seeds(sol, sd::SpeciesData;
+                             elim_small_na::Int = 0,
+                             ignore::Vector{String} = [],
+                             saveto::Union{String, Nothing} = nothing)
+    next_seeds = String[]
+    next_seed_concs = Float64[]
+    umat = reduce(vcat, sol.u')
+    for species in axes(umat, 2)
+        if !isnothing(ignore) && sd.toStr[species] in ignore
+            continue
+        end
+        spectrace = umat[:, species]
+        spec_max_conc = maximum(spectrace)
+        spec_na = sd.xyz[species]["N_atoms"]
+        if elim_small_na > 0 && spec_na < elim_small_na
+            continue
+        else
+            push!(next_seeds, sd.toStr[species])
+            push!(next_seed_concs, spec_max_conc)
+        end
+    end
+
+    if !isnothing(saveto)
+        max_smi_length = maximum(length.(next_seeds))
+        open(saveto, "w") do f
+            write(f, "$(length(next_seeds))\n")
+            write(f, "SID   $(rpad("SMILES", max_smi_length))   Max. Conc.\n")
+            for (sid, (smi, conc)) in enumerate(zip(next_seeds, next_seed_concs))
+                write(f, "$(rpad(sid, 5)) $(rpad(smi, max_smi_length))   $conc\n")
+            end
+        end
+    end
+
     return next_seeds
 end

src/exploration/methods.jl

@@ -28,8 +28,8 @@ end
 """
     IterativeExplore(rdir_head::String, reac_smiles::VEctor{String}, cde::CDE[, maxiters::Int=1000, 
                      rxn_convergence_threshold::Int=5, seed_convergence_threshold::Int=3, seed_conc=0.05,
-                     independent_blacklist::Vector{String}=[], inert_species::Vector{String}=[],
-                     modify_network_on_solve::Bool=true])
+                     n_undirected_levels::Int=0, independent_blacklist::Vector{String}=[], 
+                     inert_species::Vector{String}=[], modify_network_on_solve::Bool=true])
 
 Keyword-based container for parameters used in iterative kinetics-based CRN exploration.
 
@@ -41,6 +41,7 @@ Contains fields for:
 * Number of subspace iterations with no change in reactions to consider a subspace converged (`rxn_convergence_threshold`)
 * Number of level iterations with no change in seeds to consider the network converged (`seed_convergence_threshold`)
 * Concentration above which species will be selected as seeds each level (`seed_conc`)
+* Number of undirected levels of exploration to perform at the start (`n_undirected_levels`)
 * Blacklist of species to avoid doing independent subspace explorations on (`independent_blacklist`)
 * Inert species that should not be considered for reaction (`inert_species`)
 * Whether to allow CRN modification after solving (`modify_network_on_solve`)
@@ -53,6 +54,7 @@ Contains fields for:
     rxn_convergence_threshold::Integer = 5
     seed_convergence_threshold::Integer = 3
     seed_conc::uType = 0.05
+    n_undirected_levels::Integer = 0
     independent_blacklist::Vector{String} = String[]
     inert_species::Vector{String} = String[]
     modify_network_on_solve::Bool = true
@@ -227,9 +229,15 @@ function explore_network(exploremethod::IterativeExplore,
 
         for seed in current_seeds push!(explored_seeds, seed) end
         seeds_out_saveto = isnothing(savedir) ? nothing : joinpath(savedir, "seeds_level$(loc.level).out")
-        next_seeds = identify_next_seeds(res.sol, sd, exploremethod.seed_conc;
-                                        ignore=exploremethod.inert_species,
-                                        saveto=seeds_out_saveto)
+        if loc.level <= exploremethod.n_undirected_levels
+            next_seeds = identify_next_seeds(res.sol, sd; 
+                                             ignore=exploremethod.inert_species,
+                                             saveto=seeds_out_saveto)
+        else
+            next_seeds = identify_next_seeds(res.sol, sd, exploremethod.seed_conc;
+                                             ignore=exploremethod.inert_species,
+                                             saveto=seeds_out_saveto)
+        end
 
         if Set(current_seeds) == Set(next_seeds)
             convergence_count += 1