Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Different results when changing nb_steps #26

Open
maurosilber opened this issue Sep 7, 2024 · 2 comments
Open

Different results when changing nb_steps #26

maurosilber opened this issue Sep 7, 2024 · 2 comments

Comments

@maurosilber
Copy link
Contributor

Hi,

I get different results when running the example from the README with differing number of steps (and a fixed seed).

import rebop
from matplotlib.cm import viridis as cmap

sir = rebop.Gillespie()
sir.add_reaction(1e-4, ["S", "I"], ["I", "I"])
sir.add_reaction(0.01, ["I"], ["R"])

values = {"S": 999, "I": 1}
tmax = 250
seed = 1

ax = None
steps = [0, 10, 50, 100, 500, 1000]
for step, color in zip(steps, cmap.resampled(len(steps)).colors):
    ax = (
        sir.run(values, tmax=tmax, nb_steps=step, seed=seed)
        .to_dataframe()
        .sort_index(axis="columns")
        .plot(ax=ax, color=color, subplots=True, legend=False)
    )
ax[1].legend(steps, title="Steps", bbox_to_anchor=(1, 0.5), loc="center left")

image

I assumed that nb_steps would save a subset of the solution by interpolating at equispaced points.

I think the issue could be in using this method:

g.advance_until(t);

which can advance the simulation time without performing any reaction:

rebop/src/gillespie.rs

Lines 278 to 281 in 5ecd326

self.t += self.rng.sample::<f64, _>(Exp1) / total_rate;
if self.t > tmax {
self.t = tmax;
return;

after which the simulation is restarted from that intermediate tmax and does a different time jump, instead of doing the original jump that would have taken the simulation to a self.t > tmax.

Does this process still produce a different but correct sample of the master equation?

Thanks!

@maurosilber
Copy link
Contributor Author

It seems to produce the same distribution (at least for the SIR model).

I performed 1000 simulations of the SIR using:

  1. nb_steps=250 (blue)
  2. nb_steps=0 (orange), which was later interpolated at the same time points of 1.

and compared the distributions at the same time points:

image

where each curve corresponds to a different quantile.

Here is the script to produce the figure:

import numpy as np
import rebop
import xarray

sir = rebop.Gillespie()
sir.add_reaction(1e-4, ["S", "I"], ["I", "I"])
sir.add_reaction(0.01, ["I"], ["R"])

values = {"S": 999, "I": 1}
tmax = 250

quantiles = np.linspace(0, 1, 11)

seeds = range(1000)
x_steps = (
    xarray.merge(
        [
            sir.run(values, tmax=tmax, nb_steps=250, seed=seed).to_dataarray(
                "species", name=seed
            )
            for seed in seeds
        ]
    )
    .to_dataarray("seed")
    .quantile(quantiles, dim="seed")
)
x_interpolated = (
    xarray.merge(
        [
            sir.run(values, tmax=tmax, nb_steps=0, seed=seed)
            .interp(time=x_steps.time)
            .to_dataarray("species", name=seed)
            for seed in seeds
        ]
    )
    .to_dataarray("seed")
    .quantile(quantiles, dim="seed")
)

ax = x_steps.plot.line(x="time", col="species", add_legend=False, color="C0")
for ax, s in zip(ax.axes.flat, ax.name_dicts.flat):
    x_interpolated.sel(s).plot.line(
        x="time", ax=ax, add_legend=False, color="C1", linestyle="--"
    )

@Armavica
Copy link
Owner

Armavica commented Sep 9, 2024

Hi, thank you for taking the time to report this!

I was aware of this behavior, it is a consequence of the way I implemented the algorithm. But all samples returned are mathematically correct traces, as if you were using different seeds.

This is something that I have been wanting to change though, because it breaks the principle of least surprise. I will rewrite this part for the next release, to make the trace independent on nb_steps.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

2 participants