fixing commits

README.md

@@ -1,2 +1,46 @@
 # SynthEvo
-A tool for finding CRN topologies under user-defined constraint 
+
+A tool for finding CRN topologies under user-defined constraint
+
+### Installation
+
+To install the package, run the following command in `julia`:
+
+```julia
+using Pkg; Pkg.add("SynthEvo")
+```
+
+finally import the package in your code with:
+
+```julia
+using SynthEvo
+```
+
+### Examples
+
+check the `examples` directory for some to see some commented use-cases of the package.
+
+#### Use from source
+
+To use this package directly from source, clone the repository and in `julia` run the following commands:
+
+```julia
+using Pkg; activate("/path/to/SynthEvo")
+```
+
+Where `/path/to/SynthEvo` is the path to the SynthEvo directory inside the cloned repository.
+
+In this way, you can use the package for the current session. 
+
+### Changelog
+
+v0.0.0: Initial release with basic functionality
+
+v0.1.0: refactoring and optimization
+
+* Updataed code to imporve the performance of the symbolic computation
+* Added notebooks with examples
+* Added function documentation 
+
+
+

SynthEvo/Project.toml

@@ -0,0 +1,21 @@
+name = "SynthEvo"
+uuid = "6cd3eeb9-5a42-42a0-b734-1cbe889bcc08"
+authors = ["[email protected] <[email protected]>"]
+version = "0.9.9"
+
+[deps]
+Catalyst = "479239e8-5488-4da2-87a7-35f2df7eef83"
+Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"
+DiffEqParamEstim = "1130ab10-4a5a-5621-a13d-e4788d82bd4c"
+DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa"
+GraphViz = "f526b714-d49f-11e8-06ff-31ed36ee7ee0"
+JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
+LatinHypercubeSampling = "a5e1c1ea-c99a-51d3-a14d-a9a37257b02d"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
+Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
+RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd"
+SimpleDiffEq = "05bca326-078c-5bf0-a5bf-ce7c7982d7fd"
+StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
+Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
+Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"

src/CRNUtils.jl → SynthEvo/src/CRNUtils.jl

@@ -1,82 +1,4 @@
 
-# function V(basename, zzz)
-#     """
-#     Function to create a variable given a basename and a number/suffix
-
-#     Args:
-#     - basename: the base name of the variable
-#     - zzz: the suffix of the variable
-
-#     Returns:
-#     - the variable
-#     """
-#     return Meta.parse("$(basename)_$(zzz)")
-# end
-
-### events and integration
-
-function make_perturbation_event_random_step_function(input, intensity)
-    """
-    Define a callback function for a pertubation event.
-
-    Args:
-    - `input`: the input to the system
-    - `intensity`: the intensity of the perturbation
-
-    Returns:
-    - the callback function
-    """
-    function perturb!(integrator, input = input, intensity = intensity)
-        integrator[:U] = max(0., input + (rand()*intensity-intensity/2) )
-    end
-    return perturb!
-end
-
-
-function integrate(prob, tspan, callback, when, algorithm, reltol=1e-8, abstol=1e-8, maxiters=100)
-    """
-    Integrate the problem with a callback function 
-    that will be called only halfway through the time span.
-
-    Args:
-    - `prob`: the ODE problem
-    - `tspan`: the time span
-    - `callback`: the callback function
-    - `when`: the time when the callback should be called
-    - `algorithm`: the algorithm to use Tsit5() 
-    - `reltol`: the relative tolerance
-    - `abstol`: the absolute tolerance
-    - `maxiters`: the maximum number of iterations
-
-    Returns:
-    - the result of the numerical integration
-    """
-    condition = when #[(tspan[2] - tspan[1])/2]
-    ps_cb = PresetTimeCallback(condition, callback)
-    sol = solve(prob, algorithm, reltol=reltol, abstol=abstol, callback=ps_cb, maxiters=maxiters)
-end
-
-### parameters and variables
-
-# function assemble_opt_parameters_and_varables(p, N)
-#     """
-#     Assemble the optimization parameters and variables.
-#     It will return a 0. initial condition for all the species.
-
-#     Args:
-#     - `p`: the parameters
-#     - `N`: the number of species
-
-#     Returns:
-#     - a named tuple with the parameters and the initial conditions (p, u0)
-#     """
-#     np = count_parameters(N)
-#     p = Dict([Meta.parse(string("k_", i)) => p[i] for i in 1:np])
-#     p = push!(p, Meta.parse("U") => 1.)
-#     u0 = [Meta.parse(string("x_", i)) => 0. for i in 1:N]
-#     return (p = p, u0 = u0)
-# end
-
 function perturbation_events(input, perturbation_list)
     """
     Create a list of perturbation events.
@@ -87,46 +9,7 @@ function perturbation_events(input, perturbation_list)
     return [max.(0, [input + rand()*p - p/2]) for p in perturbation_list]
 end
 
-# function dictator(C, p)
-#     dict = Dict()
-#     for i in 1:length(C.parameters)
-#         dict[C.parameters[i]] = p[i]
-#     end
-#     return dict
-# end
 
-# using ModelingToolkit
-
-# function reorder(C, par_v)
-#     dict = Dict()
-#     for i in 1:length(ModelingToolkit.parameters(C.ext_ode))
-#         dict[ModelingToolkit.parameters(C.ext_ode)[i]] = par_v[i]
-#     end
-#     out = []
-#     for i in 1:length(C.parameters)
-#         push!(out, dict[C.parameters[i]])
-#     end
-#     return out
-# end
-
-
-# function dictator_u(C, u)
-#     dict = Dict()
-#     for i in 1:length(C.species)
-#         dict[C.species[i]] = u[i]
-#     end
-#     for i in 1:length(C.control)
-#         dict[C.control[i]] = u[i+length(C.species)]
-#     end
-
-#     for i in 1:size(C.sensitivity_variables)[1]
-#         for j in 1:size(C.sensitivity_variables)[2]
-#             dict[C.sensitivity_variables[i,j]] = u[(i-1)*size(C.sensitivity_variables)[1]+j+length(C.species)+length(C.control)]
-#         end
-#     end
-#     println(dict)
-#     return dict
-# end
 
 function make_base_problem(C, ode, state0, control0, parameters, t_shared, verbose=true, where0=nothing)
     """
@@ -258,4 +141,126 @@ function sortAndFilterReactions(C, p, threshold)
         reaction_selection = [(p[indexes[i]-1], picked_reactions[i]) for i in order],
         filtered_parameters = [ x > threshold ? x : 0. for x in p ]
     )
-end
+end
+
+### --- discarded functions --- ###
+
+# function dictator(C, p)
+#     dict = Dict()
+#     for i in 1:length(C.parameters)
+#         dict[C.parameters[i]] = p[i]
+#     end
+#     return dict
+# end
+
+# using ModelingToolkit
+
+# function reorder(C, par_v)
+#     dict = Dict()
+#     for i in 1:length(ModelingToolkit.parameters(C.ext_ode))
+#         dict[ModelingToolkit.parameters(C.ext_ode)[i]] = par_v[i]
+#     end
+#     out = []
+#     for i in 1:length(C.parameters)
+#         push!(out, dict[C.parameters[i]])
+#     end
+#     return out
+# end
+
+
+# function dictator_u(C, u)
+#     dict = Dict()
+#     for i in 1:length(C.species)
+#         dict[C.species[i]] = u[i]
+#     end
+#     for i in 1:length(C.control)
+#         dict[C.control[i]] = u[i+length(C.species)]
+#     end
+
+#     for i in 1:size(C.sensitivity_variables)[1]
+#         for j in 1:size(C.sensitivity_variables)[2]
+#             dict[C.sensitivity_variables[i,j]] = u[(i-1)*size(C.sensitivity_variables)[1]+j+length(C.species)+length(C.control)]
+#         end
+#     end
+#     println(dict)
+#     return dict
+# end
+
+
+# function V(basename, zzz)
+#     """
+#     Function to create a variable given a basename and a number/suffix
+
+#     Args:
+#     - basename: the base name of the variable
+#     - zzz: the suffix of the variable
+
+#     Returns:
+#     - the variable
+#     """
+#     return Meta.parse("$(basename)_$(zzz)")
+# end
+
+### events and integration
+
+# function make_perturbation_event_random_step_function(input, intensity)
+#     """
+#     Define a callback function for a pertubation event.
+
+#     Args:
+#     - `input`: the input to the system
+#     - `intensity`: the intensity of the perturbation
+
+#     Returns:
+#     - the callback function
+#     """
+#     function perturb!(integrator, input = input, intensity = intensity)
+#         integrator[:U] = max(0., input + (rand()*intensity-intensity/2) )
+#     end
+#     return perturb!
+# end
+
+
+# function integrate(prob, tspan, callback, when, algorithm, reltol=1e-8, abstol=1e-8, maxiters=100)
+#     """
+#     Integrate the problem with a callback function 
+#     that will be called only halfway through the time span.
+
+#     Args:
+#     - `prob`: the ODE problem
+#     - `tspan`: the time span
+#     - `callback`: the callback function
+#     - `when`: the time when the callback should be called
+#     - `algorithm`: the algorithm to use Tsit5() 
+#     - `reltol`: the relative tolerance
+#     - `abstol`: the absolute tolerance
+#     - `maxiters`: the maximum number of iterations
+
+#     Returns:
+#     - the result of the numerical integration
+#     """
+#     condition = when #[(tspan[2] - tspan[1])/2]
+#     ps_cb = PresetTimeCallback(condition, callback)
+#     sol = solve(prob, algorithm, reltol=reltol, abstol=abstol, callback=ps_cb, maxiters=maxiters)
+# end
+
+### parameters and variables
+
+# function assemble_opt_parameters_and_varables(p, N)
+#     """
+#     Assemble the optimization parameters and variables.
+#     It will return a 0. initial condition for all the species.
+
+#     Args:
+#     - `p`: the parameters
+#     - `N`: the number of species
+
+#     Returns:
+#     - a named tuple with the parameters and the initial conditions (p, u0)
+#     """
+#     np = count_parameters(N)
+#     p = Dict([Meta.parse(string("k_", i)) => p[i] for i in 1:np])
+#     p = push!(p, Meta.parse("U") => 1.)
+#     u0 = [Meta.parse(string("x_", i)) => 0. for i in 1:N]
+#     return (p = p, u0 = u0)
+# end

src/CRNevo.jl → SynthEvo/src/CRNevo.jl

@@ -2,7 +2,7 @@
 
 # is_updated means that the fitness value is up to date
 
-function initialize_state(np, ga_options)
+function initialize_state(np, ga_options, mode="random")
     """
     Initialize a pool as defined in ga_options.
     Each parameter vector will have size np.
@@ -13,6 +13,7 @@ function initialize_state(np, ga_options)
     - np: number of parameters
     - ga_options: a struct with at least the following fields:
         - genetic_pool_size: number of parameter vectors
+    -mode: a string that can be "random" or "LHC"
 
     Returns:
     - state: a struct with the following fields:
@@ -23,7 +24,13 @@ function initialize_state(np, ga_options)
             - best_loss: a list of floats
             - mean_loss: a list of floats
     """
-    parameter_pool = [[rand() for _ in 1:np] for _ in 1:ga_options.genetic_pool_size]
+    if mode == "random"
+        parameter_pool = [[rand() for _ in 1:np] for _ in 1:ga_options.genetic_pool_size]
+    elseif mode == "LHC"
+        parameter_pool = sample_parameters_from_LHC(ga_options.genetic_pool_size, np, 10)
+    else 
+        error("mode should be either 'random' or 'LHC'")
+    end
     is_updated = [false for _ in 1:ga_options.genetic_pool_size]
     fitness = [0. for _ in 1:ga_options.genetic_pool_size]
     state = (
@@ -168,7 +175,7 @@ function symbolic_evolve(loss_function_ga, state, follow_gd_fun, ga_options)
         elseif p < ga_options.crossover_rate + ga_options.duplication_rate + ga_options.death_rate + ga_options.gradient_mutation_rate + ga_options.mutation_rate
             crossover_mate = rand(1:length(others_pool))
             crossover_mask = rand(np) .> ga_options.p_cross
-            state.pool[i][:] .= ((1 .- crossover_mask)*state.pool[i]) + state.pool[others_pool[crossover_mate]]*crossover_mask
+            state.pool[i][:] .= ((1 .- crossover_mask).*state.pool[i]) .+ state.pool[others_pool[crossover_mate]].*crossover_mask
             state.is_updated[i] = false
         end
     end

src/NetworkTemplates/FullyConnectedNonExplosive.jl → SynthEvo/src/NetworkTemplates/FullyConnectedNonExplosive.jl

@@ -271,3 +271,26 @@ function create_reactions(C::FullyConnectedCRN)
 end
 
 
+function make_base_problem_for_FCNE_CRN(C, ode, t0, input=1)
+    """
+    Helper function to quickly create a base problem for the fully connected non-explosive CRN model
+    This problem can be then updated efficiently by using the remake function from the ModelingToolkit library
+
+    Args:
+    - C: the CRN model
+    - ode: the ODE system
+    - t0: the initial time
+    - input: the input
+
+    # TODO: maybe this should be part of the CRN struct eventually
+
+    Returns:
+    - the base problem (ODEProblem)
+    """
+    u0 = [0. for _ in 1:243]
+    control0 = [input]
+    p0 = [1. for i in 1:80]
+    base_problem = SynthEvo.make_base_problem(C, ode, u0, control0, p0, t0, true, 0.1)
+    return base_problem
+end
+