Custom Movement Models

If one of epymorph’s built-in Movement Models does not suit your use-case, you can of course write your own! Define an appropriate set of classes which extend the MovementModel and MovementClause parent classes, implement the required attributes and methods in those classes, and you can use it when defining a RUME just like any other movement model.

In fact, there’s no difference between how we define built-in movement models and how you would define a custom model, so you may find it helpful to examine the source code of existing models. You can even use copy/paste/modify as your strategy to get started.

We’ll use the Centroids movement model as our example; here’s the annotated source code (with a few minor details stripped out for readability):

# A. implement MovementModel
class Centroids(MovementModel):
    # B. define the tau steps for the model
    steps = (1 / 3, 2 / 3)

    # C. define the list of clauses; instances of MovementClause
    #    NOTE: this uses Python's syntax for a tuple containing one item
    clauses = (CentroidsClause(),) 


# D. implement MovementClause
class CentroidsClause(MovementClause):

    # E. data requirements for the clause
    requirements = (
        AttributeDef(
            "centroid",
            CentroidType,
            Shapes.N,
            comment="The centroids for each node as (longitude, latitude) tuples.",
        ),
        AttributeDef(
            "phi",
            float,
            Shapes.Scalar,
            default_value=40.0,
            comment="Influences the distance that movers tend to travel.",
        ),
        AttributeDef(
            "commuter_proportion",
            float,
            Shapes.Scalar,
            default_value=0.1,
            comment="The proportion of the total population which commutes.",
        ),
    )

    # F. the clause predicate
    predicate = EveryDay()

    # G. the leaves rule
    leaves = TickIndex(step=0)

    # H. the returns rule
    returns = TickDelta(step=1, days=0)

    # I. the evaluate function
    def evaluate(
        self,
        tick: Tick,
        *,
        available: NDArray[SimDType],
    ) -> NDArray[np.int64]:
        comm_prop = self.data("commuter_proportion")
        n_commuters = np.floor(available * comm_prop).astype(SimDType)
        return self.rng.multinomial(n_commuters, self.dispersal_kernel)

    # J. (not required) any additional helper functions, etc.
    @cached_property
    def dispersal_kernel(self) -> NDArray[np.float64]:
        centroid = self.data("centroid")
        phi = self.data("phi")
        distance = pairwise_haversine(centroid)
        dist_over_phi = np.clip(distance / phi, a_min=None, a_max=700.0)
        prob = np.exp(-dist_over_phi)
        return row_normalize(prob)

MovementModel

An implementation of MovementModel (A) is a very thin wrapper which simply:

  • defines the tau steps for the whole model (B), and
  • lists the clauses (C).

Tau steps are defined a sequence of lengths. You can define any number of tau steps, but they must add up to 1. The tau steps will be referenced later by index, starting from zero as is Python’s convention.

You can define any number of movement clauses, they simply have to be instances of the MovementClause class.

MovementClause

Implementations of MovementClause (D) are more involved. You must define:

  • the data parameters users will need to supply (E),
  • when the clause is applied, with a predicate and a leaves rule (F, G),
  • when movers will return home (H), and
  • the function which computes the requested number of movers for all pairs of locations in the simulation’s GeoScope (I).

Parameter values are provided by the user when they construct a RUME. Your movement clauses declare which parameters are required: their name, their data type, their maximal data shape, and a description. They can also have default values. These parameters will be available to be used by your evaluate function. epymorph will throw an error if you try to use parameters that you did not declare in this way.

We discussed the predicate, leaves, and returns attributes in the previous chapter. There are built-in MovementPredicate implementations for common use-cases (EveryDay and DayIs). And of course your leaves and returns rules must identify valid tick indices, as defined by the MovementModel that uses this clause. There is also the special TickIndex value NEVER which you can use as your returns rule to indicate that movers do not regularly return to their origin, but instead will be merged with the resident population at the destination (e.g., for migratory populations).

The evaluate function will be evaluated when this movement clause is triggered, and should return a node-by-node (NxN) array of integers describing the requested movement from each location (first axis) and to each location (second axis). The function is called with two arguments. The first argument, tick is a Tick object containing information about the current simulation time step. It includes the date, which tau step index this is, the length of this tau step, which simulation day we’re on (zero-indexed from the start of the simulation; basically how many days have elapsed since the start of the simulation), and which tau step we’re on (zero-indexed from the start of the simulation; basically how many tau steps have elapsed). The second argument, available, is the current population available for movement at each simulation node. It’s provided as an N-shaped numpy array of integers.

Besides the function arguments, typically you will make use of data parameters, the random number generator, and other properties of the simulation by accessing the Context in the body of this function. (More about the context in a moment.)

Finally, movement clauses are of course normal Python classes. You are free to define additional methods and attributes if you like (J). In this example, our clause function needs to compute a matrix called a “dispersal kernel”. The values of this matrix don’t change during the simulation, so we can save time by only computing it once. It’s declared as a property using Python’s @cached_property decorator which will cause it to be cached. This sort of thing isn’t required; we do it just for the sake of efficiency.

Simulation Context

epymorph’s modular design gives us a lot of power and flexibility, but it also adds some complexity. In the case of movement models, for instance, we want our model to work for all possible GeoScopes, but we also need to compute a matrix whose size depends on the number of nodes in the GeoScope! We don’t “know” how many nodes are going to be in the GeoScope until the RUME is assembled right before the user runs a simulation. Hence, we don’t know how many nodes there will be when we define our movement model. The answer to this problem is the “Simulation Context”. It is a collection of properties that we can use when defining modular components that need to know information about the simulation that comes from the other modules.

The GeoScope defines the number of nodes, but the Movement Model needs that for its computations? It accesses the context. The TimeFrame defines the start and end dates, but a parameter function wants to define a seasonally-varying value? It accesses the context. The context is a way to defer the usage of a value until such a time as that value is known (usually when we’re running a simulation with a complete RUME).

Every class that inherits from the BaseSimulationFunction class (which includes MovementClause) can access the simulation context as properties and methods on the object (via the self reference). Those properties are:

  • scope: the simulation geo scope,
  • time_frame: the simulation time frame,
  • ipm: the compartment model,
  • rng: the random number generator, and
  • dim: the simulation’s dimensional information (which is merely a convenient repackaging of info from the rest of the context.)

And method:

  • data(attribute): the value of the named data requirement, which is given as a RUME parameter.

evaluate

Let’s examine in detail how Centroids implements evaluate:

def evaluate(
    self,
    tick: Tick,
    *,
    available: NDArray[SimDType],
) -> NDArray[np.int64]:
    comm_prop = self.data("commuter_proportion")
    n_commuters = np.floor(available * comm_prop).astype(SimDType)
    return self.rng.multinomial(n_commuters, self.dispersal_kernel)

First, the function signature. We get the current tick and available population as a function arguments if needed. In this case we ignore the tick. It may be useful to some models, but at least be aware that this function will be evaluated every time the clause applies, which could be as often as every single tau step in the simulation. Depending on your needs, you may want to return the same or a different result every time the function is called. The return value type is a numpy array of integers.

Next we use self.data(attribute) to access the values data parameters: in this case, “commuter_proportion”. Parameter values are returned in a form determined by the AttributeDef which declared the matching requirement. If an attribute is declared as an N-shaped array of integers, its value will be a one-axis numpy array of length N (the number of geo nodes). epymorph guarantees this. Even if the user were to provide a scalar value for it in the RUME, we internally type-cast the value and broadcast its shape into the declared form for the sake of normalizing these values for internal use. Commuter proportion is declared as a scalar float, so we will get a scalar numpy float for comm_prop.

Recall our goal is to produce a matrix of requested movement, so now we combine our available data using numpy operations to suit our purpose. To start with, we have the population and we have the proportion of people we would like to have commuting, it would be helpful to get the actual number of commuters per location. np.floor(available * comm_prop) takes care of that. .astype(SimDType) is just to be careful about data types (SimDType is the “official” integer type of epymorph calculations – which is currently 64-bit integers, np.int64.)

Next we will use the random number generator to introduce stochasticity into the movement model. Real human movement has variation: people call out sick, people travel for work, not all workers have regular 9-to-5, Monday-to-Friday schedules, etc. So we use the random number generator provided by the simulation context to perform a multinomial draw to slightly randomize how n_commuters get distributed to the locations in the simulation. We can take advantage of the fact that numpy’s multinomial method is vectorized. If the first argument (the number of samples to draw) is an N-shaped array of integers, and the second argument (the probabilities of each possible outcome) is an NxN-shaped array of floats, the result will be an NxN-shaped array of integers.

If we get a little creative with the syntax, you can illustrate it like this:

multinomial(
  # arg 1: number of people:   |  # arg 2: probability of traveling:
  [                            |  [
    leaving node 0,            |    [from node 0 to node 0, 1, 2, ...]
    leaving node 1,            |    [from node 1 to node 0, 1, 2, ...]
    leaving node 2,            |    [from node 2 to node 0, 1, 2, ...]
    ...                        |    ...
  ]                            |  ]
) ->
  # result: number of people requested to travel:
  [
    [from node 0 to node 0, 1, 2, ...]
    [from node 1 to node 0, 1, 2, ...]
    [from node 2 to node 0, 1, 2, ...]
    ...
  ]

And this produces exactly what we want: a matrix describing how many people should travel from every origin to every destination. Because we rely on data which is known to be “N-shaped”, for any possible value of N, we can be certain our movement model will work for arbitrary geo scopes, as long as all required parameters are provided.

I skipped over how dispersal_kernel is computed — that is, the probabilities that any given individual who is leaving a particular location will wind up in any other location. But if you refer back to the source code you’ll see that it follows a very similar form: it uses some data requirements to compute a matrix. In this case, the data are latitudes and longitudes of a centroid for each geo node (used to compute distances between each pair of locations) and a tuning parameter “phi” (loosely describing how far people are likely to travel). The details aren’t important, what matters is that it’s also merely a product of the information available in the simulation context and some numpy array operations.

WarningOnly use the context’s random number generator!

If you’re familiar with numpy, you might be tempted to create your own random number generator (RNG) instead of using the one from the simulation context. Any time you need randomness, you should use the context’s RNG.

Although using an external RNG probably wouldn’t cause any immediate errors, part of epymorph’s design is to produce consistent results when the same simulation is run twice using the same initial random seed. This is why the context includes an RNG instance. epymorph needs to carefully maintain the state of the generator, and both the number of and order of operations it’s used for, in order to guarantee result consistency. If you use an external source of randomness or use the context RNG in unexpected ways, this guarantee will break down.

By breaking down how the Centroids movement model is implemented, we hope this gives you the insight you need to write your own.