[Python] Test coverage for Gillespie algorithm

I've written a single time step function for Gillespie algorithm for Lotka-Volterra equations.

Here is the code:

def gillespie_step(prey: int, predator: int, alpha: float, beta: float, gamma: float) -> tuple[float | None, int, int, str | None]:
"""Perform one step of the Gillespie simulation for a predator-prey system.

Args:
prey: Current population of prey, integer.
predator: Current population of predators, integer.
alpha: Prey birth rate, float.
beta: Predation rate, float.
gamma: Predator death rate, float.

Returns:
A tuple containing:
- time_step: Time duration until next event occurs (float), or None if no event occurs
- prey: Updated population of prey (integer)
- predator: Updated population of predators (integer)
- event: String describing the event ("prey_birth", "predation", or "predator_death"),
or None if no event occurs
"""
event = None
time_step = None

birth_rate = alpha * prey
predated_rate = beta * prey * predator
predator_death_rate = gamma * predator

total_rate = birth_rate + predated_rate + predator_death_rate

if total_rate == 0:
pass
else:
# Sample time
time_step = np.random.exponential(1.0 / total_rate)
q = np.random.rand()

if q <= birth_rate / total_rate:
prey += 1
event = "prey_birth"
elif q <= (birth_rate + predated_rate) / total_rate:
prey -= 1
predator += 1
event = "predation"
else:
predator -= 1
event = "predator_death"

return time_step, prey, predator, event


And here are the tests:

"""### Test Cases"""

prey, predator = 200, 200
alpha, beta, gamma = 2., 0.01, 3.
np.random.seed(0)
time_step, prey, predator, event = gillespie_step(prey, predator, alpha, beta, gamma)
a, b, c, d = 0.00056841791, 200, 199, "predator_death"
assert np.allclose(time_step, 0.00056841791, rtol= 1e-4, atol = 1e-4)
assert np.allclose(prey, 200, atol = 1e-4, rtol = 1e-4)
assert np.allclose(predator, 199, atol = 1e-4, rtol = 1e-4)
assert event == "predator_death"

prey, predator = 100, 20
alpha, beta, gamma = 3., 0.05, 1.
np.random.seed(1)
time_step, prey, predator, event = gillespie_step(prey, predator, alpha, beta, gamma)
a, b, c, d = 0.0012847758029980606, 99, 21, "predation"
assert np.allclose(time_step, 0.0012847758029980606, atol = 1e-4, rtol = 1e-4)
assert np.allclose(prey, 99, atol = 1e-4, rtol = 1e-4)
assert np.allclose(predator, 21, rtol = 1e-4, atol = 1e-4)
assert event == "predation"

prey, predator = 100, 10
alpha, beta, gamma = 1., 0.005, 5.
np.random.seed(2)
time_step, prey, predator, event = gillespie_step(prey, predator, alpha, beta, gamma)
a, b, c, d = 0.003694787024320097, 101, 10, "prey_birth"
assert np.allclose(time_step, 0.003694787024320097, atol = 1e-4, rtol = 1e-4)
assert np.allclose(prey, 101, atol = 1e-4, rtol = 1e-4)
assert np.allclose(predator, 10, atol = 1e-4, rtol = 1e-4)
assert event == "prey_birth"

prey, predator = 0, 0
alpha, beta, gamma = 1., 0.005, 1.
np.random.seed(59)
time_step, prey_out, predator_out, event = gillespie_step(prey, predator, alpha, beta, gamma)
assert time_step is None
assert prey_out == 0
assert predator_out == 0
assert event is None

prey, predator = 0, 10
alpha, beta, gamma = 2., 0.01, 1.
np.random.seed(3)
time_step, prey_out, predator_out, event = gillespie_step(prey, predator, alpha, beta, gamma)
assert event == "predator_death"

prey, predator = 10, 0
alpha, beta, gamma = 1., 0.05, 2.
np.random.seed(4)
time_step, prey_out, predator_out, event = gillespie_step(prey, predator, alpha, beta, gamma)
assert event == "prey_birth"

prey, predator = 100, 100
alpha, beta, gamma = 0., 0., 0.
time_step, prey_out, predator_out, event = gillespie_step(prey, predator, alpha, beta, gamma)
assert time_step is None
assert prey_out == 100
assert predator_out == 100
assert event is None


I have used the python coverage.py and also simple printing but QA of the project says I have only 70% test coverage: subproblem 1 test coverage is 68.3%. Please ensure all code paths are tested. Please note that I can't have more than one return statement per specification of the problem. Any thoughts?

Continue reading...