[ENH] Added read/write method for network object #651

raj1701 · 2023-05-24T20:41:06Z

GSoC 2023 Project - Develop IO routines for HNN-core outputs

The idea of the project is to add read/write capabilities for all HNN objects such as Network, Cell, Section, Dipole etc. in the HDF5 format. During the summer, I wrote IO functions for all HNN objects, the associated tests as well as refactored some existing code along the journey.

The following PR contains read/write functions and the associated tests for the Network object. The network class encapsulates objects such as Cell, Section, Cell response and Extracellular arrays. These objects are stored along the network itself. Each and every aspect of the functions is tested and documented. The IO functions are capable of handling all network models in HNN core as well any modified network instances. Capability to saving and loading up simulated networks is also present.

Some other PRs related to the project are -

[MRG] Find end pts of section and a numpy version of h.distance() #661 - This PR contains python version of NEURON functions such as define_shape() and distance(). This is done to reduce the dependency on NEURON for saving a Network. NEURON objects are only needed while simulating a Network now.
[MRG] Added read/write methods and tests for dipole #648 - This PR contains read/write functions and associated tests for the Dipole object.
[MRG] Added read/write methods and tests for cell response #644 - This PR contains read/write functions and associated tests for the Cell response object.
[MRG] Reverted hdf5 functions and added to_dict method #654 - This PR reverted the IO functionalities for the Cell response object as they are included in the IO functions for the Network object in [ENH] Added read/write method for network object #651.

Plans for the project after GSoC - Add load/save buttons in HNN GUI for the Network and Dipole object.

jasmainak · 2023-05-24T22:38:29Z

hnn_core/tests/test_network.py

+    # Check cell_types
+    # assert net_jones.cell_types == net_jones_read.cell_types
+    print(net_jones.cell_types)
+    print(net_jones_read.cell_types)


For test-driven development, you can use pytest with --pdb to directly drop into the interpreter when a test fails. Then you can use pdb commands like p net_jones.cell_types etc. Just letting you know the different options on debugging efficiently

jasmainak · 2023-05-24T22:40:10Z

hnn_core/network.py

+    def __eq__(self, other):
+        if not isinstance(other, Network):
+            return NotImplemented
+        # Check for all attributes (Discuss)


I think the right way to implement this is to first implement __eq__ for Section and Cell object and call them from here. Python builtins already have it implemented so presumably net.external_drives etc should be straightforward to check for equality

jasmainak · 2023-05-24T22:44:47Z

For partial function, I think you should expand it out. Iterate over the segments and get the gbar values. Like here:

hnn-core/hnn_core/tests/test_network.py

Lines 57 to 58 in 07374a3

    
           ca_gbar = [seg.__getattribute__('ca').gbar for 
        
                      seg in list(section.allseg())[1:-1]]

and store them. It's less compact but I think the only way to make it work

codecov-commenter · 2023-05-25T20:22:26Z

Codecov Report

Attention: 34 lines in your changes are missing coverage. Please review.

Comparison is base (2f5b15b) 91.62% compared to head (52cd670) 92.24%.
Report is 63 commits behind head on master.

❗ Current head 52cd670 differs from pull request most recent head cda152e. Consider uploading reports for the commit cda152e to get more accurate results

❗ Your organization needs to install the Codecov GitHub app to enable full functionality.

Additional details and impacted files

@@            Coverage Diff             @@
##           master     #651      +/-   ##
==========================================
+ Coverage   91.62%   92.24%   +0.61%     
==========================================
  Files          22       27       +5     
  Lines        4488     4897     +409     
==========================================
+ Hits         4112     4517     +405     
- Misses        376      380       +4

Files	Coverage Δ
hnn_core/__init__.py	`100.00% <100.00%> (ø)`
hnn_core/cell_response.py	`86.20% <ø> (+7.88%)`	⬆️
hnn_core/docs.py	`100.00% <100.00%> (ø)`
hnn_core/optimization/__init__.py	`100.00% <100.00%> (ø)`
hnn_core/optimization/optimize_evoked.py	`92.99% <100.00%> (ø)`
hnn_core/extracellular.py	`88.93% <93.33%> (+0.72%)`	⬆️
hnn_core/optimization/objective_functions.py	`92.30% <92.30%> (ø)`
hnn_core/network.py	`93.20% <90.90%> (-0.18%)`	⬇️
hnn_core/io.py	`97.14% <97.14%> (ø)`
hnn_core/cell.py	`96.82% <90.62%> (+0.28%)`	⬆️
... and 1 more

... and 3 files with indirect coverage changes

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

raj1701 · 2023-05-25T21:10:26Z

I tried rebasing but it didn't work. I followed the following steps

git fetch origin master
git checkout <feature_branch>
git rebase origin master
git push --force origin <feature_branch>

I encountered no errors but I couldn't see the changes present only in master. Where did I go wrong?

jasmainak · 2023-05-26T17:46:27Z

@raj1701 I think you just need to do:

$ git rebase master

after checking out the feature_branch if you already fetched master

raj1701 · 2023-06-02T19:37:41Z

Hey @jasmainak @ntolley, I found a bug while saving cell_response after simulation

Bug Description
The bug occurs because of the structure of vsec and isec. vsec contains a list of dict. The keys of each dict are of type int. For hdf5 keys of dict must be of type str

Solution
Convert the dict keys to the required types before writing and after reading.

I have implemented this solution in network object where I have used _vsec and _isec in cell_response. A similar solution will be required for read/write in cell_response.py (A separate PR mostly).
Since isec and vsec always remained empty when I was working explicitly on the cell_response object, the bug never occured. After running a simulation, I got vsec values and thus the bug.
LMK your thoughts on the solution and whether I should open a separate PR for cell_response object

jasmainak · 2023-06-03T01:21:05Z

Yes, you should open a new PR and rebase this branch when that PR is merged

hnn_core/network.py

jasmainak · 2023-06-03T01:25:48Z

hnn_core/cell.py

+            if not (self.sections[key] == other.sections[key]):
+                return False
+
+        return True


you should add a test in test_cell.py for this ... simple way to test is:

cell2 = cell.copy() assert cell2 == cell

jasmainak · 2023-06-03T01:26:11Z

hnn_core/network.py

@@ -412,6 +415,43 @@ def __repr__(self):
                 len(self.pos_dict['L5_basket'])))
        return '<%s | %s>' % (class_name, s)

+    def __eq__(self, other):


add a test for this method

In net.copy() method, the simulation data related to external drives and rec_arrays eg- events is deleted.
Is a similar behavior intended for read\write functionality?

jasmainak · 2023-06-08T17:12:34Z

hnn_core/network.py

@@ -1335,6 +1407,280 @@ def plot_cells(self, ax=None, show=True):
        """
        return plot_cells(net=self, ax=ax, show=show)

+    def write(self, fname):


docstring missing, also add overwrite parameter

raj1701 · 2023-06-21T08:35:16Z

I think docstrings are a better way to document as h5tree throws some random errors

jasmainak · 2023-06-21T14:53:10Z

Let's do docstrings for now

jasmainak · 2023-06-21T14:57:33Z

hnn_core/network.py

@@ -1340,6 +1415,310 @@ def plot_cells(self, ax=None, show=True):
        """
        return plot_cells(net=self, ax=ax, show=show)

+    def write(self, fname, overwrite=True, save_unsimulated=False):


This is a fair amount of code, I propose moving this to a new file ... io.py maybe? You'd basically house all the write_xxx functions and read_xxx functions there. Then here you'd just call the write_network function

There would then be a new test file test_io.py

Sure will do this

jasmainak · 2023-06-21T14:58:08Z

hnn_core/network.py

+        # Tasks
+        # 1) Remove params dict (Partial done)
+        # 2) Add overwrite warnings (Done)
+        # 3) Add boolean to read/write un-simulated network (Done)
+        # 4) Add boolean to read/write simulated network (Done)
+        # 5) Think about documentation
+        # 6) Add tests
+        # 7) Store the type of object saved and issue warning (Done)


suggest moving this to the PR description ... github has a nice checkbox feature you can use

Ooh its a pretty cool feature. I just saw it.

jasmainak · 2023-06-21T15:04:09Z

hnn_core/cell.py

@@ -443,6 +529,9 @@ def build(self, sec_name_apical=None):
            with this section. The section should belong to the apical dendrite
            of a pyramidal neuron.
        """
+        # Confirm correctness here
+        from .network_builder import load_custom_mechanisms
+        load_custom_mechanisms()


does it need to be in the other place now that you have the function call here?

I'll try removing it from the other place

rythorpe · 2023-06-22T03:53:26Z

hnn_core/cell.py

                        h.pop_section()
+                        p_mech[attr] = [seg_xs, seg_vals]


I'm assuming the purpose of this is to overwrite the callable stored in p_mech? If so, this will require calling cell.build() any time you need to write write cells attached to a Network instance, which will in turn imply that we're going to instantiate NEURON objects from within Network. I'm just thinking out loud here, but will that cause any unintended issues with the NEURON network? In the past, we've maintained the convention that NEURON objects are only instantiated from NetworkBuilder in order to silo how, when, and where such objects are created.

Yes, calling cell.build() outside the NetworkBuilder does throw an error later when Network is simulated. For resolving it, I first made a copy of the cell and called build() on the copy instead of the original.

Do we know if the NEURON objects are getting properly destroyed and de-referenced after the cell copy is destroyed by the garbage collector? The only reason I'm bringing this up is because great care needs to be taken when handling NEURON network objects which are generally global and have in the past caused caused segmentation faults that are hard to diagnose and/or make tests for.

Is there a way around this that allows us to maintain our convention of only building NEURON objects from within NetworkBuilder? CC'ing @jasmainak and @ntolley for their input.

yeah I agree we need to be careful here. See:

hnn-core/hnn_core/network_builder.py

Line 658 in c81e5b8

def _clear_last_network_objects(self):

The trade-off is between doing our computation ... which can be error-prone if we don't get the positions of the segments exactly right. Neuron does some mumbo-jumbo here:

hnn-core/hnn_core/cell.py

Lines 431 to 434 in c81e5b8

# be explicit about letting sec.L dominate over the 3d points used by

# h.pt3dadd(); see

# https://nrn.readthedocs.io/en/latest/python/modelspec/programmatic/topology/geometry.html?highlight=pt3dadd#pt3dadd # noqa

h.define_shape()

vs

potentially slowing up or causing segmentation faults in all our simulations that have run the save method

The fact that the segment coordinates prior to NEURON instantiation aren't their final values (and thus need to be re-set using length modifiers) is annoying. Maybe we should rip the bandaid off and fix this, which I think would also nullify the need to build the NEURON cells in order to write the network.

I'm all for it. Do we need the length for any purpose? I seem to remember that either you or @ntolley were experimenting with changing the length of the neurons ... or do you propose changing the coordinates of the neuron to match the length, that may lead to some tiny numerical inaccuracies

I'm proposing the latter, which I think would allow us to get rid of the length parameter. While this may introduce small numerical differences compared to what we currently have, it's a much more sustainable solution.

what would be concrete steps for @raj1701 ?

Raise warning if length is defined in param file? It should say that it has no effect. Same with network attribute : they should not be modifiable

Update the coordinates here to match the default length parameters:

hnn-core/hnn_core/cells_default.py

Lines 69 to 77 in c81e5b8

end_pts = {

'soma': [[-50, 0, 765], [-50, 0, 778]],

'apical_trunk': [[-50, 0, 778], [-50, 0, 813]],

'apical_oblique': [[-50, 0, 813], [-250, 0, 813]],

'apical_1': [[-50, 0, 813], [-50, 0, 993]],

'apical_tuft': [[-50, 0, 993], [-50, 0, 1133]],

'basal_1': [[-50, 0, 765], [-50, 0, 715]],

'basal_2': [[-50, 0, 715], [-156, 0, 609]],

'basal_3': [[-50, 0, 715], [56, 0, 609]],

Remove h.define_shape line

Anything else I missed? Does this sound doable @raj1701 ? Any inputs @ntolley ?

I think there are a quite a few params in the param files that are silently ignored already, right? Since most of this stuff is already specified in params_default.py, is the param file warning necessary?

We should also remove the lines where the NEURON length and radius attributes are set (see here).

jasmainak · 2023-06-22T13:13:36Z

hnn_core/io.py

+        key : cell type name
+        value : All co-ordintes of the cell types
+
+    params description


get rid of params

jasmainak · 2023-06-22T13:14:29Z

hnn_core/io.py

+    A hdf5 file containing the Network object
+
+    File Content Description
+    ------------------------


hnn_core/tests/test_io.py

jasmainak · 2023-06-22T13:20:23Z

hnn_core/tests/test_io.py

+                       match="The object should be of type Network."):
+        read_network(tmpdir.join('not_net.hdf5'))
+
+    # Add test to check weights are equal in connections and drives (todo)


add test to check each attribute is documented

jasmainak · 2023-06-22T13:22:17Z

hnn_core/tests/test_io.py

+    if network_model == "jones_2009":
+        net = jones_2009_model(add_drives_from_params=True)
+    elif network_model == "law_2021":
+        net = law_2021_model(add_drives_from_params=True)
+    elif network_model == "calcium":
+        net = calcium_model(add_drives_from_params=True)


Suggested change

if network_model == "jones_2009":

net = jones_2009_model(add_drives_from_params=True)

elif network_model == "law_2021":

net = law_2021_model(add_drives_from_params=True)

elif network_model == "calcium":

net = calcium_model(add_drives_from_params=True)

net = network_model(add_drives_from_params=True)

jasmainak · 2023-06-22T13:27:30Z

hnn_core/network.py

+        if not (self.cell_types == other.cell_types):
+            return False
+
+        # Check gid_ranges
+        if not (self.gid_ranges == other.gid_ranges):
+            return False
+
+        # Check pos_dict
+        if not (self.pos_dict == other.pos_dict):
+            return False
+
+        # Check cell_response
+        self.cell_response == other.cell_response
+
+        # Check external drives
+        if not (self.external_drives == other.external_drives):
+            return False
+
+        # Check external biases
+        if not (self.external_biases == other.external_biases):
+            return False


Suggested change

if not (self.cell_types == other.cell_types):

return False

# Check gid_ranges

if not (self.gid_ranges == other.gid_ranges):

return False

# Check pos_dict

if not (self.pos_dict == other.pos_dict):

return False

# Check cell_response

self.cell_response == other.cell_response

# Check external drives

if not (self.external_drives == other.external_drives):

return False

# Check external biases

if not (self.external_biases == other.external_biases):

return False

for attr in ['cell_response', 'external_drives']:

if getattr(self, attr) != getattr(other, attr):

return False

jasmainak · 2023-06-22T13:27:58Z

hnn_core/network.py

+        # Check extracellular arrays
+        if not (self.rec_arrays == other.rec_arrays):
+            return False
+
+        # Check threshold
+        if not (self.threshold == other.threshold):
+            return False
+
+        # Check delay
+        if not (self.delay == other.delay):
+            return False


include these as well

hnn_core/cell.py

doc/network_file_desc.rst

jasmainak · 2023-09-21T01:41:25Z

doc/network_file_desc.rst

+
+====================
+Network File Content
+====================


can you add some more information what the file structure is like? what is the file format etc? we need to provide some more information for users if they open the file using other tools and also for documentation

jasmainak · 2023-09-21T01:43:31Z

doc/conf.py

@@ -117,6 +117,7 @@
        ("Examples", "auto_examples/index"),
        ("GUI", "gui/index"),
        ("API", "api"),
+        ("Network File", "network_file_desc"),


don't add it here. You can add a line in the Readme linking to this file. It should say something to the effect that hnn-core uses hdf5 for IO and that the description of the file structure can be found here and the API of the functions for the IO can be found here

hnn_core/externals/mne.py

jasmainak · 2023-09-21T01:46:47Z

hnn_core/io.py

+
+
+@fill_doc
+def write_network(net, fname, overwrite=True, save_unsimulated=False):


you need to change save_unsimulated

What was the name we decided? save_output?

Then accordingly read_raw will become read_output?

yep, exactly! save_output and read_output

hnn_core/io.py

hnn_core/extracellular.py

hnn_core/docs.py

jasmainak · 2023-09-26T18:54:52Z

While reviewing Nick's PR I realized that we haven't updated any of the examples to use hdf5 files. Do you think that would be possible? In other words, we need to practice what we preach

raj1701 · 2023-09-27T16:16:40Z

While reviewing Nick's PR I realized that we haven't updated any of the examples to use hdf5 files. Do you think that would be possible? In other words, we need to practice what we preach

Surely can be done. We can have a discussion on this tomorrow

ntolley · 2023-09-28T14:57:45Z

doc/network_file_desc.rst

+
+File Content Description
+------------------------
+hdf5 is the file format used for storing the Network object. The network is stored in a layered format. The first layer consists of the network attributes.


Suggested change

hdf5 is the file format used for storing the Network object. The network is stored in a layered format. The first layer consists of the network attributes.

hdf5 is the file format used for storing the Network object. The network is stored in a multi-level format. The first level consists of the network attributes.

Using "layer" as the terminology might get confusing, as there is the biologically "layers" of the model

ntolley · 2023-09-28T14:58:18Z

doc/network_file_desc.rst

+------------------------
+hdf5 is the file format used for storing the Network object. The network is stored in a layered format. The first layer consists of the network attributes.
+The attributes of the network are then broken down until the network can be representated as a collection of key value pairs. For example - cell_types is a network
+attribute therefore in the first layer. The description of each cell type is in layer 2. Each cell has various sections. The description of a section is in layer 3.


Suggested change

attribute therefore in the first layer. The description of each cell type is in layer 2. Each cell has various sections. The description of a section is in layer 3.

attribute therefore in the first level. The description of each cell type is in level 2. Each cell has various sections. The description of a section is in level 3.

ntolley · 2023-09-28T14:59:47Z

doc/network_file_desc.rst

+N_pyr_x : int
+    Nr of cells (x).
+N_pyr_y : int
+    Nr of cells (y).


Suggested change

N_pyr_x : int

Nr of cells (x).

N_pyr_y : int

Nr of cells (y).

N_pyr_x : int

Number of cells (x).

N_pyr_y : int

Number of cells (y).

ntolley · 2023-09-28T15:24:01Z

hnn_core/tests/test_io.py

+    assert net_copy == net
+
+    # Test other not NotImplemented for Network Class
+    assert (net == "net") is False


Suggested change

assert (net == "net") is False

with pytest.raises(NotImplemented):

net == "net"

This is a more robust test for this error

ntolley · 2023-09-28T15:24:53Z

hnn_core/tests/test_extracellular.py

@@ -33,6 +33,9 @@ def test_extracellular_api():
    assert len(net.rec_arrays) == 2
    assert len(net.rec_arrays['arr1'].positions) == 2

+    # Test other not NotImplemented for ExtracellularArray Class
+    assert (net.rec_arrays['arr1'] == "extArr") is False


Suggested change

assert (net.rec_arrays['arr1'] == "extArr") is False

with pytest.raises(NotImplemented):

net.rec_arrays['arr1'] == "extArr"

ntolley · 2023-09-28T15:26:49Z

hnn_core/tests/test_cell.py

+    assert (cell1 == "cell") is False
+
+    # Test other not NotImplemented for Section Class
+    assert (cell1.sections['soma'] == "section") is False


Suggested change

assert (cell1 == "cell") is False

# Test other not NotImplemented for Section Class

assert (cell1.sections['soma'] == "section") is False

with pytest.raises(NotImplemented):

cell1 == "cell"

# Test other not NotImplemented for Section Class

with pytest.raises(NotImplemented):

cell1.sections['soma'] == "section"

ntolley · 2023-09-28T15:59:22Z

hnn_core/cell.py

@@ -220,6 +224,55 @@ def __init__(self, L, diam, Ra, cm, end_pts=None):
    def __repr__(self):
        return f'L={self.L}, diam={self.diam}, cm={self.cm}, Ra={self.Ra}'

+    def __eq__(self, other):
+        if not isinstance(other, Section):
+            return NotImplemented


Oh I think this is why the pytest.raises test didn't work. You're returning NotImplemented, but you actually need to do raise NotImplemented (check out how the other errors get raised)

jasmainak · 2023-10-04T03:29:58Z

@raj1701 can you make the CIs happy? you have flake8 errors

gtdang · 2024-01-11T22:00:17Z

@ntolley @jasmainak
I've changed the base of this PR to new branch on hnn-core so that we can more easily make edits before merging into the master. Shall we go ahead and merge this into the new branch so I can begin work on updating this PR?

ntolley · 2024-01-18T19:59:39Z

@gtdang can you open a PR with the new branch? I've figured out the problem you found in test_io.py and it'd be best to discuss on the new branch directly.

gtdang · 2024-01-18T21:47:30Z

@gtdang can you open a PR with the new branch? I've figured out the problem you found in test_io.py and it'd be best to discuss on the new branch directly.

Ok! the new PR is here #704

jasmainak reviewed May 24, 2023

View reviewed changes

raj1701 force-pushed the network_rw_test_1 branch from 0f4513e to 0053f53 Compare June 2, 2023 19:27

jasmainak reviewed Jun 3, 2023

View reviewed changes

hnn_core/network.py Outdated Show resolved Hide resolved

jasmainak reviewed Jun 3, 2023

View reviewed changes

hnn_core/network.py Outdated Show resolved Hide resolved

jasmainak reviewed Jun 3, 2023

View reviewed changes

hnn_core/network.py Outdated Show resolved Hide resolved

jasmainak reviewed Jun 3, 2023

View reviewed changes

jasmainak reviewed Jun 8, 2023

View reviewed changes

raj1701 force-pushed the network_rw_test_1 branch 2 times, most recently from d64acd6 to 3783e72 Compare June 21, 2023 07:35

jasmainak reviewed Jun 21, 2023

View reviewed changes

rythorpe reviewed Jun 22, 2023

View reviewed changes

jasmainak reviewed Jun 22, 2023

View reviewed changes

hnn_core/io.py Outdated

key : cell type name

value : All co-ordintes of the cell types

params description

Copy link

Collaborator

jasmainak Jun 22, 2023

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

get rid of params

jasmainak reviewed Jun 22, 2023

View reviewed changes

hnn_core/io.py Outdated

A hdf5 file containing the Network object

File Content Description

------------------------

Copy link

Collaborator

jasmainak Jun 22, 2023

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

jasmainak reviewed Jun 22, 2023

View reviewed changes

hnn_core/tests/test_io.py Outdated Show resolved Hide resolved

jasmainak reviewed Jun 22, 2023

View reviewed changes

raj1701 added 2 commits September 19, 2023 20:10

Minor fixes

96965c6

Updated documentation

79ab971