Expose CV boundary to Python and attach unit to CV max extent

arbor/cv_policy.cpp

@@ -120,12 +120,16 @@ struct cvp_cv_policy_max_extent {
     region domain() const { return domain_; }
 };
 
-ARB_ARBOR_API cv_policy cv_policy_max_extent(double ext, region reg, cv_policy_flag flag) {
-    return cv_policy{cvp_cv_policy_max_extent{ext, std::move(reg), flag}};
+ARB_ARBOR_API cv_policy cv_policy_max_extent(const units::quantity& ext, region reg, cv_policy_flag flag) {
+    auto ext_um = ext.value_as(units::um);
+    if (!std::isfinite(ext_um) || ext_um < 0) throw std::domain_error("Max extent requires a positive length.");
+    return cv_policy{cvp_cv_policy_max_extent{ext_um, std::move(reg), flag}};
 }
 
-ARB_ARBOR_API cv_policy cv_policy_max_extent(double ext, cv_policy_flag flag) {
-    return cv_policy{cvp_cv_policy_max_extent{ext, reg::all(), flag}};
+ARB_ARBOR_API cv_policy cv_policy_max_extent(const units::quantity& ext, cv_policy_flag flag) {
+    auto ext_um = ext.value_as(units::um);
+    if (!std::isfinite(ext_um) || ext_um < 0) throw std::domain_error("Max extent requires a positive length.");
+    return cv_policy{cvp_cv_policy_max_extent{ext_um, reg::all(), flag}};
 }
 
 struct cvp_cv_policy_explicit {

arbor/include/arbor/cv_policy.hpp

@@ -1,5 +1,6 @@
 #pragma once
 
+#include "arbor/units.hpp"
 #include <memory>
 #include <utility>
 
@@ -130,8 +131,8 @@ ARB_ARBOR_API cv_policy operator|(const cv_policy&, const cv_policy&);
 
 ARB_ARBOR_API cv_policy cv_policy_explicit(locset, region = reg::all());
 
-ARB_ARBOR_API cv_policy cv_policy_max_extent(double, region, cv_policy_flag = cv_policy_flag::none);
-ARB_ARBOR_API cv_policy cv_policy_max_extent(double, cv_policy_flag = cv_policy_flag::none);
+ARB_ARBOR_API cv_policy cv_policy_max_extent(const units::quantity& len, region, cv_policy_flag = cv_policy_flag::none);
+ARB_ARBOR_API cv_policy cv_policy_max_extent(const units::quantity& len, cv_policy_flag = cv_policy_flag::none);
 
 ARB_ARBOR_API cv_policy cv_policy_fixed_per_branch(unsigned, region, cv_policy_flag = cv_policy_flag::none);
 ARB_ARBOR_API cv_policy cv_policy_fixed_per_branch(unsigned, cv_policy_flag = cv_policy_flag::none);

arborio/cv_policy_parse.cpp

@@ -44,13 +44,13 @@ eval_map {{"default",
            make_call<int, region, cv_policy_flag>([] (int i, const region& r, cv_policy_flag f) { return arb::cv_policy_fixed_per_branch(i, r, f); },
                                                   "'fixed-per-branch' with three arguments (fixed-per-branch (count:int) (reg:region) (flags:flag))")},
           {"max-extent",
-           make_call<double>([] (double i) { return arb::cv_policy_max_extent(i); },
+           make_call<double>([] (double i) { return arb::cv_policy_max_extent(i * arb::units::um); },
                              "'max-extent' with one argument (max-extent (length:double))")},
           {"max-extent",
-           make_call<double, region>([] (double i, const region& r) { return arb::cv_policy_max_extent(i, r); },
+           make_call<double, region>([] (double i, const region& r) { return arb::cv_policy_max_extent(i * arb::units::um, r); },
                                      "'max-extent' with two arguments (max-extent (length:double) (reg:region))")},
           {"max-extent",
-           make_call<double, region, cv_policy_flag>([] (double i, const region& r, cv_policy_flag f) { return arb::cv_policy_max_extent(i, r, f); },
+           make_call<double, region, cv_policy_flag>([] (double i, const region& r, cv_policy_flag f) { return arb::cv_policy_max_extent(i * arb::units::um, r, f); },
                                                      "'max-extent' with three arguments (max-extent (length:double) (reg:region) (flags:flag))")},
           {"single",
            make_call<>([] () { return arb::cv_policy_single(); },

example/diffusion/diffusion.cpp

@@ -26,7 +26,7 @@ namespace U = arb::units;
 struct linear: public recipe {
      linear(double ext, double dx, double Xi, double beta): l{ext}, d{dx}, i{Xi}, b{beta} {
         gprop.default_parameters = neuron_parameter_defaults;
-        gprop.default_parameters.discretization = cv_policy_max_extent(d);
+        gprop.default_parameters.discretization = cv_policy_max_extent(d*U::um);
         gprop.add_ion("bla", 1, 23*U::mM, 42*U::mM, 0*U::mV, b*U::m2/U::s);
     }
 

example/ornstein_uhlenbeck/ou.cpp

@@ -33,7 +33,7 @@ class recipe: public arb::recipe {
         // single-cell tree with ncvs control volumes
         segment_tree tree;
         tree.append(mnpos, {0, 0, 0.0, 4.0}, {0, 0, ncvs*cv_size, 4.0}, 1);
-        cell_ = cable_cell(morphology(tree), dec, {}, cv_policy_max_extent(cv_size));
+        cell_ = cable_cell(morphology(tree), dec, {}, cv_policy_max_extent(cv_size * U::um));
     }
 
     arb::cell_size_type num_cells() const override { return 1; }

python/cells.cpp

@@ -111,7 +111,7 @@ arb::cv_policy make_cv_policy_fixed_per_branch(unsigned cv_per_branch, const std
     return arb::cv_policy_fixed_per_branch(cv_per_branch, arborio::parse_region_expression(reg).unwrap());
 }
 
-arb::cv_policy make_cv_policy_max_extent(double cv_length, const std::string& reg) {
+arb::cv_policy make_cv_policy_max_extent(const arb::units::quantity& cv_length, const std::string& reg) {
     return arb::cv_policy_max_extent(cv_length, arborio::parse_region_expression(reg).unwrap());
 }
 
@@ -324,6 +324,14 @@ void register_cells(py::module& m) {
             ss << p;
             return ss.str();
         })
+        .def("boundary_points",
+             [] (const arb::cv_policy& cvp,
+                 const arb::cable_cell& cell) {
+                 const auto ls = cvp.cv_boundary_points(cell);
+                 return thingify(ls, cell.provider());
+             },
+             "cell"_a, "cable cell to compute CVs for.",
+             "Compute CV boundaries for cable cell")
         .def("__str__", [](const arb::cv_policy& p) {
             std::stringstream ss;
             ss << p;
@@ -360,6 +368,14 @@ void register_cells(py::module& m) {
 
     // arb::cell_cv_data
     cell_cv_data
+        .def(py::init([](const arb::cable_cell& cell) {
+            if (auto cvd = arb::cv_data(cell); cvd) {
+                return std::move(*cvd);
+            }
+            else {
+                throw std::runtime_error("Could not construct cell_cv_data.");
+            } 
+        }))
         .def_property_readonly("num_cv", [](const arb::cell_cv_data& data){return data.size();},
                                "Return the number of CVs in the cell.")
         .def("cables",

python/example/network_two_cells_gap_junctions.py

@@ -60,7 +60,7 @@ def cell_description(self, gid):
         )
 
         if self.max_extent is not None:
-            cvp = A.cv_policy_max_extent(self.max_extent)
+            cvp = A.cv_policy_max_extent(self.max_extent * U.um)
         else:
             cvp = A.cv_policy_single()
 

python/example/single_cell_allen.py

@@ -131,7 +131,7 @@ def make_cell(base, swc, fit):
     decor.place('"midpoint"', A.threshold_detector(-40 * U.mV), "sd")
 
     # (10) discretisation strategy: max compartment length
-    cvp = A.cv_policy_max_extent(20)
+    cvp = A.cv_policy_max_extent(20 * U.um)
 
     # (11) Create cell
     return A.cable_cell(morphology, decor, labels, cvp), offset

python/example/single_cell_bluepyopt_l5pc.py

@@ -1,6 +1,7 @@
 #!/usr/bin/env python3
 
 import arbor as A
+from arbor import units as U
 import pandas
 import seaborn
 import sys
@@ -28,20 +29,22 @@
 
 # (3) Define stimulus and spike detector, adjust discretization
 decor.place(
-    '"soma_center"', A.iclamp(tstart=295, duration=5, current=1.9), "soma_iclamp"
+    '"soma_center"',
+    A.iclamp(tstart=295 * U.ms, duration=5 * U.ms, current=1.9 * U.nA),
+    "soma_iclamp",
 )
 
 # Add spike detector
-decor.place('"soma_center"', A.threshold_detector(-10), "detector")
+decor.place('"soma_center"', A.threshold_detector(-10 * U.mV), "detector")
 
 # Adjust discretization (single CV on soma, default everywhere else)
-cvp = A.cv_policy_max_extent(1.0) | A.cv_policy_single('"soma"')
+cvp = A.cv_policy_max_extent(1.0 * U.um) | A.cv_policy_single('"soma"')
 
 # (4) Create the cell.
 cell = A.cable_cell(morpho, decor, labels, cvp)
 
 # (5) Declare a probe.
-probe = A.cable_probe_membrane_voltage('"dend1"')
+probe = A.cable_probe_membrane_voltage('"dend1"', tag="Um")
 
 
 # (6) Create a class that inherits from A.recipe
@@ -92,11 +95,11 @@ def global_properties(self, gid):
 # Instruct the simulation to record the spikes and sample the probe
 sim.record(A.spike_recording.all)
 
-probe_id = A.cell_member(0, 0)
-handle = sim.sample(probe_id, A.regular_schedule(0.02))
+probe_id = A.cell_member(0, "Um")
+handle = sim.sample(probe_id, A.regular_schedule(0.02 * U.ms))
 
 # (8) Run the simulation
-sim.run(tfinal=600, dt=0.025)
+sim.run(tfinal=600 * U.ms, dt=0.025 * U.ms)
 
 # (9) Print or display the results
 spikes = sim.spikes()

python/example/single_cell_bluepyopt_simplecell.py

@@ -1,6 +1,7 @@
 #!/usr/bin/env python3
 
 import arbor as A
+from arbor import units as U
 import pandas
 import seaborn
 import sys
@@ -24,14 +25,16 @@
 
 # (3) Define stimulus and spike detector, adjust discretization
 decor.place(
-    '"soma_center"', A.iclamp(tstart=100, duration=50, current=0.05), "soma_iclamp"
+    '"soma_center"',
+    A.iclamp(tstart=100 * U.ms, duration=50 * U.ms, current=0.05 * U.nA),
+    "soma_iclamp",
 )
 
 # Add spike detector
 decor.place('"soma_center"', A.threshold_detector(-10), "detector")
 
 # Adjust discretization (single CV on soma, default everywhere else)
-cvp = A.cv_policy_max_extent(1.0) | A.cv_policy_single('"soma"')
+cvp = A.cv_policy_max_extent(1.0 * U.um) | A.cv_policy_single('"soma"')
 
 # (4) Create the cell.
 cell = A.cable_cell(morpho, decor, labels, cvp)
@@ -45,10 +48,10 @@
 m.properties.catalogue.extend(A.bbp_catalogue(), "BBP::")
 
 # (6) Attach voltage probe that samples at 50 kHz.
-m.probe("voltage", where='"soma_center"', frequency=50)
+m.probe("voltage", where='"soma_center"', tag="Um", frequency=50 * U.kHz)
 
 # (7) Simulate the cell for 200 ms.
-m.run(tfinal=200)
+m.run(tfinal=200 * U.ms)
 print("Simulation done.")
 
 # (8) Print spike times.

python/example/single_cell_cable.py

@@ -56,7 +56,7 @@ def __init__(
         self.stimulus_duration = stimulus_duration * U.ms
         self.stimulus_amplitude = stimulus_amplitude * U.nA
 
-        self.cv_policy_max_extent = cv_policy_max_extent
+        self.cv_policy_max_extent = cv_policy_max_extent * U.um
 
         self.the_props = A.neuron_cable_properties()
 

test/ubench/mech_vec.cpp

@@ -6,7 +6,6 @@
 // will need to be reworked in order to compile.
 
 #include <any>
-#include <fstream>
 
 #include <arbor/cable_cell.hpp>
 #include <arbor/morph/segment_tree.hpp>
@@ -70,7 +69,7 @@ class recipe_expsyn_1_branch: public recipe {
             decor.place(arb::mlocation{0, distribution(gen)}, arb::synapse("expsyn"), "syn");
         }
 
-        return arb::cable_cell{arb::morphology(tree), decor, {}, arb::cv_policy_max_extent((dend_length+soma_radius*2)/num_comp_)};
+        return arb::cable_cell{arb::morphology(tree), decor, {}, arb::cv_policy_max_extent((dend_length+soma_radius*2)/num_comp_ * arb::units::um)};
     }
 
     virtual cell_kind get_cell_kind(cell_gid_type) const override {
@@ -109,7 +108,7 @@ class recipe_pas_1_branch: public recipe {
 
         arb::decor decor;
         decor.paint(arb::reg::all(), arb::density("pas"));
-        return arb::cable_cell {arb::morphology(tree), decor, {}, arb::cv_policy_max_extent((dend_length+soma_radius*2)/num_comp_)};
+        return arb::cable_cell {arb::morphology(tree), decor, {}, arb::cv_policy_max_extent((dend_length+soma_radius*2)/num_comp_ * arb::units::um)};
     }
 
     virtual cell_kind get_cell_kind(cell_gid_type) const override {
@@ -151,7 +150,7 @@ class recipe_pas_3_branches: public recipe {
         arb::decor decor;
         decor.paint(arb::reg::all(), arb::density("pas"));
 
-        return arb::cable_cell{arb::morphology(tree), decor, {}, arb::cv_policy_max_extent((dend_length*3+soma_radius*2)/num_comp_)};
+        return arb::cable_cell{arb::morphology(tree), decor, {}, arb::cv_policy_max_extent((dend_length*3+soma_radius*2)/num_comp_ * arb::units::um)};
     }
 
     virtual cell_kind get_cell_kind(cell_gid_type) const override {
@@ -191,7 +190,7 @@ class recipe_hh_1_branch: public recipe {
         arb::decor decor;
         decor.paint(arb::reg::all(), arb::density("hh"));
 
-        return arb::cable_cell{arb::morphology(tree), decor, {}, arb::cv_policy_max_extent((dend_length+soma_radius*2)/num_comp_)};
+        return arb::cable_cell{arb::morphology(tree), decor, {}, arb::cv_policy_max_extent((dend_length+soma_radius*2)/num_comp_ * arb::units::um)};
     }
 
     virtual cell_kind get_cell_kind(cell_gid_type) const override {
@@ -233,7 +232,7 @@ class recipe_hh_3_branches: public recipe {
         arb::decor decor;
         decor.paint(arb::reg::all(), arb::density("hh"));
 
-        return arb::cable_cell{arb::morphology(tree), decor, {}, arb::cv_policy_max_extent((dend_length*3+soma_radius*2)/num_comp_)};
+        return arb::cable_cell{arb::morphology(tree), decor, {}, arb::cv_policy_max_extent((dend_length*3+soma_radius*2)/num_comp_ * arb::units::um)};
     }
 
     virtual cell_kind get_cell_kind(cell_gid_type) const override {

test/unit/test_cv_policy.cpp

@@ -94,8 +94,8 @@ TEST(cv_policy, empty_morphology) {
     cv_policy policies[] = {
         cv_policy_fixed_per_branch(3),
         cv_policy_fixed_per_branch(3, interior_forks),
-        cv_policy_max_extent(0.234),
-        cv_policy_max_extent(0.234, interior_forks),
+        cv_policy_max_extent(0.234 * arb::units::um),
+        cv_policy_max_extent(0.234 * arb::units::um, interior_forks),
         cv_policy_single(),
         cv_policy_single(reg::all()),
         cv_policy_explicit(ls::location(0, 0))
@@ -205,26 +205,26 @@ TEST(cv_policy, max_extent) {
         {
             // Extent of 0.25 should give exact fp calculation, giving
             // 4 CVs on the root branch.
-            cv_policy pol = cv_policy_max_extent(0.25);
+            cv_policy pol = cv_policy_max_extent(0.25 * arb::units::um);
             locset expected = as_locset(L{0, 0}, L{0, 0.25}, L{0, 0.5}, L{0, 0.75}, L{0, 1});
             EXPECT_TRUE(locset_eq(cell.provider(), expected, pol.cv_boundary_points(cell)));
         }
         {
             // Same, but applied to cable (0, 0.25, 0.75) should give 2 Cvs.
-            cv_policy pol = cv_policy_max_extent(0.25, reg::cable(0, 0.25, 0.75));
+            cv_policy pol = cv_policy_max_extent(0.25 * arb::units::um, reg::cable(0, 0.25, 0.75));
             locset expected = as_locset(L{0, 0.25}, L{0, 0.5}, L{0, 0.75});
             EXPECT_TRUE(locset_eq(cell.provider(), expected, pol.cv_boundary_points(cell)));
 
         }
         {
             // Interior forks:
-            cv_policy pol = cv_policy_max_extent(0.25, interior_forks);
+            cv_policy pol = cv_policy_max_extent(0.25 * arb::units::um, interior_forks);
             locset expected = as_locset(L{0, 0}, L{0, 0.125}, L{0, 0.375}, L{0, 0.625}, L{0, 0.875}, L{0, 1});
             EXPECT_TRUE(locset_eq(cell.provider(), expected, pol.cv_boundary_points(cell)));
         }
         {
             // Interior forks but restricted to sub-cable.
-            cv_policy pol = cv_policy_max_extent(0.25, reg::cable(0, 0.25, 0.75), interior_forks);
+            cv_policy pol = cv_policy_max_extent(0.25 * arb::units::um, reg::cable(0, 0.25, 0.75), interior_forks);
             locset expected = as_locset(L{0, 0.25}, L{0, 0.375}, L{0, 0.625}, L{0, 0.75});
             EXPECT_TRUE(locset_eq(cell.provider(), expected, pol.cv_boundary_points(cell)));
 
@@ -244,7 +244,7 @@ TEST(cv_policy, max_extent) {
         {
             // Max extent of 0.6 should give two CVs on branches of length 1,
             // four CVs on branches of length 2, and seven CVs on the branch of length 4.
-            cv_policy pol = cv_policy_max_extent(0.6);
+            cv_policy pol = cv_policy_max_extent(0.6 * arb::units::um);
             mlocation_list points = thingify(pol.cv_boundary_points(cell), cell.provider());
 
             mlocation_list points_b012 = util::assign_from(util::filter(points, [](mlocation l) { return l.branch<3; }));
@@ -321,8 +321,8 @@ TEST(cv_policy, domain) {
     EXPECT_TRUE(region_eq(cell.provider(), reg1, cv_policy_single(reg1).domain()));
     EXPECT_TRUE(region_eq(cell.provider(), reg1, cv_policy_fixed_per_branch(3, reg1).domain()));
     EXPECT_TRUE(region_eq(cell.provider(), reg1, cv_policy_fixed_per_branch(3, reg1, interior_forks).domain()));
-    EXPECT_TRUE(region_eq(cell.provider(), reg1, cv_policy_max_extent(3, reg1).domain()));
-    EXPECT_TRUE(region_eq(cell.provider(), reg1, cv_policy_max_extent(3, reg1, interior_forks).domain()));
+    EXPECT_TRUE(region_eq(cell.provider(), reg1, cv_policy_max_extent(3 * arb::units::um, reg1).domain()));
+    EXPECT_TRUE(region_eq(cell.provider(), reg1, cv_policy_max_extent(3 * arb::units::um, reg1, interior_forks).domain()));
     EXPECT_TRUE(region_eq(cell.provider(), reg1, cv_policy_every_segment(reg1).domain()));
 
     EXPECT_TRUE(region_eq(cell.provider(), join(reg1, reg2), (cv_policy_single(reg1)+cv_policy_single(reg2)).domain()));

test/unit/test_diffusion.cpp

@@ -37,7 +37,7 @@ constexpr int    with_gpu = -1;
 struct linear: public ::arb::recipe {
     linear(double x, double d, double c): extent{x}, diameter{d}, cv_length{c} {
         gprop.default_parameters = arb::neuron_parameter_defaults;
-        gprop.default_parameters.discretization = arb::cv_policy_max_extent(cv_length);
+        gprop.default_parameters.discretization = arb::cv_policy_max_extent(cv_length * arb::units::um);
         // Stick morphology
         // -----x-----
         segment_tree tree;

test/unit/test_sde.cpp

@@ -260,7 +260,7 @@ class simple_sde_recipe: public simple_recipe_base {
             segment_tree tree;
             tree.append(mnpos, {i*20., 0, 0.0, 4.0}, {i*20., 0, n1*cv_size, 4.0}, 1);
             tree.append(0, {i*20., 0, ncvs*cv_size, 4.0}, 2);
-            cells_.push_back(cable_cell(morphology(tree), dec, {}, cv_policy_max_extent(cv_size)));
+            cells_.push_back(cable_cell(morphology(tree), dec, {}, cv_policy_max_extent(cv_size * arb::units::um)));
         }
     }
 
@@ -340,7 +340,7 @@ class sde_recipe: public simple_recipe_base {
             segment_tree tree;
             tree.append(mnpos, {i*20., 0, 0.0, 4.0}, {i*20., 0, n1*cv_size, 4.0}, 1);
             tree.append(0, {i*20., 0, ncvs*cv_size, 4.0}, 2);
-            cells_.push_back(cable_cell(morphology(tree), dec, labels, cv_policy_max_extent(cv_size)));
+            cells_.push_back(cable_cell(morphology(tree), dec, labels, cv_policy_max_extent(cv_size * arb::units::um)));
         }
     }
 
@@ -972,7 +972,7 @@ class sde_recipe_gpu: public simple_recipe_base {
             segment_tree tree;
             tree.append(mnpos, {i*20., 0, 0.0, 4.0}, {i*20., 0, n1*cv_size, 4.0}, 1);
             tree.append(0, {i*20., 0, ncvs*cv_size, 4.0}, 2);
-            cells_.push_back(cable_cell(morphology(tree), dec, labels, cv_policy_max_extent(cv_size)));
+            cells_.push_back(cable_cell(morphology(tree), dec, labels, cv_policy_max_extent(cv_size * arb::units::um)));
         }
     }
 

test/unit/test_serdes.cpp

@@ -134,7 +134,7 @@ struct serdes_recipe: public arb::recipe {
     std::any get_global_properties(arb::cell_kind) const override {
         auto prop = arb::cable_cell_global_properties{};
         prop.default_parameters = arb::neuron_parameter_defaults;
-        prop.default_parameters.discretization = arb::cv_policy_max_extent(1.0);
+        prop.default_parameters.discretization = arb::cv_policy_max_extent(1.0 * arb::units::um);
         return prop;
     }