[BBPBGLIB-1036] Implement SONATA synapse reader for gap junction

## Context
Recently we have dropped the synapse-tool dependency for reading synapse parameters. SynReaderNRN is used to read the old format nrn edges files, and SonataReader is to read sonata edges. This PR is to implement a sonata reader to read gap junction synapse parameters.

## Scope
In `gap_junection.py`, a GapJunctionSonata reader derived from the SonataReader. The Sonata parameters for gap junction are `efferent_junction_id` and `afferent_junction_id`.

## Testing
Add a unit test to test GapJunctionSynapseReader for both NRN and SONATA edges files.

## Review
* [x] PR description is complete
* [x] Coding style (imports, function length, New functions, classes or files) are good
* [x] Unit/Scientific test added 
* [ ] Updated Readme, in-code, developer documentation

neurodamus/gap_junction.py

@@ -8,7 +8,8 @@
 
 from .connection_manager import ConnectionManagerBase
 from .core.configuration import ConfigurationError
-from .io.synapse_reader import SynapseReader, SynReaderNRN, SynapseParameters, FormatNotSupported
+from .io.synapse_reader import SynapseReader, SynReaderNRN, SonataReader, SynapseParameters, \
+    FormatNotSupported
 from .utils import compat
 from .utils.logging import log_verbose
 
@@ -17,40 +18,9 @@ class GapJunctionConnParameters(SynapseParameters):
     # Attribute names of synapse parameters, consistent with the normal synapses
     _synapse_fields = ("sgid", "isec", "offset", "weight", "D", "F", "ipt", "location")
 
-    # Actual fields to read from conectivity files
-    _gj_v1_fields = [
-        "connected_neurons_pre",
-        "morpho_section_id_post",
-        "morpho_offset_segment_post",
-        "conductance",
-        "junction_id_pre",
-        "junction_id_post",
-        "morpho_segment_id_post",
-
-    ]
-    _gj_v2_fields = [
-        "connected_neurons_pre",
-        "morpho_section_id_post",
-        "morpho_section_fraction_post",  # v2 field
-        "conductance",
-        "junction_id_pre",
-        "junction_id_post"
-    ]
-    # SONATA fields, see conversion map in spykfunc/schema.py and
-    # synapse-tool Naming::get_property_mapping
-    _gj_sonata_fields = [
-        "connected_neurons_pre",
-        "morpho_section_id_post",
-        "morpho_section_fraction_post",
-        "conductance",
-        "efferent_junction_id",
-        "afferent_junction_id"
-    ]
-
     @classmethod
     def create_array(cls, length):
         npa = np.recarray(length, cls.dtype)
-        npa.ipt = -1
         npa.location = 0.5
         return npa
 
@@ -65,13 +35,27 @@ class GapJunctionSynapseReader(SynapseReader):
     def create(cls, syn_src, conn_type, population=None, *args, **kw):
         """Instantiates a synapse reader, giving preference to GapJunctionSynToolReader
         """
+        if fn := cls._get_sonata_circuit(syn_src):
+            log_verbose("[SynReader] Using SonataReader.")
+            return GapJunctionSonataReader(fn, conn_type, population, **kw)
+
         # Syn2 support dropped (July 2023)
         if not ospath.isdir(syn_src) and not syn_src.endswith(".h5"):
             raise FormatNotSupported("File: {syn_src}")
         logging.info("[GapJunctionSynReader] Attempting legacy hdf5 reader.")
         return SynReaderNRN(syn_src, conn_type, None, *args, **kw)
 
 
+class GapJunctionSonataReader(SonataReader):
+    Parameters = GapJunctionConnParameters
+    parameter_mapping = {
+        "weight": "conductance",
+        "D": "efferent_junction_id",
+        "F": "afferent_junction_id",
+    }
+    # "isec", "ipt", "offset" are custom parameters as in base class
+
+
 class GapJunctionManager(ConnectionManagerBase):
     """
     The GapJunctionManager is similar to the SynapseRuleManager. It will

neurodamus/io/synapse_reader.py

@@ -12,23 +12,6 @@
 from ..utils.logging import log_verbose
 
 
-def _get_sonata_circuit(path):
-    """Returns a SONATA edge file in path if present
-    """
-    if os.path.isdir(path):
-        filename = os.path.join(path, "edges.sonata")
-        if os.path.exists(filename):
-            return filename
-    elif path.endswith(".sonata"):
-        return path
-    elif path.endswith(".h5"):
-        import h5py
-        f = h5py.File(path, 'r')
-        if "edges" in f:
-            return path
-    return None
-
-
 def _constrained_hill(K_half, y):
     K_half_fourth = K_half**4
     y_fourth = y**4
@@ -185,6 +168,23 @@ def has_property(self, field_name):
         """Checks whether source data has the given additional field.
         """
 
+    @staticmethod
+    def _get_sonata_circuit(path):
+        """Returns a SONATA edge file in path if present
+        """
+        if os.path.isdir(path):
+            filename = os.path.join(path, "edges.sonata")
+            if os.path.exists(filename):
+                return filename
+        elif path.endswith(".sonata"):
+            return path
+        elif path.endswith(".h5"):
+            import h5py
+            f = h5py.File(path, 'r')
+            if "edges" in f:
+                return path
+        return None
+
     @classmethod
     def create(cls, syn_src, conn_type=SYNAPSES, population=None, *args, **kw):
         """Instantiates a synapse reader, giving preference to SonataReader
@@ -194,7 +194,7 @@ def create(cls, syn_src, conn_type=SYNAPSES, population=None, *args, **kw):
             return cls(syn_src, conn_type, population, *args, **kw)
 
         kw["verbose"] = (MPI.rank == 0)
-        if fn := _get_sonata_circuit(syn_src):
+        if fn := cls._get_sonata_circuit(syn_src):
             log_verbose("[SynReader] Using SonataReader.")
             return SonataReader(fn, conn_type, population, **kw)
 

tests/simulations/mini_thalamus_sonata/gapjunction/README

@@ -0,0 +1,2 @@
+Extract connections with target gid = 0 from the gap junction edges file in 
+/gpfs/bbp.cscs.ch/project/proj12/SIT/thalamus_sonata/networks/edges/thalamus_neurons__thalamus_neurons__electrical_synapse/edges.h5

tests/test_synapse_reader.py

@@ -0,0 +1,48 @@
+import numpy as np
+import numpy.testing as npt
+import os
+import pytest
+from neurodamus.gap_junction import GapJunctionSynapseReader
+from pathlib import Path
+
+
+pytestmark = [
+    pytest.mark.forked,
+    pytest.mark.skipif(
+        not os.environ.get("NEURODAMUS_NEOCORTEX_ROOT"),
+        reason="Test requires loading a neocortex model to run"
+    )
+]
+
+
+def test_gapjunction_synreaderNRN():
+    nrn_file = "/gpfs/bbp.cscs.ch/project/proj12/jenkins/cellular/circuit-scx-v5-gapjunctions/gap_junctions/nrn_gj.h5" # noqa
+    nrn_reader = GapJunctionSynapseReader.create(nrn_file, 1)
+    syn_params_nrn = nrn_reader._load_synapse_parameters(100124)
+    # check reading of sgid, junction_id_pre and junction_id_post
+    ref_sgids = np.array([94669., 94723., 95634., 95823., 96581.,
+                          97338., 97455., 98139., 98432., 100725.,
+                          101360., 101506., 101696., 101696., 191567.])
+    ref_junction_id_pre = np.array([735., 736., 29., 36., 51.,
+                                    77., 744., 134., 148., 286.,
+                                    322., 337., 355., 356., 681.])
+    ref_junction_id_post = np.array([1251., 1259., 617., 1354., 1002.,
+                                     1756., 1027., 924., 709., 624.,
+                                     1050., 521., 592., 593., 590.])
+    npt.assert_allclose(syn_params_nrn.sgid, ref_sgids)
+    npt.assert_allclose(syn_params_nrn.D, ref_junction_id_pre)
+    npt.assert_allclose(syn_params_nrn.F, ref_junction_id_post)
+
+
+def test_gapjunction_sonata_reader():
+    SIM_DIR = Path(__file__).parent.absolute() / "simulations"
+    sonata_file = str(SIM_DIR / "mini_thalamus_sonata/gapjunction/edges.h5")
+    sonata_reader = GapJunctionSynapseReader.create(sonata_file, 1)
+    syn_params_sonata = sonata_reader._load_synapse_parameters(1)
+    ref_junction_id_pre = np.array([10257., 43930., 226003., 298841., 324744.,
+                                    1094745., 1167632., 1172523., 1260104.])
+    ref_junction_id_post = np.array([14., 52., 71., 76., 78., 84., 89., 90., 93.])
+    ref_weight = np.array([0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2])
+    npt.assert_allclose(syn_params_sonata.D, ref_junction_id_pre)
+    npt.assert_allclose(syn_params_sonata.F, ref_junction_id_post)
+    npt.assert_allclose(syn_params_sonata.weight, ref_weight)