BrokenProcessPool with sample_smc example

1

Hello,
I am running the example provided in pymc website for Aproximate Bayesian Computation. This is the code I run:

from scipy.integrate import odeint

# Definition of parameters
# a = 1.0  # These values are now passed as arguments to competition_model
# b = 0.1
# c = 1.5
# d = 0.75

# Lotka - Volterra equation
def dX_dt(X, t, a, b, c, d):
    """Return the growth rate of fox and rabbit populations."""
    return np.array([a * X[0] - b * X[0] * X[1], -c * X[1] + d * b * X[0] * X[1]])

# simulator function
def competition_model(rng, a, b, size=100, X0=[10.0, 5.0], time=15, c=1.5, d=0.75):
    """
    Competition model with parameters a, b, c, and d.

    Includes c and d, X0, size, and time as arguments with default values
    to avoid referencing external variables.
    This ensures that the correct values are used within the multiprocessing context.
    """
    t = np.linspace(0, time, size)
    return odeint(dX_dt, y0=X0, t=t, rtol=0.01, args=(a, b, c, d))

# function for generating noisy data to be used as observed data.
def add_noise(a, b, size=100, X0=[10.0, 5.0], time=15, c=1.5, d=0.75):
    """
    Add noise to the competition model output.

    Includes c and d, X0, size, and time as arguments with default values
    to ensure consistency with competition_model.
    """
    noise = np.random.normal(size=(size, 2))
    simulated = competition_model(None, a, b, size=size, X0=X0, time=time, c=c, d=d) + noise
    return simulated


with pm.Model() as model_lv:
    a = pm.HalfNormal("a", 1.0)
    b = pm.HalfNormal("b", 1.0)

    sim = pm.Simulator("sim", competition_model, params=(a, b), epsilon=10, observed=observed)

    idata_lv = pm.sample_smc()

Running on PyMC v5.20.1

But i get the followin error from the simulator:

BrokenProcessPool: A process in the process pool was terminated abruptly while the future was running or pending.

3

Hello,

Well the example in the page also does not work for me but it seems to be a sizing issue. The simulator class seems to resize everything to two dimensional arrays so for instance constants a,b become 1x1 arrays which breaks the dimension of odeint returning 1 x 1 x N sized arrays and an error. If I slightly modify it, it works for me and produces roughly the same results as in the page. I wonder if this resizing is an intended feature because it seems to make life generally difficult.

Figure 2025-03-25 162135
Figure 2025-03-25 1621351384×572 64.3 KB

The page needs to be updated though.

import numpy as np
import pymc as pm
import arviz as az
from scipy.integrate import odeint

# Definition of parameters
a = 1.0
b = 0.1
c = 1.5
d = 0.75

# initial population of rabbits and foxes
X0 = [10.0, 5.0]
# size of data
size = 100
# time lapse
time = 15
t = np.linspace(0, time, size)


# Lotka - Volterra equation
def dX_dt(X, t, a, b, c, d):
    """Return the growth rate of fox and rabbit populations."""

    return np.array([a * X[0] - b * X[0] * X[1], -c * X[1] + d * b * X[0] * X[1]])


# simulator function
def competition_model(rng, a, b, size=None):
  
    if isinstance(a, np.ndarray) and a.ndim==2:
      a = a[0,0]
    if isinstance(b, np.ndarray) and b.ndim==2:
      b = b[0,0]
  
    return odeint(dX_dt, y0=X0, t=t, rtol=0.01, args=(a, b, c, d))
  
# function for generating noisy data to be used as observed data.
def add_noise(a, b):
    noise = np.random.normal(size=(size, 2))
    simulated = competition_model(None, a, b) + noise
    return simulated

observed = add_noise(a, b)
with pm.Model() as model_lv:
    a = pm.HalfNormal("a", 1.0)
    b = pm.HalfNormal("b", 1.0)

    sim = pm.Simulator("sim", competition_model, params=(a, b), epsilon=10, observed=observed)
    sim.eval()
    idata_lv = pm.sample_smc()

az.plot_trace(idata_lv, kind="rank_vlines")