Implement contact process (CP)

[cstarkjp]

The next big step is to implement the contact process, and to do so in such a way that:

1. It can be at least minimally parallelized
2. The pair contact process (PCP) and two-species symbiotic pair process (2SCP) can be readily implemented as well.
3. Subsequently a granular yield version can be implemented.

[cstarkjp]

Asynchronous DP: contact process

• Zero the iteration counter: $k=0$
• Construct a $d$-dimensional lattice $\mathcal{L}$ of size $L$ along each axis
• Compute the total number of sites: $N = L^d$
• Choose an initial site occupancy probability: $p_{0}$
• Randomize the occupancy of all lattice sites ${s_i} \in \mathcal{L}$ such that, for each $i$, $s_i = \mathsf{Bern}(p_0)$
• Make a list of all the lattice sites, not just the occupied ones: $\mathcal{S}_\textsf{lattice}$
• Randomize the order of this list: $\mathcal{S}_\textsf{lattice} \mapsto \mathcal{S}\textsf{randomized}$

• Set a reporting interval: $\Delta{t}_\mathsf{report}$
• Create an empty list $\mathcal{L}_\mathsf{report}$ to record the lattice after a sequence of changes over this interval
• Create an empty list $\mathcal{T}_\mathsf{report}$ to record the report times
• Zero the timer: $t = 0$
• Zero the report interval timer: $t_\mathsf{report} = 0$
• Loop over $K$ iterations $k \in {1, \dots, K}$:
  • Consider the first site in the randomized list $\mathcal{S}_\textsf{randomized}$ and label it $i$ with occupancy $s_i$
  • Remove this site $i$ from the list $\mathcal{S}_\mathsf{randomized}$
  • Find its neighborhood and label the sites $j \in \langle i \rangle$ with occupancy $s_j$
  • Count the number of occupied sites in the neighborhood: $n_i = \sum_{j\in \langle i \rangle} s_j$
  • Compute the propagation rate: $w_i[0 \mapsto 1, n_i] = \lambda n_i/2 d$
  • Note the annihilation rate: $w_i[1 \mapsto 0, n_i] =1$
  • Compute the propagation probability: $p_\mathsf{p} = \dfrac{w_i[0 \mapsto 1, n_i]}{1+\lambda} = \dfrac{\lambda n_i}{2d(1+\lambda)}$
  • Compute the annihilation probability: $p_\mathsf{a} = \dfrac{w_i[1 \mapsto 0, n_i]}{1+\lambda} = \dfrac{1}{(1+\lambda)}$
  • Compute the time interval: $\Delta{t} = \dfrac{1}{(1+\lambda)N}$
  • Update the timer: $t \mapsto t + \Delta{t}$
    • note: the timer is incremented even if nothing happens
  • Check $s_i$
    • for an empty site $s_i==0$:
      • compute a Bernoulli sample $b_\mathsf{p} = \mathsf{Bern}(p_\mathsf{p})$
      • if $b_\mathsf{p}==\mathsf{true}$:
        • designate the site as now occupied $s_i \mapsto 1$
        • modify the lattice accordingly $\mathcal{L} \mapsto \mathcal{L}(t)$
    • for an occupied site $s_i==1$:
      • compute a Bernoulli sample $b_\mathsf{a} = \mathsf{Bern}(p_\mathsf{a})$
      • if $b_\mathsf{a}==\mathsf{true}$:
        • designate the site as now empty $s_i \mapsto 0$
        • modify the lattice accordingly $\mathcal{L} \mapsto \mathcal{L}(t)$
  • Reinsert the site $i$ at a random location in the list $\mathcal{S}_\mathsf{randomized}$
  • Update the report timer: $t_\mathsf{report} \mapsto t_\mathsf{report} + \Delta{t}$
  • Check it's time to report, $t_\mathsf{report} \geq \Delta{}t_\mathsf{report}$; if so:
    • append the updated (or not) lattice to the report list: $\mathcal{L}_\mathsf{report} \mapsto \mathcal{L}\mathsf{report} + \mathcal{L}$
    • append the current time $t$ to the report time list $\mathcal{T}_\mathsf{report}$
    • reset the report interval timer $t_\mathsf{report} \mapsto 0$
• Return the report list of lattices $\mathcal{L}_\mathsf{report}$
• Return the report time list $\mathcal{T}_\mathsf{report}$

[cstarkjp]

Simplified version

• Zero the iteration counter: $\text{ } k=0$
• Construct a $d$-dimensional lattice $\mathcal{L}$ of size $L$ along each axis
• Compute the total number of sites: $\text{ } N = L^d$
• Choose an initial site occupancy probability $p_{0}$
• Randomize the occupancy of all lattice sites ${s_i} \in \mathcal{L}$ such that, for each $i$, $s_i = \mathsf{Bern}(p_0)$
• Set a reporting interval $\Delta{t}_{\mathsf{report}}$
• Create an empty list $\mathcal{L}_{\mathsf{report}}$ to record the lattice after a sequence of changes over this interval
• Create an empty list $\mathcal{T}_{\mathsf{report}}$ to record the report times
• Zero the timer: $\text{ } t = 0$
• Zero the report interval timer: $\text{ } t_\mathsf{report} = 0$
• Loop over $K$ iterations $k \in {1, \dots, K}$:
  • Pick a site at random from the lattice, labeling its position as $i$ and its occupancy as $s_i$
  • Find its neighborhood and label the sites $j \in \langle i \rangle$ with occupancy $s_j$
  • Count the number of occupied sites in the neighborhood: $\text{ } n_i = \sum_{j\in \langle i \rangle} s_j$
  • Note the propagation rate: $\text{ } w_i[0 \mapsto 1, n_i] = \lambda n_i/2 d$
  • Note the annihilation rate: $\text{ } w_i[1 \mapsto 0, n_i] =1$
  • Check $s_i$
    • for an occupied site $s_i\text{==}1$:
      • compute the annihilation probability: $\text{ } p_\mathsf{a} = \dfrac{w_i[1 \mapsto 0, n_i]}{1+\lambda} = \dfrac{1}{(1+\lambda)}$
      • compute a Bernoulli sample $b_\mathsf{a} = \mathsf{Bern}(p_\mathsf{a})$
      • if $b_\mathsf{a}\text{==}\mathsf{true}$:
        • designate the site as now empty $s_i \mapsto 0$
        • modify (in-place) the current lattice accordingly $\mathcal{L} \rightarrow \mathcal{L}(t)$
    • for an empty site $s_i\text{==}0$ with some occupied neighbors $n_i \neq 0$:
      • compute the propagation probability: $\text{ } p_\mathsf{p} = \dfrac{w_i[0 \mapsto 1, n_i]}{1+\lambda} = \dfrac{\lambda n_i}{2d(1+\lambda)}$
      • compute a Bernoulli sample $b_\mathsf{p} = \mathsf{Bern}(p_\mathsf{p})$
      • if $b_\mathsf{p}\text{==}\mathsf{true}$:
        • designate the site as now occupied $s_i \mapsto 1$
        • modify (in-place) the current lattice accordingly $\mathcal{L} \rightarrow \mathcal{L}(t)$
    • for an empty site $s_i\text{==}0$ with no occupied neighbors $n_i = 0$:
      • do nothing
  • Compute the time interval: $\text{ } \Delta{t} = \dfrac{1}{(1+\lambda)N}$
  • Update the timer: $\text{ } t \mapsto t + \Delta{t}$
    • note: the timer is incremented even if nothing happens
  • Update the report timer: $\text{ } t_\mathsf{report} \mapsto t_\mathsf{report} + \Delta{t}$
  • Check it's time to report, $t_\mathsf{report} \geq \Delta{}t_\mathsf{report}$; if so:
    • append the updated (or not) lattice $\mathcal{L}$ to the report list: $\rightarrow \mathcal{L}_\mathsf{report} + \mathcal{L}$
    • append the current time $t$ to the report time list: $\rightarrow \mathcal{T}_{\mathsf{report}} + t$
    • reset the report interval timer $t_\mathsf{report} \mapsto 0$
• Return the list of lattices $\mathcal{L}_{\mathsf{report}}$
• Return the report-times list $\mathcal{T}_{\mathsf{report}}$