Two-Patch SIR: Spatial Correlation (Keeling & Rohani 2002)¶
Theory¶
Keeling & Rohani (2002) studied two identical SIR patches connected by symmetric migration at rate σ. The temporal correlation between the infectious time series of the two patches scales approximately as:
$$C(\sigma) \approx \frac{\sigma}{\xi + \sigma}$$
where ξ is an intrinsic desynchronisation parameter set by the patch dynamics (extinction probability, epidemic period, etc.).
Key predictions:
- At σ → 0 (no coupling), patches are independent → C ≈ 0.
- At σ → ∞ (tight coupling), patches synchronise → C → 1.
- The midpoint C = 0.5 occurs at σ = ξ.
Goal: Sweep σ over 100 log-spaced values from 10⁻⁵ to 10⁻¹, run 100 years of two-patch SIR with vital dynamics for each σ, compute the Pearson correlation between the I time series of the two patches (after a 20-year burn-in), then fit the sigmoidal C(σ) curve and recover ξ.
In [1]:
Copied!
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
from laser.core import PropertySet
from laser.core.utils import grid
from laser.generic.utils import ValuesMap
from laser.cohorts import Model
from laser.cohorts.vitaldynamics import NonDiseaseMortality, ConstantPopBirths
from laser.cohorts.migration import Migration
from laser.cohorts.utils import static_routing
import laser.cohorts.SIR as SIR
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
from laser.core import PropertySet
from laser.core.utils import grid
from laser.generic.utils import ValuesMap
from laser.cohorts import Model
from laser.cohorts.vitaldynamics import NonDiseaseMortality, ConstantPopBirths
from laser.cohorts.migration import Migration
from laser.cohorts.utils import static_routing
import laser.cohorts.SIR as SIR
In [2]:
Copied!
# Patch parameters
N_per_node = 500_000
I0 = 10 # infected individuals seeded in each patch
CDR = 20 # crude death rate per 1000 per year
beta = 0.5 # transmission rate (per day)
gamma = 0.1 # recovery rate (per day); R0 = beta/gamma = 5
nticks = 100 * 365 # 100 years of daily steps
n_patches = 2
r_mortality_val = CDR / 1000 / 365 # daily mortality rate
print(f"R0 = {beta/gamma:.1f}")
print(f"Daily mortality rate: {r_mortality_val:.6f}")
print(f"Simulation length: {nticks:,} ticks ({nticks//365} years)")
# Patch parameters
N_per_node = 500_000
I0 = 10 # infected individuals seeded in each patch
CDR = 20 # crude death rate per 1000 per year
beta = 0.5 # transmission rate (per day)
gamma = 0.1 # recovery rate (per day); R0 = beta/gamma = 5
nticks = 100 * 365 # 100 years of daily steps
n_patches = 2
r_mortality_val = CDR / 1000 / 365 # daily mortality rate
print(f"R0 = {beta/gamma:.1f}")
print(f"Daily mortality rate: {r_mortality_val:.6f}")
print(f"Simulation length: {nticks:,} ticks ({nticks//365} years)")
R0 = 5.0 Daily mortality rate: 0.000055 Simulation length: 36,500 ticks (100 years)
In [3]:
Copied!
def run_patches(sigma, nticks, N, I0, beta, gamma, r_mortality_val):
"""Run a two-patch SIR model with symmetric migration rate sigma.
Returns the Pearson correlation between the I time series of the two
patches, computed after a 20-year burn-in period.
Args:
sigma: Symmetric migration rate between the two patches.
nticks: Total simulation length in ticks (days).
N: Initial population size per patch.
I0: Initial infectious count per patch.
beta: Transmission rate (per day).
gamma: Recovery rate (per day).
r_mortality_val: Daily non-disease mortality rate.
Returns:
Pearson correlation coefficient in [-1, 1].
"""
scenario = grid(M=2, N=1)
scenario["S"] = [N - I0, N - I0]
scenario["I"] = [I0, I0]
scenario["R"] = [0, 0]
params = PropertySet({"nticks": nticks})
model = Model(scenario, params)
routing_2d = np.array([[0, sigma], [sigma, 0]], dtype=np.float64)
routing_3d = static_routing(routing_2d, nticks)
betas = ValuesMap.from_scalar(beta, nticks, 2)
r_recoveries = ValuesMap.from_scalar(gamma, nticks, 2)
model.components = [
SIR.Susceptible(model),
SIR.Infectious(model, r_recovery=r_recoveries),
SIR.Recovered(model),
NonDiseaseMortality(model, r_mortality=r_mortality_val),
ConstantPopBirths(model),
Migration(model, r_migration=sigma, routing=routing_3d),
SIR.Transmission(model, beta=betas),
]
model.run()
I0_ts = model.states.I[:, 0]
I1_ts = model.states.I[:, 1]
# Discard first 20 years as burn-in to reach quasi-stationary behaviour
burn = 20 * 365
corr = np.corrcoef(I0_ts[burn:], I1_ts[burn:])[0, 1]
return corr
def run_patches(sigma, nticks, N, I0, beta, gamma, r_mortality_val):
"""Run a two-patch SIR model with symmetric migration rate sigma.
Returns the Pearson correlation between the I time series of the two
patches, computed after a 20-year burn-in period.
Args:
sigma: Symmetric migration rate between the two patches.
nticks: Total simulation length in ticks (days).
N: Initial population size per patch.
I0: Initial infectious count per patch.
beta: Transmission rate (per day).
gamma: Recovery rate (per day).
r_mortality_val: Daily non-disease mortality rate.
Returns:
Pearson correlation coefficient in [-1, 1].
"""
scenario = grid(M=2, N=1)
scenario["S"] = [N - I0, N - I0]
scenario["I"] = [I0, I0]
scenario["R"] = [0, 0]
params = PropertySet({"nticks": nticks})
model = Model(scenario, params)
routing_2d = np.array([[0, sigma], [sigma, 0]], dtype=np.float64)
routing_3d = static_routing(routing_2d, nticks)
betas = ValuesMap.from_scalar(beta, nticks, 2)
r_recoveries = ValuesMap.from_scalar(gamma, nticks, 2)
model.components = [
SIR.Susceptible(model),
SIR.Infectious(model, r_recovery=r_recoveries),
SIR.Recovered(model),
NonDiseaseMortality(model, r_mortality=r_mortality_val),
ConstantPopBirths(model),
Migration(model, r_migration=sigma, routing=routing_3d),
SIR.Transmission(model, beta=betas),
]
model.run()
I0_ts = model.states.I[:, 0]
I1_ts = model.states.I[:, 1]
# Discard first 20 years as burn-in to reach quasi-stationary behaviour
burn = 20 * 365
corr = np.corrcoef(I0_ts[burn:], I1_ts[burn:])[0, 1]
return corr
In [4]:
Copied!
# Sweep 100 log-spaced sigma values — this may take several minutes
sigma_values = np.logspace(-5, -1, 100)
correlations = []
for i, s in enumerate(sigma_values):
c = run_patches(s, nticks, N_per_node, I0, beta, gamma, r_mortality_val)
correlations.append(c)
if (i + 1) % 10 == 0:
print(f" {i+1:3d}/100 sigma={s:.2e} corr={c:.3f}")
correlations = np.array(correlations)
print("Sweep complete.")
# Sweep 100 log-spaced sigma values — this may take several minutes
sigma_values = np.logspace(-5, -1, 100)
correlations = []
for i, s in enumerate(sigma_values):
c = run_patches(s, nticks, N_per_node, I0, beta, gamma, r_mortality_val)
correlations.append(c)
if (i + 1) % 10 == 0:
print(f" {i+1:3d}/100 sigma={s:.2e} corr={c:.3f}")
correlations = np.array(correlations)
print("Sweep complete.")
/Users/christopherlorton/projects/laser-fresh/laser.cohorts/.docs/lib/python3.12/site-packages/numpy/lib/_function_base_impl.py:3023: RuntimeWarning: invalid value encountered in divide c /= stddev[:, None] /Users/christopherlorton/projects/laser-fresh/laser.cohorts/.docs/lib/python3.12/site-packages/numpy/lib/_function_base_impl.py:3023: RuntimeWarning: invalid value encountered in divide c /= stddev[:, None] /Users/christopherlorton/projects/laser-fresh/laser.cohorts/.docs/lib/python3.12/site-packages/numpy/lib/_function_base_impl.py:3023: RuntimeWarning: invalid value encountered in divide c /= stddev[:, None] /Users/christopherlorton/projects/laser-fresh/laser.cohorts/.docs/lib/python3.12/site-packages/numpy/lib/_function_base_impl.py:3023: RuntimeWarning: invalid value encountered in divide c /= stddev[:, None] /Users/christopherlorton/projects/laser-fresh/laser.cohorts/.docs/lib/python3.12/site-packages/numpy/lib/_function_base_impl.py:3023: RuntimeWarning: invalid value encountered in divide c /= stddev[:, None] /Users/christopherlorton/projects/laser-fresh/laser.cohorts/.docs/lib/python3.12/site-packages/numpy/lib/_function_base_impl.py:3023: RuntimeWarning: invalid value encountered in divide c /= stddev[:, None] /Users/christopherlorton/projects/laser-fresh/laser.cohorts/.docs/lib/python3.12/site-packages/numpy/lib/_function_base_impl.py:3023: RuntimeWarning: invalid value encountered in divide c /= stddev[:, None] /Users/christopherlorton/projects/laser-fresh/laser.cohorts/.docs/lib/python3.12/site-packages/numpy/lib/_function_base_impl.py:3023: RuntimeWarning: invalid value encountered in divide c /= stddev[:, None] /Users/christopherlorton/projects/laser-fresh/laser.cohorts/.docs/lib/python3.12/site-packages/numpy/lib/_function_base_impl.py:3023: RuntimeWarning: invalid value encountered in divide c /= stddev[:, None] /Users/christopherlorton/projects/laser-fresh/laser.cohorts/.docs/lib/python3.12/site-packages/numpy/lib/_function_base_impl.py:3023: RuntimeWarning: invalid value encountered in divide c /= stddev[:, None]
10/100 sigma=2.31e-05 corr=nan
/Users/christopherlorton/projects/laser-fresh/laser.cohorts/.docs/lib/python3.12/site-packages/numpy/lib/_function_base_impl.py:3023: RuntimeWarning: invalid value encountered in divide c /= stddev[:, None] /Users/christopherlorton/projects/laser-fresh/laser.cohorts/.docs/lib/python3.12/site-packages/numpy/lib/_function_base_impl.py:3023: RuntimeWarning: invalid value encountered in divide c /= stddev[:, None] /Users/christopherlorton/projects/laser-fresh/laser.cohorts/.docs/lib/python3.12/site-packages/numpy/lib/_function_base_impl.py:3023: RuntimeWarning: invalid value encountered in divide c /= stddev[:, None] /Users/christopherlorton/projects/laser-fresh/laser.cohorts/.docs/lib/python3.12/site-packages/numpy/lib/_function_base_impl.py:3023: RuntimeWarning: invalid value encountered in divide c /= stddev[:, None] /Users/christopherlorton/projects/laser-fresh/laser.cohorts/.docs/lib/python3.12/site-packages/numpy/lib/_function_base_impl.py:3023: RuntimeWarning: invalid value encountered in divide c /= stddev[:, None] /Users/christopherlorton/projects/laser-fresh/laser.cohorts/.docs/lib/python3.12/site-packages/numpy/lib/_function_base_impl.py:3023: RuntimeWarning: invalid value encountered in divide c /= stddev[:, None] /Users/christopherlorton/projects/laser-fresh/laser.cohorts/.docs/lib/python3.12/site-packages/numpy/lib/_function_base_impl.py:3023: RuntimeWarning: invalid value encountered in divide c /= stddev[:, None] /Users/christopherlorton/projects/laser-fresh/laser.cohorts/.docs/lib/python3.12/site-packages/numpy/lib/_function_base_impl.py:3023: RuntimeWarning: invalid value encountered in divide c /= stddev[:, None] /Users/christopherlorton/projects/laser-fresh/laser.cohorts/.docs/lib/python3.12/site-packages/numpy/lib/_function_base_impl.py:3023: RuntimeWarning: invalid value encountered in divide c /= stddev[:, None] /Users/christopherlorton/projects/laser-fresh/laser.cohorts/.docs/lib/python3.12/site-packages/numpy/lib/_function_base_impl.py:3023: RuntimeWarning: invalid value encountered in divide c /= stddev[:, None]
20/100 sigma=5.86e-05 corr=nan
/Users/christopherlorton/projects/laser-fresh/laser.cohorts/.docs/lib/python3.12/site-packages/numpy/lib/_function_base_impl.py:3023: RuntimeWarning: invalid value encountered in divide c /= stddev[:, None] /Users/christopherlorton/projects/laser-fresh/laser.cohorts/.docs/lib/python3.12/site-packages/numpy/lib/_function_base_impl.py:3023: RuntimeWarning: invalid value encountered in divide c /= stddev[:, None] /Users/christopherlorton/projects/laser-fresh/laser.cohorts/.docs/lib/python3.12/site-packages/numpy/lib/_function_base_impl.py:3023: RuntimeWarning: invalid value encountered in divide c /= stddev[:, None] /Users/christopherlorton/projects/laser-fresh/laser.cohorts/.docs/lib/python3.12/site-packages/numpy/lib/_function_base_impl.py:3023: RuntimeWarning: invalid value encountered in divide c /= stddev[:, None] /Users/christopherlorton/projects/laser-fresh/laser.cohorts/.docs/lib/python3.12/site-packages/numpy/lib/_function_base_impl.py:3023: RuntimeWarning: invalid value encountered in divide c /= stddev[:, None] /Users/christopherlorton/projects/laser-fresh/laser.cohorts/.docs/lib/python3.12/site-packages/numpy/lib/_function_base_impl.py:3023: RuntimeWarning: invalid value encountered in divide c /= stddev[:, None] /Users/christopherlorton/projects/laser-fresh/laser.cohorts/.docs/lib/python3.12/site-packages/numpy/lib/_function_base_impl.py:3023: RuntimeWarning: invalid value encountered in divide c /= stddev[:, None] /Users/christopherlorton/projects/laser-fresh/laser.cohorts/.docs/lib/python3.12/site-packages/numpy/lib/_function_base_impl.py:3023: RuntimeWarning: invalid value encountered in divide c /= stddev[:, None] /Users/christopherlorton/projects/laser-fresh/laser.cohorts/.docs/lib/python3.12/site-packages/numpy/lib/_function_base_impl.py:3023: RuntimeWarning: invalid value encountered in divide c /= stddev[:, None] /Users/christopherlorton/projects/laser-fresh/laser.cohorts/.docs/lib/python3.12/site-packages/numpy/lib/_function_base_impl.py:3023: RuntimeWarning: invalid value encountered in divide c /= stddev[:, None]
30/100 sigma=1.48e-04 corr=nan
/Users/christopherlorton/projects/laser-fresh/laser.cohorts/.docs/lib/python3.12/site-packages/numpy/lib/_function_base_impl.py:3023: RuntimeWarning: invalid value encountered in divide c /= stddev[:, None] /Users/christopherlorton/projects/laser-fresh/laser.cohorts/.docs/lib/python3.12/site-packages/numpy/lib/_function_base_impl.py:3023: RuntimeWarning: invalid value encountered in divide c /= stddev[:, None] /Users/christopherlorton/projects/laser-fresh/laser.cohorts/.docs/lib/python3.12/site-packages/numpy/lib/_function_base_impl.py:3023: RuntimeWarning: invalid value encountered in divide c /= stddev[:, None] /Users/christopherlorton/projects/laser-fresh/laser.cohorts/.docs/lib/python3.12/site-packages/numpy/lib/_function_base_impl.py:3023: RuntimeWarning: invalid value encountered in divide c /= stddev[:, None] /Users/christopherlorton/projects/laser-fresh/laser.cohorts/.docs/lib/python3.12/site-packages/numpy/lib/_function_base_impl.py:3023: RuntimeWarning: invalid value encountered in divide c /= stddev[:, None] /Users/christopherlorton/projects/laser-fresh/laser.cohorts/.docs/lib/python3.12/site-packages/numpy/lib/_function_base_impl.py:3023: RuntimeWarning: invalid value encountered in divide c /= stddev[:, None] /Users/christopherlorton/projects/laser-fresh/laser.cohorts/.docs/lib/python3.12/site-packages/numpy/lib/_function_base_impl.py:3023: RuntimeWarning: invalid value encountered in divide c /= stddev[:, None] /Users/christopherlorton/projects/laser-fresh/laser.cohorts/.docs/lib/python3.12/site-packages/numpy/lib/_function_base_impl.py:3023: RuntimeWarning: invalid value encountered in divide c /= stddev[:, None] /Users/christopherlorton/projects/laser-fresh/laser.cohorts/.docs/lib/python3.12/site-packages/numpy/lib/_function_base_impl.py:3023: RuntimeWarning: invalid value encountered in divide c /= stddev[:, None] /Users/christopherlorton/projects/laser-fresh/laser.cohorts/.docs/lib/python3.12/site-packages/numpy/lib/_function_base_impl.py:3023: RuntimeWarning: invalid value encountered in divide c /= stddev[:, None]
40/100 sigma=3.76e-04 corr=nan
/Users/christopherlorton/projects/laser-fresh/laser.cohorts/.docs/lib/python3.12/site-packages/numpy/lib/_function_base_impl.py:3023: RuntimeWarning: invalid value encountered in divide c /= stddev[:, None] /Users/christopherlorton/projects/laser-fresh/laser.cohorts/.docs/lib/python3.12/site-packages/numpy/lib/_function_base_impl.py:3023: RuntimeWarning: invalid value encountered in divide c /= stddev[:, None] /Users/christopherlorton/projects/laser-fresh/laser.cohorts/.docs/lib/python3.12/site-packages/numpy/lib/_function_base_impl.py:3023: RuntimeWarning: invalid value encountered in divide c /= stddev[:, None] /Users/christopherlorton/projects/laser-fresh/laser.cohorts/.docs/lib/python3.12/site-packages/numpy/lib/_function_base_impl.py:3023: RuntimeWarning: invalid value encountered in divide c /= stddev[:, None] /Users/christopherlorton/projects/laser-fresh/laser.cohorts/.docs/lib/python3.12/site-packages/numpy/lib/_function_base_impl.py:3023: RuntimeWarning: invalid value encountered in divide c /= stddev[:, None] /Users/christopherlorton/projects/laser-fresh/laser.cohorts/.docs/lib/python3.12/site-packages/numpy/lib/_function_base_impl.py:3023: RuntimeWarning: invalid value encountered in divide c /= stddev[:, None] /Users/christopherlorton/projects/laser-fresh/laser.cohorts/.docs/lib/python3.12/site-packages/numpy/lib/_function_base_impl.py:3023: RuntimeWarning: invalid value encountered in divide c /= stddev[:, None] /Users/christopherlorton/projects/laser-fresh/laser.cohorts/.docs/lib/python3.12/site-packages/numpy/lib/_function_base_impl.py:3023: RuntimeWarning: invalid value encountered in divide c /= stddev[:, None] /Users/christopherlorton/projects/laser-fresh/laser.cohorts/.docs/lib/python3.12/site-packages/numpy/lib/_function_base_impl.py:3023: RuntimeWarning: invalid value encountered in divide c /= stddev[:, None] /Users/christopherlorton/projects/laser-fresh/laser.cohorts/.docs/lib/python3.12/site-packages/numpy/lib/_function_base_impl.py:3023: RuntimeWarning: invalid value encountered in divide c /= stddev[:, None]
50/100 sigma=9.55e-04 corr=nan
/Users/christopherlorton/projects/laser-fresh/laser.cohorts/.docs/lib/python3.12/site-packages/numpy/lib/_function_base_impl.py:3023: RuntimeWarning: invalid value encountered in divide c /= stddev[:, None] /Users/christopherlorton/projects/laser-fresh/laser.cohorts/.docs/lib/python3.12/site-packages/numpy/lib/_function_base_impl.py:3023: RuntimeWarning: invalid value encountered in divide c /= stddev[:, None] /Users/christopherlorton/projects/laser-fresh/laser.cohorts/.docs/lib/python3.12/site-packages/numpy/lib/_function_base_impl.py:3023: RuntimeWarning: invalid value encountered in divide c /= stddev[:, None] /Users/christopherlorton/projects/laser-fresh/laser.cohorts/.docs/lib/python3.12/site-packages/numpy/lib/_function_base_impl.py:3023: RuntimeWarning: invalid value encountered in divide c /= stddev[:, None] /Users/christopherlorton/projects/laser-fresh/laser.cohorts/.docs/lib/python3.12/site-packages/numpy/lib/_function_base_impl.py:3023: RuntimeWarning: invalid value encountered in divide c /= stddev[:, None] /Users/christopherlorton/projects/laser-fresh/laser.cohorts/.docs/lib/python3.12/site-packages/numpy/lib/_function_base_impl.py:3023: RuntimeWarning: invalid value encountered in divide c /= stddev[:, None] /Users/christopherlorton/projects/laser-fresh/laser.cohorts/.docs/lib/python3.12/site-packages/numpy/lib/_function_base_impl.py:3023: RuntimeWarning: invalid value encountered in divide c /= stddev[:, None] /Users/christopherlorton/projects/laser-fresh/laser.cohorts/.docs/lib/python3.12/site-packages/numpy/lib/_function_base_impl.py:3023: RuntimeWarning: invalid value encountered in divide c /= stddev[:, None] /Users/christopherlorton/projects/laser-fresh/laser.cohorts/.docs/lib/python3.12/site-packages/numpy/lib/_function_base_impl.py:3023: RuntimeWarning: invalid value encountered in divide c /= stddev[:, None] /Users/christopherlorton/projects/laser-fresh/laser.cohorts/.docs/lib/python3.12/site-packages/numpy/lib/_function_base_impl.py:3023: RuntimeWarning: invalid value encountered in divide c /= stddev[:, None]
60/100 sigma=2.42e-03 corr=nan
/Users/christopherlorton/projects/laser-fresh/laser.cohorts/.docs/lib/python3.12/site-packages/numpy/lib/_function_base_impl.py:3023: RuntimeWarning: invalid value encountered in divide c /= stddev[:, None] /Users/christopherlorton/projects/laser-fresh/laser.cohorts/.docs/lib/python3.12/site-packages/numpy/lib/_function_base_impl.py:3023: RuntimeWarning: invalid value encountered in divide c /= stddev[:, None] /Users/christopherlorton/projects/laser-fresh/laser.cohorts/.docs/lib/python3.12/site-packages/numpy/lib/_function_base_impl.py:3023: RuntimeWarning: invalid value encountered in divide c /= stddev[:, None] /Users/christopherlorton/projects/laser-fresh/laser.cohorts/.docs/lib/python3.12/site-packages/numpy/lib/_function_base_impl.py:3023: RuntimeWarning: invalid value encountered in divide c /= stddev[:, None] /Users/christopherlorton/projects/laser-fresh/laser.cohorts/.docs/lib/python3.12/site-packages/numpy/lib/_function_base_impl.py:3023: RuntimeWarning: invalid value encountered in divide c /= stddev[:, None] /Users/christopherlorton/projects/laser-fresh/laser.cohorts/.docs/lib/python3.12/site-packages/numpy/lib/_function_base_impl.py:3023: RuntimeWarning: invalid value encountered in divide c /= stddev[:, None] /Users/christopherlorton/projects/laser-fresh/laser.cohorts/.docs/lib/python3.12/site-packages/numpy/lib/_function_base_impl.py:3023: RuntimeWarning: invalid value encountered in divide c /= stddev[:, None] /Users/christopherlorton/projects/laser-fresh/laser.cohorts/.docs/lib/python3.12/site-packages/numpy/lib/_function_base_impl.py:3023: RuntimeWarning: invalid value encountered in divide c /= stddev[:, None] /Users/christopherlorton/projects/laser-fresh/laser.cohorts/.docs/lib/python3.12/site-packages/numpy/lib/_function_base_impl.py:3023: RuntimeWarning: invalid value encountered in divide c /= stddev[:, None] /Users/christopherlorton/projects/laser-fresh/laser.cohorts/.docs/lib/python3.12/site-packages/numpy/lib/_function_base_impl.py:3023: RuntimeWarning: invalid value encountered in divide c /= stddev[:, None]
70/100 sigma=6.14e-03 corr=nan
/Users/christopherlorton/projects/laser-fresh/laser.cohorts/.docs/lib/python3.12/site-packages/numpy/lib/_function_base_impl.py:3023: RuntimeWarning: invalid value encountered in divide c /= stddev[:, None] /Users/christopherlorton/projects/laser-fresh/laser.cohorts/.docs/lib/python3.12/site-packages/numpy/lib/_function_base_impl.py:3023: RuntimeWarning: invalid value encountered in divide c /= stddev[:, None] /Users/christopherlorton/projects/laser-fresh/laser.cohorts/.docs/lib/python3.12/site-packages/numpy/lib/_function_base_impl.py:3023: RuntimeWarning: invalid value encountered in divide c /= stddev[:, None] /Users/christopherlorton/projects/laser-fresh/laser.cohorts/.docs/lib/python3.12/site-packages/numpy/lib/_function_base_impl.py:3023: RuntimeWarning: invalid value encountered in divide c /= stddev[:, None] /Users/christopherlorton/projects/laser-fresh/laser.cohorts/.docs/lib/python3.12/site-packages/numpy/lib/_function_base_impl.py:3023: RuntimeWarning: invalid value encountered in divide c /= stddev[:, None] /Users/christopherlorton/projects/laser-fresh/laser.cohorts/.docs/lib/python3.12/site-packages/numpy/lib/_function_base_impl.py:3023: RuntimeWarning: invalid value encountered in divide c /= stddev[:, None] /Users/christopherlorton/projects/laser-fresh/laser.cohorts/.docs/lib/python3.12/site-packages/numpy/lib/_function_base_impl.py:3023: RuntimeWarning: invalid value encountered in divide c /= stddev[:, None] /Users/christopherlorton/projects/laser-fresh/laser.cohorts/.docs/lib/python3.12/site-packages/numpy/lib/_function_base_impl.py:3023: RuntimeWarning: invalid value encountered in divide c /= stddev[:, None] /Users/christopherlorton/projects/laser-fresh/laser.cohorts/.docs/lib/python3.12/site-packages/numpy/lib/_function_base_impl.py:3023: RuntimeWarning: invalid value encountered in divide c /= stddev[:, None] /Users/christopherlorton/projects/laser-fresh/laser.cohorts/.docs/lib/python3.12/site-packages/numpy/lib/_function_base_impl.py:3023: RuntimeWarning: invalid value encountered in divide c /= stddev[:, None]
80/100 sigma=1.56e-02 corr=nan
/Users/christopherlorton/projects/laser-fresh/laser.cohorts/.docs/lib/python3.12/site-packages/numpy/lib/_function_base_impl.py:3023: RuntimeWarning: invalid value encountered in divide c /= stddev[:, None] /Users/christopherlorton/projects/laser-fresh/laser.cohorts/.docs/lib/python3.12/site-packages/numpy/lib/_function_base_impl.py:3023: RuntimeWarning: invalid value encountered in divide c /= stddev[:, None] /Users/christopherlorton/projects/laser-fresh/laser.cohorts/.docs/lib/python3.12/site-packages/numpy/lib/_function_base_impl.py:3023: RuntimeWarning: invalid value encountered in divide c /= stddev[:, None] /Users/christopherlorton/projects/laser-fresh/laser.cohorts/.docs/lib/python3.12/site-packages/numpy/lib/_function_base_impl.py:3023: RuntimeWarning: invalid value encountered in divide c /= stddev[:, None] /Users/christopherlorton/projects/laser-fresh/laser.cohorts/.docs/lib/python3.12/site-packages/numpy/lib/_function_base_impl.py:3023: RuntimeWarning: invalid value encountered in divide c /= stddev[:, None] /Users/christopherlorton/projects/laser-fresh/laser.cohorts/.docs/lib/python3.12/site-packages/numpy/lib/_function_base_impl.py:3023: RuntimeWarning: invalid value encountered in divide c /= stddev[:, None] /Users/christopherlorton/projects/laser-fresh/laser.cohorts/.docs/lib/python3.12/site-packages/numpy/lib/_function_base_impl.py:3023: RuntimeWarning: invalid value encountered in divide c /= stddev[:, None] /Users/christopherlorton/projects/laser-fresh/laser.cohorts/.docs/lib/python3.12/site-packages/numpy/lib/_function_base_impl.py:3023: RuntimeWarning: invalid value encountered in divide c /= stddev[:, None] /Users/christopherlorton/projects/laser-fresh/laser.cohorts/.docs/lib/python3.12/site-packages/numpy/lib/_function_base_impl.py:3023: RuntimeWarning: invalid value encountered in divide c /= stddev[:, None] /Users/christopherlorton/projects/laser-fresh/laser.cohorts/.docs/lib/python3.12/site-packages/numpy/lib/_function_base_impl.py:3023: RuntimeWarning: invalid value encountered in divide c /= stddev[:, None]
90/100 sigma=3.94e-02 corr=nan
/Users/christopherlorton/projects/laser-fresh/laser.cohorts/.docs/lib/python3.12/site-packages/numpy/lib/_function_base_impl.py:3023: RuntimeWarning: invalid value encountered in divide c /= stddev[:, None] /Users/christopherlorton/projects/laser-fresh/laser.cohorts/.docs/lib/python3.12/site-packages/numpy/lib/_function_base_impl.py:3023: RuntimeWarning: invalid value encountered in divide c /= stddev[:, None] /Users/christopherlorton/projects/laser-fresh/laser.cohorts/.docs/lib/python3.12/site-packages/numpy/lib/_function_base_impl.py:3023: RuntimeWarning: invalid value encountered in divide c /= stddev[:, None] /Users/christopherlorton/projects/laser-fresh/laser.cohorts/.docs/lib/python3.12/site-packages/numpy/lib/_function_base_impl.py:3023: RuntimeWarning: invalid value encountered in divide c /= stddev[:, None] /Users/christopherlorton/projects/laser-fresh/laser.cohorts/.docs/lib/python3.12/site-packages/numpy/lib/_function_base_impl.py:3023: RuntimeWarning: invalid value encountered in divide c /= stddev[:, None] /Users/christopherlorton/projects/laser-fresh/laser.cohorts/.docs/lib/python3.12/site-packages/numpy/lib/_function_base_impl.py:3023: RuntimeWarning: invalid value encountered in divide c /= stddev[:, None] /Users/christopherlorton/projects/laser-fresh/laser.cohorts/.docs/lib/python3.12/site-packages/numpy/lib/_function_base_impl.py:3023: RuntimeWarning: invalid value encountered in divide c /= stddev[:, None] /Users/christopherlorton/projects/laser-fresh/laser.cohorts/.docs/lib/python3.12/site-packages/numpy/lib/_function_base_impl.py:3023: RuntimeWarning: invalid value encountered in divide c /= stddev[:, None] /Users/christopherlorton/projects/laser-fresh/laser.cohorts/.docs/lib/python3.12/site-packages/numpy/lib/_function_base_impl.py:3023: RuntimeWarning: invalid value encountered in divide c /= stddev[:, None]
100/100 sigma=1.00e-01 corr=nan Sweep complete.
/Users/christopherlorton/projects/laser-fresh/laser.cohorts/.docs/lib/python3.12/site-packages/numpy/lib/_function_base_impl.py:3023: RuntimeWarning: invalid value encountered in divide c /= stddev[:, None]
In [5]:
Copied!
# Fit C(sigma) = sigma / (xi + sigma) using scipy curve_fit
def sync_curve(sigma, xi):
"""Keeling-Rohani synchronisation curve."""
return sigma / (xi + sigma)
# Only fit over points where correlation is non-negative and finite
valid = np.isfinite(correlations) & (correlations >= 0)
popt, _ = curve_fit(sync_curve, sigma_values[valid], correlations[valid], p0=[1e-3])
xi_fit = popt[0]
print(f"Fitted desynchronisation parameter ξ = {xi_fit:.4e}")
sigma_fine = np.logspace(-5, -1, 500)
corr_fit_line = sync_curve(sigma_fine, xi_fit)
fig, ax = plt.subplots(figsize=(8, 5))
ax.scatter(np.log10(sigma_values), correlations,
s=20, color="steelblue", alpha=0.7, label="Simulated correlation")
ax.plot(np.log10(sigma_fine), corr_fit_line,
color="firebrick", linewidth=2,
label=f"Fit: C(σ) = σ/(ξ+σ), ξ = {xi_fit:.2e}")
ax.axhline(0.5, color="gray", linestyle="--", linewidth=0.8, label="C = 0.5")
ax.axvline(np.log10(xi_fit), color="orange", linestyle=":",
linewidth=1.2, label=f"σ = ξ = {xi_fit:.2e}")
ax.set_xlabel("log₁₀(σ) (migration rate)")
ax.set_ylabel("Temporal correlation C(σ)")
ax.set_title("Two-patch SIR: synchronisation vs. migration rate\n"
"(Keeling & Rohani 2002)")
ax.legend(fontsize=9)
ax.set_ylim(-0.1, 1.1)
plt.tight_layout()
plt.show()
print(f"\nMidpoint synchronisation (C=0.5) at σ ≈ ξ = {xi_fit:.4e}")
# Fit C(sigma) = sigma / (xi + sigma) using scipy curve_fit
def sync_curve(sigma, xi):
"""Keeling-Rohani synchronisation curve."""
return sigma / (xi + sigma)
# Only fit over points where correlation is non-negative and finite
valid = np.isfinite(correlations) & (correlations >= 0)
popt, _ = curve_fit(sync_curve, sigma_values[valid], correlations[valid], p0=[1e-3])
xi_fit = popt[0]
print(f"Fitted desynchronisation parameter ξ = {xi_fit:.4e}")
sigma_fine = np.logspace(-5, -1, 500)
corr_fit_line = sync_curve(sigma_fine, xi_fit)
fig, ax = plt.subplots(figsize=(8, 5))
ax.scatter(np.log10(sigma_values), correlations,
s=20, color="steelblue", alpha=0.7, label="Simulated correlation")
ax.plot(np.log10(sigma_fine), corr_fit_line,
color="firebrick", linewidth=2,
label=f"Fit: C(σ) = σ/(ξ+σ), ξ = {xi_fit:.2e}")
ax.axhline(0.5, color="gray", linestyle="--", linewidth=0.8, label="C = 0.5")
ax.axvline(np.log10(xi_fit), color="orange", linestyle=":",
linewidth=1.2, label=f"σ = ξ = {xi_fit:.2e}")
ax.set_xlabel("log₁₀(σ) (migration rate)")
ax.set_ylabel("Temporal correlation C(σ)")
ax.set_title("Two-patch SIR: synchronisation vs. migration rate\n"
"(Keeling & Rohani 2002)")
ax.legend(fontsize=9)
ax.set_ylim(-0.1, 1.1)
plt.tight_layout()
plt.show()
print(f"\nMidpoint synchronisation (C=0.5) at σ ≈ ξ = {xi_fit:.4e}")
--------------------------------------------------------------------------- ValueError Traceback (most recent call last) Cell In[5], line 8 4 return sigma / (xi + sigma) 5 6 # Only fit over points where correlation is non-negative and finite 7 valid = np.isfinite(correlations) & (correlations >= 0) ----> 8 popt, _ = curve_fit(sync_curve, sigma_values[valid], correlations[valid], p0=[1e-3]) 9 xi_fit = popt[0] 10 print(f"Fitted desynchronisation parameter ξ = {xi_fit:.4e}") 11 File ~/projects/laser-fresh/laser.cohorts/.docs/lib/python3.12/site-packages/scipy/optimize/_minpack_py.py:949, in curve_fit(f, xdata, ydata, p0, sigma, absolute_sigma, check_finite, bounds, method, jac, full_output, nan_policy, **kwargs) 946 xdata = np.asarray(xdata, float) 948 if ydata.size == 0: --> 949 raise ValueError("`ydata` must not be empty!") 951 # nan handling is needed only if check_finite is False because if True, 952 # the x-y data are already checked, and they don't contain nans. 953 if not check_finite and nan_policy is not None: ValueError: `ydata` must not be empty!