setting event_seed for poisson drive doesn't produce consistent spiking between simulations

[rythorpe]

I've tried this on two computers now (both Ubuntu 22.04) and am a bit baffled. When running the following script, roughly ~4/5 times I get the the first figure, while ~1/5 times I get the second figure. Note that both simulations are very similar, but slightly different.

import matplotlib.pyplot as plt

import hnn_core
from hnn_core import simulate_dipole, jones_2009_model, MPIBackend

net = jones_2009_model()

weights_ampa = {'L2_pyramidal': 0.001, 'L5_pyramidal': 0.01}
rate_constant = {'L2_pyramidal': 100.0, 'L5_pyramidal': 100.0}
net.add_poisson_drive(
    'poisson', rate_constant=rate_constant, weights_ampa=weights_ampa,
    location='proximal', synaptic_delays=0.1, event_seed=1)

with MPIBackend(n_procs=12):
    dpls = simulate_dipole(net, tstop=20.)
scaling_factor = 30000
dpls = [dpl.scale(scaling_factor) for dpl in dpls]  # scale in place

net.cell_response.plot_spikes_raster()

[rythorpe]

Could someone try reproducing this @jasmainak @ntolley @chenghuzi? I suspect there might be something weird going on deep in MPIBackend (possibly from nrn objects that don't get destroyed properly)?

[rythorpe]

Update: same behavior for MPIBackend(n_procs=1), so the issue isn't with MPIBackend specifically.

[rythorpe]

Here are the two different spike profiles plotted alongside the spike events of the poisson drive. The two drive sequences look identical to me.

[chenghuzi]

On my server in a notebook, it consistently outputs the second one:

[rythorpe]

Did you try running it on a few different instantiations of your jupyter kernel (i.e., restart your kernel then run, restart then run, ...)? I've only noticed the variance on different instantiations of the python/ipython shell.

[chenghuzi]

You are right. After restarting the kernel, I also got two results as you did. However, the ratio seems to be different (around 1:1), although I have not counted the frequency explicitly. Are there any other sources of randomness in this code? I tried to add an extra random seed at the beginning, but it didn’t work. Or, is there something hidden in NEURON?

[rythorpe]

The ratio might be more evenly balanced than I first indicated - it was just my hasty estimate after the first few coin tosses.

That's the thing, there aren't. At last not on the python side. There's got to be something weird going on one the NEURON/HOC side. Rounding error or leftover pointer floating around somewhere maybe?

[jasmainak]

I didn’t find time to replicate. But just to be sure, the input spike times are exactly the same in both cases?

On Tue 2 May 2023 at 22:26, Ryan Thorpe @.***> wrote:

The ratio might be more evenly balanced than I first indicated - it was just my hasty estimate after the first few coin tosses.

That's the thing, there aren't. At last not on the python side. There's got to be something weird going on one the NEURON/HOC side. Rounding error or leftover pointer floating around somewhere maybe?

— Reply to this email directly, view it on GitHub https://github.com/jonescompneurolab/hnn-core/issues/641#issuecomment-1532382652, or unsubscribe https://github.com/notifications/unsubscribe-auth/ADY6FIVI6LAOHQWGMGQSW7TXEG67DANCNFSM6AAAAAAXSD65DM . You are receiving this because you were mentioned.Message ID: @.***>

-- Sent from my iPhone

[ntolley]

Wait so does reproducing this require restarting the kernel? Or can do you get different results event when you run multiple times in the same instance

[rythorpe]

Wait so does reproducing this require restarting the kernel? Or can do you get different results event when you run multiple times in the same instance

As far as I've observed, you need to restart the kernel or use a different python shell instance.

[rythorpe]

I didn’t find time to replicate. But just to be sure, the input spike times are exactly the same in both cases?

By visual inspection of the spike time raster shown above, they're exactly the same. I'm still trying to quantify this in order to know for sure though.

[jasmainak]

If you can save the poisson drive times in double precision and check with numpy allclose or array_almost_equal that they are the same, we can exclude any issues due to random seeds

[rythorpe]

I'm investigating that right now actually! They're almost the same - it looks almost like a few different drive spike events are getting assigned to different gids. Stay tuned...

[rythorpe]

Okay, so the events stored in net.external_drives['poisson']['events'] are always the same. If I take the poisson drive spikes from net.cell_response.spike_times, however, they are only consistent after I sort them. Recall that the spike times associated with the different gids are appended to each other, so variation in the contents of net.cell_response.spike_times implies a difference in the spike event -> gid assignment.

*updated per https://github.com/jonescompneurolab/hnn-core/issues/641#issuecomment-1535485471

[jasmainak]

was this always broken? can you check if this happens, for e.g., in a commit from 6 months ago?

I don't see anything stochastic in the gid assignment ... could it be that the gid_dict is different?

[rythorpe]

Okay, so the events stored in net.external_drives['poisson']['events'] are always the same. If I take the poisson drive spikes from net.cell_response.spike_times, however, they are only consistent after I sort them. Recall that the spike times associated with the different gids are appended to each other, so variation in the contents of net.cell_response.spike_times implies a difference in the spike event -> gid assignment.

Ooops, half of what I said here was wrong. (There is apparently a slight inconsistency with how f-strings were rounding vs. numpy.) Down to the 14th decimal place, net.external_drives['poisson']['events'] is always consistent and the poisson spike times in net.cell_response.spike_times are always consistent.

I'll try this on an older commit, but this has got me pretty stumped. I honestly don't know where the variation is coming from. If I have time to try whittling down the model, I guess the next step would be see if variability is on the nrn side.