agents refactor

[thevolatilebit]

The objective of this workstream is to naturally align the framework with AlgebraicAgents.jl and further enhance the breadth of dynamics that this framework can effectively simulate.

• key simulation objects (such as ReactiveDynamicsState, Transition, or Observable) are now modeled as agents;
• the problem is no longer converted into a DiscreteProblem. Instead, an acset instance is 'materialized' into a ReactionNetworkProblem, which is simply an agent with implemented stepping;
• additionally, custom time-stepping, which was not functioning properly in the stable version, has been fixed;
• add support for structured agent types.

Minimal demo:

using ReactiveDynamics

# define the network
acs = @ReactionNetworkSchema begin
    1.0, X --> Y, name => "transition1"
end

@prob_init acs X = 10 Y = 20
@prob_params acs
@prob_meta acs tspan = 25 dt = 0.10

# convert network into an AlgAgents hierarchy
prob = ReactionNetworkProblem(acs)

# simulate
simulate(prob)

# access solution
prob.sol

# plot solution
draw(prob)

For a demo of structured agents, see this tutorial.

What's coming next?

• often, it is unnecessary to implement 'macro' forms for aspects like simulation or plotting, if a standard function call will suffice. Now, the macro forms should be reserved primarily for problem definition;
• rather than implementing naive optimization routines, it might be more prudent to implement the generic set/get parameter interface, as suggested by AlgebraicAgents.jl;
• consider whether to use Catlab.jl as the backend. Currently, changes in Catlab often disrupt acset indexing.

Along with @slwu89, we would like to concur on and outline an 'improved' schema that defines a reaction network. Additionally, we aim to establish better type names used for 1) the acset defining the reaction network and 2) the agent implementing the reaction network (i.e., materializing the definition).

[slwu89]

Thanks @thevolatilebit, I'm very excited to speak more about an improved schema. In the meantime, I'll familiarize myself with the latest changes here.

I highly agree with the decision to use more standard functions rather than macros.

I would very much like to discuss more about Catlab as a backend. In fact recently the devs have taken the parts related to the acset interface and basic acset functionality that does not directly depend on the more "category theoretic" bits of Catlab, and moved it to a new package AlgebraicJulia/ACSets.jl to enhance stability. So I think stability will improve in the future.

[slwu89]

Hi @thevolatilebit, I am updating the code locally to work with the new stable Catlab/ACSets packages. I'll be able to continue as soon as my PR AlgebraicJulia/ACSets.jl#58 is merged, as some of the attributes in the current reaction network are union types including Missing.

[thevolatilebit]

It is now possible to combine, in a single reaction network,

• "classical species" (continuous or discrete; considered as pure quantities);
• "structured agents" (possibly with custom evolutionary function; these can appear both on LHS and RHS).

Consider the following example. You can file the full script there.

@push network begin
    # With specified intensities, generate experimental resources.
    ρ1, ∅ --> R1

    # Generate "Molecule 1" (where the integer corresponds to a "state" of, e.g., experimental triage).
    ρ3, ∅ --> M1(@t(), rand(4))

    # Based on properties of particular "structured agent" assigned to the transition,
    # we can update the attributes of the instance of a transition (such as probability of success).

    # Transition "Molecule 1" into "Molecule 2."
    # Update transition probability based on properties of "M1," 
    # which was assigned as a "resource" to the transition.
    ρ4,
    R1 + M1 --> M2(@t(), rand(4)),
    preAction => update_prob_transition(state, transition)

    5.0, R2 + M1 --> @structured(M2(@t(), rand(4)), :A)
    1.0, R2 + A --> @structured(M2(@t(), rand(4)), f_species(@transition))

    1.0, R2 + M1 --> @move(M1, :M2)
    1.0, R2 + M1 --> @move(M1, :C)
end

[slwu89]

I am wondering, is there a "functional" (i.e. not macro-based) way to set up this network? I just am thinking about if this network is going to be constructed from some other method (e.g. REST API in a web service, or by some operations to automatically build large networks, like pullbacks), it would be nice to also see how to do this in that other way.

[thevolatilebit]

I am wondering, is there a "functional" (i.e. not macro-based) way to set up this network? I just am thinking about if this network is going to be constructed from some other method (e.g. REST API in a web service, or by some operations to automatically build large networks, like pullbacks), it would be nice to also see how to do this in that other way.

Great question, @slwu89! The functionality of ReactiveDynamics is built around expressions, some of which are evaluated when an instance of a reaction network is created, while others are dynamically re-evaluated during the simulation. For example, initially we parsed the LHS and RHS into a set of reactants right at the time when the network was created. Now, we store the LHS and RHS of a transition as expressions and manipulate the expressions each time we try to "fire" that transition (because we can have something like LHS -> @choose((prob, A), (1-prob, B))).

However, note that you do not have to manually write everything inside a call to @ReactionNetworkSchema. You can actually set up the defining expression programmatically or write the equations as strings (which is simpler), parse them as an expression, and call make_ReactionNetwork(expr).

Also, consider

ReactiveDynamics.jl/tutorial/agents-integration/agents_integration.jl

Lines 163 to 191 in 4791a81

	# Species and initial values.
	name = "reaction_network"
	str_init = "@prob_init $name"
	for s in data["species"]
	str_init *= " $(s["name"])=$(s["initial"])"
	end
	function get_subline(d::Vector)
	if isempty(d)
	return "∅"
	elseif haskey(first(d), "probability")
	sublines = [
	"($(sd["probability"]), " *
	join(["$(s["q"]) * $(s["name"])" for s in sd["reactants"]], " + ") *
	")" for sd in d
	]
	return "@choose(" * join(sublines, ", ") * ")"
	else
	return join(["$(s["q"]) * $(s["name"])" for s in d], " + ")
	end
	end
	# Transitions.
	str_transitions = []
	for t in data["transitions"]
	line = t["rate"] * ", " * get_subline(t["from"]) * " --> " * get_subline(t["to"])
	line *= ", name => $(t["name"]), priority => $(t["priority"]), cycletime => $(t["duration"])"
	push!(str_transitions, line)
	end

[slwu89]

Thanks @thevolatilebit for the explanation. Would it be possible to somehow use a structured Petri net representation of a process to build a ReactiveDynamics model? For example, the model built here https://algebraicjulia.github.io/AlgebraicPetri.jl/dev/generated/covid/disease_strains/#Stratify-the-SIRD-strain-vax-and-simple-trip-models which has a large number of places/transitions, but actually is quite a small number compared to other Petri net models in common use.