Assigning GIDs to netcon without cell object

[pyadav]

Hello, I am new to Neuron and I have been working on getting and existing model to work on GPU. I have been following the ringtest repository to understand how the GPU implementation works and how I should initialize the model to work with parallel context.
The model consists of k cells (motor neurons), each with one soma and four dendrites. I also want to play vectors into each of the dendrites of the cell. However, the cell I create is composed of h.SectionList() and not a neuron cell() object. As a result I am not able to use the cell().connect2target() function and I create NetCon objects separately to link with the GIDs I assign to each of the sections.

When I run simulation by using parallelcontext.psolve() I encounter the following error:

Code: Select all

~/nrn/src/coreneuron/io/phase1.cpp:108: Assertion 'neg_gid2out[nt.id].find(gid) == neg_gid2out[nt.id].end()' failed.

This made me suspect that the GIDs I have assigned to the various sections in my model are not behaving as I expect them to. When I use the gid2cell() or gid2obj() functions they return "None" for my model. For the ringtest repo these functions return the cell object from which the cell was initialized leading me to think there is an issue with how I am initializing the NetCon() or how I am linking the GIDs. The GIDs are assigned to the device and gid_exists() function does return value 2, indicating the GID was assigned and made to not communicate the spikes across machines.
These errors do not show up when I use a while loop to run simulation with h.fadvance() and integrate the simulation step by step. However, when I use psolve() then I get the aforementioned error.
Here is a minimal script to reproduce the error I am encountering and model initialization

Code: Select all

from neuron import h
import settings
from neuron import coreneuron
from neuron.units import mV
import numpy as np

pc = h.ParallelContext()
h.CVode().cache_efficient(1)
h.nrnmpi_init()

def spike_record(): # from ringtest repo
    global tvec, idvec
    tvec = h.Vector(1000000)
    idvec = h.Vector(1000000)
    pc.spike_record(-1, tvec, idvec)

def spikeout(): # from ringtest repo
    global tvec, idvec
    rank = int(pc.id())
    nhost = int(pc.nhost())
    pc.barrier()
    return tvec, idvec

class MNPool(object): # based on ringtest repo
    counter = 0
    def __init__(self, ncell, gidstart = 1, nsoma = 1, ndend = 4):
        self.gids = []
        self.delay = 1
        self.ncell = int(ncell)
        self.gidstart = gidstart
        self.mkcells(self.ncell, nsoma = 1, ndend = 4)
        MNPool.counter += 1

    def mkcells(self, ncell, nsoma = 1, ndend = 4): # creates ncell motor neuron with 4 dendrites
        self.all_MNs = []
        ctr = 0
        for i in range(self.gidstart, (ncell*(nsoma+ndend)) + self.gidstart,nsoma+ndend): # creates single MN
            
            self.gids.append([i+x for x in range(nsoma+ndend)])
            soma = h.Section(name='soma')  # Define soma section
            d1 = h.Section(name='d1')  # Define dend section
            d2 = h.Section(name='d2')  # Define dend section
            d3 = h.Section(name='d3')  # Define dend section
            d4 = h.Section(name='d4')  # Define dend section

            dend = h.SectionList()
            dend.append(sec=d1)
            dend.append(sec=d2)
            dend.append(sec=d3)
            dend.append(sec=d4)

            # Connect dendrites to soma
            d1.connect(soma(1), 0)
            d2.connect(soma(1), 0)
            d3.connect(soma(0), 0)
            d4.connect(soma(0), 0)

            MN_parts = {'soma':[soma],
                        'dend':[d1,d2,d3,d4]}
            self.all_MNs.append(MN_parts)
            for part in self.all_MNs[ctr].keys():
                if part == 'dend':
                    for section in self.all_MNs[ctr][part]:
                        section.insert('Gfluctdv') # inserts the mechanism on all 4 dendrites
            
            # This is based on ringtest repo
            cell_ctr = 0
            for part in MN_parts.keys(): # loop over both soma and dend dict
                for section in MN_parts[part]: # from each dict loop over the sections
                    settings.pc.set_gid2node(i+cell_ctr, settings.rank) # Assign GID for that section
                    nc = h.NetCon(section(0.5)._ref_v, None, sec = section) # Since the MN is not a cell objcet init a NetCon
                    settings.pc.cell(i+cell_ctr, nc, 0) # link netcon and GID
                    print("MN: ",ctr, # print statement to determine GIDs assigned to each MN and section
                    ", section: ", section, 
                    ", GID for section: ", i+cell_ctr, 
                    ", pc.gid_exists(): ", settings.pc.gid_exists(i+cell_ctr), # expected to return 2 (assigned GID but not communicating spikes)
                    ", pc.gid2obj(): ", settings.pc.gid2obj(i+cell_ctr) # expected to return the section or something similar but returns None
                    )
                    cell_ctr+=1
            ctr+=1
            
def main():
    num_MN = 2
    settings.init(num_MN)

    MN_model = MNPool(num_MN,nsoma=1, ndend= 4) # init model
    spike_record() # init spike recording

    playvect_dict = {
        "rn_vec" : [],
    }
    sim_duration = 30000
    excitatory_time = np.arange(0, sim_duration, 1)
    target_time = h.Vector(excitatory_time)

    noise_vect_w = np.random.normal(0,1,(len(MN_model.all_MNs),len(target_time))) # noise vector for the mechanism

    for ctr, MN in enumerate(MN_model.all_MNs):
        for num_dend, sec in enumerate(MN["dend"]):
            for seg in sec:
                rn_vec = h.Vector(noise_vect_w[ctr,:])
                rn_vec.play(seg._ref_rand_Gfluctdv,target_time) # play noise vector into each of the MN mechanisms
                playvect_dict['rn_vec'].append(rn_vec)

    coreneuron.enable = True
    coreneuron.gpu = True

    pc.set_maxstep(10)

    h.dt = 0.025  # unit: ms
    h.finitialize(-70 * mV)
    h.fcurrent()

    
    pc.psolve(sim_duration)

    """
    pc.psolve() raises the following error
    int nrnthread_dat2_vecplay(int, std::vector<int, std::allocator<int>> &): Assertion `0' failed.
    """


    while h.t < float(sim_duration):
        h.fadvance()
        if round(h.t, 2) % 10000.00 == 0.00:
            print(f'simulation time: {h.t:.1f}', 'ms',)
    
    """
    h.fadvance() works without error
    """
    tvec, idvec = spikeout()

if __name__ == "__main__":
    main()
    

Following is the output of the print statement in model initialization.

Code: Select all

MN:  0 , section:  soma , GID for section:  1 , pc.gid_exists():  2 , pc.gid2obj():  None
MN:  0 , section:  d1 , GID for section:  2 , pc.gid_exists():  2 , pc.gid2obj():  None
MN:  0 , section:  d2 , GID for section:  3 , pc.gid_exists():  2 , pc.gid2obj():  None
MN:  0 , section:  d3 , GID for section:  4 , pc.gid_exists():  2 , pc.gid2obj():  None
MN:  0 , section:  d4 , GID for section:  5 , pc.gid_exists():  2 , pc.gid2obj():  None
MN:  1 , section:  soma , GID for section:  6 , pc.gid_exists():  2 , pc.gid2obj():  None
MN:  1 , section:  d1 , GID for section:  7 , pc.gid_exists():  2 , pc.gid2obj():  None
MN:  1 , section:  d2 , GID for section:  8 , pc.gid_exists():  2 , pc.gid2obj():  None
MN:  1 , section:  d3 , GID for section:  9 , pc.gid_exists():  2 , pc.gid2obj():  None
MN:  1 , section:  d4 , GID for section:  10 , pc.gid_exists():  2 , pc.gid2obj():  None

I am trying to determine what could be the source for the error with the psolve() and the gid2obj() function returning "None". Please let me know if I can provide any additional information to help understand this issue.
Thank you!

[ted]

the cell I create is composed of h.SectionList() and not a neuron cell() object

First, a bit of terminological clarity. NEURON has no such thing as a "neuron cell() object". It's up to the user to decide what types of cells are to be represented in a network model, then for each of those types to create a corresponding class whose instances will have the user-defined anatomical and biophysical properties, and finally to create those instances.

SectionLists have their uses, but implementing cell instances for a network model is not one of them. Why? Many reasons, but top of my list are: SectionList objects don't encapsulate code, don't prevent collisions in name space, don't facilitate model introspection, don't allow "random access" of individual members (in fact, don't have most of the features of NEURON's own List class, let alone those of a Python list), and can't be used to write scalable code (which makes development, debugging, and performance evaluation very difficult). For example, if foo is the name of an instance of a cell class in a network, and that cell class has a soma, then foo's soma will be called foo.soma. But if foo is the name of a model cell implemented as SectionList, its soma will have a name of the form soma[j] where the value of j depends on how many other somas were already created before foo was created, and statements that involve terms like "foo.soma" or "foo.soma[j]" will generate an error that stops program execution during model setup.

Even if you have already gone to the effort of using SectionLists to create a serial implementation of your network model, I strongly advise you to abandon that and re-implement using cell classes. You'll save yourself an enormous amount of time and effort, and might even create something that others might want to reuse (the highest form of praise, as long as they cite your work . . . ). For me to advise otherwise would be a disservice to you.

[pyadav]

Thank you for the clarification!
I have updated my model to comprise of h.Section() instead of h.SectionList() to implement the cell.
I still wonder what could be the cause for the function gid2obj() or gid2cell() returning "None". Also, if the error raised by psolve() is a consequence of incorrect initialization of the GIDs in the model.

As far as I can understand the following lines of code should have:
1. Set the GID to a specific node
2. Initialized a NetCon() for the POINT_PROCESS (voltage) of that section
3. Linked the GID and NetCon()

Code: Select all

            for part in MN_parts.keys(): # loop over both soma and dend dict
                for section in MN_parts[part]: # from each dict loop over the sections
                    settings.pc.set_gid2node(i+cell_ctr, settings.rank) # Assign GID for that section
                    nc = h.NetCon(section(0.5)._ref_v, None, sec = section) # Since the MN is not a cell objcet init a NetCon
                    settings.pc.cell(i+cell_ctr, nc, 0) # link netcon and GID

My expectation was to receive the section as the output of gid2obj()/gid2cell() functions.
Is my understanding correct?

[ted]

Documentation says "The cell or artificial cell object is returned that is associated with the global id." Says "cell or artificial cell object", not "(top level) section" (that is, it doesn't say "a section that is not encapsulated in an object"). Is this the source of the problem, or am I overinterpreting the documentation? Try this simple test and see what happens, and we'll both learn something.

Phase 1.
Create an instance of a cell class.
Attach a NetCon to one of its sections and use a ParallelContext.cell(gid, NetCon) statement to associate a gid with that NetCon's spike source.
Then verify that gid2obj and gid2cell work.

Phase 2.
Restart NEURON.
Create an instance of a section.
Attach a NetCon to it and use a ParallelContext.cell(gid, NetCon) statement to associate a gid with that NetCon's spike source.
Then see if gid2obj and gid2cell work, or if they instead generate error messages.

Comment: documentation is generally terse, and usually does not include all caveats (hey, it's software documentation, not a EULA that denies any imaginable rights to the user and indemnifies the vendor from all possible damages). Sometimes disambiguation requires running a test to see if one's interpretation is or isn't correct.

[pyadav]

Thank you for the suggested test. After running the test and some variations of initializing a NetCon I observed the following.
Phase 1: After initializing a Cell class and attaching a NetCon to its soma and associating with a GID to the spike source, gid2obj and gid2cell still return "None".
I initialized the cell using the following object

Code: Select all

class Cell(object):
    def __init__(self):
        self.soma = h.Section(name='soma')
        self.d1 = h.Section(name='d1')  # Define dend section
        self.d1.connect(self.soma(1), 0)
        self.d1.insert('Gfluctdv')
        gid = 10
        settings.pc.set_gid2node(gid, settings.rank) # Assign GID for that section
        nc = h.NetCon(self.soma(0.5)._ref_v, None, sec = self.soma) # Since the MN is not a cell objcet init a NetCon
        settings.pc.cell(gid, nc, 0) # link netcon and GID

Phase 2: After creating a section (without Cell object) and following the steps, gid2obj and gid2cell still return "None".
I used the following lines to initialize the section

Code: Select all

    gid = 10
    soma = h.Section(name='soma')
    d1 = h.Section(name='d1')  # Define dend section
    d1.connect(soma(1), 0)
    d1.insert('Gfluctdv')
    settings.pc.set_gid2node(gid, settings.rank) # Assign GID for that section
    nc = h.NetCon(soma(0.5), None, sec = soma) # Since the MN is not a cell objcet init a NetCon
    settings.pc.cell(gid, nc, 0) # link netcon and GID
    

Please correct me if I am wrong, but during the initialization of a NetCon, the first argument has to be a point process. If I pass a h.Section or the Cell object the following error is raised "if arg 1 is an object it must be a point process or NULLObject".
This means that the NetCon counts spikes at the POINT_PROCESS (x._ref_v) and we link the GID to this point process. When gid2cell/gid2obj are used for a specific GID, then these functions return the POINT_PROCESS and not the "cell or artificial cell object" which is returned as "None"

[ted]

I guess I wasn't sufficiently explicit. By "simple test" I meant dirt simple--no dependence on anything (other than e.g. MPI and your operating system) that that isn't built into NEURON. That means no ion channels other than pas and HH. And no reference to 3rd party modules e.g. "settings". And don't include creation of the NetCon, or gid assignment, in the cell class--do those tasks at the top level of the interpreter (outside of any class instance).

Why? Because now the task has become
1. to discover whether ParallelContext and NetCon work in your hardware/software environment,
and, assuming that they do work,
2. to discover what is not working, and why, when you add your own elaborations

And one way to figure this out is to start with the simplest thing that works, then add complications one at a time and run tests.