diff --git a/.JuliaFormatter.toml b/.JuliaFormatter.toml
index 453925c..012b25b 100644
--- a/.JuliaFormatter.toml
+++ b/.JuliaFormatter.toml
@@ -1 +1,10 @@
-style = "sciml"
\ No newline at end of file
+remove_extra_newlines=true
+long_to_short_function_def=true
+format_docstrings=true
+normalize_line_endings="unix"
+separate_kwargs_with_semicolon=true
+trailing_comma=true
+conditional_to_if=true
+yas_style_nesting=true
+short_to_long_function_def=true
+always_use_return=true
\ No newline at end of file
diff --git a/.github/workflows/CompatHelper.yml b/.github/workflows/CompatHelper.yml
index 87ac2b2..c7f7a64 100644
--- a/.github/workflows/CompatHelper.yml
+++ b/.github/workflows/CompatHelper.yml
@@ -17,7 +17,7 @@ jobs:
       - name: Install Julia, but only if it is not already available in the PATH
         uses: julia-actions/setup-julia@v1
         with:
-          version: '1.8'
+          version: '1.9'
           arch: ${{ runner.arch }}
         if: steps.julia_in_path.outcome != 'success'
       - name: "Add the General registry via Git"
diff --git a/.github/workflows/Documenter.yml b/.github/workflows/Documenter.yml
index 2701a45..20bbd92 100644
--- a/.github/workflows/Documenter.yml
+++ b/.github/workflows/Documenter.yml
@@ -14,7 +14,7 @@ jobs:
       - uses: actions/checkout@v3
       - uses: julia-actions/setup-julia@latest
         with:
-          version: '1.8'
+          version: '1.9'
       - uses: actions/cache@v3
         env:
           cache-name: cache-artifacts
diff --git a/.github/workflows/Tests.yml b/.github/workflows/Tests.yml
index f7c7a74..ddaf9c1 100644
--- a/.github/workflows/Tests.yml
+++ b/.github/workflows/Tests.yml
@@ -11,7 +11,7 @@ jobs:
     strategy:
       fail-fast: false
       matrix:
-        julia_version: ['1.8']
+        julia_version: ['1.9']
         os: [ubuntu-latest]
     steps:
       - uses: actions/checkout@v3
diff --git a/Project.toml b/Project.toml
index 71a7f2f..b1a51d5 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,6 +1,6 @@
 name = "ReactiveDynamics"
 uuid = "c7456e7d-545a-4b79-91ea-6e93d96dd4d4"
-version = "0.2.6"
+version = "0.2.7"
 
 [deps]
 CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
@@ -33,30 +33,31 @@ Tables = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [compat]
-CSV = "0.10"
-Catlab = "0.14.13"
-ComponentArrays = "0.13"
-Crayons = "4.1"
-DataFrames = "1.4"
-DiffEqBase = "6.114"
-DifferentialEquations = "7.6"
-Distributions = "0.25"
-Documenter = "0.27"
+julia = "1.9"
+DifferentialEquations = "7.9"
+StatsFuns = "1.3"
+Catlab = "0.14"
+DataFrames = "1.6"
+PlutoUI = "0.7"
+Statistics = "1.9"
 DocumenterMarkdown = "0.2"
-GeneratedExpressions = "0.1"
-IJulia = "1.24"
+ComponentArrays = "0.14"
 JLD2 = "0.4"
+GeneratedExpressions = "0.1"
+DiffEqBase = "6.128"
 JSON = "0.21"
-MacroTools = "0.5"
-NLopt = "0.6"
-OrdinaryDiffEq = "6.38"
-Plots = "1.38"
-Pluto = "0.19"
-PlutoUI = "0.7"
+NLopt = "1.0"
+OrdinaryDiffEq = "6.55"
+Symbolics = "5.5"
+IJulia = "1.24"
+SafeTestsets = "0.1"
+CSV = "0.10"
+Plots = "1.39"
 Reexport = "1.2"
-SafeTestsets = "0.0"
-StatsFuns = "1.1"
-Symbolics = "4.14, 5"
 TOML = "1.0"
+MacroTools = "0.5"
+Crayons = "4.1"
+Documenter = "0.27"
 Tables = "1.10"
-julia = "1.8"
+Distributions = "0.25"
+Pluto = "0.19"
\ No newline at end of file
diff --git a/docs/make.jl b/docs/make.jl
index 42cdd5f..8638062 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -1,7 +1,10 @@
 using Documenter, DocumenterMarkdown
 using ReactiveDynamics
 
-makedocs(format = Documenter.HTML(prettyurls = false, edit_link = "main"),
-         sitename = "ReactiveDynamics.jl API Documentation", pages = ["index.md"])
+makedocs(;
+    format = Documenter.HTML(; prettyurls = false, edit_link = "main"),
+    sitename = "ReactiveDynamics.jl API Documentation",
+    pages = ["index.md"],
+)
 
-deploydocs(repo = "github.com/Merck/ReactiveDynamics.jl.git")
+deploydocs(; repo = "github.com/Merck/ReactiveDynamics.jl.git")
diff --git a/src/ReactiveDynamics.jl b/src/ReactiveDynamics.jl
index e74daa2..083152f 100644
--- a/src/ReactiveDynamics.jl
+++ b/src/ReactiveDynamics.jl
@@ -8,12 +8,17 @@ using ComponentArrays
 
 @reexport using GeneratedExpressions
 
-const SampleableValues = Union{Expr, Symbol, AbstractString, Float64, Int, Function}
-const ActionableValues = Union{Function, Symbol, Float64, Int}
-
-const SampleableRange = Union{Float64, Int64, AbstractString, Expr, Symbol,
-                              Tuple{Float64,
-                                    Union{Float64, Int64, AbstractString, Expr, Symbol}}}
+const SampleableValues = Union{Expr,Symbol,AbstractString,Float64,Int,Function}
+const ActionableValues = Union{Function,Symbol,Float64,Int}
+
+const SampleableRange = Union{
+    Float64,
+    Int64,
+    AbstractString,
+    Expr,
+    Symbol,
+    Tuple{Float64,Union{Float64,Int64,AbstractString,Expr,Symbol}},
+}
 
 Base.convert(::Type{SampleableRange}, x::Tuple) = (Float64(x[1]), x[2])
 
@@ -26,8 +31,14 @@ end
 @present TheoryReactionNetwork(FreeSchema) begin
     (S, T)::Ob # species, transitions
 
-    (SymbolicAttributeT, DescriptiveAttributeT, SampleableAttributeT,
-     ModalityAttributeT, PcsOptT, PrmAttributeT)::AttrType
+    (
+        SymbolicAttributeT,
+        DescriptiveAttributeT,
+        SampleableAttributeT,
+        ModalityAttributeT,
+        PcsOptT,
+        PrmAttributeT,
+    )::AttrType
 
     specName::Attr(S, SymbolicAttributeT)
     specModality::Attr(S, ModalityAttributeT)
@@ -64,112 +75,150 @@ end
 
 @acset_type FoldedReactionNetworkType(TheoryReactionNetwork)
 
-const ReactionNetwork = FoldedReactionNetworkType{Symbol, Union{String, Symbol, Missing},
-                                                  SampleableValues, Set{Symbol},
-                                                  FoldedObservable, Any}
+const ReactionNetwork = FoldedReactionNetworkType{
+    Symbol,
+    Union{String,Symbol,Missing},
+    SampleableValues,
+    Set{Symbol},
+    FoldedObservable,
+    Any,
+}
 
 Base.convert(::Type{Symbol}, ex::String) = Symbol(ex)
 
-function Base.convert(::Type{Union{String, Symbol, Missing}}, ex::String)
+Base.convert(::Type{Union{String,Symbol,Missing}}, ex::String) =
     try
         Symbol(ex)
     catch
         string(ex)
     end
-end
 
 Base.convert(::Type{SampleableValues}, ex::String) = MacroTools.striplines(Meta.parse(ex))
 Base.convert(::Type{Set{Symbol}}, ex::String) = eval(Meta.parse(ex))
 Base.convert(::Type{FoldedObservable}, ex::String) = eval(Meta.parse(ex))
 
-prettynames = Dict(:transRate => [:rate],
-                   :specInitUncertainty => [:uncertainty, :stoch, :stochasticity],
-                   :transPostAction => [:postAction, :post],
-                   :transName => [:name, :interpretation],
-                   :transPriority => [:priority],
-                   :transProbOfSuccess => [:probability, :prob, :pos],
-                   :transCapacity => [:cap, :capacity],
-                   :transCycleTime => [:ct, :cycletime],
-                   :transMaxLifeTime => [:lifetime, :maxlifetime, :maxtime, :timetolive])
-
-defargs = Dict(:T => Dict{Symbol, Any}(:transPriority => 1, :transProbOfSuccess => 1,
-                                       :transCapacity => Inf, :transCycleTime => 0.0,
-                                       :transMaxLifeTime => Inf, :transMultiplier => 1,
-                                       :transPostAction => :(), :transName => missing),
-               :S => Dict{Symbol, Any}(:specInitUncertainty => 0.0, :specInitVal => 0.0,
-                                       :specCost => 0.0, :specReward => 0.0,
-                                       :specValuation => 0.0),
-               :P => Dict{Symbol, Any}(:prmVal => missing),
-               :M => Dict{Symbol, Any}(:metaVal => missing))
-
-compilable_attrs = filter(attr -> eltype(attr) == SampleableValues,
-                          propertynames(ReactionNetwork()))
+prettynames = Dict(
+    :transRate => [:rate],
+    :specInitUncertainty => [:uncertainty, :stoch, :stochasticity],
+    :transPostAction => [:postAction, :post],
+    :transName => [:name, :interpretation],
+    :transPriority => [:priority],
+    :transProbOfSuccess => [:probability, :prob, :pos],
+    :transCapacity => [:cap, :capacity],
+    :transCycleTime => [:ct, :cycletime],
+    :transMaxLifeTime => [:lifetime, :maxlifetime, :maxtime, :timetolive],
+)
+
+defargs = Dict(
+    :T => Dict{Symbol,Any}(
+        :transPriority => 1,
+        :transProbOfSuccess => 1,
+        :transCapacity => Inf,
+        :transCycleTime => 0.0,
+        :transMaxLifeTime => Inf,
+        :transMultiplier => 1,
+        :transPostAction => :(),
+        :transName => missing,
+    ),
+    :S => Dict{Symbol,Any}(
+        :specInitUncertainty => 0.0,
+        :specInitVal => 0.0,
+        :specCost => 0.0,
+        :specReward => 0.0,
+        :specValuation => 0.0,
+    ),
+    :P => Dict{Symbol,Any}(:prmVal => missing),
+    :M => Dict{Symbol,Any}(:metaVal => missing),
+)
+
+compilable_attrs =
+    filter(attr -> eltype(attr) == SampleableValues, propertynames(ReactionNetwork()))
 
 species_modalities = [:nonblock, :conserved, :rate]
 
 function assign_defaults!(acs::ReactionNetwork)
     for (_, v_) in defargs, (k, v) in v_
-        for i in 1:length(subpart(acs, k))
+        for i = 1:length(subpart(acs, k))
             isnothing(acs[i, k]) && (subpart(acs, k)[i] = v)
         end
     end
 
-    foreach(i -> !isnothing(acs[i, :specModality]) ||
-                (subpart(acs, :specModality)[i] = Set{Symbol}()), 1:nparts(acs, :S))
+    foreach(
+        i ->
+            !isnothing(acs[i, :specModality]) ||
+                (subpart(acs, :specModality)[i] = Set{Symbol}()),
+        1:nparts(acs, :S),
+    )
     k = [:specCost, :specReward, :specValuation]
-    foreach(k -> foreach(i -> !isnothing(acs[i, k]) || (subpart(acs, k)[i] = 0.0),
-                         1:nparts(acs, :S)), k)
-
-    acs
+    foreach(
+        k -> foreach(
+            i -> !isnothing(acs[i, k]) || (subpart(acs, k)[i] = 0.0),
+            1:nparts(acs, :S),
+        ),
+        k,
+    )
+
+    return acs
 end
 
 function ReactionNetwork(transitions, reactants, obs, events)
-    merge_acs!(ReactionNetwork(), transitions, reactants, obs, events)
+    return merge_acs!(ReactionNetwork(), transitions, reactants, obs, events)
 end
 
 function ReactionNetwork(transitions, reactants, obs)
-    merge_acs!(ReactionNetwork(), transitions, reactants, obs, [])
+    return merge_acs!(ReactionNetwork(), transitions, reactants, obs, [])
 end
 
 function add_obs!(acs, obs)
     for p in obs
         sym = p.args[3].value
         i = incident(acs, sym, :obsName)
-        i = isempty(incident(acs, sym, :obsName)) ?
-            add_part!(acs, :obs; obsName = sym, obsOpts = FoldedObservable()) : i[1]
+        i = if isempty(incident(acs, sym, :obsName))
+            add_part!(acs, :obs; obsName = sym, obsOpts = FoldedObservable())
+        else
+            i[1]
+        end
         for opt in p.args[4:end]
             if isexpr(opt, :(=)) && (opt.args[1] ∈ fieldnames(FoldedObservable))
                 opt.args[1] == :every &&
                     (acs[i, :obsOpts].every = min(acs[i, :obsOpts].every, opt.args[2]))
                 opt.args[1] == :on && union!(acs[i, :obsOpts].on, [opt.args[2]])
             elseif isexpr(opt, :tuple) || opt isa SampleableValues
-                push!(acs[i, :obsOpts].range,
-                      isexpr(opt, :tuple) ? tuple(opt.args...) : opt)
+                push!(
+                    acs[i, :obsOpts].range,
+                    isexpr(opt, :tuple) ? tuple(opt.args...) : opt,
+                )
             end
         end
     end
 
-    acs
+    return acs
 end
 
 function merge_acs!(acs::ReactionNetwork, transitions, reactants, obs, events)
-    foreach(t -> add_part!(acs, :T; trans = t[1][2], transRate = t[1][1], t[2]...),
-            transitions)
+    foreach(
+        t -> add_part!(acs, :T; trans = t[1][2], transRate = t[1][1], t[2]...),
+        transitions,
+    )
     add_obs!(acs, obs)
     unique!(reactants)
-    foreach(ev -> add_part!(acs, :E; eventTrigger = ev.trigger, eventAction = ev.action),
-            events)
-    foreach(r -> isempty(incident(acs, r, :specName)) && add_part!(acs, :S; specName = r),
-            reactants)
-
-    assign_defaults!(acs)
+    foreach(
+        ev -> add_part!(acs, :E; eventTrigger = ev.trigger, eventAction = ev.action),
+        events,
+    )
+    foreach(
+        r -> isempty(incident(acs, r, :specName)) && add_part!(acs, :S; specName = r),
+        reactants,
+    )
+
+    return assign_defaults!(acs)
 end
 
 include("state.jl")
 include("compilers.jl")
-include.(readdir(joinpath(@__DIR__, "interface"), join = true))
-include.(readdir(joinpath(@__DIR__, "utils"), join = true))
-include.(readdir(joinpath(@__DIR__, "operators"), join = true))
+include.(readdir(joinpath(@__DIR__, "interface"); join = true))
+include.(readdir(joinpath(@__DIR__, "utils"); join = true))
+include.(readdir(joinpath(@__DIR__, "operators"); join = true))
 include("solvers.jl")
 include("optim.jl")
 include("loadsave.jl")
diff --git a/src/compilers.jl b/src/compilers.jl
index 8288e09..d1b5314 100644
--- a/src/compilers.jl
+++ b/src/compilers.jl
@@ -1,33 +1,50 @@
 using MacroTools: prewalk
 
-"Recursively find model variables in expressions."
+"""
+Recursively find model variables in expressions.
+"""
 function recursively_find_vars!(r, exs...)
     for ex in exs
-        ex isa Symbol ? push!(r, ex) :
-        (ex isa Expr &&
-         recursively_find_vars!(r, ex.args[(!isexpr(ex, :call) ? 1 : 2):end]...))
+        if ex isa Symbol
+            push!(r, ex)
+        else
+            (
+                ex isa Expr &&
+                recursively_find_vars!(r, ex.args[(!isexpr(ex, :call) ? 1 : 2):end]...)
+            )
+        end
     end
-    r
+    return r
 end
 
 function get_contained_params(ex, prms)
     r = []
-    recursively_find_vars!(r, ex) ∩ prms
+    return recursively_find_vars!(r, ex) ∩ prms
 end
 
-"Recursively substitute model variables. Subsitution pairs are specified in `varmap`."
+"""
+Recursively substitute model variables. Subsitution pairs are specified in `varmap`.
+"""
 function recursively_substitute_vars!(varmap, ex)
     ex isa Symbol && return (haskey(varmap, ex) ? varmap[ex] : ex)
-    ex isa Expr && for i in 1:length(ex.args)
-        ex.args[i] isa Expr ? recursively_substitute_vars!(varmap, ex.args[i]) :
-        (ex.args[i] isa Symbol && haskey(varmap, ex.args[i]) &&
-         (ex.args[i] = varmap[ex.args[i]]))
+    ex isa Expr && for i = 1:length(ex.args)
+        if ex.args[i] isa Expr
+            recursively_substitute_vars!(varmap, ex.args[i])
+        else
+            (
+                ex.args[i] isa Symbol &&
+                haskey(varmap, ex.args[i]) &&
+                (ex.args[i] = varmap[ex.args[i]])
+            )
+        end
     end
 
-    ex
+    return ex
 end
 
-"Recursively normalize dotted notation in `ex` to species names whenever a name is contained in `vars`, else left unchanged."
+"""
+Recursively normalize dotted notation in `ex` to species names whenever a name is contained in `vars`, else left unchanged.
+"""
 function recursively_expand_dots_in_ex!(ex, vars)
     if isexpr(ex, :.)
         expanded = recursively_expand_dots(ex)
@@ -38,35 +55,32 @@ function recursively_expand_dots_in_ex!(ex, vars)
             ex
         end
     end
-    ex isa Expr && for i in 1:length(ex.args)
-        ex.args[i] isa Union{Expr, Symbol} &&
+    ex isa Expr && for i = 1:length(ex.args)
+        ex.args[i] isa Union{Expr,Symbol} &&
             (ex.args[i] = recursively_expand_dots_in_ex!(ex.args[i], vars))
     end
 
-    ex
+    return ex
 end
 
-reserved_names = [
-    :t,
-    :state,
-    :obs,
-    :resample,
-    :solverarg,
-    :take,
-    :log,
-    :periodic,
-    :set_params,
-]
+reserved_names =
+    [:t, :state, :obs, :resample, :solverarg, :take, :log, :periodic, :set_params]
 push!(reserved_names, :state)
 
 function escape_ref(ex, species)
-    ex isa Symbol ? ex :
-    prewalk(ex -> isexpr(ex, :ref) && Symbol(string(ex)) ∈ species ? Symbol(string(ex)) : ex,
-            ex)
+    return if ex isa Symbol
+        ex
+    else
+        prewalk(
+            ex ->
+                isexpr(ex, :ref) && Symbol(string(ex)) ∈ species ? Symbol(string(ex)) : ex,
+            ex,
+        )
+    end
 end
 
 function wrap_expr(fex, species_names, prm_names, varmap)
-    !isa(fex, Union{Expr, Symbol}) && return fex
+    !isa(fex, Union{Expr,Symbol}) && return fex
     # escape refs in species names: A[1] -> Symbol("A[1]")
     fex = escape_ref(fex, species_names)
     # escape dots in species' names: A.B -> Symbol("A.B")
@@ -74,8 +88,10 @@ function wrap_expr(fex, species_names, prm_names, varmap)
     fex = recursively_expand_dots_in_ex!(fex, species_names)
 
     # prepare the function's body
-    letex = :(let
-              end)
+    letex = :(
+        let
+        end
+    )
     # expression walking (MacroTools): visit each expression, subsitute with the body's return value
     fex = prewalk(fex) do x
         # here we convert the query metalanguage: @t() -> time(state) etc. 
@@ -90,13 +106,15 @@ function wrap_expr(fex, species_names, prm_names, varmap)
     fex = recursively_substitute_vars!(varmap, fex)
     # substitute the params names with "pointers" into the parameter space: β -> state.p[:β]
     # params can't be mutated!
-    foreach(v -> push!(letex.args[1].args, :($v = state.p[$(QuoteNode(v))])),
-            get_contained_params(fex, prm_names))
+    foreach(
+        v -> push!(letex.args[1].args, :($v = state.p[$(QuoteNode(v))])),
+        get_contained_params(fex, prm_names),
+    )
     push!(letex.args[2].args, fex)
 
     # the function shall be a function of the dynamic ReactiveDynamicsState structure: letex -> :(state -> $letex)
     # eval the expression to a Julia function, save that function into the "compiled" acset
-    eval(:(state -> $letex))
+    return eval(:(state -> $letex))
 end
 
 function get_wrap_fun(acs::ReactionNetwork)
@@ -107,11 +125,12 @@ function get_wrap_fun(acs::ReactionNetwork)
         push!(varmap, name => :(state.p[$(QuoteNode(name))]))
     end
 
-    ex -> wrap_expr(ex, species_names, prm_names, varmap)
+    return ex -> wrap_expr(ex, species_names, prm_names, varmap)
 end
 
 function skip_compile(attr)
-    any(contains.(Ref(string(attr)), ("Name", "obs", "meta"))) || (string(attr) == "trans")
+    return any(contains.(Ref(string(attr)), ("Name", "obs", "meta"))) ||
+           (string(attr) == "trans")
 end
 
 function compile_attrs(acs::ReactionNetwork)
@@ -122,22 +141,32 @@ function compile_attrs(acs::ReactionNetwork)
         push!(varmap, name => :(state.p[$(QuoteNode(name))]))
     end
     wrap_fun = ex -> wrap_expr(ex, species_names, prm_names, varmap)
-    attrs = Dict{Symbol, Vector}()
-    transitions = Dict{Symbol, Vector}()
+    attrs = Dict{Symbol,Vector}()
+    transitions = Dict{Symbol,Vector}()
     for attr in propertynames(acs.subparts)
         attrs_ = subpart(acs, attr)
-        !contains(string(attr), "trans") ?
-        (attrs[attr] = map(i -> skip_compile(attr) ? attrs_[i] : wrap_fun(attrs_[i]),
-                           1:length(attrs_))) :
-        (transitions[attr] = map(i -> skip_compile(attr) ? attrs_[i] : wrap_fun(attrs_[i]),
-                                 1:length(attrs_)))
+        if !contains(string(attr), "trans")
+            (
+                attrs[attr] = map(
+                    i -> skip_compile(attr) ? attrs_[i] : wrap_fun(attrs_[i]),
+                    1:length(attrs_),
+                )
+            )
+        else
+            (
+                transitions[attr] = map(
+                    i -> skip_compile(attr) ? attrs_[i] : wrap_fun(attrs_[i]),
+                    1:length(attrs_),
+                )
+            )
+        end
     end
     transitions[:transActivated] = fill(true, nparts(acs, :T))
     transitions[:transToSpawn] = zeros(nparts(acs, :T))
-    transitions[:transHash] = [coalesce(acs[i, :transName], gensym())
-                               for i in 1:nparts(acs, :T)]
+    transitions[:transHash] =
+        [coalesce(acs[i, :transName], gensym()) for i = 1:nparts(acs, :T)]
 
-    attrs, transitions, wrap_fun
+    return attrs, transitions, wrap_fun
 end
 
 function remove_choose(acs::ReactionNetwork)
@@ -145,10 +174,15 @@ function remove_choose(acs::ReactionNetwork)
     pcs = []
     for attr in propertynames(acs.subparts)
         attrs_ = subpart(acs, attr)
-        foreach(i -> !isnothing(attrs_[i]) && attrs_[i] isa Expr &&
-                         (attrs_[i] = normalize_pcs!(pcs, attrs_[i])), 1:length(attrs_))
+        foreach(
+            i ->
+                !isnothing(attrs_[i]) &&
+                    attrs_[i] isa Expr &&
+                    (attrs_[i] = normalize_pcs!(pcs, attrs_[i])),
+            1:length(attrs_),
+        )
     end
 
     add_obs!(acs, pcs)
-    acs
+    return acs
 end
diff --git a/src/interface/create.jl b/src/interface/create.jl
index 796d91e..1a82374 100644
--- a/src/interface/create.jl
+++ b/src/interface/create.jl
@@ -7,22 +7,8 @@ using Symbolics: build_function, get_variables
 
 empty_set = Set{Symbol}([:∅])
 fwd_arrows = Set{Symbol}([:>, :→, :↣, :↦, :⇾, :⟶, :⟼, :⥟, :⥟, :⇀, :⇁, :⇒, :⟾])
-bwd_arrows = Set{Symbol}([
-                             :<,
-                             :←,
-                             :↢,
-                             :↤,
-                             :⇽,
-                             :⟵,
-                             :⟻,
-                             :⥚,
-                             :⥞,
-                             :↼,
-                             :↽,
-                             :⇐,
-                             :⟽,
-                             Symbol("<--"),
-                         ])
+bwd_arrows =
+    Set{Symbol}([:<, :←, :↢, :↤, :⇽, :⟵, :⟻, :⥚, :⥞, :↼, :↽, :⇐, :⟽, Symbol("<--")])
 double_arrows = Set{Symbol}([:↔, :⟷, :⇄, :⇆, :⇌, :⇋, :⇔, :⟺, Symbol("<-->")])
 
 arrows = fwd_arrows ∪ bwd_arrows ∪ double_arrows ∪ [:-->]
@@ -57,39 +43,42 @@ Most arrows accepted (both right, left, and bi-drectional arrows). Use 0 or ∅
 Custom functions and sampleable objects can be used as numeric parameters. Note that these have to be accessible from ReactiveDynamics's source code.
 
 # Examples
+
 ```julia
 acs = @ReactionNetwork begin
     1.0, X ⟶ Y
-    1.0, X ⟶ Y, priority=>6., prob=>.7, capacity=>3.
-    1.0, ∅ --> (Poisson(.3γ)X, Poisson(.5)Y)
+    1.0, X ⟶ Y, priority => 6.0, prob => 0.7, capacity => 3.0
+    1.0, ∅ --> (Poisson(0.3γ)X, Poisson(0.5)Y)
     (XY > 100) && (XY -= 1)
 end
-@push acs 1.0 X ⟶ Y 
-@prob_init acs X=1 Y=2 XY=α
-@prob_params acs γ=1 α=4
+@push acs 1.0 X ⟶ Y
+@prob_init acs X = 1 Y = 2 XY = α
+@prob_params acs γ = 1 α = 4
 @solve_and_plot acs
 ```
 """
 macro ReactionNetwork end
 
 macro ReactionNetwork()
-    make_ReactionNetwork(:())
+    return make_ReactionNetwork(:())
 end
 
 macro ReactionNetwork(ex)
-    make_ReactionNetwork(ex; eval_module = __module__)
+    return make_ReactionNetwork(ex; eval_module = __module__)
 end
 
 macro ReactionNetwork(ex, args...)
-    make_ReactionNetwork(generate(Expr(:braces, ex, args...); eval_module = __module__);
-                         eval_module = __module__)
+    return make_ReactionNetwork(
+        generate(Expr(:braces, ex, args...); eval_module = __module__);
+        eval_module = __module__,
+    )
 end
 
 function make_ReactionNetwork(ex::Expr; eval_module = @__MODULE__)
     blockex = generate(ex; eval_module)
     blockex = unblock_shallow!(blockex)
 
-    :(ReactionNetwork(get_data($(QuoteNode(blockex)))...))
+    return :(ReactionNetwork(get_data($(QuoteNode(blockex)))...))
 end
 
 ### Functions that process the input and rephrase it as a reaction system ###
@@ -98,12 +87,12 @@ function esc_dollars!(ex)
         if ex.head == :$
             return esc(:($(ex.args[1])))
         else
-            for i in 1:length(ex.args)
+            for i = 1:length(ex.args)
                 ex.args[i] = esc_dollars!(ex.args[i])
             end
         end
     end
-    ex
+    return ex
 end
 
 symbolize(pairex) = pairex isa Number ? pairex : (pairex.args[2] => pairex.args[3])
@@ -116,29 +105,33 @@ function get_data(ex)
     if isexpr(ex, :block)
         ex = striplines(ex)
         esc_dollars!(ex)
-        foreach(l -> get_data!(trans, reactants, pcs, evs, isexpr(l, :tuple) ? l.args : [l]),
-                ex.args)
+        foreach(
+            l -> get_data!(trans, reactants, pcs, evs, isexpr(l, :tuple) ? l.args : [l]),
+            ex.args,
+        )
     elseif ex != :()
         get_data!(trans, reactants, pcs, evs, ex)
     end
 
-    trans, reactants, pcs, evs
+    return trans, reactants, pcs, evs
 end
 
 function get_data!(trans, reactants, pcs, evs, exs)
-    (length(exs) == 0 && return;
-     exs[1] isa Expr && (exs[1].head ∈ ifs) ?
-     get_events!(evs, normalize_pcs!(pcs, exs[1])) :
-     get_transitions!(trans, reactants, pcs, exs))
+    length(exs) == 0 && return
+
+    if exs[1] isa Expr && (exs[1].head ∈ ifs)
+        get_events!(evs, normalize_pcs!(pcs, exs[1]))
+    else
+        get_transitions!(trans, reactants, pcs, exs)
+    end
 end
 
-function get_events!(evs, ex)
+get_events!(evs, ex) =
     if ex.head == :&&
         push!(evs, Event(ex.args...))
     else
         recursively_expand_actions!(evs, Expr(:call, :&), ex)
     end
-end
 
 function recursively_expand_actions!(evs, condex, event)
     if isexpr(event, :if)
@@ -181,49 +174,56 @@ function get_transitions!(trans, reactants, pcs, exs)
     while ix <= length(exs)
         (!isa(exs[ix], Expr) || (exs[ix].head != :call)) && (ix += 1; continue)
         karg = (xi = findfirst(k -> exs[ix].args[2] ∈ k, prettynames);
-                isnothing(xi) &&
-                    (ix += 1; continue);
-                xi)
+        isnothing(xi) && (ix += 1; continue);
+        xi)
         push!(args, karg => normalize_pcs!(pcs, exs[ix].args[3]))
         deleteat!(exs, ix)
     end
     args = merge(defargs[:T], args)
 
     append!(trans, tuple.(rxs, Ref(args)))
-    trans
+    return trans
 end
 
 function replace_in_expr(expr, pairs...)
     dict = Dict(pairs...)
 
-    prewalk(ex -> haskey(dict, ex) ? dict[ex] : ex, expr)
+    return prewalk(ex -> haskey(dict, ex) ? dict[ex] : ex, expr)
 end
 
 function normalize_pcs!(pcs, expr)
     postwalk(expr) do ex
-        isexpr(ex, :macrocall) && macroname(ex) == :register &&
-            (push!(pcs, deepcopy(ex)); ex.args[1] = Symbol("@", :take); ex.args = ex.args[1:3])
+        isexpr(ex, :macrocall) &&
+            macroname(ex) == :register &&
+            (push!(pcs, deepcopy(ex));
+            ex.args[1] = Symbol("@", :take);
+            ex.args = ex.args[1:3])
         if isexpr(ex, :macrocall) && macroname(ex) == :register
             r_sym = gensym()
             (push!(pcs, (ex_ = deepcopy(ex); insert!(ex_.args, 3, r_sym); ex_));
-             ex.args[1] = Symbol("@",
-                                 :take);
-             ex.args = [r_sym;
-                        ex.args[1:2]])
+            ex.args[1] = Symbol("@", :take);
+            ex.args = [
+                r_sym
+                ex.args[1:2]
+            ])
         end
-        ex
+        return ex
     end
 end
 
 function get_reaction_line(expr)
     biarrow = nothing
     prewalk(ex --> (ex ∈ double_arrows && (biarrow = ex); ex), expr)
-    !isnothing(biarrow) ? [expr] :
-    [replace_in_expr(expr, biarrow => :⟶), replace_in_expr(expr, biarrow => :⟵)]
+    return if !isnothing(biarrow)
+        [expr]
+    else
+        [replace_in_expr(expr, biarrow => :⟶), replace_in_expr(expr, biarrow => :⟵)]
+    end
 end
 
 function prune_reaction_line!(pcs, reactants, line)
-    line isa Expr && (line.head == :-->) &&
+    line isa Expr &&
+        (line.head == :-->) &&
         (line = Expr(:call, :→, line.args[1], line.args[2]))
     if isexpr(line, :macrocall)
         lines = copy(line.args[3:end])
@@ -231,28 +231,44 @@ function prune_reaction_line!(pcs, reactants, line)
         for l in lines
             biarrow = nothing
             prewalk(ex -> (ex ∈ double_arrows && (biarrow = ex); ex), l)
-            append!(line.args,
-                    isnothing(biarrow) ? [l] :
-                    [replace_in_expr(l, biarrow => :⟶), replace_in_expr(l, biarrow => :⟵)])
+            append!(
+                line.args,
+                if isnothing(biarrow)
+                    [l]
+                else
+                    [replace_in_expr(l, biarrow => :⟶), replace_in_expr(l, biarrow => :⟵)]
+                end,
+            )
         end
-        for i in 3:length(line.args)
-            line.args[i] = isexpr(line.args[i], :tuple) ?
-                           Expr(:tuple, line.args[i].args[1],
-                                prune_reaction_line!(pcs, reactants, line.args[i].args[2])) :
-                           prune_reaction_line!(pcs, reactants, line.args[i])
+        for i = 3:length(line.args)
+            line.args[i] = if isexpr(line.args[i], :tuple)
+                Expr(
+                    :tuple,
+                    line.args[i].args[1],
+                    prune_reaction_line!(pcs, reactants, line.args[i].args[2]),
+                )
+            else
+                prune_reaction_line!(pcs, reactants, line.args[i])
+            end
         end
     elseif line isa Expr && line.args[1] ∈ union(fwd_arrows, bwd_arrows)
-        line.args[2:3] = recursively_find_reactants!.(Ref(reactants), Ref(pcs),
-                                                      line.args[2:3])
+        line.args[2:3] =
+            recursively_find_reactants!.(Ref(reactants), Ref(pcs), line.args[2:3])
     elseif line isa Expr && line.args[1] ∈ double_arrows
         biarrow = nothing
         prewalk(ex -> (ex ∈ double_arrows && (biarrow = ex); ex), line)
-        line = prune_reaction_line!.(Ref(pcs), Ref(reactants),
-                                     (replace_in_expr(line, biarrow => :⟶),
-                                      replace_in_expr(line, biarrow => :⟵)))
+        line =
+            prune_reaction_line!.(
+                Ref(pcs),
+                Ref(reactants),
+                (
+                    replace_in_expr(line, biarrow => :⟶),
+                    replace_in_expr(line, biarrow => :⟵),
+                ),
+            )
     end
 
-    line
+    return line
 end
 
 function recursively_find_reactants!(reactants, pcs, ex::SampleableValues)
@@ -264,16 +280,18 @@ function recursively_find_reactants!(reactants, pcs, ex::SampleableValues)
         end
     elseif ex.args[1] == :*
         recursively_find_reactants!(reactants, pcs, ex.args[end])
-        foreach(i -> ex.args[i] = normalize_pcs!(pcs, ex.args[i]), 2:(length(ex.args) - 1))
+        foreach(i -> ex.args[i] = normalize_pcs!(pcs, ex.args[i]), 2:(length(ex.args)-1))
     elseif ex.args[1] == :+
-        for i in 2:length(ex.args)
+        for i = 2:length(ex.args)
             recursively_find_reactants!(reactants, pcs, ex.args[i])
         end
     elseif isexpr(ex, :macrocall) && macroname(ex) == :choose
-        for i in 3:length(ex.args)
-            recursively_find_reactants!(reactants, pcs,
-                                        isexpr(ex.args[i], :tuple) ? ex.args[i].args[2] :
-                                        ex.args[i])
+        for i = 3:length(ex.args)
+            recursively_find_reactants!(
+                reactants,
+                pcs,
+                isexpr(ex.args[i], :tuple) ? ex.args[i].args[2] : ex.args[i],
+            )
         end
     elseif isexpr(ex, :macrocall)
         recursively_find_reactants!(reactants, pcs, ex.args[3])
@@ -281,5 +299,5 @@ function recursively_find_reactants!(reactants, pcs, ex::SampleableValues)
         push!(reactants, underscorize(ex))
     end
 
-    ex
+    return ex
 end
diff --git a/src/interface/parsing_utils.jl b/src/interface/parsing_utils.jl
index 5336bfb..51d765b 100644
--- a/src/interface/parsing_utils.jl
+++ b/src/interface/parsing_utils.jl
@@ -18,20 +18,24 @@ function multiplex(mult, mults...)
     all(m -> isa(m, Number), [mult] ∪ mults) && return mult * prod(mults; init = 1.0)
     multarray = SampleableValues[]
     recursively_find_mults!(multarray, mults...)
-    mults_numeric = prod(filter(m -> isa(m, Number), multarray); init = 1.0) *
-                    (mult isa Number ? mult : 1.0)
+    mults_numeric =
+        prod(filter(m -> isa(m, Number), multarray); init = 1.0) *
+        (mult isa Number ? mult : 1.0)
     mults_expr = filter(m -> !isa(m, Number), multarray)
     mult = mult isa Expr ? deepcopy(mult) : :(*())
     mults_numeric == 1 && length(mults_expr) == 1 && return mults_expr[1]
     mults_numeric != 1.0 && push!(mult.args, mults_numeric)
     append!(mult.args, mults_expr)
 
-    mult
+    return mult
 end
 
 function recursively_find_mults!(multarray, mults...)
     for m in mults
-        isa(m, Expr) && m.args[1] == :* ?
-        recursively_find_mults!(multarray, m.args[2:end]...) : push!(multarray, m)
+        if isa(m, Expr) && m.args[1] == :*
+            recursively_find_mults!(multarray, m.args[2:end]...)
+        else
+            push!(multarray, m)
+        end
     end
 end
diff --git a/src/interface/reaction_parser.jl b/src/interface/reaction_parser.jl
index e79078e..ae92ae3 100644
--- a/src/interface/reaction_parser.jl
+++ b/src/interface/reaction_parser.jl
@@ -11,9 +11,18 @@ end
 function recursively_choose(r_line, state)
     postwalk(r_line) do ex
         if isexpr(ex, :macrocall) && (macroname(ex) == :choose)
-            sample_range([(isexpr(r, :tuple) ?
-                           (r.args[1], recursively_choose(r.args[2], state)) :
-                           recursively_choose(r, state)) for r in ex.args[3:end]], state)
+            sample_range(
+                [
+                    (
+                        if isexpr(r, :tuple)
+                            (r.args[1], recursively_choose(r.args[2], state))
+                        else
+                            recursively_choose(r, state)
+                        end
+                    ) for r in ex.args[3:end]
+                ],
+                state,
+            )
         else
             ex
         end
@@ -23,12 +32,20 @@ end
 function extract_reactants(r_line, state::ReactiveDynamicsState)
     r_line = recursively_choose(r_line, state)
 
-    recursive_find_reactants!(escape_ref(r_line, state[:, :specName]), 1.0, Set{Symbol}(),
-                              Vector{FoldedReactant}(undef, 0))
+    return recursive_find_reactants!(
+        escape_ref(r_line, state[:, :specName]),
+        1.0,
+        Set{Symbol}(),
+        Vector{FoldedReactant}(undef, 0),
+    )
 end
 
-function recursive_find_reactants!(ex::SampleableValues, mult::SampleableValues,
-                                   mods::Set{Symbol}, reactants::Vector{FoldedReactant})
+function recursive_find_reactants!(
+    ex::SampleableValues,
+    mult::SampleableValues,
+    mods::Set{Symbol},
+    reactants::Vector{FoldedReactant},
+)
     if typeof(ex) != Expr || isexpr(ex, :.) || (ex.head == :escape)
         if (ex == 0 || in(ex, empty_set))
             return reactants
@@ -36,21 +53,27 @@ function recursive_find_reactants!(ex::SampleableValues, mult::SampleableValues,
             push!(reactants, FoldedReactant(recursively_expand_dots(ex), mult, mods))
         end
     elseif ex.args[1] == :*
-        recursive_find_reactants!(ex.args[end], multiplex(mult, ex.args[2:(end - 1)]...),
-                                  mods, reactants)
+        recursive_find_reactants!(
+            ex.args[end],
+            multiplex(mult, ex.args[2:(end-1)]...),
+            mods,
+            reactants,
+        )
     elseif ex.args[1] == :+
-        for i in 2:length(ex.args)
+        for i = 2:length(ex.args)
             recursive_find_reactants!(ex.args[i], mult, mods, reactants)
         end
     elseif ex.head == :macrocall
         mods = copy(mods)
         macroname(ex) in species_modalities && push!(mods, macroname(ex))
-        foreach(i -> push!(mods, ex.args[i] isa Symbol ? ex.args[i] : ex.args[i].value),
-                4:length(ex.args))
+        foreach(
+            i -> push!(mods, ex.args[i] isa Symbol ? ex.args[i] : ex.args[i].value),
+            4:length(ex.args),
+        )
         recursive_find_reactants!(ex.args[3], mult, mods, reactants)
     else
         @error("malformed reaction")
     end
 
-    reactants
+    return reactants
 end
diff --git a/src/interface/solve.jl b/src/interface/solve.jl
index 0ce096c..a7484e9 100644
--- a/src/interface/solve.jl
+++ b/src/interface/solve.jl
@@ -9,15 +9,19 @@ import Plots
 Convert a model to a `DiscreteProblem`. If passed a problem instance, return the instance.
 
 # Examples
+
 ```julia
-@problematize acs tspan=1:100
+@problematize acs tspan = 1:100
 ```
 """
 macro problematize(acsex, args...)
     args, kwargs = args_kwargs(args)
     quote
-        $(esc(acsex)) isa DiscreteProblem ? $(esc(acsex)) :
-        DiscreteProblem($(esc(acsex)), $(args...); $(kwargs...))
+        if $(esc(acsex)) isa DiscreteProblem
+            $(esc(acsex))
+        else
+            DiscreteProblem($(esc(acsex)), $(args...); $(kwargs...))
+        end
     end
 end
 
@@ -25,10 +29,11 @@ end
 Solve the problem. Solverargs passed at the calltime take precedence.
 
 # Examples
+
 ```julia
 @solve prob
-@solve prob tspan=1:100
-@solve prob tspan=100 trajectories=20
+@solve prob tspan = 1:100
+@solve prob tspan = 100 trajectories = 20
 ```
 """
 macro solve(probex, args...)
@@ -37,11 +42,17 @@ macro solve(probex, args...)
     !isnothing(findfirst(el -> el.args[1] == :trajectories, kwargs)) && (mode = :ensemble)
 
     quote
-        prob = $(esc(probex)) isa DiscreteProblem ? $(esc(probex)) :
-               DiscreteProblem($(esc(probex)), $(args...); $(kwargs...))
+        prob = if $(esc(probex)) isa DiscreteProblem
+            $(esc(probex))
+        else
+            DiscreteProblem($(esc(probex)), $(args...); $(kwargs...))
+        end
         if $(preserve_sym(mode)) == :ensemble
-            solve(EnsembleProblem(prob; prob_func = get_prob_func(prob)), FunctionMap(),
-                  $(kwargs...))
+            solve(
+                EnsembleProblem(prob; prob_func = get_prob_func(prob)),
+                FunctionMap(),
+                $(kwargs...),
+            )
         else
             solve(prob)
         end
@@ -53,23 +64,38 @@ function plot_summary(s, labels, ixs; kwargs...)
     isempty(ixs) && return @warn "Set of species to plot must be non-empty!"
     s = EnsembleSummary(s)
     f_ix = first(ixs)
-    p = Plots.plot(s.t, s.u[f_ix, :],
-                   ribbon = (-s.qlow[f_ix, :] + s.u[f_ix, :],
-                             s.qhigh[f_ix, :] - s.u[f_ix, :]), label = labels[f_ix],
-                   fillalpha = 0.2, w = 2.0; kwargs...)
-    foreach(i -> Plots.plot!(p, s.t, s.u[i, :],
-                             ribbon = (-s.qlow[i, :] + s.u[i, :],
-                                       s.qhigh[i, :] - s.u[i, :]), label = labels[i],
-                             fillalpha = 0.2, w = 2.0; kwargs...),
-            ixs[2:end])
-
-    p
+    p = Plots.plot(
+        s.t,
+        s.u[f_ix, :];
+        ribbon = (-s.qlow[f_ix, :] + s.u[f_ix, :], s.qhigh[f_ix, :] - s.u[f_ix, :]),
+        label = labels[f_ix],
+        fillalpha = 0.2,
+        w = 2.0,
+        kwargs...,
+    )
+    foreach(
+        i -> Plots.plot!(
+            p,
+            s.t,
+            s.u[i, :];
+            ribbon = (-s.qlow[i, :] + s.u[i, :], s.qhigh[i, :] - s.u[i, :]),
+            label = labels[i],
+            fillalpha = 0.2,
+            w = 2.0,
+            kwargs...,
+        ),
+        ixs[2:end],
+    )
+
+    return p
 end
 
 function plot_ensemble_sol(sol, label, ixs; kwargs...)
-    !(sol isa EnsembleSolution) ?
-    Plots.plot(sol, idxs = ixs, label = reshape(label[ixs], 1, :); kwargs...) :
-    Plots.plot(sol, idxs = ixs, alpha = 0.7; kwargs...)
+    return if !(sol isa EnsembleSolution)
+        Plots.plot(sol; idxs = ixs, label = reshape(label[ixs], 1, :), kwargs...)
+    else
+        Plots.plot(sol; idxs = ixs, alpha = 0.7, kwargs...)
+    end
 end
 
 first_sol(sol) = sol isa EnsembleSolution ? first(sol) : sol
@@ -78,13 +104,17 @@ function match_names(selector, names)
     if isnothing(selector)
         1:length(names)
     else
-        selector = selector isa Union{AbstractString, Regex, Symbol} ? [selector] : selector
+        selector = selector isa Union{AbstractString,Regex,Symbol} ? [selector] : selector
         ixs = Int[]
         for s in selector
             s isa Symbol && (s = string(s))
-            append!(ixs,
-                    findall(name -> s isa Regex ? occursin(s, name) :
-                                    (string(s) == string(name)), names))
+            append!(
+                ixs,
+                findall(
+                    name -> s isa Regex ? occursin(s, name) : (string(s) == string(name)),
+                    names,
+                ),
+            )
         end
         unique!(ixs)
     end
@@ -94,11 +124,12 @@ end
 Plot the solution (summary).
 
 # Examples
+
 ```julia
-@plot sol plot_type=summary
-@plot sol plot_type=allocation # not supported for ensemble solutions!
-@plot sol plot_type=valuations # not supported for ensemble solutions!
-@plot sol plot_type=new_transitions # not supported for ensemble solutions!
+@plot sol plot_type = summary
+@plot sol plot_type = allocation # not supported for ensemble solutions!
+@plot sol plot_type = valuations # not supported for ensemble solutions!
+@plot sol plot_type = new_transitions # not supported for ensemble solutions!
 ```
 """
 macro plot(solex, args...)
@@ -117,8 +148,12 @@ macro plot(solex, args...)
             plot_summary(sol, names, match_names($selector, names); $(kwargs...))
         else
             names = string.(sol.prob.p[:__state__][:, :specName])
-            plot_from_log(sol.prob.p[:__state__], plot_type, match_names($selector, names);
-                          $(kwargs...))
+            plot_from_log(
+                sol.prob.p[:__state__],
+                plot_type,
+                match_names($selector, names);
+                $(kwargs...),
+            )
         end
     end
 end
@@ -130,7 +165,7 @@ function get_transitions(log)
             union!(transitions, map(r -> r[1], r[3:end]))
     end
 
-    transitions
+    return transitions
 end
 
 function complete_log_transitions(log, hashes, record_type)
@@ -147,7 +182,7 @@ function complete_log_transitions(log, hashes, record_type)
         end
     end
 
-    vals, steps
+    return vals, steps
 end
 
 function complete_log_species(log, n_species, record_type)
@@ -161,7 +196,7 @@ function complete_log_species(log, n_species, record_type)
         end
     end
 
-    vals, steps
+    return vals, steps
 end
 
 function complete_log_valuations(log)
@@ -180,19 +215,21 @@ function complete_log_valuations(log)
         end
     end
 
-    vals, steps
+    return vals, steps
 end
 
 function get_names(state, hashes)
     names = String[]
     for hash in hashes
         ix = findfirst(==(hash), state[:, :transHash])
-        push!(names,
-              assigned(state.transition_recipes[:transName], ix) ?
-              string(state[ix, :transName]) : "transition $ix")
+        push!(names, if assigned(state.transition_recipes[:transName], ix)
+            string(state[ix, :transName])
+        else
+            "transition $ix"
+        end)
     end
 
-    names
+    return names
 end
 
 function plot_from_log(state, record_type, ixs; kwargs...)
@@ -210,8 +247,13 @@ function plot_from_log(state, record_type, ixs; kwargs...)
         vals, steps = complete_log_valuations(state.log)
     end
 
-    Plots.plot(steps, vals, title = string(record_type), label = reshape(label, 1, :);
-               kwargs...)
+    return Plots.plot(
+        steps,
+        vals;
+        title = string(record_type),
+        label = reshape(label, 1, :),
+        kwargs...,
+    )
 end
 
 """
@@ -221,16 +263,18 @@ Take an acset and optimize given functional.
 
 Objective is an expression which may reference the model's variables and parameters, i.e., `A+β`.
 The values to optimized are listed using their symbolic names; unless specified, the initial value is inferred from the model.
-The vector of free variables passed to the `NLopt` solver has the form `[free_vars; free_params]`; order of vars and params, respectively, is preserved. 
+The vector of free variables passed to the `NLopt` solver has the form `[free_vars; free_params]`; order of vars and params, respectively, is preserved.
 
-By default, the functional is minimized. Specify `objective=max` to perform maximization. 
+By default, the functional is minimized. Specify `objective=max` to perform maximization.
 
 Propagates `NLopt` solver arguments; see [NLopt documentation](https://github.com/JuliaOpt/NLopt.jl).
 
 # Examples
+
 ```julia
-@optimize acs abs(A-B) A B=20. α=2. lower_bounds=0 upper_bounds=100
-@optimize acss abs(A-B) A B=20. α=2. upper_bounds=[200,300,400] maxeval=200 objective=min
+@optimize acs abs(A - B) A B = 20.0 α = 2.0 lower_bounds = 0 upper_bounds = 100
+@optimize acss abs(A - B) A B = 20.0 α = 2.0 upper_bounds = [200, 300, 400] maxeval = 200 objective =
+    min
 ```
 """
 macro optimize(acsex, obex, args...)
@@ -254,9 +298,17 @@ macro optimize(acsex, obex, args...)
             ComponentVector{Float64}(; init_p...)
         end
 
-        o = build_loss_objective($(esc(acsex)), init_vec, u0, p, $(QuoteNode(obex));
-                                 min_t = $min_t, max_t = $max_t, final_only = $final_only,
-                                 $(okwargs...))
+        o = build_loss_objective(
+            $(esc(acsex)),
+            init_vec,
+            u0,
+            p,
+            $(QuoteNode(obex));
+            min_t = $min_t,
+            max_t = $max_t,
+            final_only = $final_only,
+            $(okwargs...),
+        )
 
         optim!(o, init_vec; $(kwargs...))
     end
@@ -268,15 +320,16 @@ end
 Take an acset and fit initial values and parameters to empirical data.
 
 The values to optimized are listed using their symbolic names; unless specified, the initial value is inferred from the model.
-The vector of free variables passed to the `NLopt` solver has the form `[free_vars; free_params]`; order of vars and params, respectively, is preserved. 
+The vector of free variables passed to the `NLopt` solver has the form `[free_vars; free_params]`; order of vars and params, respectively, is preserved.
 
 Propagates `NLopt` solver arguments; see [NLopt documentation](https://github.com/JuliaOpt/NLopt.jl).
 
 # Examples
+
 ```julia
 t = [1, 50, 100]
 data = [80 30 20]
-@fit acs data t vars=A B=20 A α # fit B, A, α; empirical data is for variable A
+@fit acs data t vars = A B = 20 A α # fit B, A, α; empirical data is for variable A
 ```
 """
 macro fit(acsex, data, t, args...)
@@ -284,12 +337,9 @@ macro fit(acsex, data, t, args...)
     args, kwargs = args_kwargs(args)
     okwargs = filter(ex -> ex.args[1] in [:loss, :trajectories], kwargs)
     vars = (ix = findfirst(ex -> ex.args[1] == :vars, kwargs);
-            !isnothing(ix) ?
-            (v = kwargs[ix].args[2];
-             deleteat!(kwargs,
-                       ix);
-             v) :
-            :())
+    !isnothing(ix) ? (v = kwargs[ix].args[2];
+    deleteat!(kwargs, ix);
+    v) : :())
 
     quote
         u0, p = get_free_vars($(esc(acsex)), $(QuoteNode(args_all)))
@@ -305,8 +355,16 @@ macro fit(acsex, data, t, args...)
             ComponentVector{Float64}(; init_p...)
         end
 
-        o = build_loss_objective_datapoints($(esc(acsex)), init_vec, u0, p, $(esc(t)),
-                                            $(esc(data)), vars; $(okwargs...))
+        o = build_loss_objective_datapoints(
+            $(esc(acsex)),
+            init_vec,
+            u0,
+            p,
+            $(esc(t)),
+            $(esc(data)),
+            vars;
+            $(okwargs...),
+        )
 
         optim!(o, init_vec; $(kwargs...))
     end
@@ -318,15 +376,16 @@ end
 Take an acset, fit initial values and parameters to empirical data, and plot the result.
 
 The values to optimized are listed using their symbolic names; unless specified, the initial value is inferred from the model.
-The vector of free variables passed to the `NLopt` solver has the form `[free_vars; free_params]`; order of vars and params, respectively, is preserved. 
+The vector of free variables passed to the `NLopt` solver has the form `[free_vars; free_params]`; order of vars and params, respectively, is preserved.
 
 Propagates `NLopt` solver arguments; see [NLopt documentation](https://github.com/JuliaOpt/NLopt.jl).
 
 # Examples
+
 ```julia
 t = [1, 50, 100]
 data = [80 30 20]
-@fit acs data t vars=A B=20 A α # fit B, A, α; empirical data is for variable A
+@fit acs data t vars = A B = 20 A α # fit B, A, α; empirical data is for variable A
 ```
 """
 macro fit_and_plot(acsex, data, t, args...)
@@ -335,12 +394,9 @@ macro fit_and_plot(acsex, data, t, args...)
     args, kwargs = args_kwargs(args)
     okwargs = filter(ex -> ex.args[1] in [:loss, :trajectories], kwargs)
     vars = (ix = findfirst(ex -> ex.args[1] == :vars, kwargs);
-            !isnothing(ix) ?
-            (v = kwargs[ix].args[2];
-             deleteat!(kwargs,
-                       ix);
-             v) :
-            :())
+    !isnothing(ix) ? (v = kwargs[ix].args[2];
+    deleteat!(kwargs, ix);
+    v) : :())
 
     quote
         u0, p = get_free_vars($(esc(acsex)), $(QuoteNode(args_all)))
@@ -356,26 +412,49 @@ macro fit_and_plot(acsex, data, t, args...)
             ComponentVector{Float64}(; init_p...)
         end
 
-        o = build_loss_objective_datapoints($(esc(acsex)), init_vec, u0, p, $(esc(t)),
-                                            $(esc(data)), vars; $(okwargs...))
+        o = build_loss_objective_datapoints(
+            $(esc(acsex)),
+            init_vec,
+            u0,
+            p,
+            $(esc(t)),
+            $(esc(data)),
+            vars;
+            $(okwargs...),
+        )
 
         r = optim!(o, init_vec; $(kwargs...))
         if r[3] != :FORCED_STOP
-            s_ = build_parametrized_solver($(esc(acsex)), init_vec, u0, p;
-                                           trajectories = $trajectories)
+            s_ = build_parametrized_solver(
+                $(esc(acsex)),
+                init_vec,
+                u0,
+                p;
+                trajectories = $trajectories,
+            )
             sol = first(s_(init_vec))
             sol_ = first(s_(r[2]))
 
-            p = Plots.plot(sol; idxs = vars,
-                           label = "(initial) " .*
-                                   reshape(String.($(esc(acsex))[:, :specName])[vars], 1,
-                                           :))
-            Plots.plot!(p, $(esc(t)), transpose($(esc(data))),
-                        label = "(empirical) " .*
-                                reshape(String.($(esc(acsex))[:, :specName])[vars], 1, :))
-            Plots.plot!(p, sol_; idxs = vars,
-                        label = "(fitted) " .*
-                                reshape(String.($(esc(acsex))[:, :specName])[vars], 1, :))
+            p = Plots.plot(
+                sol;
+                idxs = vars,
+                label = "(initial) " .*
+                        reshape(String.($(esc(acsex))[:, :specName])[vars], 1, :),
+            )
+            Plots.plot!(
+                p,
+                $(esc(t)),
+                transpose($(esc(data)));
+                label = "(empirical) " .*
+                        reshape(String.($(esc(acsex))[:, :specName])[vars], 1, :),
+            )
+            Plots.plot!(
+                p,
+                sol_;
+                idxs = vars,
+                label = "(fitted) " .*
+                        reshape(String.($(esc(acsex))[:, :specName])[vars], 1, :),
+            )
             p
         else
             :FORCED_STOP
@@ -389,6 +468,7 @@ end
 Take an acset and export a solution as a function of free vars and free parameters.
 
 # Examples
+
 ```julia
 solver = @build_solver acs S α β # function of variable S and parameters α, β
 solver([S, α, β])
@@ -411,7 +491,12 @@ macro build_solver(acsex, args...)
             ComponentVector{Float64}(; init_p...)
         end
 
-        build_parametrized_solver_($(esc(acsex)), init_vec, u0, p;
-                                   trajectories = $trajectories)
+        build_parametrized_solver_(
+            $(esc(acsex)),
+            init_vec,
+            u0,
+            p;
+            trajectories = $trajectories,
+        )
     end
 end
diff --git a/src/interface/update.jl b/src/interface/update.jl
index 7483128..326c3f6 100644
--- a/src/interface/update.jl
+++ b/src/interface/update.jl
@@ -16,8 +16,11 @@ function push_to_acs!(acsex, exs...)
         ex = striplines(exs[1])
     else
         args = Any[]
-        map(el -> isexpr(el, :(=)) ? push!(args, :($(el.args[1]) => $(el.args[2]))) :
-                  push!(args, el), exs)
+        map(el -> if isexpr(el, :(=))
+            push!(args, :($(el.args[1]) => $(el.args[2])))
+        else
+            push!(args, el)
+        end, exs)
         ex = Expr(:tuple, args...)
     end
 
@@ -31,38 +34,44 @@ end
 Add reactions to an acset.
 
 # Examples
+
 ```julia
-@push sir_acs β*S*I*tdecay(@time()) S+I --> 2I name=>SI2I
-@push sir_acs begin 
-    ν*I, I --> R, name=>I2R
-    γ, R --> S, name=>R2S
+@push sir_acs β * S * I * tdecay(@time()) S + I --> 2I name => SI2I
+@push sir_acs begin
+    ν * I, I --> R, name => I2R
+    γ, R --> S, name => R2S
 end
 ```
 """
 macro push(acsex, exs...)
-    push_to_acs!(acsex, exs...)
+    return push_to_acs!(acsex, exs...)
 end
 
 """
 Set name of a transition in the model.
 
 # Examples
+
 ```julia
-@name_transition acs 1="name"
-@name_transition acs name="transition_name"
-@name_transition acs "name"="transition_name"
+@name_transition acs 1 = "name"
+@name_transition acs name = "transition_name"
+@name_transition acs "name" = "transition_name"
 ```
 """
 macro name_transition(acsex, exs...)
-    call = :(begin end)
+    call = :(
+        begin end
+    )
     for ex in exs
         call_ = if ex.args[1] isa Number
             :($(esc(acsex))[$(ex.args[1]), :transName] = $(QuoteNode(ex.args[2])))
         else
             quote
                 acs = $(esc(acsex))
-                ixs = findall(i -> string(acs[i, :transName]) == $(string(ex.args[1])),
-                              1:nparts(acs, :T))
+                ixs = findall(
+                    i -> string(acs[i, :transName]) == $(string(ex.args[1])),
+                    1:nparts(acs, :T),
+                )
                 foreach(i -> acs[i, :transName] = $(string(ex.args[2])), ixs)
             end
         end
@@ -70,19 +79,19 @@ macro name_transition(acsex, exs...)
         push!(call.args, call_)
     end
 
-    call
+    return call
 end
 
 function incident_pattern(pattern, attr)
     ix = []
-    for i in 1:length(attr)
+    for i = 1:length(attr)
         !isnothing(attr[i]) &&
             (m = match(pattern, string(attr[i]));
-             !isnothing(m) &&
-                 (string(attr[i]) == m.match)) && push!(ix, i)
+            !isnothing(m) && (string(attr[i]) == m.match)) &&
+            push!(ix, i)
     end
 
-    ix
+    return ix
 end
 
 function mode!(acs, dict)
@@ -94,8 +103,7 @@ function mode!(acs, dict)
         end
 
         for ix in i
-            isnothing(acs[ix, specModality]) &&
-                (acs[ix, specModality] = Set{Symbol}())
+            isnothing(acs[ix, specModality]) && (acs[ix, specModality] = Set{Symbol}())
             union!(acs[ix, :specModality], mods)
         end
     end
@@ -105,11 +113,13 @@ end
 Set species modality.
 
 # Supported modalities
- - nonblock
- - conserved
- - rate
+
+  - nonblock
+  - conserved
+  - rate
 
 # Examples
+
 ```julia
 @mode acs (r"proj\\w+", r"experimental\\w+") conserved
 @mode acs (S, I) conserved
@@ -127,8 +137,10 @@ macro mode(acsex, spexs, mexs)
         foreach(s -> push!(exs_, striplines(blockize(s))), $(QuoteNode(exs)))
         exs__ = []
         foreach(s -> foreach(s -> push!(exs__, s), s.args), exs_)
-        foreach(ex -> push!(dictcall, get_pattern(recursively_expand_dots(ex)) => $mods),
-                exs__)
+        foreach(
+            ex -> push!(dictcall, get_pattern(recursively_expand_dots(ex)) => $mods),
+            exs__,
+        )
 
         mode!($(esc(acsex)), dictcall)
     end
@@ -142,48 +154,64 @@ function set_valuation!(acs, dict, valuation_type)
             incident(acs, Symbol(spex), :specName)
         end
 
-        foreach(ix -> acs[ix, Symbol(:spec, Symbol(uppercasefirst(string(valuation_type))))] = eval(val),
-                i)
+        foreach(
+            ix ->
+                acs[ix, Symbol(:spec, Symbol(uppercasefirst(string(valuation_type))))] =
+                    eval(val),
+            i,
+        )
     end
 end
 
 export @cost, @reward, @valuation
 for valuation_type in (:cost, :reward, :valuation)
-    eval(quote
-             export $(Symbol(Symbol("@"), valuation_type))
-             @doc """
-             Set $($(string(valuation_type))).
-
-             # Examples
-             ```julia
-             @$($(string(valuation_type))) model experimental1=2 experimental2=3
-             ```
-             """
-             macro $valuation_type(acsex, exs...)
-                 dictcall = Dict()
-                 exs_ = []
-                 foreach(s -> push!(exs_, striplines(blockize(s))), exs)
-                 exs__ = []
-                 foreach(s -> foreach(s -> push!(exs__, s), s.args), exs_)
-                 foreach(ex -> push!(dictcall,
-                                     get_pattern(recursively_expand_dots(ex.args[1])) => ex.args[2]),
-                         exs__)
-                 :(set_valuation!($(esc(acsex)), $dictcall,
-                                  $(QuoteNode($(QuoteNode(valuation_type))))))
-             end
-         end)
+    eval(
+        quote
+            export $(Symbol(Symbol("@"), valuation_type))
+            @doc """
+            Set $($(string(valuation_type))).
+
+            # Examples
+            ```julia
+            @$($(string(valuation_type))) model experimental1=2 experimental2=3
+            ```
+            """
+            macro $valuation_type(acsex, exs...)
+                dictcall = Dict()
+                exs_ = []
+                foreach(s -> push!(exs_, striplines(blockize(s))), exs)
+                exs__ = []
+                foreach(s -> foreach(s -> push!(exs__, s), s.args), exs_)
+                foreach(
+                    ex -> push!(
+                        dictcall,
+                        get_pattern(recursively_expand_dots(ex.args[1])) => ex.args[2],
+                    ),
+                    exs__,
+                )
+                return :(set_valuation!(
+                    $(esc(acsex)),
+                    $dictcall,
+                    $(QuoteNode($(QuoteNode(valuation_type)))),
+                ))
+            end
+        end,
+    )
 end
 
 """
 Add new species to a model.
 
 # Examples
+
 ```julia
 @add_species acs S I R
 ```
 """
 macro add_species(acsex, exs...)
-    call = :(begin end)
+    call = :(
+        begin end
+    )
     spexs_ = []
     foreach(s -> push!(spexs_, s), exs)
 
@@ -192,7 +220,7 @@ macro add_species(acsex, exs...)
     end
 
     push!(call.args, :(assign_defaults!($(esc(acsex)))))
-    call
+    return call
 end
 
 get_pattern(ex) = ex isa Expr && (macroname(ex) == :r_str) ? eval(ex) : ex
@@ -201,9 +229,10 @@ get_pattern(ex) = ex isa Expr && (macroname(ex) == :r_str) ? eval(ex) : ex
 Set initial values of species in an acset.
 
 # Examples
+
 ```julia
-@prob_init acs X=1 Y=2 Z=h(α)
-@prob_init acs [1., 2., 3.]
+@prob_init acs X = 1 Y = 2 Z = h(α)
+@prob_init acs [1.0, 2.0, 3.0]
 ```
 """
 macro prob_init(acsex, exs...)
@@ -218,9 +247,13 @@ macro prob_init(acsex, exs...)
             foreach(s -> push!(exs_, striplines(blockize(s))), $(QuoteNode(exs)))
             exs__ = []
             foreach(s -> foreach(s -> push!(exs__, s), s.args), exs_)
-            foreach(ex -> push!(dictcall,
-                                get_pattern(recursively_expand_dots(ex.args[1])) => ex.args[2]),
-                    exs__)
+            foreach(
+                ex -> push!(
+                    dictcall,
+                    get_pattern(recursively_expand_dots(ex.args[1])) => ex.args[2],
+                ),
+                exs__,
+            )
 
             init!($(esc(acsex)), dictcall)
         end
@@ -229,31 +262,39 @@ end
 
 # deprecate
 macro prob_init_from_vec(acsex, vecex)
-    :(init!($(esc(acsex)), $(esc(vecex))))
+    return :(init!($(esc(acsex)), $(esc(vecex))))
 end
 
 function init!(acs, inits)
-    inits isa AbstractVector && length(inits) == nparts(acs, :S) &&
+    inits isa AbstractVector &&
+        length(inits) == nparts(acs, :S) &&
         (subpart(acs, :specInitVal) .= inits; return)
     inits isa AbstractDict && for (k, init_val) in inits
-        k isa Number ? acs[k, :specInitVal] = init_val :
-        begin
-            i = k isa Regex ? incident_pattern(k, subpart(acs, :specName)) :
-                incident(acs, k, :specName)
-            foreach(ix -> (acs[ix, :specInitVal] = init_val), i)
+        if k isa Number
+            acs[k, :specInitVal] = init_val
+        else
+            begin
+                i = if k isa Regex
+                    incident_pattern(k, subpart(acs, :specName))
+                else
+                    incident(acs, k, :specName)
+                end
+                foreach(ix -> (acs[ix, :specInitVal] = init_val), i)
+            end
         end
     end
 
-    acs
+    return acs
 end
 
 """
 Set uncertainty in initial values of species in an acset (stderr).
 
 # Examples
+
 ```julia
-@prob_uncertainty acs X=.1 Y=.2
-@prob_uncertainty acs [.1, .2,]
+@prob_uncertainty acs X = 0.1 Y = 0.2
+@prob_uncertainty acs [0.1, 0.2]
 ```
 """
 macro prob_uncertainty(acsex, exs...)
@@ -268,9 +309,13 @@ macro prob_uncertainty(acsex, exs...)
             foreach(s -> push!(exs_, striplines(blockize(s))), $(QuoteNode(exs)))
             exs__ = []
             foreach(s -> foreach(s -> push!(exs__, s), s.args), exs_)
-            foreach(ex -> push!(dictcall,
-                                get_pattern(recursively_expand_dots(ex.args[1])) => ex.args[2]),
-                    exs__)
+            foreach(
+                ex -> push!(
+                    dictcall,
+                    get_pattern(recursively_expand_dots(ex.args[1])) => ex.args[2],
+                ),
+                exs__,
+            )
 
             uncinit!($(esc(acsex)), dictcall)
         end
@@ -278,28 +323,35 @@ macro prob_uncertainty(acsex, exs...)
 end
 
 function uncinit!(acs, inits)
-    inits isa AbstractVector && length(inits) == nparts(acs, :S) &&
+    inits isa AbstractVector &&
+        length(inits) == nparts(acs, :S) &&
         (subpart(acs, :specInitUncertainty) .= inits; return)
     inits isa AbstractDict && for (k, init_val) in inits
-        k isa Number ? acs[k, :specInitUncertainty] = init_val :
-        begin
-            i = k isa Regex ? incident_pattern(k, subpart(acs, :specName)) :
-                incident(acs, k, :specName)
-            foreach(ix -> (acs[ix, :specInitUncertainty] = init_val), i)
+        if k isa Number
+            acs[k, :specInitUncertainty] = init_val
+        else
+            begin
+                i = if k isa Regex
+                    incident_pattern(k, subpart(acs, :specName))
+                else
+                    incident(acs, k, :specName)
+                end
+                foreach(ix -> (acs[ix, :specInitUncertainty] = init_val), i)
+            end
         end
     end
 
-    acs
+    return acs
 end
 
 function set_params!(acs, params)
-    params isa AbstractDict && for (k, init_val) in params
+    return params isa AbstractDict && for (k, init_val) in params
         k = get_pattern(k)
         if k isa Regex
             i = incident_pattern(k, subpart(acs, :prmName))
         else
             i = incident(acs, k, :prmName)
-            isempty(i) && (i = add_part!(acs, :P, prmName = k))
+            isempty(i) && (i = add_part!(acs, :P; prmName = k))
         end
 
         foreach(ix -> acs[ix, :prmVal] = eval(init_val), i)
@@ -310,8 +362,9 @@ end
 Set parameter values in an acset.
 
 # Examples
+
 ```julia
-@prob_params acs α=1. β=2.
+@prob_params acs α = 1.0 β = 2.0
 ```
 """
 macro prob_params(acsex, exs...)
@@ -323,117 +376,133 @@ macro prob_params(acsex, exs...)
         foreach(s -> push!(exs_, striplines(blockize(s))), $(QuoteNode(exs)))
         exs__ = []
         foreach(s -> foreach(s -> push!(exs__, s), s.args), exs_)
-        foreach(ex -> push!(dictcall,
-                            get_pattern(recursively_expand_dots(ex.args[1])) => ex.args[2]),
-                exs__)
+        foreach(
+            ex -> push!(
+                dictcall,
+                get_pattern(recursively_expand_dots(ex.args[1])) => ex.args[2],
+            ),
+            exs__,
+        )
 
         set_params!($(esc(acsex)), dictcall)
     end
 end
 
-function meta!(acs, metas)
+meta!(acs, metas) =
     for (k, metaval) in metas
         i = incident(acs, k, :metaKeyword)
-        isempty(i) && (i = add_part!(acs, :M, metaKeyword = k))
+        isempty(i) && (i = add_part!(acs, :M; metaKeyword = k))
         set_subpart!(acs, first(i), :metaVal, metaval)
     end
-end
 
 """
 Set model metadata (e.g. solver arguments)
 
 # Examples
+
 ```julia
-@prob_meta acs tspan=(0, 100.) schedule=schedule_weighted!
-@prob_meta sir_acs tspan=250 tstep=1
+@prob_meta acs tspan = (0, 100.0) schedule = schedule_weighted!
+@prob_meta sir_acs tspan = 250 tstep = 1
 ```
 """
 macro prob_meta(acsex, exs...)
     dictcall = :(Dict([]))
     foreach(ex -> push!(dictcall.args[2].args, (ex.args[1] => eval(ex.args[2]))), exs)
 
-    :(meta!($(esc(acsex)), $dictcall))
+    return :(meta!($(esc(acsex)), $dictcall))
 end
 
 """
 Alias object name in an acs.
 
 # Default names
-| name | short name |
-| :--- | :--- |
-| species | S |
-| transition | T |
-| action | A |
-| event | E |
-| param | P |
-| meta | M |
+
+| name       | short name |
+|:---------- |:---------- |
+| species    | S          |
+| transition | T          |
+| action     | A          |
+| event      | E          |
+| param      | P          |
+| meta       | M          |
 
 # Examples
+
 ```julia
-@aka acs species=resource transition=reaction
+@aka acs species = resource transition = reaction
 ```
 """
 macro aka(acsex, exs...)
     dictcall = :(Dict([]))
-    foreach(ex -> push!(dictcall.args[2].args,
-                        Symbol("alias_", findfirst(==(ex.args[1]), alias_default)) => ex.args[2]),
-            exs)
-    :(meta!($(esc(acsex)), $dictcall))
+    foreach(
+        ex -> push!(
+            dictcall.args[2].args,
+            Symbol("alias_", findfirst(==(ex.args[1]), alias_default)) => ex.args[2],
+        ),
+        exs,
+    )
+    return :(meta!($(esc(acsex)), $dictcall))
 end
 
-alias_default = Dict(:S => :species, :T => :transition, :A => :action, :E => :event,
-                     :P => :param, :M => :meta)
+alias_default = Dict(
+    :S => :species,
+    :T => :transition,
+    :A => :action,
+    :E => :event,
+    :P => :param,
+    :M => :meta,
+)
 
-function get_alias(acs, ob)
-    (i = incident(acs, Symbol(:alias_, ob), :metaKeyword);
-     !isempty(i) ?
-     acs[first(i), :metaVal] :
-     alias_default[ob])
-end
+get_alias(acs, ob) = (i = incident(acs, Symbol(:alias_, ob), :metaKeyword);
+!isempty(i) ? acs[first(i), :metaVal] : alias_default[ob])
 
-"Check model parameters have been set."
+"""
+Check model parameters have been set. # msg as return value
+"""
 macro prob_check_verbose(acsex) # msg as return value
-    :(missing_params = check_params($(esc(acsex)));
-      isempty(missing_params) ? "Params OK." :
-      "Missing params: $missing_params")
+    return :(missing_params = check_params($(esc(acsex)));
+    isempty(missing_params) ? "Params OK." : "Missing params: $missing_params")
 end
 
 """
 Add a periodic callback to a model.
 
 # Examples
+
 ```julia
-@periodic acs 1. X += 1
+@periodic acs 1.0 X += 1
 ```
 """
 macro periodic(acsex, pex, acex)
-    push_to_acs!(acsex, Expr(:&&, :(@periodic($pex)), acex))
+    return push_to_acs!(acsex, Expr(:&&, :(@periodic($pex)), acex))
 end
 
 """
 Add a jump process (with specified Poisson intensity per unit time step) to a model.
 
 # Examples
+
 ```julia
-@jump acs λ Z += rand(Poisson(1.))
+@jump acs λ Z += rand(Poisson(1.0))
 ```
 """
 macro jump(acsex, inex, acex)
-    push_to_acs!(acsex,
-                 Expr(:&&,
-                      Expr(:call, :rand, :(Poisson(max(@solverarg(:tstep) * $inex, 0)))),
-                      acex))
+    return push_to_acs!(
+        acsex,
+        Expr(:&&, Expr(:call, :rand, :(Poisson(max(@solverarg(:tstep) * $inex, 0)))), acex),
+    )
 end
 
 """
 Evaluate expression in ReactiveDynamics scope.
 
 # Examples
+
 ```julia
 @register bool_cond(t) = (100 < t < 200) || (400 < t < 500)
-@register tdecay(t) = exp(-t/10^3)
+@register tdecay(t) = exp(-t / 10^3)
 ```
 """
 macro register(ex)
-    :(@eval ReactiveDynamics $ex)
+    return :(@eval ReactiveDynamics $ex)
 end
diff --git a/src/loadsave.jl b/src/loadsave.jl
index 8263166..da3d8b4 100644
--- a/src/loadsave.jl
+++ b/src/loadsave.jl
@@ -7,22 +7,28 @@ export @export, @import
 using TOML, JLD2, CSV
 using DataFrames
 
-const objects_aliases = Dict(:S => "spec", :T => "trans", :P => "prm", :M => "meta",
-                             :E => "event", :obs => "obs")
+const objects_aliases = Dict(
+    :S => "spec",
+    :T => "trans",
+    :P => "prm",
+    :M => "meta",
+    :E => "event",
+    :obs => "obs",
+)
 
 const RN_attrs = string.(propertynames(ReactionNetwork().subparts))
 
 function get_attrs(object)
     object = object isa Symbol ? objects_aliases[object] : object
 
-    filter(x -> occursin(object, x), RN_attrs)
+    return filter(x -> occursin(object, x), RN_attrs)
 end
 
 function export_network(acs::ReactionNetwork)
     dict = Dict()
     for (key, val) in objects_aliases
         push!(dict, val => [])
-        for i in 1:nparts(acs, key)
+        for i = 1:nparts(acs, key)
             dict_ = Dict()
             for attr in get_attrs(val)
                 attr_val = acs[i, Symbol(attr)]
@@ -34,7 +40,7 @@ function export_network(acs::ReactionNetwork)
         end
     end
 
-    dict
+    return dict
 end
 
 function load_network(dict::Dict)
@@ -63,7 +69,7 @@ function load_network(dict::Dict)
         eval(Meta.parseall(row["body"]))
     end
 
-    assign_defaults!(acs)
+    return assign_defaults!(acs)
 end
 
 function import_network_csv(pathmap)
@@ -71,8 +77,15 @@ function import_network_csv(pathmap)
     for (key, paths) in pathmap
         push!(dict, key => [])
         for path in paths
-            data = DataFrame(CSV.File(path; delim = ";;", types = String,
-                                      stripwhitespace = true, comment = "#"))
+            data = DataFrame(
+                CSV.File(
+                    path;
+                    delim = ";;",
+                    types = String,
+                    stripwhitespace = true,
+                    comment = "#",
+                ),
+            )
             for row in eachrow(data)
                 object = Dict()
                 for (attr, val) in Iterators.zip(keys(row), values(row))
@@ -83,13 +96,13 @@ function import_network_csv(pathmap)
         end
     end
 
-    load_network(dict)
+    return load_network(dict)
 end
 
 function import_network(path::AbstractString)
     if splitext(path)[2] == ".csv"
-        pathmap = Dict(val => []
-                       for val in [collect(values(objects_aliases)); "registered"])
+        pathmap =
+            Dict(val => [] for val in [collect(values(objects_aliases)); "registered"])
         for row in CSV.File(path; delim = ";;", stripwhitespace = true, comment = "#")
             push!(pathmap[row.type], joinpath(dirname(path), row.path))
         end
@@ -103,19 +116,24 @@ end
 function export_network(acs::ReactionNetwork, path::AbstractString)
     if splitext(path)[2] == ".csv"
         exported_network = export_network(acs)
-        paths = DataFrame(type = [], path = [])
+        paths = DataFrame(; type = [], path = [])
         for (key, objs) in exported_network
             push!(paths, (key, "export-$key.csv"))
             objs_exported = DataFrame(Dict(attr => [] for attr in get_attrs(key)))
             for obj in objs
-                push!(objs_exported,
-                      [get(obj, key, missing) for key in names(objs_exported)])
+                push!(
+                    objs_exported,
+                    [get(obj, key, missing) for key in names(objs_exported)],
+                )
             end
 
-            CSV.write(joinpath(dirname(path), "export-$key.csv"), objs_exported,
-                      delim = ";;")
+            CSV.write(
+                joinpath(dirname(path), "export-$key.csv"),
+                objs_exported;
+                delim = ";;",
+            )
         end
-        CSV.write(path, paths, delim = ";;")
+        CSV.write(path, paths; delim = ";;")
     else
         open(io -> TOML.print(io, export_network(acs)), path, "w+")
     end
@@ -128,13 +146,14 @@ or a batch of CSV files (a root file and a number of files, each per a class of
 See `tutorials/loadsave` for an example.
 
 # Examples
+
 ```julia
 @export_network acs "acs_data.toml" # as a TOML
 @export_network acs "csv/model.csv" # as a CSV
 ```
 """
 macro export_network(acsex, pathex)
-    :(export_network($(esc(acsex)), $(string(pathex))))
+    return :(export_network($(esc(acsex)), $(string(pathex))))
 end
 
 """
@@ -144,65 +163,74 @@ or a batch of CSV files (a root file and a number of files, each per a class of
 See `tutorials/loadsave` for an example.
 
 # Examples
+
 ```julia
 @import_network "model.toml"
 @import_network "csv/model.toml"
 ```
 """
 macro import_network(pathex, name = gensym())
-    :($(esc(name)) = import_network($(string(pathex))))
+    return :($(esc(name)) = import_network($(string(pathex))))
 end
 
 macro load_models(pathex)
-    callex = :(begin end)
+    callex = :(
+        begin end
+    )
     for line in readlines(string(pathex))
         name, pathex = split(line, ';')
         name = isempty(name) ? gensym() : Symbol(name)
         push!(callex.args, :($(esc(name)) = import_network($(string(pathex)))))
     end
 
-    callex
+    return callex
 end
 
 """
     @import_solution "sol.jld2"
     @import_solution "sol.jld2" sol
+
 Import a solution from a file.
 
 # Examples
+
 ```julia
 @import_solution "sir_acs_sol/serialized/sol.jld2"
 ```
 """
 macro import_solution(pathex, varname = "sol")
-    :(JLD2.load($(string(pathex)), $(string(varname))))
+    return :(JLD2.load($(string(pathex)), $(string(varname))))
 end
 
 """
     @export_solution sol
     @export_solution sol "sol.jld2"
+
 Export a solution to a file.
 
 # Examples
+
 ```julia
 @export_solution sol "sol.jdl2"
 ```
 """
 macro export_solution(solex, pathex = "sol.jld2")
-    :(JLD2.save($(string(pathex)), $(string(solex)), $(esc(solex))))
+    return :(JLD2.save($(string(pathex)), $(string(solex)), $(esc(solex))))
 end
 
 """
     @export_solution_as_table sol
+
 Export a solution as a `DataFrame`.
 
 # Examples
+
 ```julia
 @export_solution_as_table sol
 ```
 """
 macro export_solution_as_table(solex, pathex = "sol.jld2")
-    :(DataFrame($(esc(solex))))
+    return :(DataFrame($(esc(solex))))
 end
 
 get_DataFrame(sol) = sol isa EnsembleSolution ? DataFrame(sol)[!, [:u, :t]] : DataFrame(sol)
@@ -210,13 +238,15 @@ get_DataFrame(sol) = sol isa EnsembleSolution ? DataFrame(sol)[!, [:u, :t]] : Da
 """
     @export_solution_as_csv sol
     @export_solution_as_csv sol "sol.csv"
+
 Export a solution to a file.
 
 # Examples
+
 ```julia
 @export_solution_as_csv sol "sol.csv"
 ```
 """
 macro export_solution_as_csv(solex, pathex = "sol.csv")
-    :(CSV.write($(string(pathex)), get_DataFrame($(esc(solex)))))
+    return :(CSV.write($(string(pathex)), get_DataFrame($(esc(solex)))))
 end
diff --git a/src/operators/equalize.jl b/src/operators/equalize.jl
index 1f98c9d..cbe357d 100644
--- a/src/operators/equalize.jl
+++ b/src/operators/equalize.jl
@@ -1,15 +1,21 @@
 export @equalize
 
-function expand_name_ff(ex)
-    ex isa Expr && isexpr(ex, :macrocall) ? (macroname(ex), underscorize(ex.args[end])) :
-    (nothing, underscorize(ex))
-end
+expand_name_ff(ex) =
+    if ex isa Expr && isexpr(ex, :macrocall)
+        (macroname(ex), underscorize(ex.args[end]))
+    else
+        (nothing, underscorize(ex))
+    end
 
-"Parse species equation blocks."
+"""
+Parse species equation blocks.
+"""
 function get_eqs_ff(eq)
     if eq isa Expr && isexpr(eq, :(=))
-        [expand_name_ff(eq.args[1]);
-         isexpr(eq.args[2], :(=)) ? get_eqs_ff(eq.args[2]) : expand_name_ff(eq.args[2])]
+        [
+            expand_name_ff(eq.args[1])
+            isexpr(eq.args[2], :(=)) ? get_eqs_ff(eq.args[2]) : expand_name_ff(eq.args[2])
+        ]
     else
         [expand_name_ff(eq)]
     end
@@ -21,14 +27,16 @@ function equalize!(acs::ReactionNetwork, eqs = [])
         block_alias = findfirst(e -> e[1] == :alias, block)
         block_alias = !isnothing(block_alias) ? block[block_alias][2] : first(block)[2]
         species_ixs = Int64[]
-        for e in block, i in 1:nparts(acs, :S)
-            ((i == e[2]) ||
-             (e[1] == :catchall &&
-              occursin(Regex("(__$(e[2])|$(e[2]))\$"), string(acs[i, :specName]))) ||
-             (e[2] == acs[i, :specName])) &&
-                (push!(species_ixs, i);
-                 push!(specmap,
-                       acs[i, :specName] => (acs[i, :specName] = block_alias)))
+        for e in block, i = 1:nparts(acs, :S)
+            (
+                (i == e[2]) ||
+                (
+                    e[1] == :catchall &&
+                    occursin(Regex("(__$(e[2])|$(e[2]))\$"), string(acs[i, :specName]))
+                ) ||
+                (e[2] == acs[i, :specName])
+            ) && (push!(species_ixs, i);
+            push!(specmap, acs[i, :specName] => (acs[i, :specName] = block_alias)))
         end
         isempty(species_ixs) && continue
         species_ixs = sort(unique!(species_ixs))
@@ -45,21 +53,22 @@ function equalize!(acs::ReactionNetwork, eqs = [])
     for attr in propertynames(acs.subparts)
         attr == :specName && continue
         attr_ = acs[:, attr]
-        for i in 1:length(attr_)
+        for i = 1:length(attr_)
             attr_[i] = escape_ref(attr_[i], collect(keys(specmap)))
             attr_[i] = recursively_substitute_vars!(specmap, attr_[i])
         end
     end
 
-    acs
+    return acs
 end
 
 """
 Identify (collapse) a set of species in a model.
 
 # Examples
+
 ```julia
-@join acs acs1.A=acs2.A B=C
+@join acs acs1.A = acs2.A B = C
 ```
 """
 macro equalize(acsex, exs...)
@@ -68,5 +77,5 @@ macro equalize(acsex, exs...)
     eqs = []
     foreach(ex -> ex isa Expr && merge_eqs!(eqs, get_eqs_ff(ex)), exs)
 
-    :(equalize!($(esc(acsex)), $eqs))
+    return :(equalize!($(esc(acsex)), $eqs))
 end
diff --git a/src/operators/joins.jl b/src/operators/joins.jl
index ecd3f3b..40ebc73 100644
--- a/src/operators/joins.jl
+++ b/src/operators/joins.jl
@@ -4,14 +4,16 @@ export @join
 using MacroTools
 using MacroTools: prewalk
 
-"Merge `acs2` onto `acs1`, the attributes in `acs2` taking precedence. Identify respective species given `eqs`, renaming species in `acs2`."
+"""
+Merge `acs2` onto `acs1`, the attributes in `acs2` taking precedence. Identify respective species given `eqs`, renaming species in `acs2`.
+"""
 function union_acs!(acs1, acs2, name = gensym("acs_"), eqs = [])
     acs2 = deepcopy(acs2)
     prepend!(acs2, name, eqs)
-    for i in 1:nparts(acs2, :S)
+    for i = 1:nparts(acs2, :S)
         inc = incident(acs1, acs2[i, :specName], :specName)
-        isempty(inc) &&
-            (inc = add_part!(acs1, :S; specName = acs2[i, :specName]); assign_defaults!(acs1))
+        isempty(inc) && (inc = add_part!(acs1, :S; specName = acs2[i, :specName]);
+        assign_defaults!(acs1))
         return (acs1, acs2)
         println(first(inc))
         println(acs1[first(inc), :specModality])
@@ -31,29 +33,35 @@ function union_acs!(acs1, acs2, name = gensym("acs_"), eqs = [])
         acs1[new_trans_ix, attr] .= acs2[:, attr]
     end
 
-    foreach(i -> (acs1[i, :transName] = normalize_name(Symbol(coalesce(acs1[i, :transName],
-                                                                       i)), name)),
-            new_trans_ix)
+    foreach(
+        i -> (
+            acs1[i, :transName] =
+                normalize_name(Symbol(coalesce(acs1[i, :transName], i)), name)
+        ),
+        new_trans_ix,
+    )
 
-    for i in 1:nparts(acs2, :P)
+    for i = 1:nparts(acs2, :P)
         inc = incident(acs1, acs2[i, :prmName], :prmName)
         isempty(inc) && (inc = add_part!(acs1, :P; prmName = acs2[i, :prmName]))
         !ismissing(acs2[i, :prmVal]) && (acs1[first(inc), :prmVal] = acs2[i, :prmVal])
     end
 
-    for i in 1:nparts(acs2, :M)
+    for i = 1:nparts(acs2, :M)
         inc = incident(acs1, acs2[i, :metaKeyword], :metaKeyword)
         isempty(inc) && (inc = add_part!(acs1, :M; metaKeyword = acs2[i, :metaKeyword]))
         !ismissing(acs2[i, :metaVal]) && (acs1[first(inc), :metaVal] = acs2[i, :metaVal])
     end
 
-    acs1
+    return acs1
 end
 
-"Prepend species names with a model identifier (unless a global species name)."
+"""
+Prepend species names with a model identifier (unless a global species name).
+"""
 function prepend!(acs::ReactionNetwork, name = gensym("acs"), eqs = [])
     specmap = Dict()
-    for i in 1:nparts(acs, :S)
+    for i = 1:nparts(acs, :S)
         new_name = normalize_name(name, i, acs[i, :specName], eqs)
         push!(specmap, acs[i, :specName] => (acs[i, :specName] = new_name))
     end
@@ -61,22 +69,27 @@ function prepend!(acs::ReactionNetwork, name = gensym("acs"), eqs = [])
     for attr in propertynames(acs.subparts)
         attr == :specName && continue
         attr_ = acs[:, attr]
-        for i in 1:length(attr_)
+        for i = 1:length(attr_)
             attr_[i] = escape_ref(attr_[i], collect(keys(specmap)))
             attr_[i] = recursively_substitute_vars!(specmap, attr_[i])
             attr_[i] isa Expr && (attr_[i] = prepend_obs(attr_[i], name))
         end
     end
 
-    acs
+    return acs
 end
 
-"Prepend identifier of an observable with a model identifier."
+"""
+Prepend identifier of an observable with a model identifier.
+"""
 function prepend_obs end
 
 function prepend_obs(ex::Expr, name)
-    prewalk(ex -> (isexpr(ex, :macrocall) && (macroname(ex) == :obs)) ?
-                  (ex.args[end] = Symbol(name, "__", ex.args[end]); ex) : ex, ex)
+    return prewalk(ex -> if (isexpr(ex, :macrocall) && (macroname(ex) == :obs))
+        (ex.args[end] = Symbol(name, "__", ex.args[end]); ex)
+    else
+        ex
+    end, ex)
 end
 
 prepend_obs(ex, _) = ex
@@ -89,46 +102,68 @@ normalize_name(name, parent_name) = Symbol(parent_name, "__", name)
 function normalize_name(acs_name, i::Int, name::Symbol, eqs = [])
     for (block_ix, block) in enumerate(eqs)
         block_alias = findfirst(e -> e[1] == :alias, block)
-        block_alias = !isnothing(block_alias) ? block[block_alias][2] :
-                      Symbol(:shared_species_, block_ix)
+        block_alias = if !isnothing(block_alias)
+            block[block_alias][2]
+        else
+            Symbol(:shared_species_, block_ix)
+        end
         for e in block
-            ((i == e[2]) ||
-             (e[1] == :catchall &&
-              (normalize_name(e[2], acs_name) == normalize_name(name, acs_name))) ||
-             (e[1] == acs_name &&
-              (normalize_name(e[2], acs_name) == normalize_name(name, acs_name)))) &&
-                return block_alias
+            (
+                (i == e[2]) ||
+                (
+                    e[1] == :catchall &&
+                    (normalize_name(e[2], acs_name) == normalize_name(name, acs_name))
+                ) ||
+                (
+                    e[1] == acs_name &&
+                    (normalize_name(e[2], acs_name) == normalize_name(name, acs_name))
+                )
+            ) && return block_alias
         end
     end
 
-    normalize_name(name, acs_name)
+    return normalize_name(name, acs_name)
 end
 
 matching_name(name::Symbol, parent_name) = [name, Symbol("$(parent_name)__$name")]
 
-function expand_name(ex)
-    isexpr(ex, :.) ? reconstruct(ex) :
-    isexpr(ex, :macrocall) ?
-    (macroname(ex) == :alias ? [(:alias, ex.args[3])] :
-     [(:catchall, ex.args[3]), (:alias, ex.args[3])]) :
-    (:catchall, ex)
-end
+expand_name(ex) =
+    if isexpr(ex, :.)
+        reconstruct(ex)
+    elseif isexpr(ex, :macrocall)
+        (
+            if macroname(ex) == :alias
+                [(:alias, ex.args[3])]
+            else
+                [(:catchall, ex.args[3]), (:alias, ex.args[3])]
+            end
+        )
+    else
+        (:catchall, ex)
+    end
 
 function recursively_get_syms(ex)
-    isexpr(ex, :.) ? [recursively_get_syms(ex.args[1]); ex.args[2].value] : ex
+    return isexpr(ex, :.) ? [recursively_get_syms(ex.args[1]); ex.args[2].value] : ex
 end
 function reconstruct(ex)
-    ex isa Symbol ? (ex,) :
-    (syms = recursively_get_syms(ex); (syms[1], Symbol(join(syms[2:end], "__"))))
+    return if ex isa Symbol
+        (ex,)
+    else
+        (syms = recursively_get_syms(ex); (syms[1], Symbol(join(syms[2:end], "__"))))
+    end
 end
 
-"Parse species equation blocks."
+"""
+Parse species equation blocks.
+"""
 function get_eqs(eq)
     if isexpr(eq, :macrocall)
         expand_name(eq)
     elseif eq isa Expr
-        [expand_name(eq.args[1]);
-         isexpr(eq.args[2], :(=)) ? get_eqs(eq.args[2]) : expand_name(eq.args[2])]
+        [
+            expand_name(eq.args[1])
+            isexpr(eq.args[2], :(=)) ? get_eqs(eq.args[2]) : expand_name(eq.args[2])
+        ]
     else
         [eq]
     end
@@ -145,7 +180,7 @@ function merge_eqs!(eqs, eqblock)
     foreach(i -> deleteat!(eqs, i), Iterators.reverse(sort!(eqs_)))
 
     push!(eqs, eqblock)
-    eqs_
+    return eqs_
 end
 
 """
@@ -156,12 +191,17 @@ Performs join of models and identifies model variables, as specified.
 Model variables / parameter values and metadata are propagated; the last model takes precedence.
 
 # Examples
+
 ```julia
-@join acs1 acs2 @catchall(A)=acs2.Z @catchall(XY) @catchall(B)
+@join acs1 acs2 @catchall(A) = acs2.Z @catchall(XY) @catchall(B)
 ```
 """
 macro join(exs...)
-    callex = :(begin acs_new = ReactionNetwork() end)
+    callex = :(
+        begin
+            acs_new = ReactionNetwork()
+        end
+    )
     exs = collect(exs)
     foreach(i -> (exs[i] = MacroTools.striplines(exs[i])), 1:length(exs))
     eqs = []
@@ -177,8 +217,9 @@ macro join(exs...)
 
     for acsex in exs
         (acsex, symex) = if isexpr(acsex, :macrocall)
-            str_inc = string(isexpr(acsex.args[3], :(=)) ? acsex.args[3].args[2] :
-                             acsex.args[3])
+            str_inc = string(
+                isexpr(acsex.args[3], :(=)) ? acsex.args[3].args[2] : acsex.args[3],
+            )
             if isexpr(acsex.args[3], :(=))
                 (:(include_model($str_inc)), acsex.args[3].args[1])
             else
@@ -191,5 +232,5 @@ macro join(exs...)
     end
     push!(callex.args, :(acs_new))
 
-    callex
+    return callex
 end
diff --git a/src/optim.jl b/src/optim.jl
index 714b74f..a30e179 100644
--- a/src/optim.jl
+++ b/src/optim.jl
@@ -6,12 +6,15 @@ function build_parametrized_solver(acs, init_vec, u0, params; trajectories = 1)
     function (vec)
         vec = vec isa ComponentVector ? vec : (init_vec .= vec)
         data = []
-        for _ in 1:trajectories
+        for _ = 1:trajectories
             prob.p[:__state__] = deepcopy(prob.p[:__state0__])
             for i in eachindex(prob.u0)
                 rv = randn() * vars[i]
-                prob.u0[i] = (sign(rv + prob.u0[i]) == sign(prob.u0[i])) ? rv + prob.u0[i] :
-                             prob.u0[i]
+                prob.u0[i] = if (sign(rv + prob.u0[i]) == sign(prob.u0[i]))
+                    rv + prob.u0[i]
+                else
+                    prob.u0[i]
+                end
             end
 
             for (i, k) in enumerate(wkeys(u0))
@@ -25,7 +28,7 @@ function build_parametrized_solver(acs, init_vec, u0, params; trajectories = 1)
             push!(data, solve(prob))
         end
 
-        data
+        return data
     end
 end
 
@@ -40,8 +43,11 @@ function build_parametrized_solver_(acs, init_vec, u0, params; trajectories = 1)
             prob.p[:__state__] = deepcopy(prob.p[:__state0__])
             for i in eachindex(prob.u0)
                 rv = randn() * vars[i]
-                prob.u0[i] = (sign(rv + prob.u0[i]) == sign(prob.u0[i])) ? rv + prob.u0[i] :
-                             prob.u0[i]
+                prob.u0[i] = if (sign(rv + prob.u0[i]) == sign(prob.u0[i]))
+                    rv + prob.u0[i]
+                else
+                    prob.u0[i]
+                end
             end
 
             for (i, k) in enumerate(wkeys(u0))
@@ -53,10 +59,10 @@ function build_parametrized_solver_(acs, init_vec, u0, params; trajectories = 1)
 
             sync!(prob.p[:__state__], prob.u0, prob.p)
 
-            solve(prob)
+            return solve(prob)
         end
 
-        trajectories == 1 ? data[1] : EnsembleSolution(data, 0.0, true)
+        return trajectories == 1 ? data[1] : EnsembleSolution(data, 0.0, true)
     end
 end
 
@@ -71,20 +77,34 @@ function optim!(obj, init; nlopt_kwargs...)
     opt = Opt(alg, length(init))
 
     # match to a ComponentVector
-    foreach(o -> setproperty!(opt, o...),
-            filter(x -> x[1] in propertynames(opt), nlopt_kwargs))
-    get(nlopt_kwargs, :objective, min) == min ? (opt.min_objective = obj) :
-    (opt.max_objective = obj)
+    foreach(
+        o -> setproperty!(opt, o...),
+        filter(x -> x[1] in propertynames(opt), nlopt_kwargs),
+    )
+    if get(nlopt_kwargs, :objective, min) == min
+        (opt.min_objective = obj)
+    else
+        (opt.max_objective = obj)
+    end
 
-    optimize(opt, deepcopy(init))
+    return optimize(opt, deepcopy(init))
 end
 
 const n_steps = 100
 
 # loss objective given an objective expression
-function build_loss_objective(acs, init_vec, u0, params, obex; loss = identity,
-                              trajectories = 1, min_t = -Inf, max_t = Inf,
-                              final_only = false)
+function build_loss_objective(
+    acs,
+    init_vec,
+    u0,
+    params,
+    obex;
+    loss = identity,
+    trajectories = 1,
+    min_t = -Inf,
+    max_t = Inf,
+    final_only = false,
+)
     ob = eval(get_wrap_fun(acs)(obex))
     obj_ = build_parametrized_solver(acs, init_vec, u0, params; trajectories)
 
@@ -100,29 +120,44 @@ function build_loss_objective(acs, init_vec, u0, params, obex; loss = identity,
                 range(min_t, max_t; length = n_steps)
             end
 
-            push!(ls,
-                  mean(t -> loss(ob(as_state(sol(t), t, sol.prob.p[:__state__]))),
-                       t_points))
+            push!(
+                ls,
+                mean(t -> loss(ob(as_state(sol(t), t, sol.prob.p[:__state__]))), t_points),
+            )
         end
 
-        mean(ls)
+        return mean(ls)
     end
 end
 
 # loss objective given empirical data
-function build_loss_objective_datapoints(acs, init_vec, u0, params, t, data, vars;
-                                         loss = abs2, trajectories = 1)
+function build_loss_objective_datapoints(
+    acs,
+    init_vec,
+    u0,
+    params,
+    t,
+    data,
+    vars;
+    loss = abs2,
+    trajectories = 1,
+)
     obj_ = build_parametrized_solver(acs, init_vec, u0, params; trajectories)
 
     function (vec, _)
         ls = []
         for sol in obj_(vec)
-            push!(ls,
-                  mean(t -> sum(i -> loss(sol(t[2])[i[2]] - data[i[1], t[1]]),
-                                enumerate(vars)), enumerate(t)))
+            push!(
+                ls,
+                mean(
+                    t ->
+                        sum(i -> loss(sol(t[2])[i[2]] - data[i[1], t[1]]), enumerate(vars)),
+                    enumerate(t),
+                ),
+            )
         end
 
-        mean(ls)
+        return mean(ls)
     end
 end
 
@@ -133,7 +168,7 @@ function prep_params!(params, prob)
     end
     any(p -> (p[2] === NaN) && @warn("Uninitialized prm: $p"), params)
 
-    params
+    return params
 end
 
 # set initial model variable values in an optimization problem
@@ -143,10 +178,12 @@ function prep_u0!(u0, prob)
     end
     any(u -> (u[2] === NaN) && @warn("Uninitialized prm: $(u[1])"), u0)
 
-    u0
+    return u0
 end
 
-"Extract symbolic variables referenced in `acs`, `args`."
+"""
+Extract symbolic variables referenced in `acs`, `args`.
+"""
 function get_free_vars(acs, args)
     u0_syms = collect(acs[:, :specName])
     p_syms = collect(acs[:, :prmName])
@@ -174,16 +211,18 @@ function get_free_vars(acs, args)
         if k isa Number
             push!(u0_, Int(k) => v)
         else
-            for i in 1:length(subpart(acs, :specName))
+            for i = 1:length(subpart(acs, :specName))
                 (acs[i, :specName] == k) && (push!(u0_, i => v); break)
             end
         end
     end
 
-    u0_, p
+    return u0_, p
 end
 
-"Resolve symbolic / positional model variable names to positional."
+"""
+Resolve symbolic / positional model variable names to positional.
+"""
 function get_vars(acs, args)
     (args == :()) && return args
     args_ = []
@@ -193,12 +232,13 @@ function get_vars(acs, args)
         if arg isa Number
             push!(args_, Int(arg))
         else
-            for i in 1:length(subpart(acs, :specName))
-                !isnothing(acs[i, :specName]) && (acs[i, :specName] == arg) &&
+            for i = 1:length(subpart(acs, :specName))
+                !isnothing(acs[i, :specName]) &&
+                    (acs[i, :specName] == arg) &&
                     (push!(args_, i); break)
             end
         end
     end
 
-    args_
+    return args_
 end
diff --git a/src/solvers.jl b/src/solvers.jl
index 96d09dd..e2f22f1 100644
--- a/src/solvers.jl
+++ b/src/solvers.jl
@@ -6,29 +6,33 @@ using DiffEqBase, DifferentialEquations
 using Distributions
 
 function get_sampled_transition(state, i)
-    transition = Dict{Symbol, Any}()
+    transition = Dict{Symbol,Any}()
     foreach(k -> push!(transition, k => state[i, k]), keys(state.transitions))
 
-    transition
+    return transition
 end
 
-"Compute resource requirements given transition quantities."
+"""
+Compute resource requirements given transition quantities.
+"""
 function get_reqs_init!(reqs, qs, state)
     reqs .= 0.0
-    for i in 1:size(reqs, 2)
+    for i = 1:size(reqs, 2)
         for tok in state[i, :transLHS]
             !any(m -> m in tok.modality, [:rate, :nonblock]) &&
                 (reqs[tok.index, i] += qs[i] * tok.stoich)
         end
     end
 
-    reqs
+    return reqs
 end
 
-"Compute resource requirements given transition quantities."
+"""
+Compute resource requirements given transition quantities.
+"""
 function get_reqs_ongoing!(reqs, qs, state)
     reqs .= 0.0
-    for i in 1:length(state.ongoing_transitions)
+    for i = 1:length(state.ongoing_transitions)
         for tok in state.ongoing_transitions[i][:transLHS]
             in(:rate, tok.modality) &&
                 (state.ongoing_transitions[i][:transCycleTime] > 0) &&
@@ -37,18 +41,23 @@ function get_reqs_ongoing!(reqs, qs, state)
         end
     end
 
-    reqs
+    return reqs
 end
 
-"Given requirements, return available allocation."
+"""
+Given requirements, return available allocation.
+"""
 function get_allocs!(reqs, u, state, priorities, strategy = :weighted)
-    strategy == :weighted ? alloc_weighted!(reqs, u, priorities, state) :
-    alloc_greedy!(reqs, u, priorities, state)
+    return if strategy == :weighted
+        alloc_weighted!(reqs, u, priorities, state)
+    else
+        alloc_greedy!(reqs, u, priorities, state)
+    end
 end
 
 function alloc_weighted!(reqs, u, priorities, state)
     allocs = zero(reqs)
-    for i in 1:size(reqs, 1)
+    for i = 1:size(reqs, 1)
         s = sum(reqs[i, :])
         u[i] >= s && (allocs[i, :] .= reqs[i, :]; continue)
         foreach(j -> allocs[i, j] = reqs[i, j] * priorities[j], 1:size(reqs, 2))
@@ -56,13 +65,13 @@ function alloc_weighted!(reqs, u, priorities, state)
         allocs[i, :] .*= (s == 0) ? 0.0 : (max(0, u[i]) / s)
     end
 
-    allocs
+    return allocs
 end
 
 function alloc_greedy!(reqs, u, priorities, state)
     allocs = zero(reqs)
-    sorted_trans = sort(1:size(reqs, 2), by = i -> -priorities[i])
-    for i in 1:size(reqs, 1)
+    sorted_trans = sort(1:size(reqs, 2); by = i -> -priorities[i])
+    for i = 1:size(reqs, 1)
         s = sum(reqs[i, :])
         u[i] >= s && (allocs[i, :] .= reqs[i, :]; continue)
         a = u[i]
@@ -74,10 +83,12 @@ function alloc_greedy!(reqs, u, priorities, state)
         end
     end
 
-    allocs
+    return allocs
 end
 
-"Given resource requirements and available allocations, output resulting shift size for each transition."
+"""
+Given resource requirements and available allocations, output resulting shift size for each transition.
+"""
 function get_frac_satisfied(allocs, reqs, state)
     for i in eachindex(allocs)
         allocs[i] = min(1, (reqs[i] == 0.0 ? 1 : (allocs[i] / reqs[i])))
@@ -85,13 +96,15 @@ function get_frac_satisfied(allocs, reqs, state)
     qs = vec(minimum(allocs; dims = 1))
     foreach(i -> allocs[:, i] .= reqs[:, i] * qs[i], 1:size(reqs, 2))
 
-    qs
+    return qs
 end
 
-"Given available allocations and qties of transitions requested to spawn, return number of spawned transitions. Update `alloc` to match actual allocation."
+"""
+Given available allocations and qties of transitions requested to spawn, return number of spawned transitions. Update `alloc` to match actual allocation.
+"""
 function get_init_satisfied(allocs, qs, state)
     reqs = zero(allocs)
-    for i in 1:size(allocs, 2)
+    for i = 1:size(allocs, 2)
         all(allocs[:, i] .>= 0) || (allocs[:, i] .= 0.0; qs[i] = 0)
         for tok in state[i, :transLHS]
             !any(m -> m in tok.modality, [:rate, :nonblock]) &&
@@ -104,10 +117,12 @@ function get_init_satisfied(allocs, qs, state)
     foreach(i -> qs[i] = min(qs[i], minimum(allocs[:, i])), 1:size(reqs, 2))
     foreach(i -> allocs[:, i] .= reqs[:, i] * qs[i], 1:size(reqs, 2))
 
-    qs
+    return qs
 end
 
-"Evolve transitions, spawn new transitions."
+"""
+Evolve transitions, spawn new transitions.
+"""
 function evolve!(u, state)
     update_u!(state, u)
     actual_allocs = zero(u)
@@ -116,34 +131,43 @@ function evolve!(u, state)
     reqs = zeros(nparts(state, :S), nparts(state, :T))
     qs = zeros(nparts(state, :T))
 
-    foreach(i -> qs[i] = state[i, :transRate] * state[i, :transMultiplier],
-            1:nparts(state, :T))
+    foreach(
+        i -> qs[i] = state[i, :transRate] * state[i, :transMultiplier],
+        1:nparts(state, :T),
+    )
     qs .= ceil.(Ref(Int), qs)
-    for i in 1:nparts(state, :T)
+    for i = 1:nparts(state, :T)
         new_instances = state.solverargs[:tstep] * qs[i] + state[i, :transToSpawn]
-        capacity = state[i, :transCapacity] -
-                   count(t -> t[:transHash] == state[i, :transHash],
-                         state.ongoing_transitions)
+        capacity =
+            state[i, :transCapacity] -
+            count(t -> t[:transHash] == state[i, :transHash], state.ongoing_transitions)
         (capacity < new_instances) &&
             add_to_spawn!(state, state[i, :transHash], new_instances - capacity)
         qs[i] = min(capacity, new_instances)
     end
 
     reqs = get_reqs_init!(reqs, qs, state)
-    allocs = get_allocs!(reqs, u, state, state[:, :transPriority],
-                         state.solverargs[:strategy])
+    allocs =
+        get_allocs!(reqs, u, state, state[:, :transPriority], state.solverargs[:strategy])
     qs .= get_init_satisfied(allocs, qs, state)
 
-    push!(state.log,
-          (:new_transitions, state.t,
-           [(hash, q) for (hash, q) in zip(state[:, :transHash], qs)]...))
-    u .-= sum(allocs, dims = 2)
-    actual_allocs .+= sum(allocs, dims = 2)
+    push!(
+        state.log,
+        (
+            :new_transitions,
+            state.t,
+            [(hash, q) for (hash, q) in zip(state[:, :transHash], qs)]...,
+        ),
+    )
+    u .-= sum(allocs; dims = 2)
+    actual_allocs .+= sum(allocs; dims = 2)
 
     # add spawned transitions to the heap
-    for i in 1:nparts(state, :T)
-        qs[i] != 0 && push!(state.ongoing_transitions,
-              Transition(get_sampled_transition(state, i), state.t, qs[i], 0.0))
+    for i = 1:nparts(state, :T)
+        qs[i] != 0 && push!(
+            state.ongoing_transitions,
+            Transition(get_sampled_transition(state, i), state.t, qs[i], 0.0),
+        )
     end
 
     update_u!(state, u)
@@ -152,34 +176,52 @@ function evolve!(u, state)
     qs = map(t -> t.q, state.ongoing_transitions)
 
     get_reqs_ongoing!(reqs, qs, state)
-    allocs = get_allocs!(reqs, u, state,
-                         map(t -> t[:transPriority], state.ongoing_transitions),
-                         state.solverargs[:strategy])
+    allocs = get_allocs!(
+        reqs,
+        u,
+        state,
+        map(t -> t[:transPriority], state.ongoing_transitions),
+        state.solverargs[:strategy],
+    )
     qs .= get_frac_satisfied(allocs, reqs, state)
-    push!(state.log,
-          (:saturation, state.t,
-           [(state.ongoing_transitions[i][:transHash], qs[i])
-            for i in eachindex(state.ongoing_transitions)]...))
-    u .-= sum(allocs, dims = 2)
-    actual_allocs .+= sum(allocs, dims = 2)
-
-    foreach(i -> state.ongoing_transitions[i].state += qs[i] * state.solverargs[:tstep],
-            eachindex(state.ongoing_transitions))
+    push!(
+        state.log,
+        (
+            :saturation,
+            state.t,
+            [
+                (state.ongoing_transitions[i][:transHash], qs[i]) for
+                i in eachindex(state.ongoing_transitions)
+            ]...,
+        ),
+    )
+    u .-= sum(allocs; dims = 2)
+    actual_allocs .+= sum(allocs; dims = 2)
+
+    foreach(
+        i -> state.ongoing_transitions[i].state += qs[i] * state.solverargs[:tstep],
+        eachindex(state.ongoing_transitions),
+    )
 
     push!(state.log, (:allocation, state.t, actual_allocs))
-    push!(state.log,
-          (:valuation_cost, state.t,
-           actual_allocs' * [state[i, :specCost] for i in 1:nparts(state, :S)]))
+    return push!(
+        state.log,
+        (
+            :valuation_cost,
+            state.t,
+            actual_allocs' * [state[i, :specCost] for i = 1:nparts(state, :S)],
+        ),
+    )
 end
 
 # execute callbacks
 function event_action!(state)
-    for i in 1:nparts(state, :E)
-        !isnothing(state[i, :eventTrigger]) &&
-            !isnothing(state[i, :eventAction]) || continue
+    for i = 1:nparts(state, :E)
+        !isnothing(state[i, :eventTrigger]) && !isnothing(state[i, :eventAction]) ||
+            continue
         v = state[i, :eventTrigger]
         q = v isa Bool ? (v ? 1 : 0) : (v isa Number ? rand(Poisson(v)) : 0)
-        for _ in 1:q
+        for _ = 1:q
             state[i, :eventAction]
         end
     end
@@ -189,40 +231,53 @@ end
 function finish!(u, state)
     update_u!(state, u)
     val_reward = 0
-    terminated_all = Dict{Symbol, Float64}()
-    terminated_success = Dict{Symbol, Float64}()
+    terminated_all = Dict{Symbol,Float64}()
+    terminated_success = Dict{Symbol,Float64}()
 
     ix = 1
     while ix <= length(state.ongoing_transitions)
         trans_ = state.ongoing_transitions[ix]
         ((state.t - trans_.t) < trans_.trans[:transMaxLifeTime]) &&
-            trans_.state < trans_[:transCycleTime] && (ix += 1; continue)
+            trans_.state < trans_[:transCycleTime] &&
+            (ix += 1; continue)
         toks_rhs = []
         for r in extract_reactants(trans_[:transRHS], state)
             i = find_index(r.species, state)
-            push!(toks_rhs,
-                  UnfoldedReactant(i, r.species,
-                                   context_eval(state, state.wrap_fun(r.stoich)),
-                                   r.modality ∪ state[i, :specModality]))
+            push!(
+                toks_rhs,
+                UnfoldedReactant(
+                    i,
+                    r.species,
+                    context_eval(state, state.wrap_fun(r.stoich)),
+                    r.modality ∪ state[i, :specModality],
+                ),
+            )
         end
         for tok in trans_[:transLHS]
-            in(:conserved, tok.modality) && (u[tok.index] += trans_.q * tok.stoich *
-                             (in(:rate, tok.modality) ? trans_[:transCycleTime] : 1))
+            in(:conserved, tok.modality) && (
+                u[tok.index] +=
+                    trans_.q *
+                    tok.stoich *
+                    (in(:rate, tok.modality) ? trans_[:transCycleTime] : 1)
+            )
         end
-        q = trans_.state >= trans_[:transCycleTime] ?
-            rand(Distributions.Binomial(Int(trans_.q), trans_[:transProbOfSuccess])) : 0
-        foreach(tok -> (u[tok.index] += q * tok.stoich;
-                        val_reward += state[tok.index,
-                                            :specReward] *
-                                      q * tok.stoich),
-                toks_rhs)
+        q = if trans_.state >= trans_[:transCycleTime]
+            rand(Distributions.Binomial(Int(trans_.q), trans_[:transProbOfSuccess]))
+        else
+            0
+        end
+        foreach(
+            tok -> (u[tok.index] += q * tok.stoich;
+            val_reward += state[tok.index, :specReward] * q * tok.stoich),
+            toks_rhs,
+        )
 
         update_u!(state, u)
         context_eval(state, trans_.trans[:transPostAction])
-        terminated_all[trans_[:transHash]] = get(terminated_all, trans_[:transHash], 0) +
-                                             trans_.q
-        terminated_success[trans_[:transHash]] = get(terminated_success, trans_[:transHash],
-                                                     0) + q
+        terminated_all[trans_[:transHash]] =
+            get(terminated_all, trans_[:transHash], 0) + trans_.q
+        terminated_success[trans_[:transHash]] =
+            get(terminated_success, trans_[:transHash], 0) + q
 
         ix += 1
     end
@@ -233,7 +288,7 @@ function finish!(u, state)
     push!(state.log, (:terminated_success, state.t, terminated_success...))
     push!(state.log, (:valuation_reward, state.t, val_reward))
 
-    u
+    return u
 end
 
 function free_blocked_species!(state)
@@ -246,12 +301,16 @@ end
 Transform an `acs` to a `DiscreteProblem` instance, compatible with standard solvers.
 
 # Examples
+
 ```julia
-transform(DiscreteProblem, acs, schedule=schedule_weighted!)
+transform(DiscreteProblem, acs; schedule = schedule_weighted!)
 ```
 """
-function transform(::Type{DiffEqBase.DiscreteProblem}, state::ReactiveDynamicsState;
-                   kwargs...)
+function transform(
+    ::Type{DiffEqBase.DiscreteProblem},
+    state::ReactiveDynamicsState;
+    kwargs...,
+)
     f = function (du, u, p, t)
         state = p[:__state__]
         free_blocked_species!(state)
@@ -264,21 +323,26 @@ function transform(::Type{DiffEqBase.DiscreteProblem}, state::ReactiveDynamicsSt
         event_action!(state)
 
         du .= state.u
-        push!(state.log,
-              (:valuation, t,
-               du' * [state[i, :specValuation] for i in 1:nparts(state, :S)]))
+        push!(
+            state.log,
+            (:valuation, t, du' * [state[i, :specValuation] for i = 1:nparts(state, :S)]),
+        )
 
         t = (state.t += state.solverargs[:tstep])
         update_u!(state, du)
         save!(state)
         sync_p!(p, state)
 
-        du
+        return du
     end
 
-    DiffEqBase.DiscreteProblem(f, state.u, (0.0, 2.0),
-                               Dict(state.p..., :__state__ => state,
-                                    :__state0__ => deepcopy(state)); kwargs...)
+    return DiffEqBase.DiscreteProblem(
+        f,
+        state.u,
+        (0.0, 2.0),
+        Dict(state.p..., :__state__ => state, :__state0__ => deepcopy(state));
+        kwargs...,
+    )
 end
 
 ## resolve tspan, tstep
@@ -290,7 +354,7 @@ function get_tcontrol(tspan, args)
 
     tstep = get(args, :tstep, haskey(args, :tstops) ? tspan / args[:tstops] : tunit) / tunit
 
-    ((0.0, tspan), tstep)
+    return ((0.0, tspan), tstep)
 end
 
 """
@@ -299,48 +363,75 @@ Transform an `acs` to a `DiscreteProblem` instance, compatible with standard sol
 Optionally accepts initial values and parameters, which take precedence over specifications in `acs`.
 
 # Examples
+
 ```julia
-DiscreteProblem(acs, u0, p; tspan=(.0, 100.), schedule=schedule_weighted!)
+DiscreteProblem(acs, u0, p; tspan = (0.0, 100.0), schedule = schedule_weighted!)
 ```
 """
-function DiffEqBase.DiscreteProblem(acs::ReactionNetwork, u0 = Dict(),
-                                    p = DiffEqBase.NullParameters(); kwargs...)
+function DiffEqBase.DiscreteProblem(
+    acs::ReactionNetwork,
+    u0 = Dict(),
+    p = DiffEqBase.NullParameters();
+    kwargs...,
+)
     assign_defaults!(acs)
-    keywords = Dict{Symbol, Any}([acs[i, :metaKeyword] => acs[i, :metaVal]
-                                  for i in 1:nparts(acs, :M)
-                                  if !isnothing(acs[i, :metaKeyword]) &&
-                                     !isnothing(acs[i, :metaVal])])
+    keywords = Dict{Symbol,Any}([
+        acs[i, :metaKeyword] => acs[i, :metaVal] for i = 1:nparts(acs, :M) if
+        !isnothing(acs[i, :metaKeyword]) && !isnothing(acs[i, :metaVal])
+    ])
     merge!(keywords, Dict(collect(kwargs)))
     merge!(keywords, Dict(:strategy => get(keywords, :alloc_strategy, :weighted)))
     keywords[:tspan], keywords[:tstep] = get_tcontrol(keywords[:tspan], keywords)
 
     acs = remove_choose(acs)
     attrs, transitions, wrap_fun = compile_attrs(acs)
-    state = ReactiveDynamicsState(acs, attrs, transitions, wrap_fun, keywords[:tspan][1];
-                                  keywords...)
+    state = ReactiveDynamicsState(
+        acs,
+        attrs,
+        transitions,
+        wrap_fun,
+        keywords[:tspan][1];
+        keywords...,
+    )
     init_u!(state)
     save!(state)
 
     prob = transform(DiffEqBase.DiscreteProblem, state; kwargs...)
 
-    u0 isa Dict &&
-        foreach(i -> prob.u0[i] = !isnothing(acs[i, :specName]) &&
-                                  haskey(u0, acs[i, :specName]) ? u0[acs[i, :specName]] :
-                                  prob.u0[i], 1:nparts(state, :S))
+    u0 isa Dict && foreach(
+        i ->
+            prob.u0[i] =
+                if !isnothing(acs[i, :specName]) && haskey(u0, acs[i, :specName])
+                    u0[acs[i, :specName]]
+                else
+                    prob.u0[i]
+                end,
+        1:nparts(state, :S),
+    )
     p_ = p == DiffEqBase.NullParameters() ? Dict() : Dict(k => v for (k, v) in p)
-    prob = remake(prob; u0 = prob.u0, tspan = keywords[:tspan],
-                  dt = get(keywords, :tstep, 1),
-                  p = merge(prob.p, p_,
-                            Dict(:tstep => get(keywords, :tstep, 1),
-                                 :strategy => get(keywords, :alloc_strategy, :weighted))))
-
-    prob
+    prob = remake(
+        prob;
+        u0 = prob.u0,
+        tspan = keywords[:tspan],
+        dt = get(keywords, :tstep, 1),
+        p = merge(
+            prob.p,
+            p_,
+            Dict(
+                :tstep => get(keywords, :tstep, 1),
+                :strategy => get(keywords, :alloc_strategy, :weighted),
+            ),
+        ),
+    )
+
+    return prob
 end
 
 function fetch_params(acs::ReactionNetwork)
-    Dict{Symbol, Any}((acs[i, :prmName] => acs[i, :prmVal]
-                       for i in Iterators.filter(i -> !isnothing(acs[i, :prmVal]),
-                                                 1:nparts(acs, :P))))
+    return Dict{Symbol,Any}((
+        acs[i, :prmName] => acs[i, :prmVal] for
+        i in Iterators.filter(i -> !isnothing(acs[i, :prmVal]), 1:nparts(acs, :P))
+    ))
 end
 
 # EnsembleProblem's prob_func: sample initial values
@@ -351,13 +442,16 @@ function get_prob_func(prob)
         prob.p[:__state__] = deepcopy(prob.p[:__state0__])
         for i in eachindex(prob.u0)
             rv = randn() * vars[i]
-            prob.u0[i] = (sign(rv + prob.u0[i]) == sign(prob.u0[i])) ? rv + prob.u0[i] :
-                         prob.u0[i]
+            prob.u0[i] = if (sign(rv + prob.u0[i]) == sign(prob.u0[i]))
+                rv + prob.u0[i]
+            else
+                prob.u0[i]
+            end
         end
         sync!(prob.p[:__state__], prob.u0, prob.p)
 
-        prob
+        return prob
     end
 
-    prob_func
+    return prob_func
 end
diff --git a/src/state.jl b/src/state.jl
index 1b51273..9567317 100644
--- a/src/state.jl
+++ b/src/state.jl
@@ -7,9 +7,11 @@ struct UnfoldedReactant
     modality::Set{Symbol}
 end
 
-"Ongoing transition auxiliary structure."
+"""
+Ongoing transition auxiliary structure.
+"""
 mutable struct Transition
-    trans::Dict{Symbol, Any}
+    trans::Dict{Symbol,Any}
 
     t::Float64
     q::Float64
@@ -20,7 +22,7 @@ Base.getindex(state::Transition, key) = state.trans[key]
 
 mutable struct Observable
     last::Float64 # last sampling time
-    range::Vector{Union{Tuple{Float64, SampleableValues}, SampleableValues}}
+    range::Vector{Union{Tuple{Float64,SampleableValues},SampleableValues}}
     every::Float64
     on::Vector{ActionableValues}
 
@@ -30,37 +32,58 @@ end
 mutable struct ReactiveDynamicsState
     acs::ReactionNetwork
 
-    attrs::Dict{Symbol, Vector}
-    transition_recipes::Dict{Symbol, Vector}
+    attrs::Dict{Symbol,Vector}
+    transition_recipes::Dict{Symbol,Vector}
 
     u::Vector{Float64}
     p::Any
     t::Float64
 
-    transitions::Dict{Symbol, Vector}
+    transitions::Dict{Symbol,Vector}
     ongoing_transitions::Vector{Transition}
     log::Vector{Tuple}
 
-    observables::Dict{Symbol, Observable}
+    observables::Dict{Symbol,Observable}
     solverargs::Any
 
     wrap_fun::Any
     history_u::Vector{Vector{Float64}}
     history_t::Vector{Float64}
 
-    function ReactiveDynamicsState(acs::ReactionNetwork, attrs, transition_recipes,
-                                   wrap_fun, t0 = 0; kwargs...)
+    function ReactiveDynamicsState(
+        acs::ReactionNetwork,
+        attrs,
+        transition_recipes,
+        wrap_fun,
+        t0 = 0;
+        kwargs...,
+    )
         ongoing_transitions = Transition[]
         log = NamedTuple[]
         observables = compile_observables(acs)
-        transitions_attrs = setdiff(filter(a -> contains(string(a), "trans"),
-                                           propertynames(acs.subparts)), (:trans,)) ∪
-                            [:transLHS, :transRHS, :transToSpawn, :transHash]
-        transitions = Dict{Symbol, Vector}(a => [] for a in transitions_attrs)
-
-        new(acs, attrs, transition_recipes, zeros(nparts(acs, :S)), fetch_params(acs), t0,
-            transitions, ongoing_transitions, log, observables, kwargs, wrap_fun,
-            Vector{Float64}[], Float64[])
+        transitions_attrs =
+            setdiff(
+                filter(a -> contains(string(a), "trans"), propertynames(acs.subparts)),
+                (:trans,),
+            ) ∪ [:transLHS, :transRHS, :transToSpawn, :transHash]
+        transitions = Dict{Symbol,Vector}(a => [] for a in transitions_attrs)
+
+        return new(
+            acs,
+            attrs,
+            transition_recipes,
+            zeros(nparts(acs, :S)),
+            fetch_params(acs),
+            t0,
+            transitions,
+            ongoing_transitions,
+            log,
+            observables,
+            kwargs,
+            wrap_fun,
+            Vector{Float64}[],
+            Float64[],
+        )
     end
 end
 
@@ -69,41 +92,45 @@ end
 function context_eval(state::ReactiveDynamicsState, o)
     o = o isa Function ? Base.invokelatest(o, state) : o
 
-    o isa Sampleable ? rand(o) : o
+    return o isa Sampleable ? rand(o) : o
 end
 
 function Base.getindex(state::ReactiveDynamicsState, keys...)
-    context_eval(state,
-                 (contains(string(keys[2]), "trans") ? state.transitions : state.attrs)[keys[2]][keys[1]])
+    return context_eval(
+        state,
+        (contains(string(keys[2]), "trans") ? state.transitions : state.attrs)[keys[2]][keys[1]],
+    )
 end
 
 function init_u!(state::ReactiveDynamicsState)
-    (u = fill(0.0, nparts(state, :S));
-     foreach(i -> u[i] = state[i, :specInitVal],
-             1:nparts(state, :S));
-     state.u = u)
+    return (u = fill(0.0, nparts(state, :S));
+    foreach(i -> u[i] = state[i, :specInitVal], 1:nparts(state, :S));
+    state.u = u)
 end
 function save!(state::ReactiveDynamicsState)
-    (push!(state.history_u, state.u); push!(state.history_t, state.t))
+    return (push!(state.history_u, state.u); push!(state.history_t, state.t))
 end
 
 function compile_observables(acs::ReactionNetwork)
-    observables = Dict{Symbol, Observable}()
+    observables = Dict{Symbol,Observable}()
     species_names = collect(acs[:, :specName])
     prm_names = collect(acs[:, :prmName])
     varmap = Dict([name => :(state.u[$i]) for (i, name) in enumerate(species_names)])
 
     for (name, opts) in Iterators.zip(acs[:, :obsName], acs[:, :obsOpts])
         on = map(on -> wrap_expr(on, species_names, prm_names, varmap), opts.on)
-        range = map(r -> begin
-                        r = r isa Tuple ? r : (1.0, r)
-                        (r[1], wrap_expr(r[2], species_names, prm_names, varmap))
-                    end, opts.range)
+        range = map(
+            r -> begin
+                r = r isa Tuple ? r : (1.0, r)
+                (r[1], wrap_expr(r[2], species_names, prm_names, varmap))
+            end,
+            opts.range,
+        )
 
         push!(observables, name => Observable(-Inf, range, opts.every, on, missing))
     end
 
-    observables
+    return observables
 end
 
 compileval(ex, state) = !isa(ex, Expr) ? ex : eval(state.wrap_fun(ex))
@@ -119,21 +146,23 @@ function sample_range(rng, state)
     end
 
     r = rng[ix] isa Tuple ? rng[ix][2] : rng[ix]
-    r isa Sampleable ? rand(r) : r
+    return r isa Sampleable ? rand(r) : r
 end
 
 function resample!(state::ReactiveDynamicsState, o::Observable)
     o.last = state.t
     isempty(o.range) && (return o.val = missing)
 
-    o.sampled = context_eval(state, sample_range(o.range, state))
+    return o.sampled = context_eval(state, sample_range(o.range, state))
 end
 
 resample(state::ReactiveDynamicsState, o::Symbol) = resample!(state, state.observables[o])
 
 function update_observables(state::ReactiveDynamicsState)
-    foreach(o -> (state.t - o.last) >= o.every && resample!(state, o),
-            values(state.observables))
+    return foreach(
+        o -> (state.t - o.last) >= o.every && resample!(state, o),
+        values(state.observables),
+    )
 end
 
 function prune_r_line(r_line)
@@ -142,41 +171,58 @@ function prune_r_line(r_line)
     elseif r_line isa Expr && r_line.args[1] ∈ bwd_arrows
         r_line.args[[3, 2]]
     elseif isexpr(r_line, :macrocall) && (macroname(r_line) == :choose)
-        sample_range([(isexpr(r, :tuple) ? (r.args[1], prune_r_line(r.args[2])) :
-                       prune_r_line(r)) for r in r_line.args[3:end]], state)
+        sample_range(
+            [(
+                if isexpr(r, :tuple)
+                    (r.args[1], prune_r_line(r.args[2]))
+                else
+                    prune_r_line(r)
+                end
+            ) for r in r_line.args[3:end]],
+            state,
+        )
     end
 end
 
 function find_index(species::Symbol, state::ReactiveDynamicsState)
-    findfirst(i -> state[i, :specName] == species, 1:nparts(state, :S))
+    return findfirst(i -> state[i, :specName] == species, 1:nparts(state, :S))
 end
 
 function sample_transitions!(state::ReactiveDynamicsState)
     for (_, v) in state.transitions
         empty!(v)
     end
-    for i in 1:length(state.transition_recipes[:trans])
+    for i = 1:length(state.transition_recipes[:trans])
         !(state.transition_recipes[:transActivated][i]) && continue
         l_line, r_line = prune_r_line(state.transition_recipes[:trans][i])
 
         for attr in keys(state.transition_recipes)
             attr ∈ [:trans, :transPostAction, :transActivated, :transHash] && continue
-            push!(state.transitions[attr],
-                  context_eval(state, state.transition_recipes[attr][i]))
+            push!(
+                state.transitions[attr],
+                context_eval(state, state.transition_recipes[attr][i]),
+            )
         end
 
         reactants = []
         for r in extract_reactants(l_line, state)
             j = find_index(r.species, state)
-            push!(reactants,
-                  UnfoldedReactant(j, r.species,
-                                   context_eval(state, state.wrap_fun(r.stoich)),
-                                   r.modality ∪ state[j, :specModality]))
+            push!(
+                reactants,
+                UnfoldedReactant(
+                    j,
+                    r.species,
+                    context_eval(state, state.wrap_fun(r.stoich)),
+                    r.modality ∪ state[j, :specModality],
+                ),
+            )
         end
         push!(state.transitions[:transLHS], reactants)
         push!(state.transitions[:transRHS], r_line)
-        foreach(k -> push!(state.transitions[k], state.transition_recipes[k][i]),
-                [:transPostAction, :transToSpawn, :transHash])
+        foreach(
+            k -> push!(state.transitions[k], state.transition_recipes[k][i]),
+            [:transPostAction, :transToSpawn, :transHash],
+        )
         state.transition_recipes[:transToSpawn] .= 0
     end
 end
@@ -194,11 +240,11 @@ function sync!(state::ReactiveDynamicsState, u, p)
 end
 
 function as_state(u, t, state::ReactiveDynamicsState)
-    (state = deepcopy(state); state.u .= u; state.t = t; state)
+    return (state = deepcopy(state); state.u .= u; state.t = t; state)
 end
 
 function Catlab.CategoricalAlgebra.nparts(state::ReactiveDynamicsState, obj::Symbol)
-    obj == :T ? length(state.transitions[:transLHS]) : nparts(state.acs, obj)
+    return obj == :T ? length(state.transitions[:transLHS]) : nparts(state.acs, obj)
 end
 
 ## query the state
@@ -210,18 +256,18 @@ log(state::ReactiveDynamicsState, msg) = (println(msg); push!(state.log, (:log,
 state(state::ReactiveDynamicsState) = state
 
 function periodic(state::ReactiveDynamicsState, period)
-    period == 0.0 ||
-        (length(state.history_t) > 1 &&
-         (fld(state.t, period) - fld(state.history_t[end - 1], period) > 0))
+    return period == 0.0 || (
+        length(state.history_t) > 1 &&
+        (fld(state.t, period) - fld(state.history_t[end-1], period) > 0)
+    )
 end
 
-function set_params(state::ReactiveDynamicsState, vals...)
+set_params(state::ReactiveDynamicsState, vals...) =
     for (p, v) in vals
         state.p[p] = v
     end
-end
 
 function add_to_spawn!(state::ReactiveDynamicsState, hash, n)
     ix = findfirst(ix -> state.transition_recipes[:transHash][ix] == hash)
-    !isnothing(ix) && (state.transition_recipes[:transHash][ix] += n)
+    return !isnothing(ix) && (state.transition_recipes[:transHash][ix] += n)
 end
diff --git a/src/utils/safeinclude.jl b/src/utils/safeinclude.jl
index d25e9cf..f3c50e5 100644
--- a/src/utils/safeinclude.jl
+++ b/src/utils/safeinclude.jl
@@ -10,7 +10,9 @@ macro safeinclude(args...)
         td = mktempdir()
         cp(dirname($ex), td; force = true)
         cd(td)
-        @testset $name begin include(joinpath(td, basename($ex))) end
+        @testset $name begin
+            include(joinpath(td, basename($ex)))
+        end
         cd(path)
     end
 end
diff --git a/src/utils/utils.jl b/src/utils/utils.jl
index 1fccd3c..2ed10e4 100644
--- a/src/utils/utils.jl
+++ b/src/utils/utils.jl
@@ -2,34 +2,49 @@
 
 using MacroTools: striplines
 
-"Return array of keyword expressions in `args`."
+"""
+Return array of keyword expressions in `args`.
+"""
 function kwargize(args)
     kwargs = Any[]
     map(el -> isexpr(el, :(=)) && push!(kwargs, Expr(:kw, el.args[1], el.args[2])), args)
 
-    kwargs
+    return kwargs
 end
 
-"""Split an array of expressions into arrays of "argument" expressions and keyword expressions."""
+"""
+Split an array of expressions into arrays of "argument" expressions and keyword expressions.
+"""
 function args_kwargs(args)
     args_ = Any[]
     kwargs = Any[]
-    map(el -> isexpr(el, :(=)) ? push!(kwargs, Expr(:kw, el.args[1], el.args[2])) :
-              push!(args_, esc(el)), args)
+    map(el -> if isexpr(el, :(=))
+        push!(kwargs, Expr(:kw, el.args[1], el.args[2]))
+    else
+        push!(args_, esc(el))
+    end, args)
 
-    args_, kwargs
+    return args_, kwargs
 end
 
-"Return the (expression) value stored for the given key in a collection of keyword expression, or the given default value if no mapping for the key is present."
+"""
+Return the (expression) value stored for the given key in a collection of keyword expression, or the given default value if no mapping for the key is present.
+"""
 function find_kwargex_delete!(collection, key, default = :())
     ix = findfirst(ex -> ex.args[1] == key, collection)
-    !isnothing(ix) ? (v = collection[ix].args[2]; deleteat!(collection, ix); v) : default
+    return if !isnothing(ix)
+        (v = collection[ix].args[2]; deleteat!(collection, ix); v)
+    else
+        default
+    end
 end
 
-function macroname(ex)
-    isexpr(ex, :macrocall) ? (str = string(ex.args[1]); Symbol(strip(str, '@'))) :
-    error("expr $ex is not a macrocall")
-end
+macroname(ex) =
+    if isexpr(ex, :macrocall)
+        (str = string(ex.args[1]); Symbol(strip(str, '@')))
+    else
+        error("expr $ex is not a macrocall")
+    end
 strip_sym(sym) = (str = string(sym); Symbol(strip(str, '@')))
 
 blockize(ex) = striplines(isexpr(ex, :block) ? ex : Expr(:block, ex))
@@ -45,12 +60,15 @@ function unblock_shallow!(ex)
         isexpr(ex, :block) ? append!(args, ex.args) : push!(args, ex)
     end
     ex.args = args
-    ex
+    return ex
 end
 
 function underscorize(ex)
-    ex isa Number ? ex :
-    (str = string(ex); replace(str, '.' => "__", '(' => "", ')' => "") |> Symbol)
+    return if ex isa Number
+        ex
+    else
+        (str = string(ex); replace(str, '.' => "__", '(' => "", ')' => "") |> Symbol)
+    end
 end
 
 wkeys(itr) = map(x -> x isa Pair ? x[1] : x, itr)
@@ -60,13 +78,17 @@ function wset!(col, key, val)
     ix = findfirst(==(key), wkeys(col))
     !isnothing(ix) && (col[ix] = (key => val))
 
-    col
+    return col
 end
 
-"Return the (expression) value stored for the given key in a collection of keyword expression, or the given default value if no mapping for the key is present."
+"""
+Return the (expression) value stored for the given key in a collection of keyword expression, or the given default value if no mapping for the key is present.
+"""
 function get_kwarg(collection, key, default = :())
-    ix = findfirst(ex -> ex isa Expr && hasproperty(ex, :args) && (ex.args[1] == key),
-                   collection)
+    ix = findfirst(
+        ex -> ex isa Expr && hasproperty(ex, :args) && (ex.args[1] == key),
+        collection,
+    )
 
-    !isnothing(ix) ? collection[ix].args[2] : default
+    return !isnothing(ix) ? collection[ix].args[2] : default
 end
diff --git a/test/runtests.jl b/test/runtests.jl
index 6f8d9cb..af2e288 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,5 +1,7 @@
 using SafeTestsets, BenchmarkTools
 
 @time begin
-    @time @safetestset "Tutorial tests" begin include("tutorial_tests.jl") end
+    @time @safetestset "Tutorial tests" begin
+        include("tutorial_tests.jl")
+    end
 end
diff --git a/tutorial/basics.jl b/tutorial/basics.jl
index 071f762..3b0267e 100644
--- a/tutorial/basics.jl
+++ b/tutorial/basics.jl
@@ -4,18 +4,20 @@ using ReactiveDynamics
 #     1.0, X ⟺ Y
 # end
 
-acs = @ReactionNetwork begin 1.0, X ⟺ Y, name => "transition1" end
+acs = @ReactionNetwork begin
+    1.0, X ⟺ Y, name => "transition1"
+end
 
-@prob_init acs X=10 Y=20
+@prob_init acs X = 10 Y = 20
 @prob_params acs
-@prob_meta acs tspan=250 dt=0.1
+@prob_meta acs tspan = 250 dt = 0.1
 
 prob = @problematize acs
 
 # sol = ReactiveDynamics.solve(prob)
 
-sol = @solve prob trajectories=20
+sol = @solve prob trajectories = 20
 
 using Plots
 
-@plot sol plot_type=summary
+@plot sol plot_type = summary
diff --git a/tutorial/example.jl b/tutorial/example.jl
index 06625bc..b7fb3d1 100644
--- a/tutorial/example.jl
+++ b/tutorial/example.jl
@@ -4,11 +4,12 @@ using ReactiveDynamics
 sir_acs = @ReactionNetwork
 
 # set up ontology: try ?@aka
-@aka sir_acs transition=reaction species=population_group
+@aka sir_acs transition = reaction species = population_group
 # add single reaction
-@push sir_acs β*S*I*tdecay(@t()) S +
-                                 I-->@choose((1.5, 2 * @choose(1, 2) * I),
-                                             (1.6, @register(:RATE, 1, every=2) * I)) name=>SI2I
+@push sir_acs β * S * I * tdecay(@t()) S + I --> @choose(
+    (1.5, 2 * @choose(1, 2) * I),
+    (1.6, @register(:RATE, 1, every = 2) * I)
+) name => SI2I
 
 # additional dynamics
 @push sir_acs begin
@@ -20,30 +21,30 @@ end
 @register bool_cond(t) = (100 < t < 200) || (400 < t < 500)
 @register tdecay(t) = exp(-t / 10^3)
 
-@valuation sir_acs I=0.1 # for testing purpose only
+@valuation sir_acs I = 0.1 # for testing purpose only
 
 # atomic data: initial values, parameter values
-@prob_init sir_acs S=999 I=100 R=100
+@prob_init sir_acs S = 999 I = 100 R = 100
 u0 = [999, 10, 0] # alternative specification
 @prob_init sir_acs u0
-@prob_params sir_acs β=0.0001 ν=0.01 γ=5
-@prob_meta sir_acs tspan=100
+@prob_params sir_acs β = 0.0001 ν = 0.01 γ = 5
+@prob_meta sir_acs tspan = 100
 
 #prob = @problematize sir_acs
-prob = @problematize sir_acs tspan=200
+prob = @problematize sir_acs tspan = 200
 
-sol = @solve prob trajectories=20
-@plot sol plot_type=summary
+sol = @solve prob trajectories = 20
+@plot sol plot_type = summary
 sol = @solve prob
-@plot sol plot_type=allocation
-@plot sol plot_type=valuations
-@plot sol plot_type=new_transitions
+@plot sol plot_type = allocation
+@plot sol plot_type = valuations
+@plot sol plot_type = new_transitions
 
 # acs as a model : incomplete dynamics
 sir_acs = @ReactionNetwork
 
 # set up ontology: try ?@aka
-@aka sir_acs transition=reaction species=population_group
+@aka sir_acs transition = reaction species = population_group
 # add single reaction
 
 # additional dynamics
@@ -51,17 +52,17 @@ sir_acs = @ReactionNetwork
     ν * I * tdecay(@t()), I --> @choose(E, R, S), name => I2R
     γ, R --> S, name => R2S
 end
-@jump sir_acs 3 (S > 0&&bool_cond(@t()) && (I += 1; S -= 1)) # drift term, support for event conditioning
+@jump sir_acs 3 (S > 0 && bool_cond(@t()) && (I += 1; S -= 1)) # drift term, support for event conditioning
 @periodic sir_acs 5.0 (β += 1; println(β))
 # register custom function
 @register bool_cond(t) = (100 < t < 200) || (400 < t < 500)
 @register tdecay(t) = exp(-t / 10^3)
 
 # atomic data: initial values, parameter values
-@prob_init sir_acs S=999 I=10 R=0
+@prob_init sir_acs S = 999 I = 10 R = 0
 u0 = [999, 10, 0] # alternative specification
 @prob_init sir_acs u0
-@prob_params sir_acs β=0.0001 ν=0.01 γ=5
+@prob_params sir_acs β = 0.0001 ν = 0.01 γ = 5
 
 # model dynamics
 sir_acs = @ReactionNetwork begin
@@ -70,20 +71,24 @@ sir_acs = @ReactionNetwork begin
 end
 
 # simulation parameters
-@prob_init sir_acs S=999 I=10 R=0
-@prob_uncertainty sir_acs S=10.0 I=5.0
-@prob_params sir_acs α=0.0001 β=0.01
-@prob_meta sir_acs tspan=250 dt=0.1
+@prob_init sir_acs S = 999 I = 10 R = 0
+@prob_uncertainty sir_acs S = 10.0 I = 5.0
+@prob_params sir_acs α = 0.0001 β = 0.01
+@prob_meta sir_acs tspan = 250 dt = 0.1
 
 # batch simulation
 prob = @problematize sir_acs
-sol = @solve prob trajectories=20
-@plot sol plot_type=summary
-@plot sol plot_type=summary show=:S
-@plot sol plot_type=summary c=:green xlimits=(0.0, 100.0)
+sol = @solve prob trajectories = 20
+@plot sol plot_type = summary
+@plot sol plot_type = summary show = :S
+@plot sol plot_type = summary c = :green xlimits = (0.0, 100.0)
 
-acs_1 = @ReactionNetwork begin 1.0, A --> B + C end
-acs_2 = @ReactionNetwork begin 1.0, A --> B + C end
+acs_1 = @ReactionNetwork begin
+    1.0, A --> B + C
+end
+acs_2 = @ReactionNetwork begin
+    1.0, A --> B + C
+end
 
-acs_union = @join acs_1 acs_2 acs_1.A=acs_2.A = @alias(A) @catchall(B) # @catchall: equalize species "B" from the joined submodels
-@equalize acs_union acs_1.C=acs_2.C = @alias(C) # the species can be set equal either within the call to @join the acsets, but you may always equalize a set of species of any given acset
+acs_union = @join acs_1 acs_2 acs_1.A = acs_2.A = @alias(A) @catchall(B) # @catchall: equalize species "B" from the joined submodels
+@equalize acs_union acs_1.C = acs_2.C = @alias(C) # the species can be set equal either within the call to @join the acsets, but you may always equalize a set of species of any given acset
diff --git a/tutorial/joins/joins.jl b/tutorial/joins/joins.jl
index bb0f6d5..ee87d7a 100644
--- a/tutorial/joins/joins.jl
+++ b/tutorial/joins/joins.jl
@@ -13,40 +13,46 @@ rd_models = ReactiveDynamics.ReactionNetwork[] # submodels
 end
 
 # submodels: dense interactions
-@generate {@fileval(submodel.jl, i=$i, r=r), i = 1:n_models}
+@generate {@fileval(submodel.jl, i = $i, r = r), i = 1:n_models}
 
 # batch join over the submodels
 rd_model = @generate "@join {rd_models[\$i], i=1:n_models, dlm=' '}"
 
 # identify resources
 @generate {
-           @equalize(rd_model,
-                     @alias(resource[$j])={rd_models[$i].resource[$j], i = 1:n_models,
-                                           dlm = :(=)}), j = 1:r}
+    @equalize(
+        rd_model,
+        @alias(resource[$j]) = {rd_models[$i].resource[$j], i = 1:n_models, dlm = :(=)}
+    ),
+    j = 1:r,
+}
 
 # sparse off-diagonal interactions, sparse declaration
 sparse_off_diagonal = zeros(sum(ReactiveDynamics.ns), sum(ReactiveDynamics.ns))
-for i in 1:n_models
+for i = 1:n_models
     j = rand(setdiff(1:n_models, (i,)))
     i_ix = rand(1:ReactiveDynamics.ns[i])
     j_ix = rand(1:ReactiveDynamics.ns[j])
-    sparse_off_diagonal[i_ix + sum(ReactiveDynamics.ns[1:(i - 1)]), j_ix + sum(ReactiveDynamics.ns[1:(j - 1)])] += 1
+    sparse_off_diagonal[
+        i_ix+sum(ReactiveDynamics.ns[1:(i-1)]),
+        j_ix+sum(ReactiveDynamics.ns[1:(j-1)]),
+    ] += 1
     interaction_ex = """@push rd_model begin 1., var"rd_models[$i].state[$i_ix]" --> var"rd_models[$j]__state[$j_ix]" end"""
     eval(Meta.parseall(interaction_ex))
 end
 
 sparse_off_diagonal += cat(ReactiveDynamics.M...; dims = (1, 2))
 using Plots;
-heatmap(1 .- sparse_off_diagonal, color = :greys, legend = false);
+heatmap(1 .- sparse_off_diagonal; color = :greys, legend = false);
 
 using ReactiveDynamics: nparts
 u0 = rand(1:1000, nparts(rd_model, :S))
 @prob_init rd_model u0
 
-@prob_meta rd_model tspan=10
+@prob_meta rd_model tspan = 10
 
 prob = @problematize rd_model
-sol = @solve prob trajectories=2
+sol = @solve prob trajectories = 2
 
 # plot "state" species only
-@plot sol plot_type=summary show=r"state"
+@plot sol plot_type = summary show = r"state"
diff --git a/tutorial/joins/submodel.jl b/tutorial/joins/submodel.jl
index d288d66..1e8fc43 100644
--- a/tutorial/joins/submodel.jl
+++ b/tutorial/joins/submodel.jl
@@ -11,12 +11,13 @@
 end
 
 # generate submodel dynamics
-push!(rd_models,
-      @ReactionNetwork begin M[$i][$m, $n],
-                             state[$m] +
-                             {demand[$i][$m, $n, $l] * resource[$l], l = 1:($r), dlm = +
-                              } -->
-                             state[$n] +
-                             {production[$i][$m, $n, $l] * resource[$l], l = 1:($r), dlm = +
-                              }, cycle_time => cycle_times[$i][$m, $n],
-                             probability_of_success => $m * $n / (n[$i])^2 end m=1:ReactiveDynamics.ns[$i] n=1:ReactiveDynamics.ns[$i])
+push!(
+    rd_models,
+    @ReactionNetwork begin
+        M[$i][$m, $n],
+        state[$m] + {demand[$i][$m, $n, $l] * resource[$l], l = 1:($r), dlm = +} -->
+        state[$n] + {production[$i][$m, $n, $l] * resource[$l], l = 1:($r), dlm = +},
+        cycle_time => cycle_times[$i][$m, $n],
+        probability_of_success => $m * $n / (n[$i])^2
+    end m = 1:ReactiveDynamics.ns[$i] n = 1:ReactiveDynamics.ns[$i]
+)
diff --git a/tutorial/loadsave/loadsave.jl b/tutorial/loadsave/loadsave.jl
index 1a3f51b..23d2be7 100644
--- a/tutorial/loadsave/loadsave.jl
+++ b/tutorial/loadsave/loadsave.jl
@@ -6,14 +6,14 @@ using ReactiveDynamics
 @assert isdefined(ReactiveDynamics, :tdecay)
 
 prob = @problematize sir_acs
-sol = @solve prob trajectories=20
+sol = @solve prob trajectories = 20
 
 @import_network "csv/model.csv" sir_acs_
 @assert @isdefined sir_acs_
 @assert isdefined(ReactiveDynamics, :foo)
 
 prob_ = @problematize sir_acs_
-sol_ = @solve prob_ trajectories=20
+sol_ = @solve prob_ trajectories = 20
 
 # export, import the solution
 @export_solution sol
diff --git a/tutorial/optimize/optimize.jl b/tutorial/optimize/optimize.jl
index efc81fa..0766369 100644
--- a/tutorial/optimize/optimize.jl
+++ b/tutorial/optimize/optimize.jl
@@ -8,19 +8,21 @@ acss = @ReactionNetwork begin
 end
 
 # initial values, check params
-@prob_init acss A=100.0 B=200 C=150.0
-@prob_params acss α=10
-@prob_meta acss tspan=100.0
+@prob_init acss A = 100.0 B = 200 C = 150.0
+@prob_params acss α = 10
+@prob_meta acss tspan = 100.0
 
 prob = @problematize acss
-sol = @solve prob trajectories=10
+sol = @solve prob trajectories = 10
 
 # check https://github.com/JuliaOpt/NLopt.jl for solver opts
-@optimize acss abs(A - B) A B=20.0 α=2.0 lower_bounds=0 upper_bounds=200 min_t=50 max_t=100 maxeval=20 final_only=true
+@optimize acss abs(A - B) A B = 20.0 α = 2.0 lower_bounds = 0 upper_bounds = 200 min_t = 50 max_t =
+    100 maxeval = 20 final_only = true
 
 t_ = [1, 50, 100]
 data = [70 50 50]
 
-@fit acss data t_ vars=[A] B=20 A α lower_bounds=0 upper_bounds=200 maxeval=20
+@fit acss data t_ vars = [A] B = 20 A α lower_bounds = 0 upper_bounds = 200 maxeval = 20
 
-@fit_and_plot acss data t_ vars=[A] B=20 A α lower_bounds=0 upper_bounds=200 maxeval=2 trajectories=10
+@fit_and_plot acss data t_ vars = [A] B = 20 A α lower_bounds = 0 upper_bounds = 200 maxeval =
+    2 trajectories = 10
diff --git a/tutorial/optimize/optimize_custom.jl b/tutorial/optimize/optimize_custom.jl
index db46b8d..47dc489 100644
--- a/tutorial/optimize/optimize_custom.jl
+++ b/tutorial/optimize/optimize_custom.jl
@@ -2,9 +2,11 @@ using ReactiveDynamics
 
 # fit unknown dynamics to empirical data
 ## some embedded neural network, etc.
-@register begin function function_to_learn(A, B, C, params)
-    [A, B, C]' * params # params: 3-element vector
-end end
+@register begin
+    function function_to_learn(A, B, C, params)
+        return [A, B, C]' * params # params: 3-element vector
+    end
+end
 
 acs = @ReactionNetwork begin
     function_to_learn(A, B, C, params), A --> B + C
@@ -13,9 +15,9 @@ acs = @ReactionNetwork begin
 end
 
 # initial values, check params
-@prob_init acs A=60.0 B=10.0 C=150.0
-@prob_params acs params=[0.01, 0.01, 0.01]
-@prob_meta acs tspan=100.0
+@prob_init acs A = 60.0 B = 10.0 C = 150.0
+@prob_params acs params = [0.01, 0.01, 0.01]
+@prob_meta acs tspan = 100.0
 
 prob = @problematize acs
 sol = @solve prob
@@ -23,13 +25,16 @@ sol = @solve prob
 time_points = [1, 50, 100]
 data = [60 30 5]
 
-@fit_and_plot acs data time_points vars=[A] params α maxeval=200 lower_bounds=0 upper_bounds=0.01
+@fit_and_plot acs data time_points vars = [A] params α maxeval = 200 lower_bounds = 0 upper_bounds =
+    0.01
 
 # a slightly more extended example: export solver as a function of given params
 ## some embedded neural network, etc.
-@register begin function learnt_function(A, B, C, params, α)
-    [A, B, C]' * params + α # params: 3-element vector
-end end
+@register begin
+    function learnt_function(A, B, C, params, α)
+        return [A, B, C]' * params + α # params: 3-element vector
+    end
+end
 acs = @ReactionNetwork begin
     learnt_function(A, B, C, params, α), A --> B + C, priority => 0.6
     1.0, B --> C
@@ -37,26 +42,26 @@ acs = @ReactionNetwork begin
 end
 
 # initial values, check params
-@prob_init acs A=60.0 B=10.0 C=150.0
-@prob_params acs params=[1, 2, 3] α=1
-@prob_meta acs tspan=100.0
+@prob_init acs A = 60.0 B = 10.0 C = 150.0
+@prob_params acs params = [1, 2, 3] α = 1
+@prob_meta acs tspan = 100.0
 
 prob = @problematize acs
 sol = @solve prob
-@optimize acs abs(C) params α maxeval=200 lower_bounds=0 upper_bounds=200 final_only=true
+@optimize acs abs(C) params α maxeval = 200 lower_bounds = 0 upper_bounds = 200 final_only =
+    true
 
 ## export solution as a function of params
-parametrized_solver = @build_solver acs A params α trajectories=10
+parametrized_solver = @build_solver acs A params α trajectories = 10
 parametrized_solver = @build_solver acs A params α
 parametrized_solver([1.0, 0.0, 0.0, 0.0, 1.0])
 
 using Statistics # for mean
-function my_little_objective(params)
+my_little_objective(params) =
     mean(1:10) do _ # take avg over 10 samples
         sol = parametrized_solver(params) # compute solution
-        sol[3, end] # some objective
+        return sol[3, end] # some objective
     end
-end
 
 ## now some optimizer...
 using NLopt # https://github.com/JuliaOpt/NLopt.jl
diff --git a/tutorial/rd_example.jl b/tutorial/rd_example.jl
index 5a97598..64638b3 100644
--- a/tutorial/rd_example.jl
+++ b/tutorial/rd_example.jl
@@ -9,22 +9,24 @@ rd_model = @ReactionNetwork
     n_workforce = 5
     n_equipment = 4
     # transition matrix: ones at the upper diagonal + rand noise
-    prob_success = [i / (i + 1) for i in 1:n_phase]
+    prob_success = [i / (i + 1) for i = 1:n_phase]
     transitions = zeros(n_phase, n_phase)
-    foreach(i -> transitions[i, i + 1] = 1, 1:(n_phase - 1))
+    foreach(i -> transitions[i, i+1] = 1, 1:(n_phase-1))
     transitions .+= 0.1 * rand(n_phase, n_phase)
 end
 
-rd_model = @ReactionNetwork begin $i, phase[$i] --> phase[$j], cycle_time => $i * $j end i=1:3 j=1:($i)
+rd_model = @ReactionNetwork begin
+    $i, phase[$i] --> phase[$j], cycle_time => $i * $j
+end i = 1:3 j = 1:($i)
 
 using ReactiveDynamics: nparts
 u0 = rand(1:1000, nparts(rd_model, :S))
 @prob_init rd_model u0
 
-@prob_meta rd_model tspan=100
+@prob_meta rd_model tspan = 100
 
 prob = @problematize rd_model
-sol = @solve prob trajectories=20
+sol = @solve prob trajectories = 20
 @plot sol plot_type
 
 # need to bring data into ReactiveDynamics's (solver) scope
@@ -32,7 +34,7 @@ sol = @solve prob trajectories=20
     k = 5
     r = 3
     M = zeros(k, k)
-    for i in 1:k
+    for i = 1:k
         for conn_in in rand(1:k, rand(1:4))
             M[conn_in, i] += 0.1 * rand(1:10)
         end
@@ -44,17 +46,20 @@ sol = @solve prob trajectories=20
     resource = rand(1:10, k, k, r)
 end
 
-rd_model = @ReactionNetwork begin M[$i, $j],
-                                  mod[$i] + {resource[$i, $j, $k] * resource[$k],
-                                   k = rand(1:(ReactiveDynamics.r)), dlm = +} --> mod[$j],
-                                  cycle_time => cycle_times[$i, $j] end i=1:(ReactiveDynamics.k) j=1:(ReactiveDynamics.k)
+rd_model = @ReactionNetwork begin
+    M[$i, $j],
+    mod[$i] +
+    {resource[$i, $j, $k] * resource[$k], k = rand(1:(ReactiveDynamics.r)), dlm = +} -->
+    mod[$j],
+    cycle_time => cycle_times[$i, $j]
+end i = 1:(ReactiveDynamics.k) j = 1:(ReactiveDynamics.k)
 
 using ReactiveDynamics: nparts
 u0 = rand(1:1000, nparts(rd_model, :S))
 @prob_init rd_model u0
 
-@prob_meta rd_model tspan=100
+@prob_meta rd_model tspan = 100
 
 prob = @problematize rd_model
 sol = @solve prob
-@plot sol plot_type=allocation
+@plot sol plot_type = allocation
diff --git a/tutorial/toy_pharma_model.jl b/tutorial/toy_pharma_model.jl
index 31623c9..4fc6c7e 100644
--- a/tutorial/toy_pharma_model.jl
+++ b/tutorial/toy_pharma_model.jl
@@ -3,79 +3,91 @@ using ReactiveDynamics
 # model dynamics
 toy_pharma_model = @ReactionNetwork begin
     α(candidate_compound, marketed_drug, κ),
-    3 * @conserved(scientist) + @rate(budget) --> candidate_compound, name => discovery,
-    probability => 0.3, cycletime => 10.0, priority => 0.5
+    3 * @conserved(scientist) + @rate(budget) --> candidate_compound,
+    name => discovery,
+    probability => 0.3,
+    cycletime => 10.0,
+    priority => 0.5
     β(candidate_compound, marketed_drug),
     candidate_compound + 5 * @conserved(scientist) + 2 * @rate(budget) -->
-    marketed_drug + 5 * budget, name => dx2market, probability => 0.5 + 0.001 * @t(),
+    marketed_drug + 5 * budget,
+    name => dx2market,
+    probability => 0.5 + 0.001 * @t(),
     cycletime => 4
     γ * marketed_drug, marketed_drug --> ∅, name => drug_killed
 end
 
-@periodic toy_pharma_model 1.0 budget+=11 * marketed_drug
+@periodic toy_pharma_model 1.0 budget += 11 * marketed_drug
 
 @register function α(number_candidate_compounds, number_marketed_drugs, κ)
-    κ + exp(-number_candidate_compounds) + exp(-number_marketed_drugs)
+    return κ + exp(-number_candidate_compounds) + exp(-number_marketed_drugs)
 end
 @register function β(number_candidate_compounds, number_marketed_drugs)
-    number_candidate_compounds + exp(-number_marketed_drugs)
+    return number_candidate_compounds + exp(-number_marketed_drugs)
 end
 
 # simulation parameters
 ## initial values
-@prob_init toy_pharma_model candidate_compound=5 marketed_drug=6 scientist=20 budget=100
+@prob_init toy_pharma_model candidate_compound = 5 marketed_drug = 6 scientist = 20 budget =
+    100
 ## parameters
-@prob_params toy_pharma_model κ=4 γ=0.1
+@prob_params toy_pharma_model κ = 4 γ = 0.1
 ## other arguments passed to the solver
-@prob_meta toy_pharma_model tspan=250 dt=0.1
+@prob_meta toy_pharma_model tspan = 250 dt = 0.1
 
 prob = @problematize toy_pharma_model
 
-sol = @solve prob trajectories=20
+sol = @solve prob trajectories = 20
 
 using Plots
 
-@plot sol plot_type=summary
+@plot sol plot_type = summary
 
-@plot sol plot_type=summary show=:marketed_drug
+@plot sol plot_type = summary show = :marketed_drug
 
 ## for deterministic rates 
 
 # model dynamics
 toy_pharma_model = @ReactionNetwork begin
     @per_step(α(candidate_compound, marketed_drug, κ)),
-    3 * @conserved(scientist) + @rate(budget) --> candidate_compound, name => discovery,
-    probability => 0.3, cycletime => 10.0, priority => 0.5
+    3 * @conserved(scientist) + @rate(budget) --> candidate_compound,
+    name => discovery,
+    probability => 0.3,
+    cycletime => 10.0,
+    priority => 0.5
     @per_step(β(candidate_compound, marketed_drug)),
     candidate_compound + 5 * @conserved(scientist) + 2 * @rate(budget) -->
-    marketed_drug + 5 * budget, name => dx2market, probability => 0.5 + 0.001 * @t(),
+    marketed_drug + 5 * budget,
+    name => dx2market,
+    probability => 0.5 + 0.001 * @t(),
     cycletime => 4
-    @per_step(γ*marketed_drug), marketed_drug --> ∅, name => drug_killed
+    @per_step(γ * marketed_drug), marketed_drug --> ∅, name => drug_killed
 end
 
-@periodic toy_pharma_model 0.0 budget+=11 * marketed_drug
+@periodic toy_pharma_model 0.0 budget += 11 * marketed_drug
 
 @register function α(number_candidate_compounds, number_marketed_drugs, κ)
-    κ + exp(-number_candidate_compounds) + exp(-number_marketed_drugs)
+    return κ + exp(-number_candidate_compounds) + exp(-number_marketed_drugs)
 end
 @register function β(number_candidate_compounds, number_marketed_drugs)
-    number_candidate_compounds + exp(-number_marketed_drugs)
+    return number_candidate_compounds + exp(-number_marketed_drugs)
 end
 
 # simulation parameters
 ## initial values
-@prob_init toy_pharma_model candidate_compound=5 marketed_drug=6 scientist=20 budget=100
+@prob_init toy_pharma_model candidate_compound = 5 marketed_drug = 6 scientist = 20 budget =
+    100
 ## parameters
-@prob_params toy_pharma_model κ=4 γ=0.1
+@prob_params toy_pharma_model κ = 4 γ = 0.1
 ## other arguments passed to the solver
-@prob_meta toy_pharma_model tspan=250
+@prob_meta toy_pharma_model tspan = 250
 
 prob = @problematize toy_pharma_model
 
-sol = @solve prob trajectories=20
+sol = @solve prob trajectories = 20
 
 using Plots
 
-@plot sol plot_type=summary
+@plot sol plot_type = summary
 
-@plot sol plot_type=summary show=:marketed_drug
+@plot sol plot_type = summary show = :marketed_drug