Add neuroML python example (#1386)

* add neuroML example
* no test, because GH actions arbor is built without neuroML support

.github/workflows/basic.yml

@@ -165,6 +165,6 @@ jobs:
             python python/example/single_cell_model.py
             python python/example/single_cell_recipe.py
             python python/example/single_cell_stdp.py
-            python python/example/single_cell_swc.py test/unit/swc/pyramidal.swc
+            python python/example/single_cell_swc.py python/example/morph.swc
             python python/example/single_cell_detailed.py python/example/morph.swc
             python python/example/single_cell_detailed_recipe.py python/example/morph.swc

doc/python/morphology.rst

@@ -562,7 +562,6 @@ constitute part of the CV boundary point set.
     A :class:`neuroml_morphology_data` object contains a representation of a morphology defined in
     NeuroML.
 
-
     .. py:attribute:: cell_id
        :type: optional<str>
 
@@ -578,6 +577,11 @@ constitute part of the CV boundary point set.
 
        A map from each segment group id to its corresponding collection of segments.
 
+    .. py:attribute:: morphology
+       :type: morphology
+
+       The morphology associated with the :class:`neuroml_morph_data` object.
+
     .. py:method:: segments
 
        Returns a label dictionary with a region entry for each segment, keyed by the segment id (as a string).

python/cells.cpp

@@ -225,6 +225,12 @@ void register_cells(pybind11::module& m) {
             [](const label_dict_proxy &ld) {
                 return pybind11::make_key_iterator(ld.cache.begin(), ld.cache.end());},
             pybind11::keep_alive<0, 1>())
+        .def("append", [](label_dict_proxy& l, const label_dict_proxy& other, const char* prefix) {
+                l.import(other, prefix);
+            },
+            "other"_a, "The label_dict to be imported"
+            "prefix"_a="", "optional prefix appended to the region and locset labels",
+            "Import the entries of a another label dictionary with an optional prefix.")
         .def_readonly("regions", &label_dict_proxy::regions,
              "The region definitions.")
         .def_readonly("locsets", &label_dict_proxy::locsets,

python/example/morph.nml

+<neuroml xmlns="http://www.neuroml.org/schema/neuroml2">
+<morphology id="m1">
+    <segment id="0">
+        <proximal x="0" y="0" z="0" diameter="4"/>
+        <distal  x="40" y="0" z="0" diameter="4"/>
+    </segment>
+    <segment id="1">
+        <parent segment="0"/>
+        <proximal x="40" y="0" z="0" diameter="1.6"/>
+        <distal   x="80" y="0" z="0" diameter="1.6"/>
+    </segment>
+    <segment id="2">
+        <parent segment="1"/>
+        <proximal x="80"  y="0"  z="0" diameter="1.6"/>
+        <distal   x="120" y="-5" z="0" diameter="1.6"/>
+    </segment>
+    <segment id="3">
+        <parent segment="2"/>
+        <proximal x="120" y="-5" z="0" diameter="1.6"/>
+        <distal   x="200" y="40" z="0" diameter="0.8"/>
+    </segment>
+    <segment id="4">
+        <parent segment="3"/>
+        <proximal x="200" y="40" z="0" diameter="0.8"/>
+        <distal   x="260" y="60" z="0" diameter="0.4"/>
+    </segment>
+    <segment id="5">
+        <parent segment="2"/>
+        <proximal x="120" y="-5"  z="0" diameter="1"/>
+        <distal   x="190" y="-30" z="0" diameter="1"/>
+    </segment>
+    <segment id="6">
+        <parent segment="5"/>
+        <proximal x="190" y="-30" z="0" diameter="1"/>
+        <distal   x="240" y="-70" z="0" diameter="0.4"/>
+    </segment>
+    <segment id="7">
+        <parent segment="5"/>
+        <proximal x="190" y="-30" z="0" diameter="1"/>
+        <distal   x="230" y="-10" z="0" diameter="0.4"/>
+    </segment>
+    <segment id="8">
+        <parent segment="7"/>
+        <proximal x="230" y="-10" z="0" diameter="0.4"/>
+        <distal   x="360" y="-20" z="0" diameter="0.4"/>
+    </segment>
+    <segment id="9">
+        <parent segment="0" fractionAlong="0.0"/>
+        <proximal x="0"   y="0" z="0" diameter="4"/>
+        <distal   x="-70" y="0" z="0" diameter="0.8"/>
+    </segment>
+    <segment id="10">
+        <parent segment="9"/>
+        <proximal x="-70"  y="0" z="0" diameter="0.8"/>
+        <distal   x="-100" y="0" z="0" diameter="0.8"/>
+    </segment>
+    <segmentGroup id="soma">
+        <member segment="0"/>
+    </segmentGroup>
+    <segmentGroup id="axon">
+        <member segment="9"/>
+        <member segment="10"/>
+    </segmentGroup>
+    <segmentGroup id="dend">
+        <member segment="1"/>
+        <member segment="2"/>
+        <member segment="3"/>
+        <member segment="4"/>
+        <member segment="5"/>
+    </segmentGroup>
+</morphology>
+</neuroml>

python/example/single_cell_nml.py

+#!/usr/bin/env python3
+import arbor
+from arbor import mechanism as mech
+from arbor import location as loc
+import pandas, seaborn
+import sys
+
+# Load a cell morphology from an nml file.
+# Example present here: morph.nml
+if len(sys.argv) < 2:
+    print("No NeuroML file passed to the program")
+    sys.exit(0)
+
+filename = sys.argv[1]
+
+# Read the NeuroML morphology from the file.
+morpho_nml = arbor.neuroml(filename)
+
+# Read the morphology data associated with morphology "m1".
+morpho_data = morpho_nml.morphology("m1")
+
+# Get the morphology.
+morpho = morpho_data.morphology
+
+# Get the region label dictionaries associated with the morphology.
+morpho_segments = morpho_data.segments()
+morpho_named = morpho_data.named_segments()
+morpho_groups = morpho_data.groups()
+
+# Create new label dict add to it all the NeuroML dictionaries.
+labels = arbor.label_dict()
+labels.append(morpho_segments)
+labels.append(morpho_named)
+labels.append(morpho_groups)
+
+# Add locsets to the label dictionary. 
+labels['stim_site'] = '(location 1 0.5)' # site for the stimulus, in the middle of branch 1.
+labels['axon_end']  = '(restrict (terminal) (region "axon"))' # end of the axon.
+labels['root']      = '(root)' # the start of the soma in this morphology is at the root of the cell.
+
+# Optional: print out the regions and locsets available in the label dictionary.
+print("Label dictionary regions: ", labels.regions, "\n")
+print("Label dictionary locsets: ", labels.locsets, "\n")
+
+decor = arbor.decor()
+
+# Set initial membrane potential to -55 mV
+decor.set_property(Vm=-55)
+# Use Nernst to calculate reversal potential for calcium.
+decor.set_ion('ca', method=mech('nernst/x=ca'))
+#decor.set_ion('ca', method='nernst/x=ca')
+# hh mechanism on the soma and axon.
+decor.paint('"soma"', 'hh')
+decor.paint('"axon"', 'hh')
+# pas mechanism the dendrites.
+decor.paint('"dend"', 'pas')
+# Increase resistivity on dendrites.
+decor.paint('"dend"', rL=500)
+# Attach stimuli that inject 4 nA current for 1 ms, starting at 3 and 8 ms.
+decor.place('"root"', arbor.iclamp(10, 1, current=5))
+decor.place('"stim_site"', arbor.iclamp(3, 1, current=0.5))
+decor.place('"stim_site"', arbor.iclamp(10, 1, current=0.5))
+decor.place('"stim_site"', arbor.iclamp(8, 1, current=4))
+# Detect spikes at the soma with a voltage threshold of -10 mV.
+decor.place('"axon_end"', arbor.spike_detector(-10))
+
+# Create the policy used to discretise the cell into CVs.
+# Use a single CV for the soma, and CVs of maximum length 1 μm elsewhere.
+soma_policy = arbor.cv_policy_single('"soma"')
+dflt_policy = arbor.cv_policy_max_extent(1.0)
+policy = dflt_policy | soma_policy
+decor.discretization(policy)
+
+# Combine morphology with region and locset definitions to make a cable cell.
+cell = arbor.cable_cell(morpho, labels, decor)
+
+print(cell.locations('"axon_end"'))
+
+# Make single cell model.
+m = arbor.single_cell_model(cell)
+
+# Attach voltage probes that sample at 50 kHz.
+m.probe('voltage', where='"root"',  frequency=50000)
+m.probe('voltage', where='"stim_site"',  frequency=50000)
+m.probe('voltage', where='"axon_end"', frequency=50000)
+
+# Simulate the cell for 15 ms.
+tfinal=15
+m.run(tfinal)
+print("Simulation done.")
+
+# Print spike times.
+if len(m.spikes)>0:
+    print('{} spikes:'.format(len(m.spikes)))
+    for s in m.spikes:
+        print('  {:7.4f}'.format(s))
+else:
+    print('no spikes')
+
+# Plot the recorded voltages over time.
+print("Plotting results ...")
+df_list = []
+for t in m.traces:
+    df_list.append(pandas.DataFrame({'t/ms': t.time, 'U/mV': t.value, 'Location': str(t.location), "Variable": t.variable}))
+
+df = pandas.concat(df_list)
+
+seaborn.relplot(data=df, kind="line", x="t/ms", y="U/mV",hue="Location",col="Variable",ci=None).savefig('single_cell_nml.svg')

python/example/single_cell_swc.py

@@ -17,7 +17,7 @@ import pandas, seaborn
 import sys
 
 # Load a cell morphology from an swc file.
-# Example present here: ../../test/unit/swc/pyramidal.swc
+# Example present here: morph.swc
 if len(sys.argv) < 2:
     print("No SWC file passed to the program")
     sys.exit(0)

python/proxy.hpp

@@ -23,12 +23,31 @@ struct label_dict_proxy {
         }
     }
 
-    label_dict_proxy(const arb::label_dict& label_dict): dict(label_dict) {}
+    label_dict_proxy(const arb::label_dict& label_dict): dict(label_dict) {
+        for (const auto& [l, reg]: dict.regions()) {
+            regions.push_back(l);
+        }
+        for (const auto& [l, ls]: dict.locsets()) {
+            locsets.push_back(l);
+        }
+    }
 
     std::size_t size() const  {
         return locsets.size() + regions.size();
     }
 
+    void import(const label_dict_proxy& other, std::string prefix) {
+        regions.clear();
+        locsets.clear();
+        dict.import(other.dict, prefix);
+        for (const auto& [l, reg]: dict.regions()) {
+            regions.push_back(l);
+        }
+        for (const auto& [l, ls]: dict.locsets()) {
+            locsets.push_back(l);
+        }
+    }
+
     void set(const char* name, const char* desc) {
         using namespace std::string_literals;
         // The following code takes an input name and a region or locset