I have to store and restore my simulations starting from 10k seconds. The C++ code is substantially faster for me, and I would like to use it instead of the native version. I recognize that the standard store and restore are not supported in C++ standalone. I have written a basic version of the code that, in principle, should work in C++, but I am clearly missing something. I have attached the minimal code, built around the standard example of the stdp code.
A few things that are important to me,
a. I wish to store 10 seconds of simulations (spikes, weights) every 1000 seconds or so - I am aware of tips and tricks to get this, i.e., to dump directly into a txt file. But for now, I use multiple runs.
b. Refractory period of my neurons is a dynamic variable and a neuron-specific property.
c. The random number generator seed should be the same. Potentially, my issue is in the Poisson process (restore_random_state equivalent in C++?)
Minimal code to reproduce problem
Expected output
The results I hope/expect to be the same at
a) The end of the simulation run for 200 seconds and (starting from 0)
b) When I run for 100 seconds, store and restore, run for another 100 seconds (starting from 100 seconds).
Actual output
Pasted the difference in the traces at the end of 200 seconds in the gist link.
Hi @ccluri. This is an excellent question. Getting something comparable to store/restore to work in C++ standalone mode is quite high on my priority list, but unfortunately we are not quite there yet. I did not verify things in detail, but from reading your code your approach looks generally correct to me. One thing that the get_state mechanism cannot recover is the internal state of the spike queue (spikes that are “in transit” and haven’t been delivered yet), but this shouldn’t matter in your example, since you don’t have any synaptic delays. I do not know whether this explains all the difference, but your PoissonGroup is definitely not the same between the 200s and the two 100s runs. The seed statement is included in the generated code, which means that if you run the two 100s runs, the second run will have exactly the same random numbers as the first run, in contrast to the second half of the 200s run which of course has independent random numbers.
While it will take a bit of time to implement a proper store/restore mechanism, I’ve made the internal state of the random generator accessible from the C++ code, and it now gets stored at the end of a simulation to the “results” directory. Note that this is only available in the current development version of Brian, not in the released version. I am planning to do a new release next week, but if you want to try it out before, you can install the development version as explained here: Installation — Brian 2 0.0.post128 documentation
With this version, you will get a file random_generator_state in the results directory at the end of a run. As I said before, we don’t yet have a “nice” mechanism to restore this state, but the following (inserted somewhere before device.build) should work by injecting C++ code for rng state restoration into the generated code:
device.insert_code("before_network_run","""
std::ifstream _rng_state("results/random_generator_state");
if (_rng_state.is_open()) {
_rng_state >> brian::_random_generators[0];
} else {
std::cerr << "Could not open random generator state file" << std::endl;
}""")
Thank you for your prompt reply and for your support.
before_network_run will apply this C++ code to all the times before the device.run calls. And it fails to compile.
So I tested it with a single run; it compiled and ran, but the results still don’t converge. I have updated the gist with the new version of the code, incorporating your suggestions.
I tried to trace where the error is creeping in from, but I’ve failed. Does this way of doing this break the initial part of the code?
Something I wondered,
a) Perhaps there should be a before_build, and before_first_run options that can be considered. I have no clue if that breaks other software architecture ideas you have.
b) If it is possible to have a ‘setter’ for network.t (time of the simulation)?
c) If it is helpful to include the network time in the random_generator_state, perhaps in the filename itself?
I had another look, and I noticed one small thing is missing (but in a quick test, it did not seem to be enough to get the result you want): the synapses have a lastupdate attribute, which is used to calculate the decay of the trace variables when a spike arrives (the event-driven equations are transformed as explained here: Synapses — Brian 2 0.0.post128 documentation). You should therefore restore these values as well (and do the same correction as for lastspike, i.e. subtract the previous run time), otherwise the first spikes for each synapse in the second run will not have the correct effect.
Regarding your remarks:
We do actually already have after_start slot, which will only insert the code once, after the arrays have been created and initialized. The reason why I used before_network_run is that after_start is before the seed is applied. So if you simply replaced before_network_run by after_start, restoring the random generator state wouldn’t have any effect, since it would be immediately overwritten by the seed value. You’d therefore have to make sure to remove the seed call in the run where you restore the rng state (as a side note, you can also manually set the seed in code by using brian::_random_generators[0].seed(some_value);, maybe that could be an alternative).
In general, we don’t want to manually set the time, since the semantics aren’t very clear if you do this during a run. But it might be good to support it as part of the run_args mechanism, I guess ? At the start of a new simulation, there wouldn’t be any harm in changing the time. Of course, the longer term plan is to support a “proper” store/restore mechanism which would handle the time as well.
Not sure I get this point – what would this be helpful for?