full trajectory

[dave-doty]

I'm sorry if this isn't the correct place to ask this...

Is there a way to get the full "trajectory", i.e., the entire sequence of reactions that Gillespie samples, together with the inter-reaction times (sampled exponential random variables), rather than sampling from fixed time points, in which several reactions (or none) might have occurred between two consecutive time points?

In particular, I'm hoping to use the Python interface for this, so if that is possible from the Python API, that would be great.

[dave-doty]

I looked into the source code and it doesn't seem to be supported, there's this loop that goes through many reactions without recording the intermediate configurations:

rebop/src/gillespie.rs

Lines 238 to 264 in a34e193

	loop {
	//let total_rate = make_rates(&self.reactions, &self.species, &mut rates);
	let total_rate = make_cumrates(&self.reactions, &self.species, &mut rates);
	// we don't want to use partial_cmp, for performance
	#[allow(clippy::neg_cmp_op_on_partial_ord)]
	if !(0. < total_rate) {
	self.t = tmax;
	return;
	}
	self.t += self.rng.sample::<f64, _>(Exp1) / total_rate;
	if self.t > tmax {
	self.t = tmax;
	return;
	}
	let chosen_rate = total_rate * self.rng.gen::<f64>();
	//let ireaction = choose_rate_sum(chosen_rate, &rates);
	//let ireaction = choose_rate_for(chosen_rate, &rates);
	let ireaction = choose_cumrate_sum(chosen_rate, &rates);
	//let ireaction = choose_cumrate_for(chosen_rate, &rates);
	//let ireaction = choose_cumrate_takewhile(chosen_rate, &rates);
	// here we have ireaction < self.reactions.len() because chosen_rate < total_rate
	let reaction = unsafe { self.reactions.get_unchecked(ireaction) };
	reaction.1.affect(&mut self.species);
	}

So this is a feature request to support getting this full trajectory information out of the simulation.

[Armavica]

Thank you for your message! That's exactly the right place, I'll be happy to add this to the next release.

[Armavica]

@dave-doty I added the feature that you requested, that you can activate by passing nb_steps=0 in the Python run function. It will stop at the first reaction after tmax. I am lacking time currently to fine tune the performance or the API of this function, so these might change in a future release: notably, I am wondering if I should always return tmax as a last point (which never corresponds to a reaction time), or the first reaction time after tmax (which corresponds to a reaction time when it is finite, and doesn't when it is infinite).
Please let me know if that works for you, or if you think that this could be improved! I will now release version 0.8.0 so you can install it with pip.

[dave-doty]

Thanks! Seems to work great.

One question: is the volume parameter supported directly, or should this be handled by adjusting rate constants (e.g., to implement volume $20$, divide bimolecular rate constants by $20$, and trimolecular rate constants by $20^2$, etc.).

[dave-doty]

BTW I think returning only the actual reaction times (not tmax) makes the most sense, I think code would expect only actual reaction times to appear. I also like your convention of adding one more time at "inf" to represent that the system reached a terminal state before time tmax.

[dave-doty]

Can I also say how awesome it is that you implemented this in one day? We found your package after trying to get Gillespy2 to do the same thing, with the discussion there suggesting that the problem is uncomputable: StochSS/GillesPy2#427

[Armavica]

One question: is the volume parameter supported directly, or should this be handled by adjusting rate constants (e.g., to implement volume 20, divide bimolecular rate constants by 20, and trimolecular rate constants by 202, etc.).

Volume is not currently supported, but that could be the next feature, if there is interest for it! If you want to open a new issue for this, it will remind me to work on it at some point.

Thank you for your feedback on the API! I am glad that it works for you.

Can I also say how awesome it is that you implemented this in one day? We found your package after trying to get Gillespy2 to do the same thing, with the discussion there suggesting that the problem is uncomputable

Thank you 😊 Most of the functionality was already there, I just had to make it return all the events generated. They are right that this might be slower for simulations with a huge number of steps happening (mainly because of the need to allocate a large amount of memory to store all this data), and computationally it's anyway always going to be a bit faster to go through steps in a loop without worrying to return all the results.

Finally, something that I am not quite satisfied about currently is that even with a fixed seed, the returned trajectory will not be the same if you choose different nb_steps values. They are all correct trajectories, but this is a consequence of the current implementation of the simulation loops with nb_steps > 0. I would like to revisit this one day, to make all discretizations return the exact same trajectory.

[dave-doty]

Finally, something that I am not quite satisfied about currently is that even with a fixed seed, the returned trajectory will not be the same if you choose different nb_steps values. They are all correct trajectories, but this is a consequence of the current implementation of the simulation loops with nb_steps > 0. I would like to revisit this one day, to make all discretizations return the exact same trajectory.

It's also the case that if you leave nb_steps alone, but change tmax, then it changes the random outcome with the same seed.

It would be very helpful for reproducing bugs to be able to increase tmax, to see the CRN a bit further out, while having it behave exactly the same up to the original tmax.

[Armavica]

Yes, that makes a lot of sense. Someone else mentioned the same thing recently (#26), so it's definitely the next thing that I am going to do in rebop.