From ec1f7520e3b9f5354f7a325249272c08efcbeb18 Mon Sep 17 00:00:00 2001
From: Ben Cumming <bcumming@cscs.ch>
Date: Tue, 11 Feb 2020 14:28:28 +0100
Subject: [PATCH] enable locset names for probe sites in single cell model.
 (#955)

Use locset names to set probe sites in single cell model.

1. Extend C++ cable_cell interface to provide a concrete list of locations associated with a locset on the cell.
2. Extend Python single_cell_model to allow placement of probes on locsets.
3. Remove "debug" output in Python ring example.

The single cell model maintains a trace callback for each probe location, and needs to know the explicit locations before starting the simulation. The first step is required to enumerate and record the locations associated with a locset when attaching probes.

Open issues:

1. The extension to `cable_cell` won't sufficiently support locsets that rely on random number generators with global seeds (it would require passing additional meta data required to fully instantiate the locset)
2. Enhancing the probe recipe interface to be able to associate one `probe_id` with all locations in a locset (the number of which is not specified at the time the `probe_id` is created) would support the interface implemented for the `single_cell_model`.

Note: The richer interface on the C++ side will be implemented in time (see #957).
---
 arbor/cable_cell.cpp                  |  8 ++++++++
 arbor/include/arbor/cable_cell.hpp    |  3 +++
 python/cells.cpp                      |  4 ++++
 python/example/ring.py                |  6 +++---
 python/example/single_cell_builder.py |  9 +++++----
 python/example/single_cell_swc.py     |  2 +-
 python/morphology.cpp                 |  8 ++------
 python/single_cell_model.cpp          | 23 +++++++++++++++++------
 8 files changed, 43 insertions(+), 20 deletions(-)

diff --git a/arbor/cable_cell.cpp b/arbor/cable_cell.cpp
index f90e5bff..462c5602 100644
--- a/arbor/cable_cell.cpp
+++ b/arbor/cable_cell.cpp
@@ -103,6 +103,10 @@ struct cable_cell_impl {
             }
         }
     }
+
+    mlocation_list locations(const locset& l) const {
+        return thingify(l, provider);
+    }
 };
 
 using impl_ptr = std::unique_ptr<cable_cell_impl, void (*)(cable_cell_impl*)>;
@@ -133,6 +137,10 @@ const mprovider& cable_cell::provider() const {
     return impl_->provider;
 }
 
+mlocation_list cable_cell::locations(const locset& l) const {
+    return impl_->locations(l);
+}
+
 const cable_cell_location_map& cable_cell::location_assignments() const {
     return impl_->location_map;
 }
diff --git a/arbor/include/arbor/cable_cell.hpp b/arbor/include/arbor/cable_cell.hpp
index dc742b5a..e268062a 100644
--- a/arbor/include/arbor/cable_cell.hpp
+++ b/arbor/include/arbor/cable_cell.hpp
@@ -178,6 +178,9 @@ public:
         return location_assignments().get<i_clamp>();
     }
 
+    // Access to a concrete list of locations for a locset.
+    mlocation_list locations(const locset&) const;
+
     // Generic access to painted and placed items.
     const cable_cell_region_map& region_assignments() const;
     const cable_cell_location_map& location_assignments() const;
diff --git a/python/cells.cpp b/python/cells.cpp
index f9807142..10d98efe 100644
--- a/python/cells.cpp
+++ b/python/cells.cpp
@@ -525,6 +525,10 @@ void register_cells(pybind11::module& m) {
             },
             "locations"_a, "detector"_a,
             "Add a voltage spike detector at each location in locations.")
+        // Get locations associated with a locset label.
+        .def("locations",
+            [](arb::cable_cell& c, const char* label) {return c.locations(label);},
+            "label"_a, "The locations of the cell morphology for a locset label.")
         // Discretization control.
         .def("compartments_on_samples",
             [](arb::cable_cell& c) {c.default_parameters.discretization = arb::cv_policy_every_sample{};},
diff --git a/python/example/ring.py b/python/example/ring.py
index 94136094..ee58e1e9 100644
--- a/python/example/ring.py
+++ b/python/example/ring.py
@@ -92,7 +92,7 @@ class ring_recipe (arbor.recipe):
         loc = arbor.location(0, 0) # at the soma
         return arbor.cable_probe('voltage', id, loc)
 
-context = arbor.context(threads=4, gpu_id=None)
+context = arbor.context(threads=12, gpu_id=None)
 print(context)
 
 meters = arbor.meter_manager()
@@ -128,7 +128,7 @@ spike_recorder = arbor.attach_spike_recorder(sim)
 # Sample rate of 10 sample every ms.
 samplers = [arbor.attach_sampler(sim, 0.1, arbor.cell_member(gid,0)) for gid in range(ncells)]
 
-tfinal=1
+tfinal=100
 sim.run(tfinal)
 print(f'{sim} finished')
 
@@ -147,7 +147,7 @@ fig, ax = plt.subplots()
 for gid in range(ncells):
     times = [s.time  for s in samplers[gid].samples(arbor.cell_member(gid,0))]
     volts = [s.value for s in samplers[gid].samples(arbor.cell_member(gid,0))]
-    ax.plot(times, volts, '.')
+    ax.plot(times, volts)
 
 legends = ['cell {}'.format(gid) for gid in range(ncells)]
 ax.legend(legends)
diff --git a/python/example/single_cell_builder.py b/python/example/single_cell_builder.py
index 7c7cd7a2..486b328a 100644
--- a/python/example/single_cell_builder.py
+++ b/python/example/single_cell_builder.py
@@ -43,6 +43,8 @@ b.add_label('dend', '(join (region "dendn") (region "dendx"))')
 b.add_label('stim_site', '(location 2 0.5)')
 # The root of the tree (equivalent to '(location 0 0)')
 b.add_label('root', '(root)')
+# The tips of the dendrites (3 locations at b4, b3, b5).
+b.add_label('dtips', '(terminal)')
 
 # Extract the cable cell from the builder.
 cell = b.build()
@@ -63,13 +65,12 @@ cell.place('stim_site', arbor.iclamp( 80, 2, 0.8))
 # Add a spike detector with threshold of -10 mV.
 cell.place('root', arbor.spike_detector(-10))
 
-# make single cell model
+# Make single cell model.
 m = arbor.single_cell_model(cell)
 
 # Attach voltage probes, sampling at 10 kHz.
-m.probe('voltage', loc(0,0), 10000) # at the soma
-m.probe('voltage', loc(3,1), 10000) # at the end of branch 3
-m.probe('voltage', loc(4,1), 10000) # at the end of branch 4
+m.probe('voltage', loc(0,0), 10000) # at the soma.
+m.probe('voltage', 'dtips',   10000) # at the tips of the dendrites.
 
 # Run simulation for 100 ms of simulated activity.
 tfinal=100
diff --git a/python/example/single_cell_swc.py b/python/example/single_cell_swc.py
index 951f502f..4cdbb9c4 100644
--- a/python/example/single_cell_swc.py
+++ b/python/example/single_cell_swc.py
@@ -42,7 +42,7 @@ cell.compartments_on_samples()
 m = arbor.single_cell_model(cell)
 
 # Attach voltage probes that sample at 50 kHz.
-m.probe('voltage', where=loc(0,0),  frequency=50000)
+m.probe('voltage', where='root',  frequency=50000)
 m.probe('voltage', where=loc(2,1),  frequency=50000)
 m.probe('voltage', where=loc(4,1),  frequency=50000)
 m.probe('voltage', where=loc(30,1), frequency=50000)
diff --git a/python/morphology.cpp b/python/morphology.cpp
index 0e503bf5..6f991901 100644
--- a/python/morphology.cpp
+++ b/python/morphology.cpp
@@ -41,13 +41,9 @@ void register_morphology(pybind11::module& m) {
         .def_readonly("position", &arb::mlocation::pos,
             "The relative position on the branch (∈ [0.,1.], where 0. means proximal and 1. distal).")
         .def("__str__",
-            [](arb::mlocation l) {
-                return util::pprintf("(location {} {})", l.branch, l.pos);
-            })
+            [](arb::mlocation l) { return util::pprintf("(location {} {})", l.branch, l.pos); })
         .def("__repr__",
-            [](arb::mlocation l) {
-                return util::pprintf("<arbor.location: branch {}, position {}>", l.branch, l.pos);
-            });
+            [](arb::mlocation l) { return util::pprintf("(location {} {})", l.branch, l.pos); });
 
     // arb::mpoint
     pybind11::class_<arb::mpoint> mpoint(m, "mpoint");
diff --git a/python/single_cell_model.cpp b/python/single_cell_model.cpp
index 6d5084c4..29d5df5d 100644
--- a/python/single_cell_model.cpp
+++ b/python/single_cell_model.cpp
@@ -157,7 +157,10 @@ public:
 
     // example use:
     //      m.probe('voltage', arbor.location(2,0.5))
-    void probe(const std::string& what, const arb::mlocation& where, double frequency) {
+    //      m.probe('voltage', '(location 2 0.5)')
+    //      m.probe('voltage', 'term')
+
+    void probe(const std::string& what, const arb::locset& where, double frequency) {
         if (what != "voltage") {
             throw pyarb_error(
                 util::pprintf("{} does not name a valid variable to trace (currently only 'voltage' is supported)", what));
@@ -166,11 +169,9 @@ public:
             throw pyarb_error(
                 util::pprintf("sampling frequency is not greater than zero", what));
         }
-        if (where.branch>=cell_.morphology().num_branches()) {
-            throw pyarb_error(
-                util::pprintf("invalid location", what));
+        for (auto& l: cell_.locations(where)) {
+            probes_.push_back({l, frequency});
         }
-        probes_.push_back({where, frequency});
     }
 
     void add_ion(const std::string& ion, double valence, double int_con, double ext_con, double rev_pot) {
@@ -238,7 +239,17 @@ void register_single_cell(pybind11::module& m) {
         .def(pybind11::init<arb::cable_cell>(),
             "cell"_a, "Initialise a single cell model for a cable cell.")
         .def("run", &single_cell_model::run, "tfinal"_a, "Run model from t=0 to t=tfinal ms.")
-        .def("probe", &single_cell_model::probe,
+        .def("probe",
+            [](single_cell_model& m, const char* what, const char* where, double frequency) {
+                m.probe(what, where, frequency);},
+            "what"_a, "where"_a, "frequency"_a,
+            "Sample a variable on the cell.\n"
+            " what:      Name of the variable to record (currently only 'voltage').\n"
+            " where:     Location on cell morphology at which to sample the variable.\n"
+            " frequency: The target frequency at which to sample [Hz].")
+        .def("probe",
+            [](single_cell_model& m, const char* what, const arb::mlocation& where, double frequency) {
+                m.probe(what, where, frequency);},
             "what"_a, "where"_a, "frequency"_a,
             "Sample a variable on the cell.\n"
             " what:      Name of the variable to record (currently only 'voltage').\n"
-- 
GitLab