Open nathanielvirgo opened 1 year ago
Does this example help?
https://docs.sciml.ai/Catalyst/stable/example_networks/smoluchowski_coagulation_equation/
It does help, but in terms of documentation it's not very clear - it's a very complicated example and that particular aspect of it isn't really explained. I guess the relevant line is
@species k[1:nr] (X(t))[1:N]
and I guess the k[1:nr]
part is optional, since that seems to be declaring reaction rates and not species. (Which seems kind of weird - why aren't they declared with @parameters
instead?) I can extrapolate from this by guesswork but it would be nice to have it documented properly.
It would also be helpful to know how to do the following things, if it's possible to do them:
X[i]
. In the example I gave calling them X[i]
makes sense, but it's easy to imagine cases where you'd want to name them differently, with the names also generated programmatically in the form of strings.Yeah, k
should be a parameter. It is a relic from a much older version of Catalyst where there was less of a distinction, and looks like we missed updating it.
You can create a species from a string as follows:
str = "AB"
strsym = Symbol(str)
@variables t
spcs = @species $(strsym)(t)
giving
julia> spcs[1]
AB(t)
We should indeed update the programmatic tutorial to show this.
If you want to add new species and/or reactions I'd suggest creating a new network using extend
and/or compose
as shown at https://docs.sciml.ai/Catalyst/stable/catalyst_functionality/compositional_modeling/
These take one or more systems and merge them together / create a tree-like system from them.
Alternatively, we do provide addspecies!
, addparam!
and addreaction!
to mutate a system, see
https://docs.sciml.ai/Catalyst/stable/api/catalyst_api/#Functions-to-extend-or-modify-a-network
however, I'd recommend not using these unless really needed as they may be deprecated in the near future. (The reason being that ModelingToolkit wants systems to be immutable, with adding/removing components done via generating new systems. So we have debated removing these features to maintain compatibility with that philosophy, and don't really advertise them in the docs.)
Great, that's really helpful, thank you
I'm having trouble using species created this way in reactions:
@variables t
@parameters α
# species created normally
@species x(t)
@species y(t)
# species created programmatically
strsym = Symbol("A")
a = @species $(strsym)(t)
strsym = Symbol("B")
b = @species $(strsym)(t)
# 'normal' species work ok:
Reaction(α,[x],[y])
# but programmatically created ones do not:
Reaction(α,[a],[b])
gives the error
MethodError: no method matching getmetadata(::Vector{Num}, ::Type{Catalyst.VariableSpecies}, ::Bool)
Closest candidates are:
getmetadata(::SymbolicUtils.Symbolic, ::Any, ::Any)
@ SymbolicUtils ~/.julia/packages/SymbolicUtils/Oyu8Z/src/types.jl:629
getmetadata(::SymbolicUtils.Symbolic, ::Any)
@ SymbolicUtils ~/.julia/packages/SymbolicUtils/Oyu8Z/src/types.jl:620
getmetadata(::Complex{Num}, ::Any)
@ Symbolics ~/.julia/packages/Symbolics/BQlmn/src/Symbolics.jl:156
...
Stacktrace:
[1] isspecies(s::Vector{Num})
@ Catalyst ~/.julia/packages/Catalyst/zTvWJ/src/reactionsystem.jl:36
[2] isvalidreactant(s::Vector{Num})
@ Catalyst ~/.julia/packages/Catalyst/zTvWJ/src/reactionsystem.jl:126
[3] _all
@ ./reduce.jl:1283 [inlined]
[4] #all#830
@ ./reducedim.jl:1007 [inlined]
[5] all
@ ./reducedim.jl:1007 [inlined]
[6] Reaction(rate::Num, subs::Vector{Vector{Num}}, prods::Vector{Vector{Num}}, substoich::Vector{Int64}, prodstoich::Vector{Int64}; netstoich::Nothing, only_use_rate::Bool, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
@ Catalyst ~/.julia/packages/Catalyst/zTvWJ/src/reactionsystem.jl:172
[7] Reaction
@ ~/.julia/packages/Catalyst/zTvWJ/src/reactionsystem.jl:128 [inlined]
[8] Reaction(rate::Num, subs::Vector{Vector{Num}}, prods::Vector{Vector{Num}}; kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
@ Catalyst ~/.julia/packages/Catalyst/zTvWJ/src/reactionsystem.jl:191
[9] Reaction(rate::Num, subs::Vector{Vector{Num}}, prods::Vector{Vector{Num}})
@ Catalyst ~/.julia/packages/Catalyst/zTvWJ/src/reactionsystem.jl:188
[10] top-level scope
@ In[63]:18
See my example above; @species
returns a vector with the first entry being your symbolic variable.
The documentation at https://docs.sciml.ai/Catalyst/stable/catalyst_functionality/programmatic_CRN_construction/ explains how to generate reactions programmatically, but it doesn't mention how to do this for species. It would obviously be a huge drawback if you had to enter the list of species manually each time!
I'm submitting this as a bug report because I think this should be fixed in the documentation, but in addition I would also just like to know how to do this.
for example, what should I do if I want to generate this reaction network?
A1 <=> A2 A2 <=> A3 ... A999 <=> A1000