python-neo icon indicating copy to clipboard operation
python-neo copied to clipboard

NestIO compatible with Nest 3.2

Open morales-gregorio opened this issue 2 years ago • 14 comments

Hi all,

@rgutzen @jasperalbers and I were looking into the NestIO today and we came up with this edits that make the current NestIO compatible with both Nest 2.x and 3.x versions.

We incorporated all the changes suggested by @essink into this PR.

Essentially we force the IO to ignore all the headers (which raises a warning that it did so) as suggested by @rgutzen, which I also wrote into PR #739

Additionally, we found that the current NestIO decides whether to read neo.AnalogSignals or neo.SpikeTrain based on the file name alone. This is problematic since Nest 3.2 (maybe others too) saves spikes to the .dat format by default, which the IO expected to be analogsignals. I fixed this issue by adding a new keyword target_object which then determines what kind of object will be read. This does not allow to read analogsignals and spikes into the same segment (two different calls are needed), but this seems a rather good practice in my opinion, since spiketrains and analogsignales are rarely combined into the same segment.

I tested the code to load a .dat file generated with Nest 3.2 and it works, I can add a proper test somewhere if you want.

Best, Aitor

morales-gregorio avatar Mar 09 '22 12:03 morales-gregorio

Looking into the capability to deal with lists it only seems to complicate things and it is not capable of dealing with multiple files of the same type, which is something I would expect. I would suggest getting rid of the capability to deal with lists altogether and stick to single files as most other IOs do. This then allows to remove code further and if someone really wants to have the analogsignals and spikes together it is just a matter of running this two times.

morales-gregorio avatar Mar 09 '22 14:03 morales-gregorio

Hi all, thanks for taking care of this. We created a NEST-elephant tutorial and we ran into this issue with the nestio and NEST 3.3.

Since this was tested with NEST 3.2, maybe add this info to the docstring(?):

https://github.com/NeuralEnsemble/python-neo/blob/50abfbce832edb1c0f8364ac34c8a6442e1a3edf/neo/io/nestio.py#L4

Moritz-Alexander-Kern avatar Aug 12 '22 10:08 Moritz-Alexander-Kern

Hi @morales-gregorio Do you have a small test file for this format that could be added to gin?

JuliaSprenger avatar Feb 16 '23 14:02 JuliaSprenger

Just rebased the branch and some chaos has followed. My local history looks as one would expect, with this branch now splitting from the current master. I am not sure why Github interprets that all the previous commits are part of this PR :/

morales-gregorio avatar Feb 17 '23 14:02 morales-gregorio

@JuliaSprenger I created a small NEST (3.3) simulation where I record both spikes and analogsignals. The results and simulation script are here: test_NEST_simulation_and_results.zip.

The current intended use is therefore something like this:

import neo
import quantities as pq
import matplotlib.pyplot as plt
# Load spikes
io = neo.NestIO(['sr1-142-0.dat', 'sr2-154-0.dat'], target_object='SpikeTrain')
bl = io.read_block([], t_start=0*pq.s, t_stop=20*pq.s)

# Plot spikes
plt.figure(figsize=(4, 3))
for i, st in enumerate(bl.segments[0].spiketrains):
    if st.annotations['source_file'] == 'sr1-142-0.dat':
        color = 'b'
    elif st.annotations['source_file'] == 'sr2-154-0.dat':
        color = 'r'
    plt.scatter(st.times, [i]*len(st.times), marker='|', c=color)
plt.xlabel('Time (ms)')
plt.ylabel('Neurons')
plt.savefig('nest_spikes.png', dpi=300)

# Load anasigs
io = neo.NestIO(['vm1-141-0.dat', 'vm2-153-0.dat'], target_object='AnalogSignal')
bl = io.read_block([], t_start=0*pq.s, t_stop=20*pq.s)

# Plot anasigs
plt.figure(figsize=(4, 3))
for i, ana in enumerate(bl.segments[0].analogsignals):
    if ana.annotations['source_file'] == 'vm1-141-0.dat':
        color = 'b'
    elif ana.annotations['source_file'] == 'vm2-153-0.dat':
        color = 'r'

    plt.plot(ana.times, 8*i + ana.magnitude, c=color, alpha=0.5)

plt.xlabel('Time (ms)')
plt.ylabel('Neurons')
plt.savefig('nest_anasigs.png', dpi=300)

Which creates these two plots in this case: nest_anasigs nest_spikes

morales-gregorio avatar Feb 17 '23 17:02 morales-gregorio

@JuliaSprenger I just fixed the history, are there any open issues with the current version of the NestIO?

morales-gregorio avatar Aug 03 '23 10:08 morales-gregorio

@JuliaSprenger I just fixed the history, are there any open issues with the current version of the NestIO?

or @apdavison

Moritz-Alexander-Kern avatar Aug 03 '23 11:08 Moritz-Alexander-Kern

Hi @heplesser , thanks for the comments and the detailed description of NEST outputs. I think this information is super helpful in driving this IO forwards. Also, thanks for your insights in problems that may arise at the "full scale" level.

mdenker avatar Feb 01 '24 13:02 mdenker

Hello everyone!

After discussion with @mdenker and @apdavison yesterday, I provide here a little script that creates a network consisting of

  • aeif_cond_alpha neurons (to show multimeter recording)
  • iaf_psc_alpha_ps neurons (to show precise time recording)
  • weight_recorder (to show weight recording)

Running the script will generate a number of spike recorder, voltmeter, multimeter and weight recorder output files with different settings. Note in particular the effect of "time_in_steps": If false (default), time is recorded in ms under heading "time_ms", if true time will be recorded in columns "time_step" and "time_offset", where the actual spike time is given by time_step * resolution - time_offset. Unfortunately, the resolution is not stored in the file, so I suggest to not support these files for now.

Do not rely on ordering of columns or rows. In parallel simulations, there will be one file per virtual process, enumerated by the last part of the file name.

NEST 3 files can be safely identified by the

# NEST version: 3.7.0-rc1
# RecordingBackendASCII version: 2

header. NEST 2 files do not have these headers. I would suggest to use the header for identification, also to be prepared for future versions of the ASCII backend.

nest.ResetKernel()
nest.overwrite_files = True

g = nest.Create('poisson_generator', params={'rate': 100000})
n = nest.Create('aeif_cond_alpha', 10)
p = nest.Create('iaf_psc_alpha_ps', 5)

rnt = nest.Create('spike_recorder', params={'record_to': 'ascii',
                                            'label': 'aeif_spikes_times'})
rns = nest.Create('spike_recorder', params={'record_to': 'ascii',
                                           'time_in_steps': True,
                                            'label': 'aeif_spikes_steps'})

rpt = nest.Create('spike_recorder', params={'record_to': 'ascii',
                                            'label': 'precise_spikes_times'})
rps = nest.Create('spike_recorder', params={'record_to': 'ascii',
                                           'time_in_steps': True,
                                            'label': 'precise_spikes_steps'})

vt = nest.Create('voltmeter', params={'record_to': 'ascii',
                                      'label': 'voltmeter_times'})
vs = nest.Create('voltmeter', params={'time_in_steps': True,
                                      'label': 'voltmeter_steps'})

mm1t = nest.Create('multimeter', params={'record_from': ('g_ex', 'g_in', 'V_m', 'w'),
                                         'record_to': 'ascii',
                                            'label': 'multimeter_1ms'})
mmrt = nest.Create('multimeter', params={'record_from': ('g_ex', 'g_in', 'V_m', 'w'),
                                         'record_to': 'ascii',
                                         'interval': nest.resolution,
                                            'label': 'multimeter_res'})

w = nest.Create('weight_recorder', params={'record_to': 'ascii',
                                            'label': 'weight_rec'})

nest.SetDefaults('stdp_synapse', {'weight_recorder': w})

nest.Connect(g, n)
nest.Connect(n, n, syn_spec={'synapse_model': 'stdp_synapse', 'weight': 100.})
nest.Connect(n, p, syn_spec={'weight': 1000.})
nest.Connect(n, rnt + rns)
nest.Connect(p, rpt + rps)
nest.Connect(vt + vs + mm1t + mmrt, n)

nest.Simulate(100)

heplesser avatar Apr 05 '24 10:04 heplesser

Hi @heplesser ! Thanks for pushing this forward and coming up with good test examples.

I think we should not introduce a NEST dependency to the neo tests, I would instead suggest to create several (small) files for multiple NEST versions. The files should have descriptive names like "nest3.7_voltmeter_ascii.dat". What do you think?

morales-gregorio avatar Apr 05 '24 10:04 morales-gregorio

Hi @heplesser ! Thanks for pushing this forward and coming up with good test examples.

I think we should not introduce a NEST dependency to the neo tests, I would instead suggest to create several (small) files for multiple NEST versions. The files should have descriptive names like "nest3.7_voltmeter_ascii.dat". What do you think?

@morales-gregorio I agree, otherwise we could easily end up with circular dependencies and all-to-heavy test setups. Frozen example files for testing should work nicely. It should suffice to base the versioning on the ASCII Recording Backend version, which changes more slowly than the NEST version.

heplesser avatar Apr 05 '24 11:04 heplesser

Hi @heplesser ! Thanks for pushing this forward and coming up with good test examples. I think we should not introduce a NEST dependency to the neo tests, I would instead suggest to create several (small) files for multiple NEST versions. The files should have descriptive names like "nest3.7_voltmeter_ascii.dat". What do you think?

@morales-gregorio I agree, otherwise we could easily end up with circular dependencies and all-to-heavy test setups. Frozen example files for testing should work nicely. It should suffice to base the versioning on the ASCII Recording Backend version, which changes more slowly than the NEST version.

@heplesser could you create those files? 😇

morales-gregorio avatar Apr 05 '24 11:04 morales-gregorio

Hey @morales-gregorio ,

the current test files for testing NEST-IO are located on NeuralEnsemble / ephy_testing_data in nest here: https://gin.g-node.org/NeuralEnsemble/ephy_testing_data/src/master/nest

The contents are:

0gid-1time-1256-0.gdf
0gid-1time-2gex-1262-0.dat
0gid-1time-2gex-3Vm-1261-0.dat
0gid-1time-2Vm-1259-0.dat
0gid-1time-2Vm-3gex-4gin-1260-0.dat
0gid-1time-2Vm-3Iex-4Iin-1264-0.dat
0gid-1time_in_steps-1258-0.gdf
0gid-1time_in_steps-2Vm-1263-0.dat
0time-1255-0.gdf
0time_in_steps-1257-0.gdf
brunel-delta-nest_mod.py  # <--- script that generated the files
N1-0gid-1time-2Vm-1265-0.dat
N1-0time-1Vm-1266-0.dat
N1-0Vm-1267-0.dat

I have copied and formatted the above script into generate_data.py:

"""
Script to generate test data for the nestIO.

This script creates a network simulation using the NEST simulator (https://www.nest-simulator.org/).
It generates output files for spike recorders, voltmeters, multimeters and weight recorders with
different settings. Note the effect of "time_in_steps": If false (default), time is recorded in ms
under heading "time_ms"; if true, time will be recorded in columns "time_step" and "time_offset",
where the actual spike time is given by time_step * resolution - time_offset. The resolution is not
stored in the file. In parallel simulations, there will be one file per virtual process, enumerated
by the last part of the file name.

Usage:
    python generate_data.py


"""

import nest


def main():
    """
    Executes the main logic of the script.

    This function sets up a simulation using the NEST simulator (https://www.nest-simulator.org/)
    for simulating neural networks. It creates a network consisting of aeif_cond_alpha neurons
    (to show multimeter recording), iaf_psc_alpha_ps neurons (to show precise time recording),
    and weight_recorder (to show weight recording). The simulation is run for a specified duration.

    Note: Do not rely on ordering of columns or rows. In parallel simulations, there will be
    one file per virtual process, enumerated by the last part of the file name.

    Parameters:
        None

    Returns:
        None

    """
    nest.ResetKernel()
    nest.overwrite_files = True

    g = nest.Create("poisson_generator", params={"rate": 100000})
    n = nest.Create("aeif_cond_alpha", 10)
    p = nest.Create("iaf_psc_alpha_ps", 5)

    rnt = nest.Create(
        "spike_recorder",
        params={"record_to": "ascii", "label": "aeif_spikes_times"},
    )
    rns = nest.Create(
        "spike_recorder",
        params={
            "record_to": "ascii",
            "time_in_steps": True,
            "label": "aeif_spikes_steps",
        },
    )

    rpt = nest.Create(
        "spike_recorder",
        params={"record_to": "ascii", "label": "precise_spikes_times"},
    )
    rps = nest.Create(
        "spike_recorder",
        params={
            "record_to": "ascii",
            "time_in_steps": True,
            "label": "precise_spikes_steps",
        },
    )

    vt = nest.Create("voltmeter", params={"record_to": "ascii", "label": "voltmeter_times"})
    vs = nest.Create("voltmeter", params={"time_in_steps": True, "label": "voltmeter_steps"})

    mm1t = nest.Create(
        "multimeter",
        params={
            "record_from": ("g_ex", "g_in", "V_m", "w"),
            "record_to": "ascii",
            "label": "multimeter_1ms",
        },
    )
    mmrt = nest.Create(
        "multimeter",
        params={
            "record_from": ("g_ex", "g_in", "V_m", "w"),
            "record_to": "ascii",
            "interval": nest.resolution,
            "label": "multimeter_res",
        },
    )

    w = nest.Create("weight_recorder", params={"record_to": "ascii", "label": "weight_rec"})

    nest.SetDefaults("stdp_synapse", {"weight_recorder": w})

    nest.Connect(g, n)
    nest.Connect(n, n, syn_spec={"synapse_model": "stdp_synapse", "weight": 100.0})
    nest.Connect(n, p, syn_spec={"weight": 1000.0})
    nest.Connect(n, rnt + rns)
    nest.Connect(p, rpt + rps)
    nest.Connect(vt + vs + mm1t + mmrt, n)

    nest.Simulate(100)


if __name__ == "__main__":
    main()

which generates: (.zip attached)

aeif_spikes_steps-18-0.dat
aeif_spikes_times-17-0.dat
generate_data.py
multimeter_1ms-23-0.dat
multimeter_res-24-0.dat
precise_spikes_steps-20-0.dat
precise_spikes_times-19-0.dat
voltmeter_times-21-0.dat
weight_rec-25-0.dat

nestIO_test_data.zip

Refer to this readme here: https://gin.g-node.org/NeuralEnsemble/ephy_testing_data to add the files.

Moritz-Alexander-Kern avatar Apr 05 '24 11:04 Moritz-Alexander-Kern

There are 8 failing tests. I've looked into the first of these, and it seems to be that the case of id_column=None is not being correctly handled.

apdavison avatar May 06 '24 14:05 apdavison