[WIP] Fix bugs in read_params workflow

[katduecker]

I have been experiencing issues when loading in network parameters from a json file, somehow the simulation always looked different. There were a couple of small bugs due to hard-coded synaptic delays in params.py. There is also an update of the event_seeds that is only applied when reading params from a file, but not when defining an event seed. I have commented it out for now but we'd have to decide if it's needed.

[asoplata]

Thanks for finding these Katharina, these are both quite the pickle(s).

In both cases, it appears that both of these issues only occur when using the add_drives_from_params argument set to True in jones_2009_model(add_drives_from_params=True,...). In order for this PR to not block your work, one possible but annoying workaround would be that you could reproduce the evdist* and evprox* drives you need that are made by this option, but "by hand" using Network.add_evoked_drive. You could then save your Network to a "hierarchical JSON" file using hnn_io.write_network_configuration and read it back in for new simulations using hnn_io.read_network_configuration. This would allow you to avoid relying on the problematic code in drives.py::_add_drives_from_params and params.py::create_pext. If you are using custom "flat JSON" files like param/default.json instead of "hierarchical JSON" files like param/jones2009_base.json, then this may not work and we will need to find another workaround until this issue is fully resolved. I apologize since this workaround kind of reads as "do it [the drives] yourself", but that may be our best option depending on if we're willing to change this code or not (we could definitely help you with this if you wanted).

CC @ntolley and @dylansdaniels : The bugs/unexpected behavior you've identified are legitimate. However, if we merged both of these changes, that would have at least small or possibly large impacts on others' simulation outputs, since this would in essence be changing the default model in a scientific sense, and would affect reproducibility. We absolutely need to discuss this at the next HNN Dev meeting, since it is not a change to be made lightly. For reproducibility reasons, and to not suddenly force scientific output of simulations by others to change, we may want to not merge your changes.

Depending on future discussion on this issue, one alternative could include: we change jones_2009_model's add_drives_from_params argument to more explicitly be a "legacy" argument, and also add a new argument that does similar functionality, but includes your fixes.

For more treatment of how we want to handle changes to the kinda-legacy params code in general, see my recent comment here #962 (comment) .

[katduecker]

Hi Austin, I was actually using this option for the new model, which I currently can't distribute with the hierarchical json, because the read_network_configuration function builds a jones model...

At the very least we should change the hard-coded synaptic delays, because there were some typos (some of the L5 delays were set to 0.1, but they should be 1). I don't know if this will change things dramatically. The seeding, however, might, especially if a developer has been using a specific seed to tune their model... Also if the seed is overridden, then why define it in the json file in the first place...

[asoplata]

Argh, right, I forgot about the read_network_configuration hard-coding we discussed before. I also completely agree with your comment on the seed; like you said before, for reproducibility purposes, there should always be a user-settable seed, and what the user provides should not be overridden. That in particular seems to be "technical debt" catching up with us.

[asoplata]

I have made a 95% complete PR that applies one possible solution at katduecker#2 .

[asoplata]

This comment first documents how the implicit event_seed's (called "seeds") of various external drives work currently (currently as master 90f4d63 ), in the case of passing add_drives_from_params=True to jones_2009_model(), both with and without legacy_mode=True. This comment then later describes what (few) changes we are making to enable legacy_mode to work as it did before, while enabling non-legacy_mode to use user-provided seeds without modification, enabling reproducibility. This comment is also intended to act as documentation for these changes, and to be referred to by comments in the code. You'll want to sit down for this one.

It is easiest if we begin by explaining how things currently work on master, as of 90f4d63 :

Current add-drives, no legacy

Let's assume you are doing the following code example:

net = jones_2009_model(add_drives_from_params=True)  # NO legacy_mode
dpl = simulate_dipole(net, tstop=170.0)

As far as I can tell, all code for loading, changing, and setting the unique event_seed of each drive is contained entirely in the appropriately-named drives.py::_add_drives_from_params() here. Let's go through this function.

1. First, we call params.py::_extract_drive_specs_from_hnn_params. This loads our parameters from param/default.json, where the initial seed for all drives is 2:
   

   hnn-core/hnn_core/param/default.json

   Lines 18 to 20 in 90f4d63

   	"prng_seedcore_evprox_1": 2,
   	"prng_seedcore_evdist_1": 2,
   	"prng_seedcore_evprox_2": 2,

   
   Because legacy_mode is off, we only load the three typical drives, along with their initial seeds:
   • evdist1: 2
   • evprox1: 2
   • evprox2: 2

2. Next, each drive is then added to the Network in the main for-loop. I will note, however, that they are added with the above event_seed, but these seeds will be overwritten later.

3. Finally, we get to the bottom drive-related block. This adds net.gid_ranges[drive_name][0] to the seeds of each of these three drives. In the default case, this changes the seeds such that they become the following:

   • evdist1: 272 = 2 (initial event_seed) + 270 (initial gid)
   • evprox1: 507 = 2 (initial event_seed) + 505 (initial gid)
   • evprox2: 777 = 2 (initial event_seed) + 775 (initial gid)

4. As far as I can tell, there are no more changes to the drives' seeds, and these become the final seed values used in this case. These are probably the seeds that many people are actually using, unless they add their drives manually or specify their own params dictionary or "flat JSON" file.

Current add-drives, yes legacy

legacy_mode=True makes things much more complicated. Let's assume you are doing the following code example, and go through drives.py::_add_drives_from_params() again.

net = jones_2009_model(add_drives_from_params=True, legacy_mode=True)
dpl = simulate_dipole(net, tstop=170.0)

1. First, we call params.py::_extract_drive_specs_from_hnn_params. This loads our parameters from param/default.json, where the initial seed for all drives is 2:
   

   hnn-core/hnn_core/param/default.json

   Lines 18 to 20 in 90f4d63

   	"prng_seedcore_evprox_1": 2,
   	"prng_seedcore_evdist_1": 2,
   	"prng_seedcore_evprox_2": 2,

   
   However, legacy_mode makes further changes to our three "typical" drives and decrements them by 18 here. legacy_mode also adds more drives than just the typical three (this will be important later). The drives' seeds at this point are as follows:
   • bursty1: 2
   • bursty2: 2
   • evdist1: -16 = 2 (initial event seed) - 18 (legacy offset)
   • evprox1: -16 = 2 (initial event seed) - 18 (legacy offset)
   • evprox2: -16 = 2 (initial event seed) - 18 (legacy offset)
   • extgauss: 2
   • extpois: 2

2. Next, each drive is then added to the Network in the main for-loop like before. Again, its seeds at this point will be overwritten later.

3. Finally, we get to the bottom drive-related block. This adds net.gid_ranges[drive_name][0] to the seeds of every drives. However, since bursty1 and bursty2 are drives specific to legacy_mode, they are created first. This has a downstream effect where it changes the gid_range-related addition to the seeds of ALL other drives, including the typical three. Yet another difference is that in the legacy_mode=True case, evdist1 has a gid_range length of 270 instead of the typical 235, which will cause evprox1 and evprox2's initial gid ranges to begin at a higher number than in the non-legacy_mode case. In summary, what this means is that the seeds of the typical three drives are very different from the legacy_mode=False case. The drives' seeds at this point are as follows:

   • bursty1: 272 = 2 (initial event seed) + 270 (initial gid)
   • bursty2: 282 = 2 (initial event seed) + 280 (initial gid)
   • evdist1: 274 = 2 (initial event seed) + 290 (initial gid) - 18 (legacy offset)
   • evprox1: 544 = 2 (initial event seed) + 560 (initial gid) - 18 (legacy offset)
   • evprox2: 814 = 2 (initial event seed) + 830 (initial gid) - 18 (legacy offset)
   • extgauss: 1102 = 2 (initial event seed) + 1100 (initial gid)
   • extpois: 1372 = 2 (initial event seed) + 1370 (initial gid)

4. As far as I can tell, there are no more changes to the drives' seeds, and these become the final seed values used in this case.

New event_seed code

In the new commit katduecker@1931f58 , we make the following changes, which will allow for per-drive seeds to be faithfully used, except in the case when both legacy_mode and add_drives_from_params are True. We:

1. Update the param/default.json file to use the true seeds in the non-legacy_mode case, i.e.:

   "prng_seedcore_evprox_1": 507,
   "prng_seedcore_evdist_1": 272,
   "prng_seedcore_evprox_2": 777,
   

2. Change the gid_range addition at

   hnn-core/hnn_core/drives.py

   Line 241 in 90f4d63

   for drive_name, drive in net.external_drives.items():

   so that only uses the old "gid_range" addition method when legacy_mode=True, AND so it uses the original values of the seeds of the three typical drives at that point in the code. To do this, we change the code from

   for drive_name, drive in net.external_drives.items():
       drive["event_seed"] += net.gid_ranges[drive_name][0]

   to:

   if net._legacy_mode:
       for drive_name, drive in net.external_drives.items():
           if drive_name in ("evdist1", "evprox1", "evprox2"):
               drive["event_seed"] = -16
           drive["event_seed"] += net.gid_ranges[drive_name][0]

3. Regenerate our default hierarchical-JSON networks accordingly.

[asoplata]

Note that so far, only the event_seed case has been dealt with. The synaptic delay changes have been reverted temporarily since they break valid tests.

[asoplata]

For the synaptic delays...we may want to split that into its own PR/issue. As of the most recent commit b94978a , the only remaining test that fails is one where we test that the dipole output is identical. Specifically, that test is test_parallel_backends.py::TestParallelBackends::test_compare_hnn_core

hnn-core/hnn_core/tests/test_parallel_backends.py

Line 371 in 90f4d63

def test_compare_hnn_core(self, run_hnn_core_fixture, backend, n_jobs=1):

In that test, we compare the dipole output of the typical sim to existing dipole output data here, on the test_data branch https://github.com/jonescompneurolab/hnn-core/blob/test_data/dpl.txt . Unsurprisingly, since the network is ever-so-slightly different, the dipoles are not an exact match. However, when I plotted them, they were more different than I was expecting, as shown in the below figure. In this figure, dpl is the existing dipole test data (i.e. excluding the synaptic delay changes), while dpl2 is the new dipole (i.e. including the synaptic delay changes).

To my untrained eye, these are different enough for concern, but I'm not an expert at interpretation of this. Possibilities for what to do next:

1. We don't make this synaptic delay change at all: probably a bad idea.
2. We decide this is close enough, and upload new dipole data to be the new default that we test against. In this case, we would add the output dpl2.txt data to a new file named something like dpl-updated-synaptic-delays.txt or whatever, on the test_data branch.
3. We "kick the can down the road": we separate the synaptic delays out from this issue/PR into its own, and investigate the differences more thoroughly. The first thing I would do in this case is reproduce everything in the ERP tutorials and check that they are qualitatively similar. Such an investigation could be as detailed as the investigator decides it should be.

[asoplata]

Decision Point: @ntolley @katduecker @dylansdaniels we need to decide what to do about the changes to the dipole from the synaptic delays for this PR. (However, this can definitely wait until after SFN...)