From a164b3fadf4bd59b7a6f28bba5852afc85fb12da Mon Sep 17 00:00:00 2001 From: Claude Date: Thu, 31 Jul 2025 18:54:34 -0400 Subject: [PATCH] Update code formatting with JuliaFormatter MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Applied latest JuliaFormatter standards to improve code consistency and readability across the codebase. 🤖 Generated with [Claude Code](https://claude.ai/code) Co-Authored-By: Claude --- docs/make.jl | 7 +- docs/pages.jl | 159 +++++++++--------- .../bifurcation_kit_extension.jl | 3 +- ...cairo_makie_extension_spatial_modelling.jl | 9 +- ...graph_makie_extension_spatial_modelling.jl | 3 +- .../rn_graph_plot.jl | 124 +++++++------- .../homotopy_continuation_extension.jl | 1 + src/Catalyst.jl | 3 +- src/chemistry_functionality.jl | 4 +- src/dsl.jl | 34 ++-- src/expression_utils.jl | 3 +- src/network_analysis.jl | 46 ++--- src/plotting.jl | 3 +- src/reaction.jl | 7 +- src/reactionsystem.jl | 36 ++-- src/reactionsystem_conversions.jl | 114 +++++++------ .../serialise_reactionsystem.jl | 15 +- .../lattice_jump_systems.jl | 4 +- .../lattice_reaction_systems.jl | 8 +- .../lattice_sim_struct_interfacing.jl | 3 +- .../spatial_ODE_systems.jl | 7 +- .../spatial_reactions.jl | 9 +- src/spatial_reaction_systems/utility.jl | 7 +- 23 files changed, 328 insertions(+), 281 deletions(-) diff --git a/docs/make.jl b/docs/make.jl index f732463793..9d0b7c0837 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -39,9 +39,10 @@ makedocs(sitename = "Catalyst.jl", collapselevel = 1, assets = ["assets/favicon.ico"], canonical = "https://docs.sciml.ai/Catalyst/stable/"), - modules = [Catalyst, ModelingToolkit, - isdefined(Base, :get_extension) ? Base.get_extension(Catalyst, :CatalystGraphMakieExtension) : - Catalyst.CatalystGraphMakieExtension], + modules = [Catalyst, ModelingToolkit, + isdefined(Base, :get_extension) ? + Base.get_extension(Catalyst, :CatalystGraphMakieExtension) : + Catalyst.CatalystGraphMakieExtension], doctest = false, clean = true, pages = pages, diff --git a/docs/pages.jl b/docs/pages.jl index cfe9c4078e..86a437e2b2 100644 --- a/docs/pages.jl +++ b/docs/pages.jl @@ -1,83 +1,80 @@ pages = Any[ - "Home" => "index.md", - "Introduction to Catalyst" => Any[ - "introduction_to_catalyst/catalyst_for_new_julia_users.md", - "introduction_to_catalyst/introduction_to_catalyst.md", - "introduction_to_catalyst/math_models_intro.md" - ], - "Model creation and properties" => Any[ - "model_creation/dsl_basics.md", - "model_creation/dsl_advanced.md", - "model_creation/programmatic_CRN_construction.md", - "model_creation/compositional_modeling.md", - "model_creation/constraint_equations.md", - "model_creation/conservation_laws.md", - "model_creation/parametric_stoichiometry.md", - "model_creation/functional_parameters.md", - "model_creation/model_file_loading_and_export.md", - "model_creation/model_visualisation.md", - "model_creation/reactionsystem_content_accessing.md", - "model_creation/chemistry_related_functionality.md", - "Examples" => Any[ - "model_creation/examples/basic_CRN_library.md", - "model_creation/examples/programmatic_generative_linear_pathway.md", - "model_creation/examples/hodgkin_huxley_equation.md", - "model_creation/examples/smoluchowski_coagulation_equation.md", - "model_creation/examples/noise_modelling_approaches.md", - ] - ], - "Model simulation and visualization" => Any[ - "model_simulation/simulation_introduction.md", - "model_simulation/simulation_plotting.md", - "model_simulation/simulation_structure_interfacing.md", - "model_simulation/ensemble_simulations.md", - "model_simulation/ode_simulation_performance.md", - "model_simulation/sde_simulation_performance.md", - "model_simulation/finite_state_projection_simulation.md", - "Examples" => Any[ - "model_simulation/examples/periodic_events_simulation.md", - "model_simulation/examples/activation_time_distribution_measurement.md", - "model_simulation/examples/interactive_brusselator_simulation.md" - ] - ], - "Network Analysis" => Any[ - "network_analysis/odes.md", - "network_analysis/crn_theory.md", - "network_analysis/network_properties.md" - ], - "Steady state analysis" => Any[ - "steady_state_functionality/homotopy_continuation.md", - "steady_state_functionality/nonlinear_solve.md", - "steady_state_functionality/steady_state_stability_computation.md", - "steady_state_functionality/bifurcation_diagrams.md", - "steady_state_functionality/dynamical_systems.md", - "Examples" => Any[ - "steady_state_functionality/examples/nullcline_plotting.md", - "steady_state_functionality/examples/bifurcationkit_periodic_orbits.md", - "steady_state_functionality/examples/bifurcationkit_codim2.md" - ] - ], - "Inverse problems" => Any[ - "inverse_problems/petab_ode_param_fitting.md", - "inverse_problems/optimization_ode_param_fitting.md", - "inverse_problems/behaviour_optimisation.md", - "inverse_problems/structural_identifiability.md", - "inverse_problems/global_sensitivity_analysis.md", - #"Examples" => Any[ - # "inverse_problems/examples/ode_fitting_oscillation.md" - #] - ], - "Spatial modelling" => Any[ - "spatial_modelling/lattice_reaction_systems.md", - "spatial_modelling/lattice_simulation_structure_ interaction.md", - "spatial_modelling/lattice_simulation_plotting.md", - "spatial_modelling/spatial_ode_simulations.md", - "spatial_modelling/spatial_jump_simulations.md" - ], - "FAQs" => "faqs.md", - "API" => Any[ - "api/core_api.md", - "api/network_analysis_api.md" - ], - "Developer Documentation" => "devdocs/dev_guide.md" +"Home" => "index.md", +"Introduction to Catalyst" => Any[ + "introduction_to_catalyst/catalyst_for_new_julia_users.md", + "introduction_to_catalyst/introduction_to_catalyst.md", + "introduction_to_catalyst/math_models_intro.md" +], +"Model creation and properties" => Any[ + "model_creation/dsl_basics.md", + "model_creation/dsl_advanced.md", + "model_creation/programmatic_CRN_construction.md", + "model_creation/compositional_modeling.md", + "model_creation/constraint_equations.md", + "model_creation/conservation_laws.md", + "model_creation/parametric_stoichiometry.md", + "model_creation/functional_parameters.md", + "model_creation/model_file_loading_and_export.md", + "model_creation/model_visualisation.md", + "model_creation/reactionsystem_content_accessing.md", + "model_creation/chemistry_related_functionality.md", + "Examples" => Any[ + "model_creation/examples/basic_CRN_library.md", + "model_creation/examples/programmatic_generative_linear_pathway.md", + "model_creation/examples/hodgkin_huxley_equation.md", + "model_creation/examples/smoluchowski_coagulation_equation.md", + "model_creation/examples/noise_modelling_approaches.md" + ] +], +"Model simulation and visualization" => Any[ + "model_simulation/simulation_introduction.md", + "model_simulation/simulation_plotting.md", + "model_simulation/simulation_structure_interfacing.md", + "model_simulation/ensemble_simulations.md", + "model_simulation/ode_simulation_performance.md", + "model_simulation/sde_simulation_performance.md", + "model_simulation/finite_state_projection_simulation.md", + "Examples" => Any[ + "model_simulation/examples/periodic_events_simulation.md", + "model_simulation/examples/activation_time_distribution_measurement.md", + "model_simulation/examples/interactive_brusselator_simulation.md" + ] +], +"Network Analysis" => Any[ + "network_analysis/odes.md", + "network_analysis/crn_theory.md", + "network_analysis/network_properties.md" +], +"Steady state analysis" => Any[ + "steady_state_functionality/homotopy_continuation.md", + "steady_state_functionality/nonlinear_solve.md", + "steady_state_functionality/steady_state_stability_computation.md", + "steady_state_functionality/bifurcation_diagrams.md", + "steady_state_functionality/dynamical_systems.md", + "Examples" => Any[ + "steady_state_functionality/examples/nullcline_plotting.md", + "steady_state_functionality/examples/bifurcationkit_periodic_orbits.md", + "steady_state_functionality/examples/bifurcationkit_codim2.md" + ] +], +"Inverse problems" => Any[ + "inverse_problems/petab_ode_param_fitting.md", + "inverse_problems/optimization_ode_param_fitting.md", + "inverse_problems/behaviour_optimisation.md", + "inverse_problems/structural_identifiability.md", + "inverse_problems/global_sensitivity_analysis.md" #"Examples" => Any[ # "inverse_problems/examples/ode_fitting_oscillation.md" #] +], +"Spatial modelling" => Any[ + "spatial_modelling/lattice_reaction_systems.md", + "spatial_modelling/lattice_simulation_structure_ interaction.md", + "spatial_modelling/lattice_simulation_plotting.md", + "spatial_modelling/spatial_ode_simulations.md", + "spatial_modelling/spatial_jump_simulations.md" +], +"FAQs" => "faqs.md", +"API" => Any[ + "api/core_api.md", + "api/network_analysis_api.md" +], +"Developer Documentation" => "devdocs/dev_guide.md" ] diff --git a/ext/CatalystBifurcationKitExtension/bifurcation_kit_extension.jl b/ext/CatalystBifurcationKitExtension/bifurcation_kit_extension.jl index 9b207bc2e8..d88d8d1b3c 100644 --- a/ext/CatalystBifurcationKitExtension/bifurcation_kit_extension.jl +++ b/ext/CatalystBifurcationKitExtension/bifurcation_kit_extension.jl @@ -37,6 +37,7 @@ function bkext_make_nsys(rs, u0) cons_default = [cons_eq.rhs for cons_eq in cons_eqs] cons_default = Catalyst.get_networkproperties(rs).conservedconst => cons_default defaults = Dict([u0; cons_default]) - nsys = convert(NonlinearSystem, rs; defaults, remove_conserved = true, conseqs_remake_warn = false) + nsys = convert( + NonlinearSystem, rs; defaults, remove_conserved = true, conseqs_remake_warn = false) return complete(nsys) end diff --git a/ext/CatalystCairoMakieExtension/cairo_makie_extension_spatial_modelling.jl b/ext/CatalystCairoMakieExtension/cairo_makie_extension_spatial_modelling.jl index 5b719558f8..8bd342c21a 100644 --- a/ext/CatalystCairoMakieExtension/cairo_makie_extension_spatial_modelling.jl +++ b/ext/CatalystCairoMakieExtension/cairo_makie_extension_spatial_modelling.jl @@ -24,7 +24,8 @@ function lattice_animation( vals, plot_min, plot_max = Catalyst.extract_vals(sol, sp, lrs, plot_min, plot_max, t) # Creates the base figure (which is modified in the animation). - fig, ax, plt = scatterlines(vals[1]; + fig, ax, + plt = scatterlines(vals[1]; axis = (xlabel = "Compartment", ylabel = "$(sp)", limits = (nothing, nothing, plot_min, plot_max)), markersize = markersize, kwargs...) @@ -100,8 +101,10 @@ function lattice_animation( x_vals, y_vals = Catalyst.extract_grid_axes(lrs) # Creates the base figure (which is modified in the animation). - fig, ax, hm = heatmap(x_vals, y_vals, vals[1]; - axis = (xgridvisible = false, ygridvisible = false, xlabel = "Compartment", ylabel = "Compartment"), + fig, ax, + hm = heatmap(x_vals, y_vals, vals[1]; + axis = (xgridvisible = false, ygridvisible = false, + xlabel = "Compartment", ylabel = "Compartment"), colormap, colorrange = (plot_min, plot_max), kwargs...) ttitle && (ax.title = "Time: $(round(t[1]; sigdigits = 3))") diff --git a/ext/CatalystGraphMakieExtension/graph_makie_extension_spatial_modelling.jl b/ext/CatalystGraphMakieExtension/graph_makie_extension_spatial_modelling.jl index 5d61a7a55e..cfb6a441ec 100644 --- a/ext/CatalystGraphMakieExtension/graph_makie_extension_spatial_modelling.jl +++ b/ext/CatalystGraphMakieExtension/graph_makie_extension_spatial_modelling.jl @@ -27,7 +27,8 @@ function lattice_animation( vals, plot_min, plot_max = Catalyst.extract_vals(sol, sp, lrs, plot_min, plot_max, t) # Creates the base figure (which is modified in the animation). - fig, ax, plt = graphplot(plot_graph; node_color = vals[1], + fig, ax, + plt = graphplot(plot_graph; node_color = vals[1], node_attr = (colorrange = (plot_min, plot_max), colormap), node_size, kwargs... ) ttitle && (ax.title = "Time: $(round(t[1]; sigdigits = 3))") diff --git a/ext/CatalystGraphMakieExtension/rn_graph_plot.jl b/ext/CatalystGraphMakieExtension/rn_graph_plot.jl index be83f6f8e6..eef6a2cf75 100644 --- a/ext/CatalystGraphMakieExtension/rn_graph_plot.jl +++ b/ext/CatalystGraphMakieExtension/rn_graph_plot.jl @@ -12,29 +12,29 @@ Wrapper intended to allow plotting of multiple edges. This is needed in the foll `gen_distances` sets the distances between the edges that allows multiple to be visible on the plot at the same time. """ struct MultiGraphWrap{T} <: Graphs.AbstractGraph{T} - g::SimpleDiGraph{T} - multiedges::Vector{Graphs.SimpleEdge{T}} - """Sets the drawing order of the edges. Needed because multiedges need to be consecutive to be drawn properly.""" - edgeorder::Vector{Int64} + g::SimpleDiGraph{T} + multiedges::Vector{Graphs.SimpleEdge{T}} + """Sets the drawing order of the edges. Needed because multiedges need to be consecutive to be drawn properly.""" + edgeorder::Vector{Int64} end # Create the SimpleDiGraph corresponding to the species and reactions, the species-reaction graph -function SRGraphWrap(rn::ReactionSystem) +function SRGraphWrap(rn::ReactionSystem) srg = species_reaction_graph(rn) multiedges = Vector{Graphs.SimpleEdge{Int}}() sm = speciesmap(rn) specs = species(rn) deps = Set() - for (i, rx) in enumerate(reactions(rn)) + for (i, rx) in enumerate(reactions(rn)) empty!(deps) get_variables!(deps, rx.rate, specs) if !isempty(deps) - for spec in deps + for spec in deps specidx = sm[spec] - has_edge(srg, specidx, i + length(specs)) ? - push!(multiedges, Graphs.SimpleEdge(specidx, i + length(specs))) : - add_edge!(srg, Graphs.SimpleEdge(specidx, i + length(specs))) + has_edge(srg, specidx, i + length(specs)) ? + push!(multiedges, Graphs.SimpleEdge(specidx, i + length(specs))) : + add_edge!(srg, Graphs.SimpleEdge(specidx, i + length(specs))) end end end @@ -44,16 +44,16 @@ function SRGraphWrap(rn::ReactionSystem) end # Automatically set edge drawing order if not supplied -function MultiGraphWrap(g::SimpleDiGraph{T}, multiedges::Vector{Graphs.SimpleEdge{T}}) where T +function MultiGraphWrap(g::SimpleDiGraph{T}, multiedges::Vector{Graphs.SimpleEdge{T}}) where {T} edgelist = vcat(collect(Graphs.edges(g)), multiedges) edgeorder = sortperm(edgelist) MultiGraphWrap(g, multiedges, edgeorder) end # Return the multigraph and reaction order corresponding to the complex graph. The reaction order is the order of reactions(rn) that would match the edge order given by g.edgeorder. -function ComplexGraphWrap(rn::ReactionSystem) +function ComplexGraphWrap(rn::ReactionSystem) img = incidencematgraph(rn) - D = incidencemat(rn; sparse=true) + D = incidencemat(rn; sparse = true) specs = species(rn) rxs = reactions(rn) @@ -78,7 +78,7 @@ function ComplexGraphWrap(rn::ReactionSystem) edgelist = edgelist[rxorder] multiedges = Vector{Graphs.SimpleEdge{Int}}() for i in 2:length(edgelist) - isequal(edgelist[i], edgelist[i-1]) && push!(multiedges, edgelist[i]) + isequal(edgelist[i], edgelist[i - 1]) && push!(multiedges, edgelist[i]) end MultiGraphWrap(img, multiedges), rxorder end @@ -87,8 +87,8 @@ Base.eltype(g::MultiGraphWrap) = eltype(g.g) Graphs.edgetype(g::MultiGraphWrap) = edgetype(g.g) Graphs.has_edge(g::MultiGraphWrap, s, d) = has_edge(g.g, s, d) Graphs.has_vertex(g::MultiGraphWrap, i) = has_vertex(g.g, i) -Graphs.inneighbors(g::MultiGraphWrap{T}, i) where T = inneighbors(g.g, i) -Graphs.outneighbors(g::MultiGraphWrap{T}, i) where T = outneighbors(g.g, i) +Graphs.inneighbors(g::MultiGraphWrap{T}, i) where {T} = inneighbors(g.g, i) +Graphs.outneighbors(g::MultiGraphWrap{T}, i) where {T} = outneighbors(g.g, i) Graphs.ne(g::MultiGraphWrap) = length(g.multiedges) + length(Graphs.edges(g.g)) Graphs.nv(g::MultiGraphWrap) = nv(g.g) Graphs.vertices(g::MultiGraphWrap) = vertices(g.g) @@ -96,7 +96,7 @@ Graphs.is_directed(::Type{<:MultiGraphWrap}) = true Graphs.is_directed(g::MultiGraphWrap) = is_directed(g.g) Graphs.is_connected(g::MultiGraphWrap) = is_connected(g.g) -function Graphs.adjacency_matrix(g::MultiGraphWrap) +function Graphs.adjacency_matrix(g::MultiGraphWrap) adj = Graphs.adjacency_matrix(g.g) for e in g.multiedges adj[src(e), dst(e)] = 1 @@ -108,12 +108,12 @@ function Graphs.edges(g::MultiGraphWrap) edgelist = vcat(collect(Graphs.edges(g.g)), g.multiedges)[g.edgeorder] end -function gen_distances(g::MultiGraphWrap; inc = 0.2) +function gen_distances(g::MultiGraphWrap; inc = 0.2) edgelist = edges(g) distances = zeros(length(edgelist)) edgedict = Dict(edgelist[1] => [1]) for (i, e) in enumerate(@view edgelist[2:end]) - if edgelist[i] != edgelist[i+1] + if edgelist[i] != edgelist[i + 1] edgedict[e] = [i+1] else push!(edgedict[e], i+1) @@ -122,7 +122,7 @@ function gen_distances(g::MultiGraphWrap; inc = 0.2) for (edge, inds) in edgedict if haskey(edgedict, Edge(dst(edge), src(edge))) - distances[inds[1]] != 0. && continue + distances[inds[1]] != 0.0 && continue inds_ = edgedict[Edge(dst(edge), src(edge))] len = length(inds) + length(inds_) @@ -130,7 +130,7 @@ function gen_distances(g::MultiGraphWrap; inc = 0.2) ep = sp + inc*(len-1) dists = collect(sp:inc:ep) distances[inds] = dists[1:length(inds)] - distances[inds_] = -dists[length(inds)+1:end] + distances[inds_] = -dists[(length(inds) + 1):end] else sp = -inc/2*(length(inds)-1) ep = sp + inc*(length(inds)-1) @@ -158,30 +158,30 @@ Notes: red and black arrows from `A` to the reaction node. For a list of accepted keyword arguments to the graph plot, please see the [GraphMakie documentation](https://graph.makie.org/stable/#The-graphplot-Recipe). -""" +""" function Catalyst.plot_network(rn::ReactionSystem; kwargs...) srg = SRGraphWrap(rn) ns = length(species(rn)) - nodecolors = vcat([:skyblue3 for i in 1:ns], - [:green for i in ns+1:nv(srg)]) - ilabels = vcat(map(s -> String(tosymbol(s, escape=false)), species(rn)), - ["R$i" for i in 1:nv(srg)-ns]) + nodecolors = vcat([:skyblue3 for i in 1:ns], + [:green for i in (ns + 1):nv(srg)]) + ilabels = vcat(map(s -> String(tosymbol(s, escape = false)), species(rn)), + ["R$i" for i in 1:(nv(srg) - ns)]) ssm = substoichmat(rn) psm = prodstoichmat(rn) # Get stoichiometry of reaction edgelabels = map(Graphs.edges(srg.g)) do e - string(src(e) > ns ? - psm[dst(e), src(e)-ns] : - ssm[src(e), dst(e)-ns]) - end + string(src(e) > ns ? + psm[dst(e), src(e) - ns] : + ssm[src(e), dst(e) - ns]) + end edgecolors = [:black for i in 1:ne(srg)] num_e = ne(srg.g) # Handle the rate edges for i in 1:length(srg.edgeorder) # If there are stoichiometry and rate edges from the same species to reaction - if srg.edgeorder[i] > num_e + if srg.edgeorder[i] > num_e edgecolors[i] = :red insert!(edgelabels, i, "") elseif edgelabels[i] == "0" @@ -190,29 +190,29 @@ function Catalyst.plot_network(rn::ReactionSystem; kwargs...) end end - layout = if !haskey(kwargs, :layout) + layout = if !haskey(kwargs, :layout) Stress() end - f = graphplot(srg; - layout, - edge_color = edgecolors, - elabels = edgelabels, - elabels_rotation = 0, - ilabels = ilabels, - node_color = nodecolors, - arrow_shift = :end, - arrow_size = 20, - curve_distance_usage = true, - curve_distance = gen_distances(srg), - kwargs... - ) + f = graphplot(srg; + layout, + edge_color = edgecolors, + elabels = edgelabels, + elabels_rotation = 0, + ilabels = ilabels, + node_color = nodecolors, + arrow_shift = :end, + arrow_size = 20, + curve_distance_usage = true, + curve_distance = gen_distances(srg), + kwargs... + ) f.axis.xautolimitmargin = (0.15, 0.15) f.axis.yautolimitmargin = (0.15, 0.15) hidedecorations!(f.axis) hidespines!(f.axis) f.axis.aspect = DataAspect() - + f end @@ -248,21 +248,21 @@ function Catalyst.plot_complexes(rn::ReactionSystem; show_rate_labels = false, k # Get complex graph and reaction order for edgecolors and edgelabels. rxorder gives the order of reactions(rn) that would match the edge order in edges(cg). cg, rxorder = ComplexGraphWrap(rn) - layout = if !haskey(kwargs, :layout) + layout = if !haskey(kwargs, :layout) Stress() end f = graphplot(cg; - layout, - edge_color = edgecolors[rxorder], - elabels = show_rate_labels ? edgelabels[rxorder] : [], - ilabels = complexlabels(rn), - node_color = :skyblue3, - elabels_rotation = 0, - arrow_shift = :end, - curve_distance_usage = true, - curve_distance = gen_distances(cg), - kwargs... - ) + layout, + edge_color = edgecolors[rxorder], + elabels = show_rate_labels ? edgelabels[rxorder] : [], + ilabels = complexlabels(rn), + node_color = :skyblue3, + elabels_rotation = 0, + arrow_shift = :end, + curve_distance_usage = true, + curve_distance = gen_distances(cg), + kwargs... + ) f.axis.xautolimitmargin = (0.15, 0.15) f.axis.yautolimitmargin = (0.15, 0.15) @@ -273,9 +273,9 @@ function Catalyst.plot_complexes(rn::ReactionSystem; show_rate_labels = false, k f end -function complexelem_tostr(e::Catalyst.ReactionComplexElement, specstrs) +function complexelem_tostr(e::Catalyst.ReactionComplexElement, specstrs) if e.speciesstoich == 1 - return "$(specstrs[e.speciesid])" + return "$(specstrs[e.speciesid])" else return "$(e.speciesstoich)$(specstrs[e.speciesid])" end @@ -285,11 +285,11 @@ end function complexlabels(rn::ReactionSystem) labels = String[] - specstrs = map(s -> String(tosymbol(s, escape=false)), species(rn)) + specstrs = map(s -> String(tosymbol(s, escape = false)), species(rn)) complexes, B = reactioncomplexes(rn) for complex in complexes - if isempty(complex) + if isempty(complex) push!(labels, "∅") elseif length(complex) == 1 push!(labels, complexelem_tostr(complex[1], specstrs)) diff --git a/ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl b/ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl index 2b353d420f..3128b3992d 100644 --- a/ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl +++ b/ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl @@ -166,6 +166,7 @@ end # Filters away solutions with negative species concentrations (and for neg_thres < val < 0.0, sets val=0.0). function filter_negative_f(sols; neg_thres = -1e-15) for sol in sols, idx in 1:length(sol) + (neg_thres < sol[idx] < 0) && (sol[idx] = 0) end return filter(sol -> all(>=(0), sol), sols) diff --git a/src/Catalyst.jl b/src/Catalyst.jl index b1f78cc469..a079816d87 100644 --- a/src/Catalyst.jl +++ b/src/Catalyst.jl @@ -113,7 +113,8 @@ export @reaction_network, @network_component, @reaction, @species # Network analysis functionality. include("network_analysis.jl") export reactioncomplexmap, reactioncomplexes, incidencemat -export complexstoichmat, laplacianmat, fluxmat, massactionvector, complexoutgoingmat, adjacencymat +export complexstoichmat, laplacianmat, fluxmat, massactionvector, complexoutgoingmat, + adjacencymat export incidencematgraph, linkageclasses, stronglinkageclasses, terminallinkageclasses, deficiency, subnetworks export linkagedeficiencies, isreversible, isweaklyreversible diff --git a/src/chemistry_functionality.jl b/src/chemistry_functionality.jl index 4ad17a95e6..ea69c68479 100644 --- a/src/chemistry_functionality.jl +++ b/src/chemistry_functionality.jl @@ -103,7 +103,8 @@ function make_compound(expr) # If no ivs were given, inserts an expression which evaluates to the union of the ivs # for the species the compound depends on. - ivs_get_expr = :(unique(reduce( vcat, (sorted_arguments(ModelingToolkit.unwrap(comp)) + ivs_get_expr = :(unique(reduce( + vcat, (sorted_arguments(ModelingToolkit.unwrap(comp)) for comp in $components)))) if isempty(ivs) species_expr = Catalyst.insert_independent_variable(species_expr, :($ivs_get_expr...)) @@ -184,6 +185,7 @@ function make_compounds(expr) # Next, loops through all 7*[Number of compounds] lines and add them to compound_declarations. compound_calls = [Catalyst.make_compound(line) for line in expr.args] for compound_call in compound_calls, line in striplines(compound_call).args + push!(compound_declarations.args, line) end diff --git a/src/dsl.jl b/src/dsl.jl index f4dbe140fa..d358150448 100644 --- a/src/dsl.jl +++ b/src/dsl.jl @@ -298,12 +298,14 @@ function make_reaction_system(ex::Expr, name) requiredec = haskey(options, :require_declaration) reactions = get_reactions(reaction_lines) sps_inferred, ps_pre_inferred = extract_sps_and_ps(reactions, syms_declared; requiredec) - vs_inferred, diffs_inferred, equations = read_equations_option!(diffsexpr, options, + vs_inferred, diffs_inferred, + equations = read_equations_option!(diffsexpr, options, union(syms_declared, sps_inferred), tiv; requiredec) ps_inferred = setdiff(ps_pre_inferred, vs_inferred, diffs_inferred) syms_inferred = union(sps_inferred, ps_inferred, vs_inferred, diffs_inferred) all_syms = union(syms_declared, syms_inferred) - obsexpr, obs_eqs, obs_syms = read_observables_option(options, ivs, + obsexpr, obs_eqs, + obs_syms = read_observables_option(options, ivs, union(sps_declared, vs_declared), all_syms; requiredec) # Read options not related to the declaration or inference of symbols. @@ -368,7 +370,7 @@ function get_reactions(exprs::Vector{Expr}) # Currently, reaction bundling where rates (but neither substrates nor products) are # bundled, is disabled. See discussion in https://github.com/SciML/Catalyst.jl/issues/1219. if !in(arrow, double_arrows) && Meta.isexpr(rate, :tuple) && - !Meta.isexpr(reaction.args[2], :tuple) && !Meta.isexpr(reaction.args[3], :tuple) + !Meta.isexpr(reaction.args[2], :tuple) && !Meta.isexpr(reaction.args[3], :tuple) error("Bundling of reactions with multiple rates but singular substrates and product sets is disallowed. This error is potentially due to a bidirectional (`<-->`) reaction being incorrectly typed as `-->`.") end @@ -650,7 +652,7 @@ function read_compounds_option(options) cmpexpr_init = options[:compounds] cmpexpr_init.args[3] = option_block_form(get_block_option(cmpexpr_init)) cmps_declared = [find_varinfo_in_declaration(arg.args[2])[1] - for arg in cmpexpr_init.args[3].args] + for arg in cmpexpr_init.args[3].args] (length(cmps_declared) == 1) && (cmpexpr_init.args[1] = Symbol("@compound")) else # If option is not used, return empty vectors and expressions. cmpexpr_init = :() @@ -666,7 +668,7 @@ function read_differentials_option(options) # If differentials were provided as options, this is used as the initial expression. # If the default differential (D(...)) was used in equations, this is added to the expression. diffsexpr = (haskey(options, :differentials) ? - get_block_option(options[:differentials]) : striplines(:(begin end))) + get_block_option(options[:differentials]) : striplines(:(begin end))) diffsexpr = option_block_form(diffsexpr) # Goes through all differentials, checking that they are correctly formatted. Adds their @@ -688,10 +690,12 @@ end # Reads the variables options. Outputs a list of the variables inferred from the equations, # as well as the equation vector. If the default differential was used, update the `diffsexpr` # expression so that this declares this as well. -function read_equations_option!(diffsexpr, options, syms_unavailable, tiv; requiredec = false) +function read_equations_option!( + diffsexpr, options, syms_unavailable, tiv; requiredec = false) # Prepares the equations. First, extract equations from the provided option (converting to block form if required). # Next, uses MTK's `parse_equations!` function to split input into a vector with the equations. - eqs_input = haskey(options, :equations) ? get_block_option(options[:equations]) : :(begin end) + eqs_input = haskey(options, :equations) ? get_block_option(options[:equations]) : + :(begin end) eqs_input = option_block_form(eqs_input) equations = Expr[] ModelingToolkit.parse_equations!(Expr(:block), equations, @@ -746,7 +750,8 @@ end # Reads the observables options. Outputs an expression for creating the observable variables, # a vector containing the observable equations, and a list of all observable symbols (this # list contains both those declared separately or inferred from the `@observables` option` input`). -function read_observables_option(options, all_ivs, us_declared, all_syms; requiredec = false) +function read_observables_option( + options, all_ivs, us_declared, all_syms; requiredec = false) syms_unavailable = setdiff(all_syms, us_declared) if haskey(options, :observables) # Gets list of observable equations and prepares variable declaration expression. @@ -757,7 +762,8 @@ function read_observables_option(options, all_ivs, us_declared, all_syms; requir for (idx, obs_eq) in enumerate(obs_eqs.args) # Extract the observable, checks for errors. - obs_name, ivs, _, defaults, metadata = find_varinfo_in_declaration(obs_eq.args[2]) + obs_name, ivs, _, defaults, + metadata = find_varinfo_in_declaration(obs_eq.args[2]) # Error checks. (requiredec && !in(obs_name, us_declared)) && @@ -772,7 +778,8 @@ function read_observables_option(options, all_ivs, us_declared, all_syms; requir error("An observable ($obs_name) uses a name that already have been already been declared or inferred as another model property.") (obs_name in us_declared) && is_escaped_expr(obs_eq.args[2]) && error("An interpolated observable have been used, which has also been explicitly declared within the system using either @species or @variables. This is not permitted.") - ((obs_name in us_declared) || is_escaped_expr(obs_eq.args[2])) && !isnothing(metadata) && + ((obs_name in us_declared) || is_escaped_expr(obs_eq.args[2])) && + !isnothing(metadata) && error("Metadata was provided to observable $obs_name in the `@observables` macro. However, the observable was also declared separately (using either @species or @variables). When this is done, metadata should instead be provided within the original @species or @variable declaration.") # This bit adds the observables to the @variables vector which is given as output. @@ -879,7 +886,7 @@ end # be used or not. If not provided, use the default (true). function read_combinatoric_ratelaws_option(options) return haskey(options, :combinatoric_ratelaws) ? - get_block_option(options[:combinatoric_ratelaws]) : true + get_block_option(options[:combinatoric_ratelaws]) : true end ### `@reaction` Macro & its Internals ### @@ -984,8 +991,9 @@ function recursive_escape_functions!(expr::ExprValues, syms_skip = []) (typeof(expr) != Expr) && (return expr) foreach(i -> expr.args[i] = recursive_escape_functions!(expr.args[i], syms_skip), 1:length(expr.args)) - if (expr.head == :call) && (expr.args[1] isa Symbol) &&!isdefined(Catalyst, expr.args[1]) && - expr.args[1] ∉ syms_skip + if (expr.head == :call) && (expr.args[1] isa Symbol) && + !isdefined(Catalyst, expr.args[1]) && + expr.args[1] ∉ syms_skip expr.args[1] = esc(expr.args[1]) end expr diff --git a/src/expression_utils.jl b/src/expression_utils.jl index 055e337137..fd997b18f3 100644 --- a/src/expression_utils.jl +++ b/src/expression_utils.jl @@ -30,7 +30,6 @@ function forbidden_symbol_check(syms) error("The following symbol(s) are used as species or parameters: $used_forbidden_syms, this is not permitted.") end - ### Catalyst-specific Expressions Manipulation ### # Many option inputs can be on a form `@option input` or `@option begin ... end`. In both these @@ -93,7 +92,7 @@ function find_varinfo_in_declaration(expr::ExprValues) Meta.isexpr(expr, :escape) && (expr = expr.args[1]) (expr isa Symbol) || error("Erroneous expression encountered in `find_varinfo_in_declaration` (got `$expr` after processing, this should be a symbol).") - return (;sym = expr, ivs, idxs, default, metadata) + return (; sym = expr, ivs, idxs, default, metadata) end # Converts an expression of the forms: diff --git a/src/network_analysis.jl b/src/network_analysis.jl index ee78f7eced..ff2c74f442 100644 --- a/src/network_analysis.jl +++ b/src/network_analysis.jl @@ -96,7 +96,8 @@ function reactioncomplexes(rn::ReactionSystem; sparse = false) nps = get_networkproperties(rn) if isempty(nps.complexes) || (sparse != issparse(nps.complexes)) complextorxsmap = reactioncomplexmap(rn) - nps.complexes, nps.incidencemat = if sparse + nps.complexes, + nps.incidencemat = if sparse reactioncomplexes(SparseMatrixCSC{Int, Int}, rn, complextorxsmap) else reactioncomplexes(Matrix{Int}, rn, complextorxsmap) @@ -215,12 +216,13 @@ Return an r×c matrix ``K`` such that, if complex ``j`` is the substrate complex **Warning**: Unlike other Catalyst functions, the `fluxmat` function will return a `Matrix{Num}` in the symbolic case. This is to allow easier computation of the matrix decomposition of the ODEs, and to ensure that multiplying the sparse form of the matrix will work. """ -function fluxmat(rn::ReactionSystem, pmap::Dict = Dict(); sparse=false) +function fluxmat(rn::ReactionSystem, pmap::Dict = Dict(); sparse = false) deps = Set() for (i, rx) in enumerate(reactions(rn)) empty!(deps) get_variables!(deps, rx.rate, species(rn)) - (!isempty(deps)) && (error("Reaction $rx's rate constant depends on species $(join(deps, ", ")). `fluxmat` cannot support rate constants of this form.")) + (!isempty(deps)) && + (error("Reaction $rx's rate constant depends on species $(join(deps, ", ")). `fluxmat` cannot support rate constants of this form.")) end rates = if isempty(pmap) @@ -238,7 +240,7 @@ function fluxmat(rn::ReactionSystem, pmap::Dict = Dict(); sparse=false) end end -function fluxmat(::Type{SparseMatrixCSC{T, Int}}, rcmap, rates) where T +function fluxmat(::Type{SparseMatrixCSC{T, Int}}, rcmap, rates) where {T} Is = Int[] Js = Int[] Vs = T[] @@ -254,7 +256,7 @@ function fluxmat(::Type{SparseMatrixCSC{T, Int}}, rcmap, rates) where T Z = sparse(Is, Js, Vs, length(rates), length(rcmap)) end -function fluxmat(::Type{Matrix{T}}, rcmap, rates) where T +function fluxmat(::Type{Matrix{T}}, rcmap, rates) where {T} nr = length(rates) nc = length(rcmap) K = zeros(T, nr, nc) @@ -278,7 +280,8 @@ end # Helper to substitute values into a (vector of) symbolic expressions. The syms are the symbols to substitute and the symexprs are the expressions to substitute into. function substitutevals(rn::ReactionSystem, map::Dict, syms, symexprs) - length(map) != length(syms) && error("Incorrect number of parameter-value pairs were specified.") + length(map) != length(syms) && + error("Incorrect number of parameter-value pairs were specified.") map = symmap_to_varmap(rn, map) map = Dict(ModelingToolkit.value(k) => v for (k, v) in map) vals = [substitute(expr, map) for expr in symexprs] @@ -295,13 +298,15 @@ If the `combinatoric_ratelaws` option is set, will include prefactors for that ( **Warning**: Unlike other Catalyst functions, the `massactionvector` function will return a `Vector{Num}` in the symbolic case. This is to allow easier computation of the matrix decomposition of the ODEs. """ -function massactionvector(rn::ReactionSystem, scmap::Dict = Dict(); combinatoric_ratelaws = Catalyst.get_combinatoric_ratelaws(rn)) +function massactionvector(rn::ReactionSystem, scmap::Dict = Dict(); + combinatoric_ratelaws = Catalyst.get_combinatoric_ratelaws(rn)) r = numreactions(rn) rxs = reactions(rn) sm = speciesmap(rn) for rx in rxs - !ismassaction(rx, rn) && error("The supplied ReactionSystem has non-mass action reaction $rx. The `massactionvector` can only be constructed for mass action networks.") + !ismassaction(rx, rn) && + error("The supplied ReactionSystem has non-mass action reaction $rx. The `massactionvector` can only be constructed for mass action networks.") end specs = if isempty(scmap) @@ -328,12 +333,14 @@ function massactionvector(rn::ReactionSystem, scmap::Dict = Dict(); combinatoric Φ end -function massactionvector(rn::ReactionSystem, scmap::Tuple; combinatoric_ratelaws = Catalyst.get_combinatoric_ratelaws(rn)) +function massactionvector(rn::ReactionSystem, scmap::Tuple; + combinatoric_ratelaws = Catalyst.get_combinatoric_ratelaws(rn)) sdict = Dict(scmap) massactionvector(rn, sdict; combinatoric_ratelaws) end -function massactionvector(rn::ReactionSystem, scmap::Vector; combinatoric_ratelaws = Catalyst.get_combinatoric_ratelaws(rn)) +function massactionvector(rn::ReactionSystem, scmap::Vector; + combinatoric_ratelaws = Catalyst.get_combinatoric_ratelaws(rn)) sdict = Dict(scmap) massactionvector(rn, sdict; combinatoric_ratelaws) end @@ -460,7 +467,6 @@ function incidencematgraph(incidencemat::SparseMatrixCSC{Int, Int}) return graph end - """ species_reaction_graph(rn::ReactionSystem) @@ -856,7 +862,8 @@ function cache_conservationlaw_eqs!(rn::ReactionSystem, N::AbstractMatrix, col_o depspecs = sts[depidxs] missingvec = [missing for _ in 1:nullity] constants = MT.unwrap(only( - @parameters $(CONSERVED_CONSTANT_SYMBOL)[1:nullity] = missing [conserved = true, guess = ones(nullity)])) + @parameters $(CONSERVED_CONSTANT_SYMBOL)[1:nullity] = missing [ + conserved = true, guess = ones(nullity)])) conservedeqs = Equation[] constantdefs = Equation[] @@ -929,7 +936,7 @@ end Constructively compute whether a kinetic system (a reaction network with a set of rate constants) will admit detailed-balanced equilibrium solutions, using the Wegscheider conditions, [Feinberg, 1989](https://www.sciencedirect.com/science/article/pii/0009250989851243). A detailed-balanced solution is one for which the rate of every forward reaction exactly equals its reverse reaction. Accepts a dictionary, vector, or tuple of variable-to-value mappings, e.g. [k1 => 1.0, k2 => 2.0,...]. """ -function isdetailedbalanced(rs::ReactionSystem, parametermap::Dict; abstol=0, reltol=1e-9) +function isdetailedbalanced(rs::ReactionSystem, parametermap::Dict; abstol = 0, reltol = 1e-9) if length(parametermap) != numparams(rs) error("Incorrect number of parameters specified.") elseif !isreversible(rs) @@ -982,7 +989,7 @@ function isdetailedbalanced(rs::ReactionSystem, parametermap::Dict; abstol=0, re end # Helper to find the index of the reaction with a given reactant and product complex. -function edgeindex(imat, src::T, dst::T) where T <: Int +function edgeindex(imat, src::T, dst::T) where {T <: Int} for i in 1:size(imat, 2) (imat[src, i] == -1) && (imat[dst, i] == 1) && return i end @@ -1072,7 +1079,8 @@ function adjacencymat(rn::ReactionSystem, pmap::Dict = Dict(); sparse = false) for (i, rx) in enumerate(reactions(rn)) empty!(deps) get_variables!(deps, rx.rate, species(rn)) - (!isempty(deps)) && (error("Reaction $rx's rate constant depends on species $(join(deps, ", ")). `adjacencymat` cannot support rate constants of this form.")) + (!isempty(deps)) && + (error("Reaction $rx's rate constant depends on species $(join(deps, ", ")). `adjacencymat` cannot support rate constants of this form.")) end rates = if isempty(pmap) @@ -1089,7 +1097,7 @@ function adjacencymat(rn::ReactionSystem, pmap::Dict = Dict(); sparse = false) end end -function adjacencymat(::Type{SparseMatrixCSC{T, Int}}, D, rates) where T +function adjacencymat(::Type{SparseMatrixCSC{T, Int}}, D, rates) where {T} Is = Int[] Js = Int[] Vs = T[] @@ -1105,7 +1113,7 @@ function adjacencymat(::Type{SparseMatrixCSC{T, Int}}, D, rates) where T A = sparse(Is, Js, Vs, nc, nc) end -function adjacencymat(::Type{Matrix{T}}, D, rates) where T +function adjacencymat(::Type{Matrix{T}}, D, rates) where {T} nc = size(D, 1) A = zeros(T, nc, nc) @@ -1117,12 +1125,12 @@ function adjacencymat(::Type{Matrix{T}}, D, rates) where T A end -function adjacencymat(rn::ReactionSystem, pmap::Vector{<:Pair}; sparse=false) +function adjacencymat(rn::ReactionSystem, pmap::Vector{<:Pair}; sparse = false) pdict = Dict(pmap) adjacencymat(rn, pdict; sparse) end -function adjacencymat(rn::ReactionSystem, pmap::Tuple; sparse=false) +function adjacencymat(rn::ReactionSystem, pmap::Tuple; sparse = false) pdict = Dict(pmap) adjacencymat(rn, pdict; sparse) end diff --git a/src/plotting.jl b/src/plotting.jl index 331a6b4f75..585b892a67 100644 --- a/src/plotting.jl +++ b/src/plotting.jl @@ -58,7 +58,8 @@ end for (_, mime) in _showables MIMEType = typeof(MIME(mime)) - @eval Base.show(io::IO, ::$MIMEType, s::Showable{>:$MIMEType}; options...) = show( + @eval Base.show(io::IO, ::$MIMEType, s::Showable{>:$MIMEType}; + options...) = show( io, $MIMEType(), s.content; s.options..., options...) end diff --git a/src/reaction.jl b/src/reaction.jl index 2d0346230e..74cd2cf444 100644 --- a/src/reaction.jl +++ b/src/reaction.jl @@ -348,8 +348,7 @@ MT.is_alg_equation(rx::Reaction) = false # MTK functions for extracting variables within equation type object MT.eqtype_supports_collect_vars(rx::Reaction) = true function MT.collect_vars!(unknowns, parameters, rx::Reaction, iv; depth = 0, - op = MT.Differential) - + op = MT.Differential) MT.collect_vars!(unknowns, parameters, rx.rate, iv; depth, op) for items in (rx.substrates, rx.products, rx.substoich, rx.prodstoich) @@ -394,6 +393,7 @@ function ModelingToolkit.get_variables!(set, rx::Reaction) foreach(sub -> push!(set, sub), rx.substrates) foreach(prod -> push!(set, prod), rx.products) for stoichs in (rx.substoich, rx.prodstoich), stoich in stoichs + (stoich isa BasicSymbolic) && get_variables!(set, stoich) end if hasnoisescaling(rx) @@ -702,7 +702,8 @@ Notes: The following values are possible: end const JUMP_SCALES = (PhysicalScale.Jump, PhysicalScale.VariableRateJump) -const NON_CONSTANT_JUMP_SCALES = (PhysicalScale.ODE, PhysicalScale.SDE, PhysicalScale.VariableRateJump) +const NON_CONSTANT_JUMP_SCALES = ( + PhysicalScale.ODE, PhysicalScale.SDE, PhysicalScale.VariableRateJump) """ has_physical_scale(rx::Reaction) diff --git a/src/reactionsystem.jl b/src/reactionsystem.jl index 6bd02f54c1..ad7f107961 100644 --- a/src/reactionsystem.jl +++ b/src/reactionsystem.jl @@ -608,7 +608,8 @@ function debug_comparer(fun, prop1, prop2, propname; debug = false) if fun(prop1, prop2) return true else - debug && println("Comparison was false for property: ", propname, "\n Found: ", prop1, " vs ", prop2) + debug && println("Comparison was false for property: ", propname, + "\n Found: ", prop1, " vs ", prop2) return false end end @@ -637,28 +638,38 @@ function isequivalent(rn1::ReactionSystem, rn2::ReactionSystem; ignorenames = tr end debug_comparer(==, get_combinatoric_ratelaws(rn1), get_combinatoric_ratelaws(rn2), "combinatoric_ratelaws"; debug) || return false - debug_comparer(==, MT.iscomplete(rn1), MT.iscomplete(rn2), "complete"; debug) || return false + debug_comparer(==, MT.iscomplete(rn1), MT.iscomplete(rn2), "complete"; debug) || + return false # symbolic variables and parameters debug_comparer(isequal, get_iv(rn1), get_iv(rn2), "ivs"; debug) || return false debug_comparer(issetequal, get_sivs(rn1), get_sivs(rn2), "sivs"; debug) || return false - debug_comparer(issetequal, get_unknowns(rn1), get_unknowns(rn2), "unknowns"; debug) || return false - debug_comparer(issetequal, get_species(rn1), get_species(rn2), "species"; debug) || return false + debug_comparer(issetequal, get_unknowns(rn1), get_unknowns(rn2), "unknowns"; debug) || + return false + debug_comparer(issetequal, get_species(rn1), get_species(rn2), "species"; debug) || + return false debug_comparer(issetequal, get_ps(rn1), get_ps(rn2), "ps"; debug) || return false - debug_comparer(issetequal, MT.get_defaults(rn1), MT.get_defaults(rn2), "defaults"; debug) || return false + debug_comparer( + issetequal, MT.get_defaults(rn1), MT.get_defaults(rn2), "defaults"; debug) || + return false # equations and reactions - debug_comparer(issetequal, MT.get_observed(rn1), MT.get_observed(rn2), "observed"; debug) || return false + debug_comparer( + issetequal, MT.get_observed(rn1), MT.get_observed(rn2), "observed"; debug) || + return false debug_comparer(issetequal, get_eqs(rn1), get_eqs(rn2), "eqs"; debug) || return false - debug_comparer(issetequal, MT.get_continuous_events(rn1), MT.get_continuous_events(rn2), "cevents"; debug) || return false - debug_comparer(issetequal, MT.get_discrete_events(rn1), MT.get_discrete_events(rn2), "devents"; debug) || return false + debug_comparer(issetequal, MT.get_continuous_events(rn1), + MT.get_continuous_events(rn2), "cevents"; debug) || return false + debug_comparer(issetequal, MT.get_discrete_events(rn1), + MT.get_discrete_events(rn2), "devents"; debug) || return false # coupled systems if (length(get_systems(rn1)) != length(get_systems(rn2))) println("Systems have different numbers of subsystems.") return false end - debug_comparer(issetequal, get_systems(rn1), get_systems(rn2), "systems"; debug) || return false + debug_comparer(issetequal, get_systems(rn1), get_systems(rn2), "systems"; debug) || + return false true end @@ -1316,7 +1327,7 @@ function set_default_metadata(rs::ReactionSystem; default_reaction_metadata = [] ns_sps = Iterators.filter(Catalyst.isspecies, ns_syms) ns_vs = Iterators.filter( sym -> !Catalyst.isspecies(sym) && - !ModelingToolkit.isparameter(sym), ns_syms) + !ModelingToolkit.isparameter(sym), ns_syms) # Adds parameters, species, and variables to the `ReactionSystem`. @set! rs.ps = union(get_ps(rs), ns_ps) sps_new = union(get_species(rs), ns_sps) @@ -1428,7 +1439,7 @@ function MT.flatten(rs::ReactionSystem; name = nameof(rs)) end function complete_check(sys, method) - if MT.iscomplete(sys) + if MT.iscomplete(sys) error("$method with one or more `ReactionSystem`s requires systems to not be marked complete, but system: $(MT.get_name(sys)) is marked complete.") end nothing @@ -1486,10 +1497,9 @@ Notes: """ function ModelingToolkit.extend(sys::MT.AbstractSystem, rs::ReactionSystem; name::Symbol = nameof(sys)) - complete_check(sys, "ModelingToolkit.extend") complete_check(rs, "ModelingToolkit.extend") - + any(T -> sys isa T, (ReactionSystem, ODESystem, NonlinearSystem)) || error("ReactionSystems can only be extended with ReactionSystems, ODESystems and NonlinearSystems currently. Received a $(typeof(sys)) system.") diff --git a/src/reactionsystem_conversions.jl b/src/reactionsystem_conversions.jl index f55a6770d9..84d6060391 100644 --- a/src/reactionsystem_conversions.jl +++ b/src/reactionsystem_conversions.jl @@ -47,7 +47,7 @@ end # including non-species variables. drop_dynamics(s) = isconstant(s) || isbc(s) || (!isspecies(s)) -function assemble_oderhs(rs, ispcs; combinatoric_ratelaws = true, remove_conserved = false, +function assemble_oderhs(rs, ispcs; combinatoric_ratelaws = true, remove_conserved = false, physical_scales = nothing, expand_catalyst_funs = true) nps = get_networkproperties(rs) species_to_idx = Dict(x => i for (i, x) in enumerate(ispcs)) @@ -58,12 +58,12 @@ function assemble_oderhs(rs, ispcs; combinatoric_ratelaws = true, remove_conserv Dict() end - for (rxidx,rx) in enumerate(get_rxs(rs)) + for (rxidx, rx) in enumerate(get_rxs(rs)) # check this reaction should be treated as an ODE - !((physical_scales === nothing) || - (physical_scales[rxidx] == PhysicalScale.ODE)) && continue + !((physical_scales === nothing) || + (physical_scales[rxidx] == PhysicalScale.ODE)) && continue - rl = oderatelaw(rx; combinatoric_ratelaw = combinatoric_ratelaws, + rl = oderatelaw(rx; combinatoric_ratelaw = combinatoric_ratelaws, expand_catalyst_funs) remove_conserved && (rl = substitute(rl, depspec_submap)) for (spec, stoich) in rx.netstoich @@ -99,7 +99,7 @@ end function assemble_drift(rs, ispcs; combinatoric_ratelaws = true, as_odes = true, include_zero_odes = true, remove_conserved = false, physical_scales = nothing, expand_catalyst_funs = true) - rhsvec = assemble_oderhs(rs, ispcs; combinatoric_ratelaws, remove_conserved, + rhsvec = assemble_oderhs(rs, ispcs; combinatoric_ratelaws, remove_conserved, physical_scales, expand_catalyst_funs) if as_odes D = Differential(get_iv(rs)) @@ -129,7 +129,7 @@ function assemble_diffusion(rs, sts, ispcs; combinatoric_ratelaws = true, end for (j, rx) in enumerate(get_rxs(rs)) - rl = oderatelaw(rx; combinatoric_ratelaw = combinatoric_ratelaws, + rl = oderatelaw(rx; combinatoric_ratelaw = combinatoric_ratelaws, expand_catalyst_funs) rlsqrt = sqrt(abs(rl)) hasnoisescaling(rx) && (rlsqrt *= getnoisescaling(rx)) @@ -182,7 +182,7 @@ Notes: """ function jumpratelaw(rx; combinatoric_ratelaw = true, expand_catalyst_funs = true) @unpack rate, substrates, substoich, only_use_rate = rx - + rl = rate expand_catalyst_funs && (rl = expand_registered_functions(rl)) @@ -368,7 +368,7 @@ function classify_vrjs(rs, physcales) isvrjvec end -function assemble_jumps(rs; combinatoric_ratelaws = true, physical_scales = nothing, +function assemble_jumps(rs; combinatoric_ratelaws = true, physical_scales = nothing, expand_catalyst_funs = true, save_positions = (true, true)) meqs = MassActionJump[] ceqs = ConstantRateJump[] @@ -376,7 +376,7 @@ function assemble_jumps(rs; combinatoric_ratelaws = true, physical_scales = noth unknownset = Set(get_unknowns(rs)) rxs = get_rxs(rs) - if physical_scales === nothing + if physical_scales === nothing physcales = [PhysicalScale.Jump for _ in enumerate(rxs)] else physcales = physical_scales @@ -401,14 +401,14 @@ function assemble_jumps(rs; combinatoric_ratelaws = true, physical_scales = noth if (!isvrj) && ismassaction(rx, rs; rxvars, haveivdep = false, unknownset) push!(meqs, makemajump(rx; combinatoric_ratelaw = combinatoric_ratelaws)) else - rl = jumpratelaw(rx; combinatoric_ratelaw = combinatoric_ratelaws, + rl = jumpratelaw(rx; combinatoric_ratelaw = combinatoric_ratelaws, expand_catalyst_funs) affect = Vector{Equation}() for (spec, stoich) in rx.netstoich # don't change species that are constant or BCs (!drop_dynamics(spec)) && push!(affect, spec ~ spec + stoich) end - if isvrj + if isvrj push!(veqs, VariableRateJump(rl, affect; save_positions)) else push!(ceqs, ConstantRateJump(rl, affect)) @@ -431,7 +431,7 @@ function addconstraints!(eqs, rs::ReactionSystem, ists, ispcs; remove_conserved initeqs = Equation[] defs = MT.defaults(rs) obs = MT.observed(rs) - + # make dependent species observables and add conservation constants as parameters if remove_conserved nps = get_networkproperties(rs) @@ -442,7 +442,7 @@ function addconstraints!(eqs, rs::ReactionSystem, ists, ispcs; remove_conserved if treat_conserved_as_eqs # add back previously removed dependent species - sts = union(sts, nps.depspecs) + sts = union(sts, nps.depspecs) # treat conserved eqs as normal eqs append!(eqs, conservedequations(rs)) @@ -450,7 +450,7 @@ function addconstraints!(eqs, rs::ReactionSystem, ists, ispcs; remove_conserved # add initialization equations for conserved parameters initialmap = Dict(u => Initial(u) for u in species(rs)) conseqs = conservationlaw_constants(rs) - initeqs = [Symbolics.substitute(conseq, initialmap) for conseq in conseqs] + initeqs = [Symbolics.substitute(conseq, initialmap) for conseq in conseqs] else # add the dependent species as observed obs = copy(obs) @@ -476,7 +476,6 @@ function addconstraints!(eqs, rs::ReactionSystem, ists, ispcs; remove_conserved eqs, sts, ps, obs, defs, initeqs end - # used by flattened systems that don't support constraint equations currently function error_if_constraints(::Type{T}, sys::ReactionSystem) where {T <: MT.AbstractSystem} any(eq -> eq isa Equation, get_eqs(sys)) && @@ -504,7 +503,6 @@ diff_2_zero(expr) = (Symbolics.is_derivative(expr) ? 0 : expr) COMPLETENESS_ERROR = "A ReactionSystem must be complete before it can be converted to other system types. A ReactionSystem can be marked as complete using the `complete` function." - ### System Conversions ### """ @@ -528,9 +526,9 @@ Keyword args and default values: """ function Base.convert(::Type{<:ODESystem}, rs::ReactionSystem; name = nameof(rs), combinatoric_ratelaws = get_combinatoric_ratelaws(rs), - include_zero_odes = true, remove_conserved = false, checks = false, - default_u0 = Dict(), default_p = Dict(), - defaults = _merge(Dict(default_u0), Dict(default_p)), expand_catalyst_funs = true, + include_zero_odes = true, remove_conserved = false, checks = false, + default_u0 = Dict(), default_p = Dict(), + defaults = _merge(Dict(default_u0), Dict(default_p)), expand_catalyst_funs = true, kwargs...) # Error checks. iscomplete(rs) || error(COMPLETENESS_ERROR) @@ -553,7 +551,6 @@ function Base.convert(::Type{<:ODESystem}, rs::ReactionSystem; name = nameof(rs) kwargs...) end - const NONLIN_PROB_REMAKE_WARNING = """ Note, when constructing `NonlinearSystem`s with `remove_conserved = true`, possibly via \ calling `NonlinearProblem`, `remake(::NonlinearProblem)` has some \ @@ -563,7 +560,7 @@ const NONLIN_PROB_REMAKE_WARNING = """ entry](https://docs.sciml.ai/Catalyst/stable/faqs/#faq_remake_nonlinprob) for more \ details. This warning can be disabled by passing `conseqs_remake_warn = false`.""" -function is_autonomous_error(iv) +function is_autonomous_error(iv) return """ Attempting to convert a non-autonomous `ReactionSystem` (e.g. where some rate depends \ on $(iv)) to a `NonlinearSystem`. This is not possible. if you are intending to \ @@ -598,7 +595,7 @@ Keyword args and default values: function Base.convert(::Type{<:NonlinearSystem}, rs::ReactionSystem; name = nameof(rs), combinatoric_ratelaws = get_combinatoric_ratelaws(rs), remove_conserved = false, conseqs_remake_warn = true, checks = false, - default_u0 = Dict(), default_p = Dict(), + default_u0 = Dict(), default_p = Dict(), defaults = _merge(Dict(default_u0), Dict(default_p)), all_differentials_permitted = false, expand_catalyst_funs = true, kwargs...) # Error checks. @@ -613,7 +610,11 @@ function Base.convert(::Type{<:NonlinearSystem}, rs::ReactionSystem; name = name ists, ispcs = get_indep_sts(fullrs, remove_conserved) eqs = assemble_drift(fullrs, ispcs; combinatoric_ratelaws, remove_conserved, as_odes = false, include_zero_odes = false, expand_catalyst_funs) - eqs, us, ps, obs, defs, initeqs = addconstraints!(eqs, fullrs, ists, ispcs; + eqs, us, + ps, + obs, + defs, + initeqs = addconstraints!(eqs, fullrs, ists, ispcs; remove_conserved, treat_conserved_as_eqs = true) # Throws a warning if there are differential equations in non-standard format. @@ -681,9 +682,9 @@ Notes: function Base.convert(::Type{<:SDESystem}, rs::ReactionSystem; name = nameof(rs), combinatoric_ratelaws = get_combinatoric_ratelaws(rs), include_zero_odes = true, checks = false, remove_conserved = false, - default_u0 = Dict(), default_p = Dict(), + default_u0 = Dict(), default_p = Dict(), defaults = _merge(Dict(default_u0), Dict(default_p)), - expand_catalyst_funs = true, + expand_catalyst_funs = true, kwargs...) # Error checks. iscomplete(rs) || error(COMPLETENESS_ERROR) @@ -695,7 +696,7 @@ function Base.convert(::Type{<:SDESystem}, rs::ReactionSystem; ists, ispcs = get_indep_sts(flatrs, remove_conserved) eqs = assemble_drift(flatrs, ispcs; combinatoric_ratelaws, include_zero_odes, remove_conserved, expand_catalyst_funs) - noiseeqs = assemble_diffusion(flatrs, ists, ispcs; combinatoric_ratelaws, + noiseeqs = assemble_diffusion(flatrs, ists, ispcs; combinatoric_ratelaws, remove_conserved, expand_catalyst_funs) eqs, us, ps, obs, defs = addconstraints!(eqs, flatrs, ists, ispcs; remove_conserved) @@ -740,7 +741,7 @@ function merge_physical_scales(rxs, physical_scales, default) scales[idx] = default end end - + scales end @@ -771,7 +772,7 @@ Notes: function Base.convert(::Type{<:JumpSystem}, rs::ReactionSystem; name = nameof(rs), combinatoric_ratelaws = get_combinatoric_ratelaws(rs), remove_conserved = nothing, checks = false, default_u0 = Dict(), default_p = Dict(), - defaults = _merge(Dict(default_u0), Dict(default_p)), expand_catalyst_funs = true, + defaults = _merge(Dict(default_u0), Dict(default_p)), expand_catalyst_funs = true, save_positions = (true, true), physical_scales = nothing, kwargs...) iscomplete(rs) || error(COMPLETENESS_ERROR) spatial_convert_err(rs::ReactionSystem, JumpSystem) @@ -780,28 +781,29 @@ function Base.convert(::Type{<:JumpSystem}, rs::ReactionSystem; name = nameof(rs flatrs = Catalyst.flatten(rs) - physical_scales = merge_physical_scales(reactions(flatrs), physical_scales, + physical_scales = merge_physical_scales(reactions(flatrs), physical_scales, PhysicalScale.Jump) - admissible_scales = (PhysicalScale.ODE, PhysicalScale.Jump, + admissible_scales = (PhysicalScale.ODE, PhysicalScale.Jump, PhysicalScale.VariableRateJump) unique_scales = unique(physical_scales) - (unique_scales ⊆ admissible_scales) || + (unique_scales ⊆ admissible_scales) || error("Physical scales must currently be one of $admissible_scales for hybrid systems.") # basic jump states and equations - eqs = assemble_jumps(flatrs; combinatoric_ratelaws, expand_catalyst_funs, + eqs = assemble_jumps(flatrs; combinatoric_ratelaws, expand_catalyst_funs, physical_scales, save_positions) ists, ispcs = get_indep_sts(flatrs) # handle coupled ODEs and BC species - if (PhysicalScale.ODE in unique_scales) || has_nonreactions(flatrs) - odeeqs = assemble_drift(flatrs, ispcs; combinatoric_ratelaws, + if (PhysicalScale.ODE in unique_scales) || has_nonreactions(flatrs) + odeeqs = assemble_drift(flatrs, ispcs; combinatoric_ratelaws, remove_conserved = false, physical_scales, include_zero_odes = true) append!(eqs, odeeqs) - eqs, us, ps, obs, defs = addconstraints!(eqs, flatrs, ists, ispcs; + eqs, us, ps, + obs, defs = addconstraints!(eqs, flatrs, ists, ispcs; remove_conserved = false) else - any(isbc, get_unknowns(flatrs)) && + any(isbc, get_unknowns(flatrs)) && (ists = vcat(ists, filter(isbc, get_unknowns(flatrs)))) us = ists ps = get_ps(flatrs) @@ -826,7 +828,7 @@ function DiffEqBase.ODEProblem(rs::ReactionSystem, u0, tspan, p = DiffEqBase.NullParameters(), args...; check_length = false, name = nameof(rs), combinatoric_ratelaws = get_combinatoric_ratelaws(rs), - include_zero_odes = true, remove_conserved = false, checks = false, + include_zero_odes = true, remove_conserved = false, checks = false, expand_catalyst_funs = true, structural_simplify = false, kwargs...) osys = convert(ODESystem, rs; name, combinatoric_ratelaws, include_zero_odes, checks, remove_conserved, expand_catalyst_funs) @@ -878,11 +880,11 @@ Keyword args and default values: function DiffEqBase.NonlinearProblem(rs::ReactionSystem, u0, p = DiffEqBase.NullParameters(), args...; name = nameof(rs), combinatoric_ratelaws = get_combinatoric_ratelaws(rs), - remove_conserved = false, conseqs_remake_warn = true, checks = false, - check_length = false, expand_catalyst_funs = true, + remove_conserved = false, conseqs_remake_warn = true, checks = false, + check_length = false, expand_catalyst_funs = true, structural_simplify = false, all_differentials_permitted = false, kwargs...) - nlsys = convert(NonlinearSystem, rs; name, combinatoric_ratelaws, checks, - all_differentials_permitted, remove_conserved, conseqs_remake_warn, + nlsys = convert(NonlinearSystem, rs; name, combinatoric_ratelaws, checks, + all_differentials_permitted, remove_conserved, conseqs_remake_warn, expand_catalyst_funs) nlsys = structural_simplify ? MT.structural_simplify(nlsys) : complete(nlsys) return NonlinearProblem(nlsys, u0, p, args...; check_length, @@ -893,10 +895,10 @@ end function DiffEqBase.SDEProblem(rs::ReactionSystem, u0, tspan, p = DiffEqBase.NullParameters(), args...; name = nameof(rs), combinatoric_ratelaws = get_combinatoric_ratelaws(rs), - include_zero_odes = true, checks = false, check_length = false, - remove_conserved = false, structural_simplify = false, + include_zero_odes = true, checks = false, check_length = false, + remove_conserved = false, structural_simplify = false, expand_catalyst_funs = true, kwargs...) - sde_sys = convert(SDESystem, rs; name, combinatoric_ratelaws, expand_catalyst_funs, + sde_sys = convert(SDESystem, rs; name, combinatoric_ratelaws, expand_catalyst_funs, include_zero_odes, checks, remove_conserved) # Handles potential differential algebraic equations (which requires `structural_simplify`). @@ -984,16 +986,16 @@ plot(sol, idxs = :A) """ function JumpInputs(rs::ReactionSystem, u0, tspan, p = DiffEqBase.NullParameters(); name = nameof(rs), combinatoric_ratelaws = get_combinatoric_ratelaws(rs), - checks = false, physical_scales = nothing, expand_catalyst_funs = true, + checks = false, physical_scales = nothing, expand_catalyst_funs = true, save_positions = (true, true), remake_warn = true, kwargs...) - jsys = complete(convert(JumpSystem, rs; name, combinatoric_ratelaws, checks, + jsys = complete(convert(JumpSystem, rs; name, combinatoric_ratelaws, checks, physical_scales, expand_catalyst_funs, save_positions)) - if MT.has_variableratejumps(jsys) || MT.has_equations(jsys) || - !isempty(MT.continuous_events(jsys)) + if MT.has_variableratejumps(jsys) || MT.has_equations(jsys) || + !isempty(MT.continuous_events(jsys)) prob = ODEProblem(jsys, u0, tspan, p; kwargs...) if remake_warn - @warn "JumpInputs has detected the system includes ODEs, variable rate jumps, or continuous events. Please note that currently remake does not work for such problems, and both JumpInputs and then JumpProblem must be called again if you wish to change any parameter or initial condition values. This warning can be disabled by passing JumpInputs the keyword argument `remake_warn = false`." + @warn "JumpInputs has detected the system includes ODEs, variable rate jumps, or continuous events. Please note that currently remake does not work for such problems, and both JumpInputs and then JumpProblem must be called again if you wish to change any parameter or initial condition values. This warning can be disabled by passing JumpInputs the keyword argument `remake_warn = false`." end else prob = DiscreteProblem(jsys, u0, tspan, p; kwargs...) @@ -1017,13 +1019,13 @@ end # DROP IN CATALYST 16 # DiscreteProblem from AbstractReactionNetwork function DiffEqBase.DiscreteProblem(rs::ReactionSystem, u0, tspan::Tuple, - p = DiffEqBase.NullParameters(), args...; name = nameof(rs), - combinatoric_ratelaws = get_combinatoric_ratelaws(rs), checks = false, + p = DiffEqBase.NullParameters(), args...; name = nameof(rs), + combinatoric_ratelaws = get_combinatoric_ratelaws(rs), checks = false, expand_catalyst_funs = true, kwargs...) Base.depwarn("DiscreteProblem(rn::ReactionSystem, ...) is deprecated and will be \ - removed in Catalyst 16. Use JumpInputs(rn, ...) instead.", + removed in Catalyst 16. Use JumpInputs(rn, ...) instead.", :DiscreteProblem) - jsys = convert(JumpSystem, rs; name, combinatoric_ratelaws, checks, + jsys = convert(JumpSystem, rs; name, combinatoric_ratelaws, checks, expand_catalyst_funs) jsys = complete(jsys) return DiscreteProblem(jsys, u0, tspan, p, args...; kwargs...) @@ -1038,7 +1040,7 @@ function JumpProcesses.JumpProblem(rs::ReactionSystem, prob::SciMLBase.AbstractD Base.depwarn("JumpProblem(rn::ReactionSystem, prob, ...) is \ deprecated and will be removed in Catalyst 16. Use \ JumpProblem(JumpInputs(rn, ...), ...) instead.", :JumpProblem) - jsys = convert(JumpSystem, rs; name, combinatoric_ratelaws, + jsys = convert(JumpSystem, rs; name, combinatoric_ratelaws, expand_catalyst_funs, checks) jsys = complete(jsys) return JumpProblem(jsys, prob, aggregator; kwargs...) @@ -1056,7 +1058,7 @@ function DiffEqBase.SteadyStateProblem(rs::ReactionSystem, u0, p = DiffEqBase.NullParameters(), args...; check_length = false, name = nameof(rs), combinatoric_ratelaws = get_combinatoric_ratelaws(rs), - remove_conserved = false, include_zero_odes = true, checks = false, + remove_conserved = false, include_zero_odes = true, checks = false, expand_catalyst_funs = true, structural_simplify = false, kwargs...) osys = convert(ODESystem, rs; name, combinatoric_ratelaws, include_zero_odes, checks, remove_conserved, expand_catalyst_funs) diff --git a/src/reactionsystem_serialisation/serialise_reactionsystem.jl b/src/reactionsystem_serialisation/serialise_reactionsystem.jl index 8d9ab5dffe..0401b207e3 100644 --- a/src/reactionsystem_serialisation/serialise_reactionsystem.jl +++ b/src/reactionsystem_serialisation/serialise_reactionsystem.jl @@ -70,18 +70,22 @@ function get_full_system_string(rn::ReactionSystem, annotate::Bool, top_level::B # to the function that creates the next sub-system declarations. file_text, _ = push_field(file_text, rn, annotate, top_level, IV_FS) file_text, has_sivs = push_field(file_text, rn, annotate, top_level, SIVS_FS) - file_text, has_parameters, has_species, has_variables = handle_us_n_ps( + file_text, has_parameters, + has_species, has_variables = handle_us_n_ps( file_text, rn, annotate, top_level) file_text, has_reactions = push_field(file_text, rn, annotate, top_level, REACTIONS_FS) file_text, has_equations = push_field(file_text, rn, annotate, top_level, EQUATIONS_FS) file_text, has_observed = push_field(file_text, rn, annotate, top_level, OBSERVED_FS) file_text, has_defaults = push_field(file_text, rn, annotate, top_level, DEFAULTS_FS) - file_text, has_continuous_events = push_field(file_text, rn, annotate, + file_text, + has_continuous_events = push_field(file_text, rn, annotate, top_level, CONTINUOUS_EVENTS_FS) - file_text, has_discrete_events = push_field(file_text, rn, annotate, + file_text, + has_discrete_events = push_field(file_text, rn, annotate, top_level, DISCRETE_EVENTS_FS) file_text, has_systems = push_systems_field(file_text, rn, annotate, top_level) - file_text, has_connection_type = push_field(file_text, rn, annotate, + file_text, + has_connection_type = push_field(file_text, rn, annotate, top_level, CONNECTION_TYPE_FS) # Finalise the system. Creates the final `ReactionSystem` call. @@ -92,7 +96,8 @@ function get_full_system_string(rn::ReactionSystem, annotate::Bool, top_level::B has_equations, has_observed, has_defaults, has_continuous_events, has_discrete_events, has_systems, has_connection_type) annotate || (@string_prepend! "\n" file_text) - annotate && top_level && @string_prepend! "\n# Serialised using Catalyst version v$(Catalyst.VERSION)." file_text + annotate && top_level && + @string_prepend! "\n# Serialised using Catalyst version v$(Catalyst.VERSION)." file_text @string_prepend! "let" file_text @string_append! file_text "\n\n" rs_creation_code "\n\nend" diff --git a/src/spatial_reaction_systems/lattice_jump_systems.jl b/src/spatial_reaction_systems/lattice_jump_systems.jl index 7f9cf3e0e4..e776ddee7e 100644 --- a/src/spatial_reaction_systems/lattice_jump_systems.jl +++ b/src/spatial_reaction_systems/lattice_jump_systems.jl @@ -19,7 +19,8 @@ function DiffEqBase.DiscreteProblem(lrs::LatticeReactionSystem, u0_in, tspan, # vert_ps values are vectors. Here, index (i) is a parameter's value in vertex i. # edge_ps values are sparse matrices. Here, index (i,j) is a parameter's value in the edge from vertex i to vertex j. # Uniform vertex/edge parameters store only a single value (a length 1 vector, or size 1x1 sparse matrix). - vert_ps, edge_ps = lattice_process_p(p_in, vertex_parameters(lrs), + vert_ps, + edge_ps = lattice_process_p(p_in, vertex_parameters(lrs), edge_parameters(lrs), lrs) # Returns a DiscreteProblem (which basically just stores the processed input). @@ -64,6 +65,7 @@ function make_hopping_constants(dprob::DiscreteProblem, lrs::LatticeReactionSyst # vector containing all edges leading out from that vertex (sorted by destination index). edge_array = [Pair{Int64, Int64}[] for _1 in 1:num_species(lrs), _2 in 1:num_verts(lrs)] for e in edge_iterator(lrs), s_idx in 1:num_species(lrs) + push!(edge_array[s_idx, e[1]], e) end foreach(e_vec -> sort!(e_vec; by = e -> e[2]), edge_array) diff --git a/src/spatial_reaction_systems/lattice_reaction_systems.jl b/src/spatial_reaction_systems/lattice_reaction_systems.jl index 46a2eafec6..e95c902176 100644 --- a/src/spatial_reaction_systems/lattice_reaction_systems.jl +++ b/src/spatial_reaction_systems/lattice_reaction_systems.jl @@ -155,7 +155,7 @@ struct LatticeReactionSystem{Q, R, S, T} <: MT.AbstractTimeDependentSystem # Additional error checks. if any(haskey(Symbolics.unwrap(symvar).metadata, Symbolics.ArrayShapeCtx) - for symvar in [ps; species(rs)]) + for symvar in [ps; species(rs)]) throw(ArgumentError("Some species and/or parameters used to create the `LatticeReactionSystem` are array variables ($(filter(symvar -> haskey(Symbolics.unwrap(symvar).metadata, Symbolics.ArrayShapeCtx), [ps; species(rs)]))). This is currently not supported.")) end @@ -401,8 +401,10 @@ end Returns `true` if `lrs` was created using a cartesian grid lattice (e.g. created via `CartesianGrid(5,5)`). Otherwise, returns `false`. """ -has_cartesian_lattice(lrs::LatticeReactionSystem) = lattice(lrs) isa - CartesianGridRej{N, T} where {N, T} +function has_cartesian_lattice(lrs::LatticeReactionSystem) + lattice(lrs) isa + CartesianGridRej{N, T} where {N, T} +end """ has_masked_lattice(lrs::LatticeReactionSystem) diff --git a/src/spatial_reaction_systems/lattice_sim_struct_interfacing.jl b/src/spatial_reaction_systems/lattice_sim_struct_interfacing.jl index 4d176bae59..6faccc5fc8 100644 --- a/src/spatial_reaction_systems/lattice_sim_struct_interfacing.jl +++ b/src/spatial_reaction_systems/lattice_sim_struct_interfacing.jl @@ -430,7 +430,8 @@ function rebuild_lat_internals!(lt_ofun::LatticeTransportODEFunction, ps_new, # Recreates the new parameters on the requisite form. ps_new = [(length(p) == 1) ? p[1] : p for p in deepcopy(ps_new)] ps_new = [p => p_val for (p, p_val) in zip(parameters(lrs), deepcopy(ps_new))] - vert_ps, edge_ps = lattice_process_p(ps_new, vertex_parameters(lrs), + vert_ps, + edge_ps = lattice_process_p(ps_new, vertex_parameters(lrs), edge_parameters(lrs), lrs) ps_new = [vert_ps; edge_ps] diff --git a/src/spatial_reaction_systems/spatial_ODE_systems.jl b/src/spatial_reaction_systems/spatial_ODE_systems.jl index 2e3a6d46f2..0f6d919a4a 100644 --- a/src/spatial_reaction_systems/spatial_ODE_systems.jl +++ b/src/spatial_reaction_systems/spatial_ODE_systems.jl @@ -71,7 +71,8 @@ struct LatticeTransportODEFunction{P, Q, R, S, T} # Computes `LatticeTransportODEFunction` functor fields. heterogeneous_vert_p_idxs = make_heterogeneous_vert_p_idxs(ps, lrs) mtk_ps, p_setters = make_mtk_ps_structs(ps, lrs, heterogeneous_vert_p_idxs) - t_rate_idx_types, leaving_rates = make_t_types_and_leaving_rates(transport_rates, + t_rate_idx_types, + leaving_rates = make_t_types_and_leaving_rates(transport_rates, lrs) # Creates and returns the `LatticeTransportODEFunction` functor. @@ -196,7 +197,8 @@ function DiffEqBase.ODEProblem(lrs::LatticeReactionSystem, u0_in, tspan, # edge_ps values are sparse matrices. Here, index (i,j) is a parameter's value in the edge from vertex i to vertex j. # Uniform vertex/edge parameters store only a single value (a length 1 vector, or size 1x1 sparse matrix). # In the `ODEProblem` vert_ps and edge_ps are merged (but for building the ODEFunction, they are separate). - vert_ps, edge_ps = lattice_process_p(p_in, vertex_parameters(lrs), + vert_ps, + edge_ps = lattice_process_p(p_in, vertex_parameters(lrs), edge_parameters(lrs), lrs) # Creates the ODEFunction. @@ -311,6 +313,7 @@ end # transport reaction values. function set_jac_transport_values!(jac_prototype, transport_rates, lrs) for (s, rates) in transport_rates, e in edge_iterator(lrs) + idx_src = get_index(e[1], s, num_species(lrs)) idx_dst = get_index(e[2], s, num_species(lrs)) val = get_transport_rate(rates, e, size(rates) == (1, 1)) diff --git a/src/spatial_reaction_systems/spatial_reactions.jl b/src/spatial_reaction_systems/spatial_reactions.jl index ca08724882..455afc46ec 100644 --- a/src/spatial_reaction_systems/spatial_reactions.jl +++ b/src/spatial_reaction_systems/spatial_reactions.jl @@ -57,7 +57,6 @@ function make_transport_reaction(rateex, species) # Checks for input errors. forbidden_symbol_check(union([species], parameters)) - # Creates expressions corresponding to actual code from the internal DSL representation. sexprs = get_usexpr([species], Dict{Symbol, Expr}()) pexprs = get_psexpr(parameters, Dict{Symbol, Expr}()) @@ -104,13 +103,13 @@ function check_spatial_reaction_validity(rs::ReactionSystem, tr::TransportReacti error("A transport reaction used a species, $(tr.species), with metadata not matching its lattice reaction system. Please fetch this species from the reaction system and use it during transport reaction creation.") end if any(isequal(rs_p, tr_p) && !isequivalent(rs_p, tr_p) - for rs_p in parameters(rs), tr_p in Symbolics.get_variables(tr.rate)) + for rs_p in parameters(rs), tr_p in Symbolics.get_variables(tr.rate)) error("A transport reaction used a parameter with metadata not matching its lattice reaction system. Please fetch this parameter from the reaction system and use it during transport reaction creation.") end # Checks that no edge parameter occurs among rates of non-spatial reactions. if any(!isempty(intersect(Symbolics.get_variables(r.rate), edge_parameters)) - for r in reactions(rs)) + for r in reactions(rs)) error("Edge parameter(s) were found as a rate of a non-spatial reaction.") end end @@ -125,10 +124,10 @@ const ep_metadata = Catalyst.EdgeParameter => true function isequivalent(sym1, sym2; ignored_metadata = [MT.SymScope]) isequal(sym1, sym2) || (return false) if any((md1 != ep_metadata) && (md1[1] ∉ ignored_metadata) && (md1 ∉ sym2.metadata) - for md1 in sym1.metadata) + for md1 in sym1.metadata) return false elseif any((md2 != ep_metadata) && (md2[1] ∉ ignored_metadata) && (md2 ∉ sym1.metadata) - for md2 in sym2.metadata) + for md2 in sym2.metadata) return false elseif typeof(sym1) != typeof(sym2) return false diff --git a/src/spatial_reaction_systems/utility.jl b/src/spatial_reaction_systems/utility.jl index f0a30e331d..b895b90107 100644 --- a/src/spatial_reaction_systems/utility.jl +++ b/src/spatial_reaction_systems/utility.jl @@ -182,10 +182,9 @@ function make_sidxs_to_transrate_map(vert_ps::Vector{Pair{R, Vector{T}}}, # First, compute a map from species in their symbolics form to their values. # Next, convert to map from species index to values. transport_rates_speciesmap = compute_all_transport_rates(p_val_dict, lrs) - return Pair{Int64, SparseMatrixCSC{T, Int64}}[ - speciesmap(reactionsystem(lrs))[spat_rates[1]] => spat_rates[2] - for spat_rates in transport_rates_speciesmap - ] + return Pair{Int64, SparseMatrixCSC{T, Int64}}[speciesmap(reactionsystem(lrs))[spat_rates[1]] => spat_rates[2] + for spat_rates in + transport_rates_speciesmap] end # Computes the transport rates for all species with transportation rates. Output is a map