alanderos91 / BioSimulator.jl

A stochastic simulation framework in Julia.
https://alanderos91.github.io/BioSimulator.jl/stable/
Other
47 stars 7 forks source link

Extract transition matrix #5

Closed rveltz closed 6 years ago

rveltz commented 6 years ago

Hi,

I am sorry for posting this here because this is not really an issue. I would like to extract the (dense) transition matrix from a model, say for example, from your Birth-Death-Immigration Process.

Is it possible?

Further, I can see that there is a macro in DiffEqJump that extract the process from a list of reactions. Would your solution work in this case too?

Thank you a lot for your help,

Best regards

alanderos91 commented 6 years ago

I would like to extract the (dense) transition matrix from a model, say for example, from your Birth-Death-Immigration Process.

I'm assuming you're referring to the matrix describing the net change due to each reaction? You can't do this out-of-the-box, but all the information is there. Here's a quick an dirty solution:

function stoichmatrix(model::Network)
    d = n_species(model)
    c = n_reactions(model)
    v = zeros(Int, d, c)

    stoichmatrix!(v, model)
end

function stoichmatrix!(v, model::Network)
    species = species_list(model)
    reactions = reaction_list(model)

    for (j, r) in enumerate(values(reactions))
        for (i, s) in enumerate(keys(species))
            vm = get(r.reactants, s, 0)
            vp = get(r.products, s, 0)

            v[i, j] = vp - vm
        end
    end
    return v
end

As an example (assuming using BioSimulator):

include(joinpath(Pkg.dir("BioSimulator"), "test", "test_models.jl"))

m = kendall() # birth-death-immigration

reaction_list(m)
# DataStructures.OrderedDict{Symbol,BioSimulator.Reaction} with 3 entries:
#   :birth       => X --> 2 X
#   :death       => X --> ∅
#   :immigration => ∅ --> X

stoichmatrix(kendall())
# 1×3 Array{Int64,2}:
#   1  -1  1

Further, I can see that there is a macro in DiffEqJump that extract the process from a list of reactions. Would your solution work in this case too?

@reaction_network, right? Any model constructed with that macro has the same information (and more) as the Network construct in BioSimulator. However, I'm not sure I understand the question.

rveltz commented 6 years ago

Thank you! You did understand the question...