Generation of stochastic models

[sdwfrost]

Hi!

Are there any plans to implement the construction of stochastic models (SDEs, jump systems for continuous systems, Markov models for discrete time systems)?

[mhlr]

I also curious how AlgbraicJulia could be used to enhance JuMP and JuliaConstraints

[slwu89]

@sdwfrost I'm interesting in returning to think about this issue, albeit starting with only jump systems for continuous time discrete event (CTDE) simulation. I am going to write down some thoughts I have here regarding using DWD or UWD composition for CTDE systems, which may not come up when dealing with deterministic differential dynamics.

To provide concreteness I'll use the Lotka-Volterra equations from the docs https://algebraicjulia.github.io/AlgebraicDynamics.jl/dev/examples/Lotka-Volterra/ but let $\beta=\gamma$ throughout. When thinking of the equations as a fluid approximation of a CTDE system, this is more reasonable, as the predation event is a discrete event that occurs at a point in time and instantaneously decreases the prey population by 1 and increases the predator population by 1.

The LV equations are then:

$$ \begin{align*} \dot{r} &= \alpha r - \beta r f \\ \dot{f} &= \beta r f - \delta f \end{align*} $$

Interpreting this as a CTDE system, there are 3 events which change the 2 states. Let's assume our CTDE system is a continuous time Markov chain, then we have 3 Poisson processes whose intensities are functions of the state and jump when things happen which update the state, we list them below:

• birth of prey, occurs with intensity $\alpha r$ and updates state as $r += 1$
• death of predator, occurs with intensity $\delta f$ and updates state as $f -= 1$
• predation, occurs with intensity $\beta r f$ and updates state as $r -= 1$ and $f += 1$

It would seem that for this specific example, UWD composition would be a nicer way to generate this stochastic CTDE model. Each of the 3 events corresponds to a box/resource sharer. This is nice because each source of randomness in the model gets its own box. Composition means identifying variables that meet at a junction.

Machines would not be as nice for this. The way I understand it, when using DWDs as the composition syntax, the number of states for the composed model is the coproduct of the states in each machine. But that won't do if we want to have a rule that says "one box, one Poisson process", because clearly the rabbits are affected by birth and predation, and the foxes are affected by death and predation. And, in the CTDE world, we'd need predation to have a shared source of randomness for the model to be consistent, which I find hard to imagine how to do with machines, although I suppose it could be done with instantaneous machines, somehow.

But there are some nice things about machines. One thing is that the input ports are "read only". This means that intensity functions can depend on parts of the state that the event doesn't change. This isn't apparent in the simple LV example, but consider disease transmission models where the force of infection has some term $1/N$ where $N$ is the total population (i.e. frequency dependent transmission). In this case, the intensity depends on the total state (or at least all those considered part of $N$), but the event would only affect 2 at most (moves a susceptible to infected). This could be done with UWDs in a somewhat "hacky" way by just setting the change of those wires to be 0, but that seems somewhat unsatisfying.

[slwu89]

Going to tag @adolgert here, as he might be interested. Especially given that the type of CTDE I'm describing is basically what he describes in https://arxiv.org/abs/1610.03939

[slwu89]

@sdwfrost I implemented a very silly quick example of "Poisson resource sharers" here: https://github.com/slwu89/AlgebraicDynamics.jl/blob/stochastic/stochastic.jl which simulates a birth-death process. The only real changes are that the interface struct now has a field affects which tells us how the jump affects the state that box is allowed to see when it fires. In the example it's just a vector of integers to which can add or subtract from the state, but it could be more general, including random variables to allow simulation of compound Poisson processes.

The composition pattern looks like this and each box is a jump, where the stochasticity comes from a time-changed Poisson process. It produces a realization of a typical birth death process.

@slibkind when looking through things I was checking against the oapply and induced_states methods from uwd_dynam.jl. I just had a few questions about things (due to my own lack of clarity). In the oapply method, the 3rd argument is S′::Pushout, which comes from induced_states. But within that function we again take the coproduct in the line S = coproduct((FinSet∘nstates).(xs)). Can you explain why we need to do this operation twice?

Also in induced_dynamics, can you check my understanding of arguments? state_map is a mapping from the coproduct of all states in each box, to the reduced set of states after identifying via junctions? And states is a mapping from boxes into the domain of state_map?

[slibkind]

This is a great example @slwu89!

To answer your questions about oapply and induced_dynamics. In both examples there are two things to keep track of: the states before the composition and the states after the composition. In your example, imagine that the birth and death boxes both contain a resource sharer with 2 states. Then there are 4 in the coproduct of the two systems (before composition) but only 3 states in the composite system (because two of them get identified). In the code the total states of the coproduct are indicated by S and the states of the composite system are indicated by S'. This is why we still need to compute S in oapply even though S' is passed in as an argument.

In induced_dynamics, states is a function which takes a box and returns a list of the states in the coproduct that correspond to that box. For example, in the example states(1) would return [1,2] and states(2) would return [3,4] because the first box corresponds to the first two states and the second box corresponds to the remaining two states. The state_map function is a way of relating the states before and after composition. It essentially is a function from S to S' which says what happens to the states after the composition. In this example state_map.func is the list [1,2,2,3] because the second state of the first box (which is the second state in S) and the first state of the second box (which is the third state in S) are identified becoming the second state in S'.

[jpfairbanks]

I think that a jump process simulator for resource sharers is a great feature. Doing just the vector of increments for each process makes a lot of sense. We could take models from AlgebraicPetri and give them a poisson process semantics for each transition and then use AlgebraicDynamics as a semantics for that. Really cool way of blending different frameworks.

[slwu89]

Hi @slibkind, thank you very much for the clear explanation!

That's good to hear @jpfairbanks; the connection to AlgPetri is nice too, that would be a nice way to wrangle stochastic models out of the PNs you can build in that package. Since it seems sufficiently interesting, I'll poke around with reasonable ways to go about it on my fork.

In general I quite like resource sharers for building up stochastic systems for the reasons listed above (each event that can change a part of state gets its own box). I could see AlgPetri and AlgDynamics working really nicely together to build complex stochastic models, the PNs in AlgPetri work really nicely with pullbacks to stratify things, and then using resource sharers to handle nasty coupling bits, like force of infection terms with non-mass action equations, covers a huge amount of models of practical interest.

[jpfairbanks]

That would be great. I think Catalyst has some similar capabilities that would be good to build on / learn from

[slwu89]

I haven't looked at Catalyst before in detail, thanks for the suggestion @jpfairbanks. One interesting observation, at the end of the first section in this manual file https://docs.sciml.ai/Catalyst/stable/introduction_to_catalyst/introduction_to_catalyst/ there is a Petri-net like representation of a ReactionSystem (apparently using Catlab's graphviz code). The red arrows in the graph refer to a "read only" relationship between that reaction and that species, that is, that species is needed to compute the rate function but is neither consumed nor produced by the firing of that reaction. This is basically what I was discussing in my wall of text earlier; I can replicate this "read only" idea with the existing resource sharers by simply saying the effect of these ports on their junctions is zero. I'll do this for initial prototyping of course. But perhaps in the future 2 types of ports may be considered, read only and effecting? @slibkind what do you think about this idea?

[jpfairbanks]

That would be an interesting Operad. It would be cool to show that we could make that work. I think @olynch has thought about UWDs with two kinds of ports.

[slibkind]

I can replicate this "read only" idea with the existing resource sharers by simply saying the effect of these ports on their junctions is zero. I'll do this for initial prototyping of course. But perhaps in the future 2 types of ports may be considered, read only and effecting?

@slwu89 That's an interesting idea. My initial reaction is that is exactly what resource sharing machines are for. Resource sharing machines can have input ports, which --- as you noted in an earlier comment --- are read-only. I haven't implemented these yet but if you have a good use case that's even more reason to do so!

[slwu89]

Thanks for the reminder re; resource sharing machines @slibkind! I'll review your work on those within the next week. In the meantime, can I ask what are the main differences between using an operad based on RSMs as opposed to the proposed UWDs with 2 kinds of ports? Just to summarize my thoughts from all this, some things I really like about UWDs with 2 ports for these kinds of models are (I'm aware this is repeating what we both are clear on...just want to list it again for my own notes in this conversation!):

1. state is "outside" of the boxes, in the junctions.
2. each box represents a single event (crucial for the jump process setting)
3. each box can depend on multiple inputs (read only); crucial for non-mass action terms
4. each box can affect multiple states (e.g. someone sneezes on a crowded subway and infects a lot of people, etc etc)

[slibkind]

Oh! I didn't understand that the state is entirely in the junctions. That would be a big difference between RSMs and the proposed algebra on UWDs with two ports. I think another difference might be that in the operad for RSMs the input ports and the junction ports have are wired together in different ways. Whereas in your framework my understanding is that they should both connect at junctions.

[olynch]

@slibkind wrote a paper on unifying parameter setting and resource sharing: https://arxiv.org/abs/2007.14442. Did we ever implement this in AlgebraicJulia? Might be good to revisit.

[slibkind]

We have not! This is a bit tangential, but one thing I was wondering about that is the morphisms of the operad can be specified by a $\mathcal C$-set where $\mathcal C$ is a pushout of the schema for UWDs and DWDs. I was wondering if it's possible to do that automatically yet?

[jpfairbanks]

@olynch, can we do schema pushouts in GATLab yet?

[olynch]

Yes.

[jpfairbanks]

@slwu89, are you thinking of the junctions as representing individual agents in the model or are they counts of agents in a state? Your item (4) makes me think that you are using one junction per agent, which would be pretty different from the current ODE approach, and the ABM approach.

[slwu89]

@jpfairbanks so I am thinking of building systems of interacting counting processes of the sort proposed in this paper https://arxiv.org/abs/1610.03939 (author @adolgert is a former colleague and we built lots of big stochastic models of this type before, albeit without the categorical machinery AlgJulia can bring to bear!).

What I am thinking is that each box in the UWD contains a source of randomness (a counting process). The intensity function for the driving process can read some number of junctions in order to calculate the intensity (the read only ports). When the counting process fires it may affect some number of junctions (the effecting ports).

I'm being purposely vague about if the junctions are agents or counts of agents, the counting process formalism should work for either. Like for an SIR model we could have 3 junctions and 2 boxes (infection process and recovery process). Or we could do for N persons in an SIR model, N junctions (perhaps coding 0,1,2 as S,I,R) and 2N boxes (infection and recovery for each person).

I'm planning to deeply review @slibkind's paper soon to be able to talk about this in context of that work.

[jpfairbanks]

Ok I think that the counts of agents approach should be a better fit for this project, of course if you have counts of agents that are either 0 or 1, then that recovers the individual scale model, but I think that the ABM world is sufficiently different from the jump/counting process view that it deserves its own package. Counting processes of RSMs would be a great fit here.

[slwu89]

Yeah, I think we are in agreement. To be clear, I don't mean "agent" to refer to the very complicated AnyLogic style models of certain ABMs, where each agent possibly has a little internal model of the world it can use to make decisions. I'm not as familiar with such models.

But "individual" based models where we deal with counts rather than densities, where state jumps at the points of counting processes, is definitely something I'm excited about! The expanded representation could sometimes be useful if for example you wanted to "run forwards" a survival analysis type setup, where perhaps each individual in the SIR example had their recovery rate drawn from a frailty distribution, to model individual heterogeneity.

[sdwfrost]

Although the distinction with ABMs is a bit blurry, I think that the "individual" models are more along the line of discrete event simulations, like this SimJulia example.

[slwu89]

After taking a break from thinking about this for awhile and coming back to it, I don't think its feasible to do right now. The birth-death composition example above would work nicely for a "flat" composition where one already had the entire UWD for the whole system and wanted to insert different counting processes in each resource sharer. This honestly is already a nice improvement over existing methods (are there any?) for structuring general models based off competing counting processes.

The issue is that each box in the UWD is associated with an independent Poisson process, and they need to remain that way under composition. A full theory of how to do this probably has to wait for more work on the math side...maybe @epatters work on semi-markov processes will eventually help. I'd be curious to see if an eventual framework could handle the birth-death example above composed with other processes.