Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add generation of functions to pass to ODEProblem, SDEProblem, JumpProblem #178

Open
sdwfrost opened this issue Sep 27, 2024 · 1 comment

Comments

@sdwfrost
Copy link

Functionality to generate ODEProblem, SDEProblem and JumpProblem present in Petri.jl are missing in AlgebraicPetri. Can this functionality be added? This has to be done carefully in order to accommodate time-varying parameters correctly. @slwu89 suggested using CompTime.jl for this. @slmontes - could you take this on? @slwu89 - would you mind giving some pointers?

@slwu89
Copy link
Member

slwu89 commented Sep 27, 2024

No problem @sdwfrost happy to be tagged for any questions. It wouldn't be that hard to modify this function

""" vectorfield_expr(pn::AbstractPetriNet)
Generates a Julia expression which is then evaluated that calculates the
vectorfield of the Petri net being simulated under the law of mass action.
The resulting function has a signature of the form `f!(du, u, p, t)` and can be
passed to the DifferentialEquations.jl solver package.
"""
vectorfield_expr(pn::AbstractPetriNet) = begin
fquote = Expr(:function, Expr(:tuple, :du, :u, :p, :t))
fcode = Expr[]
num_t = nt(pn)
# generate vector of rate constants for each transition
p_ix = [tname(pn, i) for i in 1:nt(pn)]
push!(fcode, :(
p_m = Vector{Union{Float64,Function}}(undef,$(num_t))
))
for i in 1:num_t
if eltype(p_ix) <: Symbol
push!(fcode, :(
p_m[$(i)] = p[$(Meta.quot(p_ix[i]))]
))
else
push!(fcode, :(
p_m[$(i)] = p[$(p_ix[i])]
))
end
end
# for each transition, calculate its firing intensity
for i in 1:num_t
is_ix = subpart(pn, incident(pn, i, :it), :is) # input places
is_pn_ix = [sname(pn, j) for j in is_ix]
os_ix = subpart(pn, incident(pn, i, :ot), :os) # output places
os_pn_ix = [sname(pn, j) for j in os_ix]
# generate vector of markings just for t's inputs and calc rate
n_input = length(is_pn_ix)
push!(fcode, :(
inputs = zeros($(n_input))
))
for j in 1:n_input
if eltype(is_pn_ix) <: Symbol
push!(fcode, :(
inputs[$(j)] = u[$(Meta.quot(is_pn_ix[j]))]
))
else
push!(fcode, :(
inputs[$(j)] = u[$(is_pn_ix[j])]
))
end
end
push!(fcode, :(
rate = valueat(p_m[$(i)], u, t) * prod(inputs)
))
# transition decreases inputs and increases output marking
if eltype(os_pn_ix) <: Symbol
for j in os_pn_ix
push!(fcode, :(
du[$(Meta.quot(j))] += rate
))
end
for j in is_pn_ix
push!(fcode, :(
du[$(Meta.quot(j))] -= rate
))
end
else
# dont need quote nodes
for j in os_pn_ix
push!(fcode, :(
du[$(j)] += rate
))
end
for j in is_pn_ix
push!(fcode, :(
du[$(j)] -= rate
))
end
end
end
push!(fcode, :(
return du
))
push!(fquote.args, Expr(:block, fcode...))
return mk_function(AlgebraicPetri, fquote)
end
to output something to make a CTMC (JumpProblem) simulator. That function is coded rather "manually", pushing to the Julia AST (Expr). CompTime.jl could make the whole thing look much nicer, but for a first pass, probably easier to just do it without generated code. Especially regarding the time-varying parameters which make the whole thing depend on distributions more complex than simple exponentials.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

2 participants