Different results when changing nb_steps

[maurosilber]

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")

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: https://github.com/Armavica/rebop/blob/5ecd326afec91846c28320fcce208645eb38c804/src/lib.rs#L347 which can advance the simulation time without performing any reaction: https://github.com/Armavica/rebop/blob/5ecd326afec91846c28320fcce208645eb38c804/src/gillespie.rs#L278-L281 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!