3.2.11.5. causality/Stochsim.hpp

Added in version 1.1.

The causality::Simulator class implements the Stochastic Simulation Algorithm, also known as the Gillespie algorithm, but with the ability to generate the underlying network as needed, similarly to how a derivation graph, dg::DG, can be created, and in particular the generation through graph transformation rules and strategies, Derivation Graph Strategies.

class causality::Simulator

The main class for performing stochastic simulations of chemical systems.

using ExpandNetwork = std::function<bool(std::shared_ptr<dg::Builder>, const std::vector<std::shared_ptr<graph::Graph>>&, const std::vector<std::shared_ptr<graph::Graph>>&)>

Simulator(LabelSettings labelSettings, const std::vector<std::shared_ptr<graph::Graph>> &graphDatabase, IsomorphismPolicy graphPolicy, ExpandNetwork expandNetwork, const std::vector<std::pair<std::shared_ptr<graph::Graph>, int>> &initialState, std::function<std::shared_ptr<DrawFunction>(const Marking&)> draw, std::function<double(double)> drawTime, bool withSetCompare)

Parameters:

• LabelSettings labelSettings – a settings object that will be passed to dg::DG::make().

• list[Graph] graphDatabase – a list of graphs that will be passed to dg::DG::make().

• IsomorphismPolicy graphPolicy – the policy that will be passed to dg::DG::make().

• expandNetwork –

  a callback which will be called when the simulator needs new hyperedge/reaction information. The callback will be given:

  1. the DG builder object for the underlying derivation graph.

  2. a list of graphs that were newly discovered in the last iteration.

  3. a list of all graphs in the current state.

  The callback must return a boolean which indicates whether it should be called again.

• initialState –

  the initial simulation state in terms of the number of copies of each graph/molecule. The graphs/molecules not mentioned are not considered part of the simulation yet, and are thus implicitly assumed to be 0.

  Caution

  If you use input flow through DrawMassAction then the input rate is only queried once a graph/molecule is known to the simulation, and you must thus mention them in initialState if they should be queried immediately in the simulation.

• draw – The simulator will initially create a dg::DG which is given to this function. It must then return a callable that is used in each simulation step to draw the next hyperedge/reaction.

• drawTime – The function to use for drawing the time increment in each simulation step. The function will be given the reactivity, and must return the time increment. If a null function is given, it will default to an instance of DrawTimeExponential.

• bool withSetCompare – Whether to skip network expansion if a larger subset of graphs has been used for expansion before.

Throws:

LogicError if draw is a null function.

Throws:

LogicError if draw(dg) returns a null function for the internally created derivation graph dg.

~Simulator()

std::shared_ptr<dg::DG> getDG() const

Returns:

the internal dg::DG underlying the simulation.

int getIteration() const

Returns:

the current iteration number. It starts at 0 and is incremented in the beginning of each iteration.

double getTime() const

Returns:

the current simulation time.

int state(dg::DG::Vertex v) const

int state(std::shared_ptr<graph::Graph> g) const

Returns:

the number of occurrences of the given vertex/graph in the current state.

const EventTrace &getTrace() const

Returns:

a copy of the event trace for the entire simulation.

bool isNetworkOpen() const

Returns:

whether the underlying network is still open for expansion. See also the keepNetworkOpen parameter of simulate().

void closeNetwork()

Manually close the network for expansion, instead of letting the simulate() method do it. The method does nothing if the network is already closed. If the user saved the dg::Builder object in the expansion callback, this method also does nothing.

void setOnIterationBegin(std::function<void(Simulator&)> f, int interval)

Set/remove a callback invoked in the very beginning of iterations of the simulation. just after the iteration count has been incremented. It is called at every interval th iteration, The callback is invoked with the simulator object as argument. Give a null function to remove the current callback.

Throws:

LogicError if interval is non-positive.

void setOnIterationEnd(std::function<bool(Simulator&)> f)

Set/remove the callback invoked in the very end of each iteration of the simulation, after the drawn action has been carried out and the time advanced. It is called with the simulator object, and it must return a boolean indicating whether to continue the simulation, i.e., true means continue.

void setOnDeadlock(std::function<void(Simulator&)> f)

Set/remove the callback invoked if there are no events out of the current state. It is called with the simulator object as argument.

void setOnExpand(std::function<void(Simulator&)> f)

Set/remove the callback invoked when the simulator is about to request events out of the current state. It is called with the simulator object as argument.

void setOnExpandAvoided(std::function<void(Simulator&)> f)

Set/remove the callback invoked when the simulator detected is already had all events out of the current state. It is called with the simulator object as argument.

const EventTrace &simulate(std::optional<double> time, bool advanceToEndTime, std::optional<int> iterations, bool keepNetworkOpen)

Start/continue the simulation.

Simulate an additional amount of time or number of iterations, whichever is reached first, or until no further events are possible (a deadlock).

Parameters:

• time – the additional amount of time to simulate, or std::nullopt for unbounded.

• advanceToEndTime – if a time bound is given and the simulation stops due to this bound, advance the current time to the time bound, instead of staying at the time of the last event.

• iterations – the additional number of iterations to simulate, or std::nullopt for unbounded.

• keepNetworkOpen – if false the internal dg::Builder object will be deleted before returning.

Returns:

getTrace()

class causality::Simulator::DrawTimeExponential

A shorthand for drawing time from an exponential distribution.

double operator()(double activitySum) const

Parameters:

activitySum – the total sum of activity in the system.

Returns:

\(\frac{-\ln r}{activitySum}\), where \(r\) is a random number in \([0, 1)\) drawn with rngUniformReal().

class causality::Simulator::DrawFunction

The base class that event drawing functions must inherit from.

In each iteration a causality::Simulator must draw the next event that should happen. How to do this drawing can be customized, but such a drawing function needs detailed information about the underlying network and as the network expands, the drawing function must be kept in sync with the network- The drawing function customization is therefore done when a slightly indirect manner. Instead of giving the drawing function directly, you give a function that can create a drawing function. That is, the causality::Simulator will create an internal causality::Marking and give it as argument to the function you give. Your function must then create an actual drawing function, which must inherit from causality::Simulator::DrawFunction and implement the appropriate methods.

For an example of a drawing function, see causality::Simulator::DrawMassAction::Function, and its creator function causality::Simulator::DrawMassAction, which is the one users interact with.

class Choice

The class used in the return type of causality::Simulator::DrawFunction::draw(). It is a more efficient representation of a causality::Action.

Choice() = default

Choice(Action a)

Choice(InputAction a)

Choice(OutputAction a)

Choice(EdgeAction a)

Construct either a null choice or an actual choice from an action.

std::optional<Action> asAction() const

Returns:

a conversion the choice into an action.

virtual ~DrawFunction()

virtual void syncSize() = 0

Called whenever the underlying derivation graph has changed size. If the drawing function has internal data structures, this method is where such data structures can be resized. The derivation graph must be given to this object by its creator.

virtual std::pair<Choice, double> draw() = 0

Called in order to draw the next event.

Returns:

the drawn action to take and a number indicating the activity of the system.

virtual void stateUpdated() = 0

Updates internal caches after the state has been updated. This function must be invoked inbetween calls to draw.

class causality::Simulator::DrawMassAction

A creator for a drawing function implementing the law of mass action.

It supports assigning a rate for input actions, output actions, and reactions, by taking a callback (or constant) for each type. To avoid the overhead of calling these callbacks in each iteration, a returned rate can be cached. Therefore, the return value of each callback (or the constant of each type) is a pair with the first entry being the rate, and the second entry a boolean indicating whether rate should be cached.

Each of the rate function can also be set to a null function, which means a default rate is used: input rate 0.0, reaction rate 1.0, output rate 0.0.

DrawMassAction(std::function<std::pair<double, bool>(dg::DG::Vertex)> inputRate, std::function<std::pair<double, bool>(dg::DG::HyperEdge)> reactionRate, std::function<std::pair<double, bool>(dg::DG::Vertex)> outputRate, const std::string &implementationName)

Parameters:

• inputRate – the rate used for pseudo-reactions for creating molecules.

• reactionRate – the rate used for each reaction in the system.

• outputRate – the rate used for pseudo-reactions for destroying molecules.

• implementationName – the name of the internal implementation to use.

~DrawMassAction()

std::shared_ptr<Function> operator()(const Marking &state)

Parameters:

state – the state of the underlying simulation. It must be kept alive as long as the returned drawing function is kept alive.

Returns:

a drawing function implementing the law of mass action.

class causality::Simulator::DrawMassAction::Function

A helper class for performing stochastic simulations where events are drawn according to the law of mass action. Importantly, if the underlying derivation graph is enlarged then syncSize() must be called before calling draw().

Function(const Marking &state, std::function<std::pair<double, bool>(dg::DG::Vertex)> inputRate, std::function<std::pair<double, bool>(dg::DG::HyperEdge)> reactionRate, std::function<std::pair<double, bool>(dg::DG::Vertex)> outputRate, const std::string &implementation)

Construct a new instance, that will draw based on the given state. The return value of the rate callbacks must be 1) the rate and 2) a boolean telling whether the library should cache the rate. If true no more calls with the same argument will be made. Each of the rate functions may be an empty std::function (i.e., default constructed), which means a default rate is used: input rate 0.0, reaction rate 1.0, output rate 0.0.

~Function()

virtual void syncSize() override

virtual std::pair<DrawFunction::Choice, double> draw() override

virtual void stateUpdated() override