From 4505986d5a7de14173ccb01811c57279a75f03d7 Mon Sep 17 00:00:00 2001
From: Thorsten Hater <24411438+thorstenhater@users.noreply.github.com>
Date: Tue, 23 Aug 2022 14:09:53 +0200
Subject: [PATCH] Add more visibility to decor. (#1953)

Functionality to
- delete ions from global properties
- peek into global properties
- expose painted, placed, and defaulted items
- expose voltage limit to python
- make voltage limit an optional to allow for arbitrary values.
- Add ion_settings view to python = (valence, ion_data, rev_pot_method)

Documentation
- cable cell properties in python

Closes #1578
Closes #1577
Closes #1607
---
 arbor/fvm_lowered_cell_impl.hpp          |   8 +-
 arbor/include/arbor/cable_cell_param.hpp |   4 +-
 doc/cpp/cable_cell.rst                   |   4 +-
 doc/python/cable_cell.rst                |  63 +++++++++-
 doc/python/decor.rst                     |  14 ++-
 python/cells.cpp                         | 139 +++++++++++++++++++++--
 6 files changed, 214 insertions(+), 18 deletions(-)

diff --git a/arbor/fvm_lowered_cell_impl.hpp b/arbor/fvm_lowered_cell_impl.hpp
index 9b10a0e4..bfb0a1ff 100644
--- a/arbor/fvm_lowered_cell_impl.hpp
+++ b/arbor/fvm_lowered_cell_impl.hpp
@@ -95,8 +95,8 @@ private:
     std::vector<mechanism_ptr> mechanisms_; // excludes reversal potential calculators.
     std::vector<mechanism_ptr> revpot_mechanisms_;
 
-    // Non-physical voltage check threshold, 0 => no check.
-    value_type check_voltage_mV_ = 0;
+    // Optional non-physical voltage check threshold
+    std::optional<double> check_voltage_mV_;
 
     // Flag indicating that at least one of the mechanisms implements the post_events procedure
     bool post_events_ = false;
@@ -311,9 +311,9 @@ fvm_integration_result fvm_lowered_cell_impl<Backend>::integrate(
 
         // Check for non-physical solutions:
 
-        if (check_voltage_mV_>0) {
+        if (check_voltage_mV_) {
             PE(advance:integrate:physicalcheck);
-            assert_voltage_bounded(check_voltage_mV_);
+            assert_voltage_bounded(check_voltage_mV_.value());
             PL();
         }
 
diff --git a/arbor/include/arbor/cable_cell_param.hpp b/arbor/include/arbor/cable_cell_param.hpp
index 83a7e3de..b4307d0a 100644
--- a/arbor/include/arbor/cable_cell_param.hpp
+++ b/arbor/include/arbor/cable_cell_param.hpp
@@ -323,9 +323,9 @@ ARB_ARBOR_API extern cable_cell_parameter_set neuron_parameter_defaults;
 struct ARB_SYMBOL_VISIBLE cable_cell_global_properties {
     mechanism_catalogue catalogue = global_default_catalogue();
 
-    // If >0, check membrane voltage magnitude is less than limit
+    // Optional check if membrane voltage magnitude is less than limit
     // during integration.
-    double membrane_voltage_limit_mV = 0;
+    std::optional<double> membrane_voltage_limit_mV;
 
     // True => combine linear synapses for performance.
     bool coalesce_synapses = true;
diff --git a/doc/cpp/cable_cell.rst b/doc/cpp/cable_cell.rst
index ea800af3..b44e36dc 100644
--- a/doc/cpp/cable_cell.rst
+++ b/doc/cpp/cable_cell.rst
@@ -214,9 +214,9 @@ Global properties
    by default, this is set to point to `global_default_catalogue()`, the catalogue
    that contains all mechanisms bundled with arbor.
 
-   .. cpp:member:: double membrane_voltage_limit_mv
+   .. cpp:member:: optional<double> membrane_voltage_limit_mv
 
-   if non-zero, check to see if the membrane voltage ever exceeds this value
+   if set, check to see if the membrane voltage ever exceeds this value
    in magnitude during the course of a simulation. if so, throw an exception
    and abort the simulation.
 
diff --git a/doc/python/cable_cell.rst b/doc/python/cable_cell.rst
index 20018d83..03603b98 100644
--- a/doc/python/cable_cell.rst
+++ b/doc/python/cable_cell.rst
@@ -74,7 +74,68 @@ Cable cells
         :type index: int
         :rtype: tuple(int, int)
 
-
 .. py:class:: ion
 
     properties of an ionic species.
+
+
+.. py:class:: cable_cell_global_properties
+
+   .. property:: catalogue
+
+       All mechanism names refer to mechanism instances in this mechanism
+       catalogue. by default, this is set to point to ``default_catalogue()``.
+
+   .. property:: membrane_voltage_limit
+
+       Set a limiter ``U_max`` (mV) on the membrane potential; if ``U > U_max``
+       at any point and location the simulation is aborted with an error.
+       Defaults to ``None``, if set to a numeric value the limiter is armed.
+
+   .. property:: ion_data
+
+     Return a read-only view onto concentrations, diffusivity, and reversal potential settings.
+
+  .. property:: ion_valence
+
+     Return a read-only view onto species and charge.
+
+  .. property:: ion_reversal_potential
+
+     Return a read-only view onto reversal potential methods (if set).
+
+   .. property:: ions
+
+     Return a read-only view onto all settings.
+
+   .. function:: set_ion(name,
+                         charge,
+                         internal_concentration, external_concentration,
+                         reversal_potential, reversal_potential_method,
+                         diffusivty)
+
+      Add a new ion to the global set of known species.
+
+   .. function:: unset_ion(name)
+
+      Remove the named ion.
+
+   .. property:: membrane_potential
+
+    Set the default value for the membrane potential. (``mV``)
+
+   .. property:: membrane_capacitance
+
+    Set the default value for the membrane potential. (``F/m²``)
+
+    .. property:: temperature
+
+    Set the default value for the temperature (``K``).
+
+    .. property:: axial_resisitivity
+
+    Set the default value for the membrane axial resisitivity. (``Ω·cm``)
+
+
+For convenience, ``neuron_cable_properties`` is a predefined value that holds
+values that correspond to NEURON defaults.
diff --git a/doc/python/decor.rst b/doc/python/decor.rst
index 6eb79d24..f8fbc43e 100644
--- a/doc/python/decor.rst
+++ b/doc/python/decor.rst
@@ -172,4 +172,16 @@ Cable cell decoration
 
         Set the cv_policy used to discretise the cell into control volumes for simulation.
 
-        :param str policy: :ref:`string representation <morph-cv-sexpr>` of a cv_policy.
\ No newline at end of file
+        :param str policy: :ref:`string representation <morph-cv-sexpr>` of a cv_policy.
+
+    .. method:: paintings()
+
+        Returns a list of tuples ``(region, painted_object)`` for inspection.
+
+    .. method:: placements()
+
+        Returns a list of tuples ``(locset, placed_object)`` for inspection.
+
+    .. method:: defaults()
+
+        Returns a list of all set defaults for inspection.
diff --git a/python/cells.cpp b/python/cells.cpp
index ca4716f2..9f1c1237 100644
--- a/python/cells.cpp
+++ b/python/cells.cpp
@@ -34,14 +34,23 @@
 #include "error.hpp"
 #include "proxy.hpp"
 #include "pybind11/cast.h"
+#include "pybind11/stl.h"
 #include "pybind11/pytypes.h"
 #include "schedule.hpp"
 #include "strprintf.hpp"
 
 namespace pyarb {
 
+template <typename T>
+std::string to_string(const T& t) {
+    std::stringstream ss;
+    ss << t;
+    return ss.str();
+}
+
 // This isn't pretty. Partly because the information in the global parameters
 // is all over the place.
+template <>
 std::string to_string(const arb::cable_cell_global_properties& props) {
     std::string s = "{arbor.cable_global_properties";
 
@@ -417,6 +426,46 @@ void register_cells(pybind11::module& m) {
           "`proportion` is the proportion of the CV (itegrated by area or length) included in the region."
     );
 
+    pybind11::class_<arb::init_membrane_potential> membrane_potential(m, "membrane_potential", "Setting the initial membrane voltage.");
+    membrane_potential
+        .def(pybind11::init([](double v) -> arb::init_membrane_potential { return {v}; }))
+        .def("__repr__", [](const arb::init_membrane_potential& d){return "Vm=" + std::to_string(d.value);});
+
+    pybind11::class_<arb::membrane_capacitance> membrane_capacitance(m, "membrane_capacitance", "Setting the membrane capacitance.");
+    membrane_capacitance
+        .def(pybind11::init([](double v) -> arb::membrane_capacitance { return {v}; }))
+        .def("__repr__", [](const arb::membrane_capacitance& d){return "Cm=" + std::to_string(d.value);});
+
+    pybind11::class_<arb::temperature_K> temperature_K(m, "temperature_K", "Setting the temperature.");
+    temperature_K
+        .def(pybind11::init([](double v) -> arb::temperature_K { return {v}; }))
+        .def("__repr__", [](const arb::temperature_K& d){return "T=" + std::to_string(d.value);});
+
+    pybind11::class_<arb::axial_resistivity> axial_resistivity(m, "axial_resistivity", "Setting the axial resistivity.");
+    axial_resistivity
+        .def(pybind11::init([](double v) -> arb::axial_resistivity { return {v}; }))
+        .def("__repr__", [](const arb::axial_resistivity& d){return "Ra" + std::to_string(d.value);});
+
+    pybind11::class_<arb::init_reversal_potential> reversal_potential(m, "reversal_potential", "Setting the initial reversal potential.");
+    reversal_potential
+        .def(pybind11::init([](const std::string& i, double v) -> arb::init_reversal_potential { return {i, v}; }))
+        .def("__repr__", [](const arb::init_reversal_potential& d){return "e" + d.ion + "=" + std::to_string(d.value);});
+
+    pybind11::class_<arb::init_int_concentration> int_concentration(m, "int_concentration", "Setting the initial internal ion concentration.");
+    int_concentration
+        .def(pybind11::init([](const std::string& i, double v) -> arb::init_int_concentration { return {i, v}; }))
+        .def("__repr__", [](const arb::init_int_concentration& d){return d.ion + "i" + "=" + std::to_string(d.value);});
+
+    pybind11::class_<arb::init_ext_concentration> ext_concentration(m, "ext_concentration", "Setting the initial external ion concentration.");
+    ext_concentration
+        .def(pybind11::init([](const std::string& i, double v) -> arb::init_ext_concentration { return {i, v}; }))
+        .def("__repr__", [](const arb::init_ext_concentration& d){return d.ion + "o" + "=" + std::to_string(d.value);});
+
+    pybind11::class_<arb::ion_diffusivity> ion_diffusivity(m, "ion_diffusivity", "Setting the ion diffusivity.");
+    ion_diffusivity
+        .def(pybind11::init([](const std::string& i, double v) -> arb::ion_diffusivity { return {i, v}; }))
+        .def("__repr__", [](const arb::ion_diffusivity& d){return "D" + d.ion + "=" + std::to_string(d.value);});
+
     // arb::density
 
     pybind11::class_<arb::density> density(m, "density", "For painting a density mechanism on a region.");
@@ -544,10 +593,29 @@ void register_cells(pybind11::module& m) {
     // arb::cable_cell_global_properties
     pybind11::class_<arb::cable_cell_ion_data> ion_data(m, "ion_data");
     ion_data
-        .def_readonly("internal_concentration", &arb::cable_cell_ion_data::init_int_concentration, "Internal concentration.")
-        .def_readonly("external_concentration", &arb::cable_cell_ion_data::init_ext_concentration, "External concentration.")
+        .def_readonly("internal_concentration", &arb::cable_cell_ion_data::init_int_concentration,  "Internal concentration.")
+        .def_readonly("external_concentration", &arb::cable_cell_ion_data::init_ext_concentration,  "External concentration.")
+        .def_readonly("diffusivity",            &arb::cable_cell_ion_data::diffusivity,             "Diffusivity.")
         .def_readonly("reversal_concentration", &arb::cable_cell_ion_data::init_reversal_potential, "Reversal potential.");
 
+    struct ion_settings {
+        int charge = 0;
+        std::optional<double> internal_concentration;
+        std::optional<double> external_concentration;
+        std::optional<double> diffusivity;
+        std::optional<double> reversal_potential;
+        std::string reversal_potential_method = "const";
+    };
+
+    pybind11::class_<ion_settings> py_ion_data(m, "ion_settings");
+    ion_data
+        .def_property_readonly("charge",                    [](const ion_settings& s) { return s.charge; },                    "Valence.")
+        .def_property_readonly("internal_concentration",    [](const ion_settings& s) { return s.internal_concentration; },    "Internal concentration.")
+        .def_property_readonly("external_concentration",    [](const ion_settings& s) { return s.external_concentration; },    "External concentration.")
+        .def_property_readonly("diffusivity",               [](const ion_settings& s) { return s.diffusivity; },               "Diffusivity.")
+        .def_property_readonly("reversal_potential",        [](const ion_settings& s) { return s.reversal_potential; },        "Reversal potential.")
+        .def_property_readonly("reversal_potential_method", [](const ion_settings& s) { return s.reversal_potential_method; }, "Reversal potential method.");
+
     pybind11::class_<arb::cable_cell_global_properties> gprop(m, "cable_global_properties");
     gprop
         .def(pybind11::init<>())
@@ -559,6 +627,9 @@ void register_cells(pybind11::module& m) {
         .def_property("membrane_potential",
                       [](const arb::cable_cell_global_properties& props) { return props.default_parameters.init_membrane_potential; },
                       [](arb::cable_cell_global_properties& props, double u) { props.default_parameters.init_membrane_potential = u; })
+        .def_property("membrane_voltage_limit",
+                      [](const arb::cable_cell_global_properties& props) { return props.membrane_voltage_limit_mV; },
+                      [](arb::cable_cell_global_properties& props, std::optional<double> u) { props.membrane_voltage_limit_mV = u; })
         .def_property("membrane_capacitance",
                       [](const arb::cable_cell_global_properties& props) { return props.default_parameters.membrane_capacitance; },
                       [](arb::cable_cell_global_properties& props, double u) { props.default_parameters.membrane_capacitance = u; })
@@ -568,12 +639,6 @@ void register_cells(pybind11::module& m) {
         .def_property("axial_resisitivity",
                       [](const arb::cable_cell_global_properties& props) { return props.default_parameters.axial_resistivity; },
                       [](arb::cable_cell_global_properties& props, double u) { props.default_parameters.axial_resistivity = u; })
-        .def_property_readonly("ion_data",
-                      [](const arb::cable_cell_global_properties& props) { return props.default_parameters.ion_data; })
-        .def_property_readonly("ion_valence",
-                      [](const arb::cable_cell_global_properties& props) { return props.ion_species; })
-        .def_property_readonly("ion_reversal_potential",
-                      [](const arb::cable_cell_global_properties& props) { return props.default_parameters.reversal_potential_method; })
         .def("set_property",
             [](arb::cable_cell_global_properties& props,
                optional<double> Vm, optional<double> cm,
@@ -590,6 +655,13 @@ void register_cells(pybind11::module& m) {
             pybind11::arg_v("tempK", pybind11::none(), "temperature [Kelvin]."),
             "Set global default values for cable and cell properties.")
         // add/modify ion species
+        .def("unset_ion",
+             [](arb::cable_cell_global_properties& props, const char* ion) {
+                 props.ion_species.erase(ion);
+                 props.default_parameters.ion_data.erase(ion);
+                 props.default_parameters.reversal_potential_method.erase(ion);
+             },
+             "Remove ion species from properties.")
         .def("set_ion",
             [](arb::cable_cell_global_properties& props, const char* ion,
                optional<double> valence, optional<double> int_con,
@@ -627,6 +699,34 @@ void register_cells(pybind11::module& m) {
             "specific regions using the paint interface, while the method for calculating\n"
             "reversal potential is global for all compartments in the cell, and can't be\n"
              "overriden locally.")
+        .def_property_readonly("ion_data",
+                      [](const arb::cable_cell_global_properties& props) { return props.default_parameters.ion_data; })
+        .def_property_readonly("ion_valence",
+                      [](const arb::cable_cell_global_properties& props) { return props.ion_species; })
+        .def_property_readonly("ion_reversal_potential",
+                      [](const arb::cable_cell_global_properties& props) { return props.default_parameters.reversal_potential_method; })
+        .def_property_readonly("ions",
+                               [](arb::cable_cell_global_properties& g) {
+                                   std::unordered_map<std::string, ion_settings> result;
+                                   for (const auto& [k, v]: g.ion_species) {
+                                       auto& ion = result[k];
+                                       ion.charge = v;
+                                       auto& data = g.default_parameters.ion_data;
+                                       if (data.count(k)) {
+                                           auto& i = data.at(k);
+                                           ion.diffusivity = i.diffusivity;
+                                           ion.external_concentration = i.init_ext_concentration;
+                                           ion.internal_concentration = i.init_int_concentration;
+                                           ion.reversal_potential     = i.init_reversal_potential;
+                                       }
+                                       auto& revpot = g.default_parameters.reversal_potential_method;
+                                       if (revpot.count(k)) {
+                                           ion.reversal_potential_method = revpot.at(k).name();
+                                       }
+                                   }
+                                   return result;
+                               },
+                               "Return a view of all ion settings.")
         .def_readwrite("catalogue",
                        &arb::cable_cell_global_properties::catalogue,
                        "The mechanism catalogue.")
@@ -689,6 +789,29 @@ void register_cells(pybind11::module& m) {
             "potential can be overridden on specific regions using the paint interface, \n"
             "while the method for calculating reversal potential is global for all\n"
             "compartments in the cell, and can't be overriden locally.")
+        .def("paintings",
+            [](arb::decor& dec) {
+                std::vector<std::tuple<std::string, arb::paintable>> result;
+                for (const auto& [k, v]: dec.paintings()) {
+                    result.emplace_back(to_string(k), v);
+                }
+                return result;
+            },
+            "Return a view of all painted items.")
+        .def("placements",
+            [](arb::decor& dec) {
+                std::vector<std::tuple<std::string, arb::placeable, std::string>> result;
+                for (const auto& [k, v, t]: dec.placements()) {
+                    result.emplace_back(to_string(k), v, t);
+                }
+                return result;
+            },
+            "Return a view of all placed items.")
+        .def("defaults",
+            [](arb::decor& dec) {
+                return dec.defaults().serialize();
+            },
+            "Return a view of all defaults.")
         // Paint mechanisms.
         .def("paint",
             [](arb::decor& dec, const char* region, const arb::density& mechanism) {
-- 
GitLab