graph_tool.dynamics - Dynamical processes

This module contains implementations of some often-studied dynamical processes that take place on networks.

Summary

Discrete-time dynamics

DiscreteStateBase

Base state for discrete-time dynamics.

EpidemicStateBase

Base state for epidemic dynamics.

SIState

SI compartmental epidemic model.

SISState

SIS compartmental epidemic model.

SIRState

SIR compartmental epidemic model.

SIRSState

SIRS compartmental epidemic model.

VoterState

Generalized q-state voter model dynamics.

MajorityVoterState

Generalized q-state majority voter model dynamics.

BinaryThresholdState

Generalized binary threshold dynamics.

IsingGlauberState

Glauber dynamics of the Ising model.

CIsingGlauberState

Glauber dynamics of the continuous Ising model.

IsingMetropolisState

Metropolis-Hastings dynamics of the Ising model.

PottsGlauberState

Glauber dynamics of the Potts model.

PottsMetropolisState

Metropolis-Hastings dynamics of the Potts model.

AxelrodState

Axelrod’s culture dissemination model.

BooleanState

Boolean network dynamics.

KirmanState

Kirman’s “ant colony” model.

Continuous-time dynamics

ContinuousStateBase

Base state for continuous-time dynamics.

KuramotoState

The Kuramoto model.

Contents

class graph_tool.dynamics.DiscreteStateBase(g, make_state, params, s=None, stype='int32_t')[source]

Bases: object

Base state for discrete-time dynamics. This class it not meant to be instantiated directly.

copy()[source]

Return a copy of the state.

get_state()[source]

Returns the internal VertexPropertyMap with the current state.

get_active()[source]

Returns list of “active” nodes, for states where this concept is used.

reset_active()[source]

Resets list of “active” nodes, for states where this concept is used.

iterate_sync(niter=1)[source]

Updates nodes synchronously (i.e. a full “sweep” of all nodes in parallel), niter number of times. This function returns the number of nodes that changed state.

If enabled during compilation, this algorithm runs in parallel (i.e. using more than one thread.)

iterate_async(niter=1)[source]

Updates nodes asynchronously (i.e. single vertex chosen randomly), niter number of times. This function returns the number of nodes that changed state.

class graph_tool.dynamics.EpidemicStateBase(g, v0=None, s=None)[source]

Bases: graph_tool.dynamics.DiscreteStateBase

Base state for epidemic dynamics. This class it not meant to be instantiated directly.

class graph_tool.dynamics.SIState(g, beta=1.0, r=0, exposed=False, epsilon=0.1, v0=None, s=None)[source]

Bases: graph_tool.dynamics.EpidemicStateBase

SI compartmental epidemic model.

Parameters
gGraph

Graph to be used for the dynamics

betafloat (optional, default: 1.)

Transmission probability.

rfloat (optional, default: 0.)

Spontaneous infection probability.

exposedboolean (optional, default: False)

If True, an SEI model is simulated, with an additional “exposed” state.

epsilonfloat (optional, default: .1)

Susceptible to exposed transition probability. This only has an effect if exposed=True.

v0int or Vertex (optional, default: None)

Initial infected vertex. If not provided, and if the global state is also not provided via paramter s, a random vertex will be chosen.

sVertexPropertyMap (optional, default: None)

Initial global state. If not provided, all vertices will be initialized to the susceptible state.

Notes

This implements an SI epidemic process [pastor-satorras-epidemic-2015], where nodes in the susceptible state (value 0) are infected by neighbours in the infected state (value 1).

If a node \(i\) is updated at time \(t\), the transition probabilities from state \(s_i(t)\) to state \(s_i(t+1)\) are given as follows:

  1. If \(s_i(t) = 0\), we have \(s_i(t+1) = 1\) with probability
    \[(1-r)\left(1-\prod_j(1-\beta)^{A_{ij}\delta_{s_j(t),1}}\right) + r,\]

    otherwise \(s_i(t+1) = 0\).

  2. If \(s_i(t) = 1\), we have \(s_i(t+1) = 1\) with probability 1.

If the option exposed == True is given, then the states transit first from 0 to -1 (exposed) with probability given by 1. above, and then finally from -1 to 1 with probability \(\epsilon\).

References

pastor-satorras-epidemic-2015(1,2)

Romualdo Pastor-Satorras, Claudio Castellano, Piet Van Mieghem, and Alessandro Vespignani, “Epidemic processes in complex networks”, Rev. Mod. Phys. 87, 925 (2015) DOI: 10.1103/RevModPhys.87.925 [sci-hub, @tor], arXiv: 1408.2701

Examples

>>> g = gt.collection.data["pgp-strong-2009"]
>>> state = gt.SIState(g, beta=0.01)
>>> X = []
>>> for t in range(1000):
...     ret = state.iterate_sync()
...     X.append(state.get_state().fa.sum())

>>> figure(figsize=(6, 4))
<...>
>>> plot(X)
[...]
>>> xlabel(r"Time")
Text(...)
>>> ylabel(r"Infected nodes")
Text(...)
>>> tight_layout()
>>> savefig("SI.pdf")
_images/SI.svg

Number of infected nodes vs. time for an SI dynamics.

class graph_tool.dynamics.SISState(g, beta=1.0, gamma=0.1, r=0, exposed=False, epsilon=0.1, v0=None, s=None)[source]

Bases: graph_tool.dynamics.DiscreteStateBase

SIS compartmental epidemic model.

Parameters
gGraph

Graph to be used for the dynamics

betafloat (optional, default: 1.)

Transmission probability.

gammafloat (optional, default: .1)

Recovery probability.

rfloat (optional, default: 0.)

Spontaneous infection probability.

exposedboolean (optional, default: False)

If True, an SEIS model is simulated, with an additional “exposed” state.

epsilonfloat (optional, default: .1)

Susceptible to exposed transition probability. This only has an effect if exposed=True.

v0int or Vertex (optional, default: None)

Initial infected vertex. If not provided, and if the global state is also not provided via paramter s, a random vertex will be chosen.

sVertexPropertyMap (optional, default: None)

Initial global state. If not provided, all vertices will be initialized to the susceptible state.

Notes

This implements an SIS epidemic process [pastor-satorras-epidemic-2015], where nodes in the susceptible state (value 0) are infected by neighbours in the infected state (value 1), which can then eventually recover to a susceptible state.

If a node \(i\) is updated at time \(t\), the transition probabilities from state \(s_i(t)\) to state \(s_i(t+1)\) are given as follows:

  1. If \(s_i(t) = 0\), we have \(s_i(t+1) = 1\) with probability
    \[(1-r)\left(1-\prod_j(1-\beta)^{A_{ij}\delta_{s_j(t),1}}\right) + r,\]

    otherwise \(s_i(t+1) = 0\).

  2. If \(s_i(t) = 1\), we have \(s_i(t+1) = 0\) with probability \(\gamma\), or \(s_i(t+1) = 1\) with probability \(1-\gamma\).

If the option exposed == True is given, then the states transit first from 0 to -1 (exposed) with probability given by 1. above, and then finally from -1 to 1 with probability \(\epsilon\).

References

pastor-satorras-epidemic-2015(1,2)

Romualdo Pastor-Satorras, Claudio Castellano, Piet Van Mieghem, and Alessandro Vespignani, “Epidemic processes in complex networks”, Rev. Mod. Phys. 87, 925 (2015) DOI: 10.1103/RevModPhys.87.925 [sci-hub, @tor], arXiv: 1408.2701

Examples

>>> g = gt.collection.data["pgp-strong-2009"]
>>> state = gt.SISState(g, beta=0.01, gamma=0.007)
>>> X = []
>>> for t in range(1000):
...     ret = state.iterate_sync()
...     X.append(state.get_state().fa.sum())

>>> figure(figsize=(6, 4))
<...>
>>> plot(X)
[...]
>>> xlabel(r"Time")
Text(...)
>>> ylabel(r"Infected nodes")
Text(...)
>>> tight_layout()
>>> savefig("SI.pdf")
_images/SIS.svg

Number of infected nodes vs. time for an SIS dynamics.

class graph_tool.dynamics.SIRState(g, beta=1.0, gamma=0.1, r=0, exposed=False, epsilon=0.1, v0=None, s=None)[source]

Bases: graph_tool.dynamics.DiscreteStateBase

SIR compartmental epidemic model.

Parameters
gGraph

Graph to be used for the dynamics

betafloat (optional, default: 1.)

Transmission probability.

gammafloat (optional, default: .1)

Recovery probability.

rfloat (optional, default: 0.)

Spontaneous infection probability.

exposedboolean (optional, default: False)

If True, an SEIR model is simulated, with an additional “exposed” state.

epsilonfloat (optional, default: .1)

Susceptible to exposed transition probability. This only has an effect if exposed=True.

v0int or Vertex (optional, default: None)

Initial infected vertex. If not provided, and if the global state is also not provided via paramter s, a random vertex will be chosen.

sVertexPropertyMap (optional, default: None)

Initial global state. If not provided, all vertices will be initialized to the susceptible state.

Notes

This implements an SIR epidemic process [pastor-satorras-epidemic-2015], where nodes in the susceptible state (value 0) are infected by neighbours in the infected state (value 1), which can then eventually recover to a recovered state (value 2).

If a node \(i\) is updated at time \(t\), the transition probabilities from state \(s_i(t)\) to state \(s_i(t+1)\) are given as follows:

  1. If \(s_i(t) = 0\), we have \(s_i(t+1) = 1\) with probability
    \[(1-r)\left(1-\prod_j(1-\beta)^{A_{ij}\delta_{s_j(t),1}}\right) + r,\]

    otherwise \(s_i(t+1) = 0\).

  2. If \(s_i(t) = 1\), we have \(s_i(t+1) = 2\) with probability \(\gamma\), or \(s_i(t+1) = 1\) with probability \(1-\gamma\).

If the option exposed == True is given, then the states transit first from 0 to -1 (exposed) with probability given by 1. above, and then finally from -1 to 1 with probability \(\epsilon\).

References

pastor-satorras-epidemic-2015(1,2)

Romualdo Pastor-Satorras, Claudio Castellano, Piet Van Mieghem, and Alessandro Vespignani, “Epidemic processes in complex networks”, Rev. Mod. Phys. 87, 925 (2015) DOI: 10.1103/RevModPhys.87.925 [sci-hub, @tor], arXiv: 1408.2701

Examples

>>> g = gt.collection.data["pgp-strong-2009"]
>>> state = gt.SIRState(g, beta=0.01, gamma=0.0025)
>>> S, X, R = [], [], []
>>> for t in range(2000):
...     ret = state.iterate_sync()
...     s = state.get_state().fa
...     S.append((s == 0).sum())
...     X.append((s == 1).sum())
...     R.append((s == 2).sum())

>>> figure(figsize=(6, 4))
<...>
>>> plot(S, label="Susceptible")
[...]
>>> plot(X, label="Infected")
[...]
>>> plot(R, label="Recovered")
[...]
>>> xlabel(r"Time")
Text(...)
>>> ylabel(r"Number of nodes")
Text(...)
>>> legend(loc="best")
<...>
>>> tight_layout()
>>> savefig("SIR.pdf")
_images/SIR.svg

Number of susceptible, infected, and recovered nodes vs. time for an SIR dynamics.

class graph_tool.dynamics.SIRSState(g, beta=1, gamma=0.1, mu=0.1, r=0, exposed=False, epsilon=0.1, v0=None, s=None)[source]

Bases: graph_tool.dynamics.DiscreteStateBase

SIRS compartmental epidemic model.

Parameters
gGraph

Graph to be used for the dynamics

betafloat (optional, default: 1.)

Transmission probability.

gammafloat (optional, default: .1)

I to R recovery probability.

mufloat (optional, default: .1)

R to S recovery probability.

rfloat (optional, default: 0.)

Spontaneous infection probability.

exposedboolean (optional, default: False)

If True, an SEIRS model is simulated, with an additional “exposed” state.

epsilonfloat (optional, default: .1)

Susceptible to exposed transition probability. This only has an effect if exposed=True.

v0int or Vertex (optional, default: None)

Initial infected vertex. If not provided, and if the global state is also not provided via paramter s, a random vertex will be chosen.

sVertexPropertyMap (optional, default: None)

Initial global state. If not provided, all vertices will be initialized to the susceptible state.

Notes

This implements an SIRS epidemic process [pastor-satorras-epidemic-2015], where nodes in the susceptible state (value 0) are infected by neighbours in the infected state (value 1), which can then eventually recover to a recovered state (value 2), and finally back to the susceptible state.

If a node \(i\) is updated at time \(t\), the transition probabilities from state \(s_i(t)\) to state \(s_i(t+1)\) are given as follows:

  1. If \(s_i(t) = 0\), we have \(s_i(t+1) = 1\) with probability
    \[(1-r)\left(1-\prod_j(1-\beta)^{A_{ij}\delta_{s_j(t),1}}\right) + r,\]

    otherwise \(s_i(t+1) = 0\).

  2. If \(s_i(t) = 1\), we have \(s_i(t+1) = 2\) with probability \(\gamma\), or \(s_i(t+1) = 1\) with probability \(1-\gamma\).

  3. If \(s_i(t) = 2\), we have \(s_i(t+1) = 1\) with probability \(\mu\), or \(s_i(t+1) = 2\) with probability \(1-\mu\).

If the option exposed == True is given, then the states transit first from 0 to -1 (exposed) with probability given by 1. above, and then finally from -1 to 1 with probability \(\epsilon\).

References

pastor-satorras-epidemic-2015(1,2)

Romualdo Pastor-Satorras, Claudio Castellano, Piet Van Mieghem, and Alessandro Vespignani, “Epidemic processes in complex networks”, Rev. Mod. Phys. 87, 925 (2015) DOI: 10.1103/RevModPhys.87.925 [sci-hub, @tor], arXiv: 1408.2701

Examples

>>> g = gt.collection.data["pgp-strong-2009"]
>>> state = gt.SIRSState(g, beta=0.2, gamma=0.025, mu=0.02)
>>> S, X, R = [], [], []
>>> for t in range(2000):
...     ret = state.iterate_sync()
...     s = state.get_state().fa
...     S.append((s == 0).sum())
...     X.append((s == 1).sum())
...     R.append((s == 2).sum())

>>> figure(figsize=(6, 4))
<...>
>>> plot(S, label="Susceptible")
[...]
>>> plot(X, label="Infected")
[...]
>>> plot(R, label="Recovered")
[...]
>>> xlabel(r"Time")
Text(...)
>>> ylabel(r"Number of nodes")
Text(...)
>>> legend(loc="best")
<...>
>>> tight_layout()
>>> savefig("SIRS.pdf")
_images/SIRS.svg

Number of susceptible, infected, and recovered nodes vs. time for an SIRS dynamics.

class graph_tool.dynamics.VoterState(g, q=2, r=0.0, s=None)[source]

Bases: graph_tool.dynamics.DiscreteStateBase

Generalized q-state voter model dynamics.

Parameters
gGraph

Graph to be used for the dynamics

qint (optional, default: 2)

Number of opinions.

rfloat (optional, default: 0.)

Random opinion probability.

sVertexPropertyMap (optional, default: None)

Initial global state. If not provided, a random state will be chosen.

Notes

This implements the voter model dynamics [clifford-model-1973] [holley-ergodic-1075] on a network.

If a node \(i\) is updated at time \(t\), the transition probabilities from state \(s_i(t)\) to state \(s_i(t+1)\) are given as follows:

1. With a probability \(r\) one of the \(q\) opinions, \(x\), is chosen uniformly at random, and assigned to \(i\), i.e. \(s_i(t+1) = x\).

2. Otherwise, a random (in-)neighbour \(j\) is chosen. and its opinion is copied, i.e. \(s_i(t+1) = s_j(t)\).

References

clifford-model-1973(1,2)

Clifford, P., Sudbury, A., “A model for spatial conflict”, Biometrika 60, 581–588 (1973). DOI: 10.1093/biomet/60.3.581 [sci-hub, @tor].

holley-ergodic-1075(1,2)

Holley, R. A., Liggett, T. M., “Ergodic Theorems for Weakly Interacting Infinite Systems and the Voter Model”, Ann. Probab. 3, 643–663 (1975). DOI: 10.1214/aop/1176996306 [sci-hub, @tor].

Examples

>>> g = gt.collection.data["pgp-strong-2009"]
>>> state = gt.VoterState(g, q=4)
>>> x = [[] for r in range(4)]
>>> for t in range(2000):
...     ret = state.iterate_sync()
...     s = state.get_state().fa
...     for r in range(4):
...         x[r].append((s == r).sum())
>>> figure(figsize=(6, 4))
<...>
>>> for r in range(4):
...     plot(x[r], label="Opinion %d" % r)
[...]
>>> xlabel(r"Time")
Text(...)
>>> ylabel(r"Number of nodes")
Text(...)
>>> legend(loc="best")
<...>
>>> tight_layout()
>>> savefig("voter.pdf")
_images/voter.svg

Number of nodes with a given opinion vs. time for a voter model dynamics with \(q=4\) opinions.

class graph_tool.dynamics.MajorityVoterState(g, q=2, r=0, s=None)[source]

Bases: graph_tool.dynamics.DiscreteStateBase

Generalized q-state majority voter model dynamics.

Parameters
gGraph

Graph to be used for the dynamics

qint (optional, default: 2)

Number of opinions.

rfloat (optional, default: 0.)

Random opinion probability.

sVertexPropertyMap (optional, default: None)

Initial global state. If not provided, a random state will be chosen.

Notes

This implements the majority voter model dynamics [oliveira-isotropic-1992] on a network.

If a node \(i\) is updated at time \(t\), the transition probabilities from state \(s_i(t)\) to state \(s_i(t+1)\) are given as follows:

1. With a probability \(r\) one of the \(q\) opinions, \(x\), is chosen uniformly at random, and assigned to \(i\), i.e. \(s_i(t+1) = x\).

2. Otherwise, the majority opinion \(x\) held by all (in-)neighbours of \(i\) is chosen. In case of a tie between two or more opinions, a random choice between them is made. The chosen opinion is then copied, i.e. \(s_i(t+1) = x\).

References

oliveira-isotropic-1992(1,2)

de Oliveira, M.J., “Isotropic majority-vote model on a square lattice”, J Stat Phys 66: 273 (1992). DOI: 10.1007/BF01060069 [sci-hub, @tor].

Examples

>>> g = gt.collection.data["pgp-strong-2009"]
>>> state = gt.MajorityVoterState(g, q=4)
>>> x = [[] for r in range(4)]
>>> for t in range(2000):
...     ret = state.iterate_async(niter=g.num_vertices())
...     s = state.get_state().fa
...     for r in range(4):
...         x[r].append((s == r).sum())
>>> figure(figsize=(6, 4))
<...>
>>> for r in range(4):
...     plot(x[r], label="Opinion %d" % r)
[...]
>>> xlabel(r"Time")
Text(...)
>>> ylabel(r"Number of nodes")
Text(...)
>>> legend(loc="best")
<...>
>>> tight_layout()
>>> savefig("majority-voter.pdf")
_images/majority-voter.svg

Number of nodes with a given opinion vs. time for a majority voter model dynamics with \(q=4\) opinions.

class graph_tool.dynamics.BinaryThresholdState(g, w=1.0, h=0.5, r=0.0, s=None)[source]

Bases: graph_tool.dynamics.DiscreteStateBase

Generalized binary threshold dynamics.

Parameters
gGraph

Graph to be used for the dynamics

wEdgePropertyMap or float (optional, default: 1.)

Edge weights. If a scalar is provided, it’s used for all edges.

hfloat (optional, default: .5)

Relative threshold value.

rfloat (optional, default: 0.)

Input random flip probability.

sVertexPropertyMap (optional, default: None)

Initial global state. If not provided, a random state will be chosen.

Notes

This implements a Boolean threshold model on a network.

If a node \(i\) is updated at time \(t\), the transition to state \(s_i(t+1)\) is given by

\[\begin{split}s_i(t+1) = \begin{cases} 1, & \text{ if } \sum_jA_{ij}w_{ij}\hat s_j(t) > h k_i,\\ 0, & \text{ otherwise.} \end{cases}\end{split}\]

where \(k_i=\sum_jA_{ij}\) and \(\hat s_i(t)\) are the flipped inputs sampled with probability

\[P(\hat s_i(t)|s_i(t)) = r^{1-\delta_{\hat s_i(t),s_i(t)}}(1-r)^{\delta_{\hat s_i(t),s_i(t)}}.\]

Examples

>>> g = gt.GraphView(gt.collection.data["polblogs"], directed=False)
>>> gt.remove_parallel_edges(g)
>>> g = gt.extract_largest_component(g, prune=True)
>>> state = gt.BinaryThresholdState(g, r=0.25)
>>> ret = state.iterate_sync(niter=1000)
>>> gt.graph_draw(g, g.vp.pos, vertex_fill_color=state.s, vcmap=cm.bone,
...               output_size=(700,400), output="binary-threshold.pdf")
<...>
_images/binary-threshold.svg

State of a binary threshold dynamics on a political blog network.

class graph_tool.dynamics.IsingGlauberState(g, beta=1.0, w=1.0, h=0.0, s=None)[source]

Bases: graph_tool.dynamics.DiscreteStateBase

Glauber dynamics of the Ising model.

Parameters
gGraph

Graph to be used for the dynamics

betafloat (optional, default: 1.)

Inverse temperature.

wEdgePropertyMap or float (optional, default: 1.)

Edge interaction strength. If a scalar is provided, it’s used for all edges.

hVertexPropertyMap or float (optional, default: 0.)

Vertex local field. If a scalar is provided, it’s used for all vertices.

sVertexPropertyMap (optional, default: None)

Initial global state. If not provided, a random state will be chosen.

Notes

This implements the Glauber dynamics of the Ising model [ising-model] on a network.

If a node \(i\) is updated at time \(t\), the transition to state \(s_i(t+1) \in \{-1,+1\}\) is done with probability

\[P(s_i(t+1)|\boldsymbol s(t)) = \frac{\exp(\beta s_i(t+1)\sum_jA_{ij}w_{ij}s_j(t) + h_is_i(t+1))} {2\cosh(\beta\sum_jA_{ij}w_{ij}s_j(t) + h_i)}.\]

References

ising-model(1,2)

https://en.wikipedia.org/wiki/Ising_model

Examples

>>> g = gt.GraphView(gt.collection.data["polblogs"], directed=False)
>>> gt.remove_parallel_edges(g)
>>> g = gt.extract_largest_component(g, prune=True)
>>> state = gt.IsingGlauberState(g, beta=.05)
>>> ret = state.iterate_async(niter=1000 * g.num_vertices())
>>> gt.graph_draw(g, g.vp.pos, vertex_fill_color=state.s, vcmap=cm.bone,
...               output_size=(700,400), output="glauber-ising.pdf")
<...>
_images/glauber-ising.svg

State of a Glauber Ising dynamics on a political blog network.

class graph_tool.dynamics.CIsingGlauberState(g, beta=1.0, w=1.0, h=0.0, s=None)[source]

Bases: graph_tool.dynamics.DiscreteStateBase

Glauber dynamics of the continuous Ising model.

Parameters
gGraph

Graph to be used for the dynamics

betafloat (optional, default: 1.)

Inverse temperature.

wEdgePropertyMap or float (optional, default: 1.)

Edge interaction strength. If a scalar is provided, it’s used for all edges.

hVertexPropertyMap or float (optional, default: 0.)

Vertex local field. If a scalar is provided, it’s used for all vertices.

sVertexPropertyMap (optional, default: None)

Initial global state. If not provided, a random state will be chosen.

Notes

This implements the Glauber dynamics of the continuous Ising model [ising-model] on a network.

If a node \(i\) is updated at time \(t\), the transition to state \(s_i(t+1) \in [-1,+1]\) is done with probability density

\[P(s_i(t+1)|\boldsymbol s(t)) = \frac{\exp(\beta s_i(t+1)\sum_jA_{ij}w_{ij}s_j(t) + h_is_i(t+1))} {Z(\beta\sum_jA_{ij}w_{ij}s_j(t) + h_i)},\]

with \(Z(x) = 2\sinh(x)/x\).

References

ising-model(1,2)

https://en.wikipedia.org/wiki/Ising_model

Examples

>>> g = gt.GraphView(gt.collection.data["polblogs"], directed=False)
>>> gt.remove_parallel_edges(g)
>>> g = gt.extract_largest_component(g, prune=True)
>>> state = gt.CIsingGlauberState(g, beta=.2)
>>> ret = state.iterate_async(niter=1000 * g.num_vertices())
>>> gt.graph_draw(g, g.vp.pos, vertex_fill_color=state.s, vcmap=cm.magma,
...               output_size=(700,400), output="glauber-cising.pdf")
<...>
_images/glauber-cising.svg

State of a continuous Glauber Ising dynamics on a political blog network.

class graph_tool.dynamics.IsingMetropolisState(g, beta=1, w=1, h=0, s=None)[source]

Bases: graph_tool.dynamics.DiscreteStateBase

Metropolis-Hastings dynamics of the Ising model.

Parameters
gGraph

Graph to be used for the dynamics

betafloat (optional, default: 1.)

Inverse temperature.

wEdgePropertyMap or float (optional, default: 1.)

Edge interaction strength. If a scalar is provided, it’s used for all edges.

hVertexPropertyMap or float (optional, default: 0.)

Vertex local field. If a scalar is provided, it’s used for all vertices.

sVertexPropertyMap (optional, default: None)

Initial global state. If not provided, a random state will be chosen.

Notes

This implements the Metropolis-Hastings dynamics [metropolis-equations-1953] [hastings-monte-carlo-1970] of the Ising model [ising-model] on a network.

If a node \(i\) is updated at time \(t\), the transition to state \(s_i(t+1) = -s_i(t)\) is done with probability

\[\min\left\{1, \exp\left[-2s_i(t)\left(h_i + \beta\sum_jA_{ij}w_{ij}s_j(t)\right)\right]\right\}\]

otherwise we have \(s_i(t+1) = s_i(t)\).

References

ising-model(1,2)

https://en.wikipedia.org/wiki/Ising_model

metropolis-equations-1953(1,2)

Metropolis, N., A.W. Rosenbluth, M.N. Rosenbluth, A.H. Teller, and E. Teller, “Equations of State Calculations by Fast Computing Machines,” Journal of Chemical Physics, 21, 1087–1092 (1953). DOI: 10.1063/1.1699114 [sci-hub, @tor]

hastings-monte-carlo-1970(1,2)

Hastings, W.K., “Monte Carlo Sampling Methods Using Markov Chains and Their Applications,” Biometrika, 57, 97–109, (1970). DOI: 10.1093/biomet/57.1.97 [sci-hub, @tor]

Examples

>>> g = gt.GraphView(gt.collection.data["polblogs"], directed=False)
>>> gt.remove_parallel_edges(g)
>>> g = gt.extract_largest_component(g, prune=True)
>>> state = gt.IsingMetropolisState(g, beta=.1)
>>> ret = state.iterate_async(niter=1000 * g.num_vertices())
>>> gt.graph_draw(g, g.vp.pos, vertex_fill_color=state.s, vcmap=cm.bone,
...               output_size=(700,400), output="metropolis-ising.pdf")
<...>
_images/metropolis-ising.svg

State of a Metropolis-Hastings Ising dynamics on a political blog network.

class graph_tool.dynamics.PottsGlauberState(g, f, w=1, h=0, s=None)[source]

Bases: graph_tool.dynamics.DiscreteStateBase

Glauber dynamics of the Potts model.

Parameters
gGraph

Graph to be used for the dynamics

flist of lists or two-dimensional ndarray

Matrix of interactions between spin values, of dimension \(q\times q\), where \(q\) is the number of spins.

wEdgePropertyMap or float (optional, default: 1.)

Edge interaction strength. If a scalar is provided, it’s used for all edges.

hVertexPropertyMap or iterable or float (optional, default: 0.)

Vertex local field. If an iterable is provided, it will be used as the field for all vertices. If a scalar is provided, it will be used for all spins values and vertices.

sVertexPropertyMap (optional, default: None)

Initial global state. If not provided, a random state will be chosen.

Notes

This implements the Glauber dynamics of the Potts model [potts-model] on a network.

If a node \(i\) is updated at time \(t\), the transition to state \(s_i(t+1) \in \{0,\dots,q-1\}\) is done with probability

\[P(s_i(t+1)|\boldsymbol s(t)) \propto \exp\left(\sum_jA_{ij}w_{ij}f_{s_i(t+1), s_j(t)} + h^{(i)}_{s_i(t+1)}\right)\]

References

potts-model(1,2)

https://en.wikipedia.org/wiki/Potts_model

Examples

>>> g = gt.GraphView(gt.collection.data["polblogs"], directed=False)
>>> gt.remove_parallel_edges(g)
>>> g = gt.extract_largest_component(g, prune=True)
>>> f = np.eye(4) * 0.1
>>> state = gt.PottsGlauberState(g, f)
>>> ret = state.iterate_async(niter=1000 * g.num_vertices())
>>> gt.graph_draw(g, g.vp.pos, vertex_fill_color=state.s,
...               output_size=(700,400), output="glauber-potts.pdf")
<...>
_images/glauber-potts.svg

State of a Glauber Potts dynamics with \(q=4\) on a political blog network.

class graph_tool.dynamics.PottsMetropolisState(g, f, w=1, h=0, s=None)[source]

Bases: graph_tool.dynamics.DiscreteStateBase

Metropolis-Hastings dynamics of the Potts model.

Parameters
gGraph

Graph to be used for the dynamics

flist of lists or two-dimensional ndarray

Matrix of interactions between spin values, of dimension \(q\times q\), where \(q\) is the number of spins.

wEdgePropertyMap or float (optional, default: 1.)

Edge interaction strength. If a scalar is provided, it’s used for all edges.

hVertexPropertyMap or iterable or float (optional, default: 0.)

Vertex local field. If an iterable is provided, it will be used as the field for all vertices. If a scalar is provided, it will be used for all spins values and vertices.

sVertexPropertyMap (optional, default: None)

Initial global state. If not provided, a random state will be chosen.

Notes

This implements the Metropolis-Hastings dynamics [metropolis-equations-1953] [hastings-monte-carlo-1970] of the Potts model [potts-model] on a network.

If a node \(i\) is updated at time \(t\), the transition to state \(s_i(t+1) \in \{0,\dots,q-1\}\) is done as follows:

  1. A spin value \(r\) is sampled uniformly at random from the set \(\{0,\dots,q-1\}\).

  2. The transition \(s_i(t+1)=r\) is made with probability

    \[\min\left[1, \exp\left(\sum_jA_{ij}w_{ij}(f_{r, s_j(t)}-f_{s_i(t), s_j(t)}) + h^{(i)}_{r} - h^{(i)}_{s_i(t)}\right)\right]\]

    otherwise we have \(s_i(t+1)=s_i(t)\).

References

potts-model(1,2)

https://en.wikipedia.org/wiki/Potts_model

metropolis-equations-1953(1,2)

Metropolis, N., A.W. Rosenbluth, M.N. Rosenbluth, A.H. Teller, and E. Teller, “Equations of State Calculations by Fast Computing Machines,” Journal of Chemical Physics, 21, 1087–1092 (1953). DOI: 10.1063/1.1699114 [sci-hub, @tor]

hastings-monte-carlo-1970(1,2)

Hastings, W.K., “Monte Carlo Sampling Methods Using Markov Chains and Their Applications,” Biometrika, 57, 97–109, (1970) DOI: 10.1093/biomet/57.1.97 [sci-hub, @tor]

Examples

>>> g = gt.GraphView(gt.collection.data["polblogs"], directed=False)
>>> gt.remove_parallel_edges(g)
>>> g = gt.extract_largest_component(g, prune=True)
>>> f = np.eye(4) * 0.1
>>> state = gt.PottsGlauberState(g, f)
>>> ret = state.iterate_async(niter=1000 * g.num_vertices())
>>> gt.graph_draw(g, g.vp.pos, vertex_fill_color=state.s,
...               output_size=(700,400), output="metropolis-potts.pdf")
<...>
_images/metropolis-potts.svg

State of a Metropolis-Hastings Potts dynamics with \(q=4\) on a political blog network.

class graph_tool.dynamics.KirmanState(g, d=0.1, c1=0.001, c2=0.001, s=None)[source]

Bases: graph_tool.dynamics.DiscreteStateBase

Kirman’s “ant colony” model.

Parameters
gGraph

Graph to be used for the dynamics

dfloat (optional, default: .1)

Strategy infection probability.

c1float (optional, default: .001)

Spontaneous transition probability to first strategy.

c2float (optional, default: .001)

Spontaneous transition probability to second strategy.

sVertexPropertyMap (optional, default: None)

Initial global state. If not provided, a random state will be chosen.

Notes

This implements Kirman’s “ant colony” model [kirman_ants_1993] on a network.

If a node \(i\) is updated at time \(t\), the transition to state \(s_i(t+1) \in \{0,1\}\) is done as follows:

  1. If \(s_i(t) = 0\), we have \(s_i(t) = 1\) with probability \(c_1\).

  2. Otherwise if \(s_i(t) = 1\), we have \(s_i(t) = 0\) with probability \(c_2\).

  3. Otherwise we have \(s_i(t+1) = 1 - s_i(t)\) with probability

    \[1 - (1-d)^{\sum_jA_{ij}(1-\delta_{s_i(t), s_j(t)})}\]
  4. Otherwise we have \(s_i(t+1) = s_i(t)\).

References

kirman_ants_1993(1,2)

A. Kirman, “Ants, Rationality, and Recruitment”, The Quarterly Journal of Economics 108, 137 (1993), DOI: 10.2307/2118498 [sci-hub, @tor].

Examples

>>> g = gt.GraphView(gt.collection.data["polblogs"], directed=False)
>>> gt.remove_parallel_edges(g)
>>> g = gt.extract_largest_component(g, prune=True)
>>> state = gt.KirmanState(g)
>>> ret = state.iterate_sync(niter=1000)
>>> gt.graph_draw(g, g.vp.pos, vertex_fill_color=state.s, vcmap=cm.bone,
...               output_size=(700,400), output="kirman.pdf")
<...>
_images/kirman.svg

State of Kirman’s model on a political blog network.

class graph_tool.dynamics.AxelrodState(g, f=10, q=2, r=0, s=None)[source]

Bases: graph_tool.dynamics.DiscreteStateBase

Axelrod’s culture dissemination model.

Parameters
gGraph

Graph to be used for the dynamics

fint (optional, default: 10)

Number of features.

qint (optional, default: 2)

Number of traits for each feature.

rfloat (optional, default: .0)

Spontaneous trait change probability.

sVertexPropertyMap (optional, default: None)

Initial global state. If not provided, a random state will be chosen.

Notes

This implements Axelrod’s model for culture dissemination [axelrod-dissemination-1997].

Each node has a vector state \(\boldsymbol s^{(i)} \in \{0,\dots,q-1\}^f\).

If a node \(i\) is updated at time \(t\), the transition to state \(\boldsymbol s^{(i)}(t+1)\) is done as follows:

  1. With probability \(r\) a feature \(l\) is chosen uniformly at random from the interval \(\{0,\dots,f-1\}\), and a trait \(r\) is chosen uniformly at random from the interval \(\{0,\dots,q-1\}\), and the new state is set as \(s^{(i)}_l(t+1)=r\).

  2. Otherwise, a neighbour \(j\) is chosen uniformly at random, and we let \(d\) be the number of equal traits across features between \(i\) and \(j\),

    \[d = \sum_{l=0}^{f-1} \delta_{s^{(i)}_l(t), s^{(j)}_l(t)}.\]

    Then with probability \(d/f\) a trait \(l\) is chosen uniformly at random from the set of differing features of size \(f-d\), i.e. \(\{l|s^{(i)}_l(t) \ne s^{(j)}_l(t)\}\), and the corresponding trait of \(j\) is copied to \(i\): \(s^{(i)}_l(t+1) = s^{(j)}_l(t)\).

  3. Otherwise we have \(\boldsymbol s_i(t+1) = \boldsymbol s_i(t)\).

References

axelrod-dissemination-1997(1,2)

Axelrod, R., “The Dissemination of Culture: A Model with Local Convergence and Global Polarization”, Journal of Conflict Resolution, 41(2), 203–226 (1997). DOI: 10.1177/0022002797041002001 [sci-hub, @tor]

Examples

>>> g = gt.GraphView(gt.collection.data["polblogs"], directed=False)
>>> gt.remove_parallel_edges(g)
>>> g = gt.extract_largest_component(g, prune=True)
>>> state = gt.AxelrodState(g, f=10, q=30, r=0.005)
>>> ret = state.iterate_async(niter=10000000)
>>> gt.graph_draw(g, g.vp.pos,
...               vertex_fill_color=gt.perfect_prop_hash([state.s])[0],
...               vcmap=cm.magma, output_size=(700,400), output="axelrod.pdf")
<...>
_images/axelrod.svg

State of Axelrod’s model on a political blog network.

class graph_tool.dynamics.BooleanState(g, f=None, p=0.5, r=0, s=None)[source]

Bases: graph_tool.dynamics.DiscreteStateBase

Boolean network dynamics.

Parameters
gGraph

Graph to be used for the dynamics

fVertexPropertyMap (optional, default: None)

Vertex property map of type vector<bool> containing the Boolean functions. If not provided, the functions will be randomly chosen.

pfloat (optional, default: .5)

Output probability of random functions. This only has an effect if f is None.

rfloat (optional, default: 0.)

Input random flip probability.

sVertexPropertyMap (optional, default: None)

Initial global state. If not provided, a random state will be chosen.

Notes

This implements a Boolean network model.

If a node \(i\) is updated at time \(t\), the transition to state \(s_i(t+1)\) is given by

\[s_i(t+1) = f^{(i)}_{\sum_{j\in \partial i}2^{\hat s_j(t)}}\]

where \(\partial i\) are the (in-)neighbors of \(i\), indexed from \(0\) to \(k-1\), and \(\hat s_i(t)\) are the flipped inputs sampled with probability

\[P(\hat s_i(t)|s_i(t)) = r^{1-\delta_{\hat s_i(t),s_i(t)}}(1-r)^{\delta_{\hat s_i(t),s_i(t)}}.\]

Examples

>>> g = gt.random_graph(50, lambda: (2,2))
>>> state = gt.BooleanState(g)
>>> ret = state.iterate_sync(niter=1000)
>>> s0 = state.s.copy()
>>> ret = state.iterate_sync(niter=1)
>>> l = 1
>>> while any(state.s.a != s0.a):
...     ret = state.iterate_sync(niter=1)
...     l += 1
>>> print("Period length:", l)
Period length: 3
class graph_tool.dynamics.ContinuousStateBase(g, make_state, params, t0=0, s=None, stype='double')[source]

Bases: object

Base state for continuous-time dynamics. This class it not meant to be instantiated directly.

copy()[source]

Copy state.

get_state()[source]

Returns the internal VertexPropertyMap with the current state.

get_diff(dt)[source]

Returns the current time derivative for all the nodes. The parameter dt is the time interval in consideration, which is used only if the ODE has a stochastic component.

If enabled during compilation, this algorithm runs in parallel.

solve(t, *args, **kwargs)[source]

Integrate the system up to time t. The remaining parameters are passed to scipy.integrate.solve_ivp(). This solver is not suitable for stochastic ODEs.

solve_euler(t, dt=0.001)[source]

Integrate the system up o time t using a simple Euler’s method with step size dt. This solver is suitable for stochastic ODEs.

class graph_tool.dynamics.KuramotoState(g, omega=1, w=1, sigma=0, t0=0, s=None)[source]

Bases: graph_tool.dynamics.ContinuousStateBase

The Kuramoto model.

Parameters
gGraph

Graph to be used for the dynamics

omegaVertexPropertyMap or float (optional, default: 1)

Intrinsic frequencies for each node. If a scalar is given, it will be used for all nodes.

wEdgePropertyMap or float (optional, default: 1)

Coupling strength of each edge. If a scalar is given, it will be used for all edges.

sigmafloat (optional, default: .0)

Stochastic noise magnitude.

sVertexPropertyMap (optional, default: None)

Initial global state. If not provided, a random state will be chosen.

Notes

This implements Kuramoto’s model for synchronization [kuramoto_self-entrainment_1975] [rodrigues_kuramoto_2016].

Each node has an angle \(\theta_i\), which evolves in time obeying the differential equation:

\[\frac{\mathrm{d}\theta_i}{\mathrm{d}t} = \omega_i + \sum_{j}A_{ij}w_{ij}\sin(\theta_j-\theta_i) + \sigma\xi_i(t),\]

where \(\xi_i(t)\) is a Gaussian noise with zero mean and unit variance.

References

kuramoto_self-entrainment_1975(1,2)

Y. Kuramoto, “Self-entrainment of a population of coupled non-linear oscillators”, International Symposium on Mathematical Problems in Theoretical Physics. Lecture Notes in Physics, vol 39. Springer, Berlin, Heidelberg (1975), DOI: 10.1007/BFb0013365 [sci-hub, @tor]

rodrigues_kuramoto_2016(1,2)

Francisco A. Rodrigues, Thomas K. DM.Peron, Peng Ji, Jürgen Kurth, “The Kuramoto model in complex networks”, Physics Reports 610 1-98 (2016) DOI: 10.1016/j.physrep.2015.10.008 [sci-hub, @tor], arXiv: 1511.07139

Examples

>>> g = gt.collection.data["karate"]
>>> omega = g.new_vp("double", np.random.normal(0, 1, g.num_vertices())) 
>>> state = gt.KuramotoState(g, omega=omega, w=1.5)
>>> thetas = []
>>> ts = linspace(0, 40, 1000)
>>> for t in ts:
...     ret = state.solve(t, first_step=0.0001)
...     thetas.append(state.get_state().fa % (2 * pi))

>>> figure(figsize=(6, 4))
<...>
>>> for v in g.vertices():
...    plot(ts, [t[int(v)] - t.mean() for t in thetas])
[...]
>>> xlabel(r"Time")
Text(...)
>>> ylabel(r"$\theta_i - \left<\theta\right>$")
Text(...)
>>> tight_layout()
>>> savefig("karate-kuramoto.pdf")
_images/karate-kuramoto.svg

Kuramoto oscilator dynamics on the Karate Club network.

get_r_phi()[source]

Returns the phase coherence \(r\) and average phase \(\phi\), defined as

\[re^{i\phi} = \frac{1}{N}\sum_j e^{i\theta_j}.\]