compartment sets 2: reports

include/bbp/sonata/compartment_sets.h

@@ -30,26 +30,45 @@ class CompartmentSets;
 struct CompartmentLocation {
   public:
     uint64_t nodeId = 0;
-    uint64_t sectionIndex = 0;
+    uint64_t sectionId = 0;
     double offset = 0.0;
 
     /// Comparator. Used to compare vectors in CompartmentSet. More idiomatic than defining a
     /// comaprator on the fly
     bool operator==(const CompartmentLocation& other) const {
-        return nodeId == other.nodeId && sectionIndex == other.sectionIndex &&
-               offset == other.offset;
+        return nodeId == other.nodeId && sectionId == other.sectionId && offset == other.offset;
     }
 
     bool operator!=(const CompartmentLocation& other) const {
         return !(*this == other);
     }
+
+    bool operator<(const CompartmentLocation& other) const {
+        if (nodeId != other.nodeId)
+            return nodeId < other.nodeId;
+        if (sectionId != other.sectionId)
+            return sectionId < other.sectionId;
+        return offset < other.offset;
+    }
+
+    bool operator>(const CompartmentLocation& other) const {
+        return other < *this;
+    }
+
+    bool operator<=(const CompartmentLocation& other) const {
+        return !(other < *this);
+    }
+
+    bool operator>=(const CompartmentLocation& other) const {
+        return !(*this < other);
+    }
 };
 
 /// Ostream << operator used by catch2 when there are problems for example
 inline std::ostream& operator<<(std::ostream& os, const CompartmentLocation& cl) {
     os << "CompartmentLocation("
        << "nodeId: " << cl.nodeId << ", "
-       << "sectionIndex: " << cl.sectionIndex << ", "
+       << "sectionId: " << cl.sectionId << ", "
        << "offset: " << cl.offset << ")";
     return os;
 }
@@ -186,4 +205,4 @@ class SONATA_API CompartmentSets
 };
 
 }  // namespace sonata
-}  // namespace bbp
+}  // namespace bbp

include/bbp/sonata/config.h

@@ -382,10 +382,11 @@ class SONATA_API SimulationConfig
      */
     struct Report {
         enum class Sections { invalid = -1, soma, axon, dend, apic, all };
-        enum class Type { invalid = -1, compartment, lfp, summation, synapse };
+        enum class Type { invalid = -1, compartment, lfp, summation, synapse, compartment_set };
         enum class Scaling { invalid = -1, none, area };
         enum class Compartments { invalid = -1, center, all };
 
+
         /// Node sets on which to report
         std::string cells;
         /// Sections on which to report. Default value: "soma"
@@ -398,6 +399,8 @@ class SONATA_API SimulationConfig
         /// For compartment type, select compartments to report.
         /// Default value: "center"(for sections: soma), "all"(for other sections)
         Compartments compartments;
+        /// Name of the compartment set (from compartment_set.json) used for generating the report.
+        std::string compartment_set;
         /// The simulation variable to access. The variables available are model dependent. For
         /// summation type, it supports multiple variables by comma separated strings. E.g. “ina”,
         /// "AdEx.V_M, v", "i_membrane, IClamp".

python/bindings.cpp

@@ -551,19 +551,23 @@ PYBIND11_MODULE(_libsonata, m) {
     py::class_<CompartmentLocation>(m, "CompartmentLocation")
         .def("__eq__", &CompartmentLocation::operator==)
         .def("__ne__", &CompartmentLocation::operator!=)
+        .def("__lt__", &CompartmentLocation::operator<)
+        .def("__le__", &CompartmentLocation::operator<=)
+        .def("__gt__", &CompartmentLocation::operator>)
+        .def("__ge__", &CompartmentLocation::operator>=)
         .def("__repr__",
              [](const CompartmentLocation& self) {
                  return fmt::format("CompartmentLocation({}, {}, {})",
                                     self.nodeId,
-                                    self.sectionIndex,
+                                    self.sectionId,
                                     self.offset);
              })
         .def("__str__",
              [](const CompartmentLocation& self) { return py::str(py::repr(py::cast(self))); })
         .def_readonly("node_id", &CompartmentLocation::nodeId, DOC_COMPARTMENTLOCATION(nodeId))
-        .def_readonly("section_index",
-                      &CompartmentLocation::sectionIndex,
-                      DOC_COMPARTMENTLOCATION(sectionIndex))
+        .def_readonly("section_id",
+                      &CompartmentLocation::sectionId,
+                      DOC_COMPARTMENTLOCATION(sectionId))
         .def_readonly("offset", &CompartmentLocation::offset, DOC_COMPARTMENTLOCATION(offset));
 
     py::class_<CompartmentSet>(m, "CompartmentSet")
@@ -632,7 +636,7 @@ PYBIND11_MODULE(_libsonata, m) {
 
     py::class_<CompartmentSets>(m, "CompartmentSets")
         .def(py::init<const std::string&>())
-        .def_static("fromFile", &CompartmentSets::fromFile, py::arg("path"))
+        .def_static("from_file", &CompartmentSets::fromFile, py::arg("path"))
         .def("__contains__",
              &CompartmentSets::contains,
              py::arg("key"),
@@ -869,6 +873,9 @@ PYBIND11_MODULE(_libsonata, m) {
         .def_readonly("cells",
                       &SimulationConfig::Report::cells,
                       DOC_SIMULATIONCONFIG(Report, cells))
+        .def_readonly("compartment_set",
+                      &SimulationConfig::Report::compartment_set,
+                      DOC_SIMULATIONCONFIG(Report, compartmentSet))
         .def_readonly("sections",
                       &SimulationConfig::Report::sections,
                       DOC_SIMULATIONCONFIG(Report, sections))
@@ -898,6 +905,7 @@ PYBIND11_MODULE(_libsonata, m) {
                       DOC_SIMULATIONCONFIG(Report, enabled));
 
     py::enum_<SimulationConfig::Report::Sections>(report, "Sections")
+        .value("invalid", SimulationConfig::Report::Sections::invalid)
         .value("soma",
                SimulationConfig::Report::Sections::soma,
                DOC_SIMULATIONCONFIG(Report, Sections, soma))
@@ -918,13 +926,15 @@ PYBIND11_MODULE(_libsonata, m) {
         .value("compartment", SimulationConfig::Report::Type::compartment)
         .value("lfp", SimulationConfig::Report::Type::lfp)
         .value("summation", SimulationConfig::Report::Type::summation)
-        .value("synapse", SimulationConfig::Report::Type::synapse);
+        .value("synapse", SimulationConfig::Report::Type::synapse)
+        .value("compartment_set", SimulationConfig::Report::Type::compartment_set);
 
     py::enum_<SimulationConfig::Report::Scaling>(report, "Scaling")
         .value("none", SimulationConfig::Report::Scaling::none)
         .value("area", SimulationConfig::Report::Scaling::area);
 
     py::enum_<SimulationConfig::Report::Compartments>(report, "Compartments")
+        .value("invalid", SimulationConfig::Report::Compartments::invalid)
         .value("center", SimulationConfig::Report::Compartments::center)
         .value("all", SimulationConfig::Report::Compartments::all);
 

python/generated/docstrings.h

@@ -383,7 +383,7 @@ static const char *__doc_bbp_sonata_NodeSets_toJSON = R"doc(Return the nodesets
 
 static const char *__doc_bbp_sonata_CompartmentLocation_nodeId = R"doc(Id of the node.)doc";
 
-static const char *__doc_bbp_sonata_CompartmentLocation_sectionIndex = R"doc(Absolute section index. Progressive index that uniquely identifies the section. There is a mapping between neuron section names (i.e. dend[10]) and this index.)doc";
+static const char *__doc_bbp_sonata_CompartmentLocation_sectionId = R"doc(Absolute section id. Progressive index that uniquely identifies the section. It is different from the relative section index. There is a mapping between neuron section names (i.e. dend[10]) and this index.)doc";
 
 static const char *__doc_bbp_sonata_CompartmentLocation_offset = R"doc(Offset of the compartment along the section)doc";
 
@@ -1156,6 +1156,8 @@ static const char *__doc_bbp_sonata_SimulationConfig_Report_Type_synapse = R"doc
 
 static const char *__doc_bbp_sonata_SimulationConfig_Report_cells = R"doc(Node sets on which to report)doc";
 
+static const char *__doc_bbp_sonata_SimulationConfig_Report_compartmentSet = R"doc(Compartment set on which to report)doc";
+
 static const char *__doc_bbp_sonata_SimulationConfig_Report_compartments =
 R"doc(For compartment type, select compartments to report. Default value:
 "center"(for sections: soma), "all"(for other sections))doc";

python/tests/test_compartment_sets.py

@@ -17,20 +17,50 @@ def setUp(self):
             "population": "pop0",
             "compartment_set": [
                 [4, 40, 0.9],
-                [4, 40, 0.9],
-                [3, 30, 0.75]
+                [4, 41, 0.9],
+                [5, 30, 0.75]
             ]
         }'''
         self.cs = CompartmentSet(self.json)
     def test_constructor_from_values(self):
         loc = self.cs[0]
         self.assertEqual(loc.node_id, 4)
-        self.assertEqual(loc.section_index, 40)
+        self.assertEqual(loc.section_id, 40)
         self.assertAlmostEqual(loc.offset, 0.9)
 
-    def test_equality(self):
-        self.assertEqual(self.cs[0], self.cs[1])
-        self.assertNotEqual(self.cs[0], self.cs[2])
+    def test_comparison_operators(self):
+        cs2 = CompartmentSet(self.json)
+
+        # For == and != use self.cs and cs2
+        self.assertTrue(self.cs[0] == cs2[0])
+        self.assertFalse(self.cs[0] != cs2[0])
+        self.assertFalse(self.cs[0] == cs2[1])
+        self.assertTrue(self.cs[0] != cs2[1])
+
+        # For other comparisons use self.cs only
+        loc1 = self.cs[0]
+        loc2 = self.cs[1]
+        loc3 = self.cs[0]
+
+        # Less than
+        self.assertTrue(loc1 < loc2 or loc2 < loc1)
+        self.assertFalse(loc1 < loc3)
+
+        # Less than or equal
+        self.assertTrue(loc1 <= loc3)
+        self.assertTrue(loc1 <= loc2 or loc2 <= loc1)
+
+        # Greater than
+        self.assertTrue(loc2 > loc1 or loc1 > loc2)
+        self.assertFalse(loc1 > loc3)
+
+        # Greater than or equal
+        self.assertTrue(loc1 >= loc3)
+        self.assertTrue(loc1 >= loc2 or loc2 >= loc1)
+
+        # Consistency check for self.cs pairs
+        for a, b in [(loc1, loc2), (loc1, loc3)]:
+            self.assertTrue((a < b) + (a == b) + (a > b) == 1)
 
     def test_repr_and_str(self):
         loc = self.cs[0]
@@ -44,8 +74,8 @@ def setUp(self):
             "compartment_set": [
                 [1, 10, 0.5],
                 [2, 20, 0.25],
-                [3, 30, 0.75],
-                [2, 20, 0.25]
+                [2, 20, 0.26],
+                [4, 20, 0.25]
             ]
         }'''
         self.cs = CompartmentSet(self.json)
@@ -64,11 +94,11 @@ def test_len_dunder(self):
 
     def test_getitem(self):
         loc = self.cs[0]
-        self.assertEqual((loc.node_id, loc.section_index, loc.offset), (1, 10, 0.5))
+        self.assertEqual((loc.node_id, loc.section_id, loc.offset), (1, 10, 0.5))
 
     def test_getitem_negative_index(self):
         loc = self.cs[-1]
-        self.assertEqual((loc.node_id, loc.section_index, loc.offset), (2, 20, 0.25))
+        self.assertEqual((loc.node_id, loc.section_id, loc.offset), (4, 20, 0.25))
 
     def test_getitem_out_of_bounds_raises(self):
         with self.assertRaises(IndexError):
@@ -78,15 +108,15 @@ def test_getitem_out_of_bounds_raises(self):
 
     def test_iterators(self):
         node_ids = [loc.node_id for loc in self.cs]
-        self.assertEqual(node_ids, [1, 2, 3, 2])
+        self.assertEqual(node_ids, [1, 2, 2, 4])
         node_ids = [loc.node_id for loc in self.cs.filtered_iter([2, 3])]
-        self.assertEqual(node_ids, [2, 3, 2])
+        self.assertEqual(node_ids, [2, 2])
         conv_to_list = list(self.cs.filtered_iter([2, 3]))
         self.assertTrue(all(isinstance(i, CompartmentLocation) for i in conv_to_list))
 
     def test_node_ids(self):
         node_ids = self.cs.node_ids()
-        self.assertEqual(node_ids, Selection([1, 2, 3]))
+        self.assertEqual(node_ids, Selection([1, 2, 4]))
 
     def test_filter_identity(self):
         filtered = self.cs.filter()
@@ -183,7 +213,7 @@ def test_toJSON_roundtrip(self):
         self.assertEqual(self.cs, cs2)
 
     def test_static_fromFile(self):
-        cs_file = CompartmentSets.fromFile(os.path.join(PATH, 'compartment_sets.json'))
+        cs_file = CompartmentSets.from_file(os.path.join(PATH, 'compartment_sets.json'))
         self.assertEqual(cs_file, self.cs)
 
     def test_repr_and_str(self):

python/tests/test_config.py

@@ -446,7 +446,7 @@ def test_basic(self):
         self.assertEqual(modifications["no_SK_E2"].section_configure, "%s.gSK_E2bar_SK_E2 = 0")
 
         self.assertEqual(self.config.list_report_names,
-                         { "axonal_comp_centers", "cell_imembrane", "compartment", "soma", "lfp" })
+                         { "axonal_comp_centers", "cell_imembrane", "compartment", "soma", "lfp", "compartment_set_v" })
 
         Report = SimulationConfig.Report
         self.assertEqual(self.config.report('soma').cells, 'Column')
@@ -468,6 +468,10 @@ def test_basic(self):
         self.assertEqual(self.config.report('cell_imembrane').type.name, 'summation')
         self.assertEqual(self.config.report('cell_imembrane').variable_name, 'i_membrane, IClamp')
         self.assertEqual(self.config.report('lfp').type, Report.Type.lfp)
+        self.assertEqual(self.config.report('compartment_set_v').type, Report.Type.compartment_set)
+        self.assertEqual(self.config.report('compartment_set_v').sections, Report.Sections.invalid)
+        self.assertEqual(self.config.report('compartment_set_v').compartments, Report.Compartments.invalid)
+        self.assertEqual(self.config.report('compartment_set_v').compartment_set, "cs0")
 
         self.assertEqual(self.config.network,
                          os.path.abspath(os.path.join(PATH, 'config/circuit_config.json')))

src/compartment_sets.cpp

@@ -151,6 +151,22 @@ class CompartmentSet {
         compartment_locations_.reserve(comp_it->size());
         for (auto&& el : *comp_it) {
             compartment_locations_.emplace_back(CompartmentSet::_parseCompartmentLocation(el));
+            if (compartment_locations_.size() >= 2) {
+                const auto& prev = compartment_locations_[compartment_locations_.size() - 2];
+                const auto& curr = compartment_locations_.back();
+                if (curr <= prev) {
+                    throw SonataError(
+                        fmt::format("CompartmentSet 'compartment_set' must be strictly sorted "
+                                    "(no duplicates). Found CompartmentLocation({}, {}, {}) before "
+                                    "CompartmentLocation({}, {}, {})",
+                                    prev.nodeId,
+                                    prev.sectionId,
+                                    prev.offset,
+                                    curr.nodeId,
+                                    curr.sectionId,
+                                    curr.offset));
+                }
+            }
         }
         compartment_locations_.shrink_to_fit();
     }
@@ -216,7 +232,7 @@ class CompartmentSet {
         j["compartment_set"] = nlohmann::json::array();
         for (const auto& elem : compartment_locations_) {
             j["compartment_set"].push_back(
-                nlohmann::json::array({elem.nodeId, elem.sectionIndex, elem.offset}));
+                nlohmann::json::array({elem.nodeId, elem.sectionId, elem.offset}));
         }
 
         return j;

src/config.cpp

@@ -83,7 +83,8 @@ NLOHMANN_JSON_SERIALIZE_ENUM(SimulationConfig::Report::Type,
                               {SimulationConfig::Report::Type::compartment, "compartment"},
                               {SimulationConfig::Report::Type::lfp, "lfp"},
                               {SimulationConfig::Report::Type::summation, "summation"},
-                              {SimulationConfig::Report::Type::synapse, "synapse"}})
+                              {SimulationConfig::Report::Type::synapse, "synapse"},
+                              {SimulationConfig::Report::Type::compartment_set, "compartment_set"}})
 NLOHMANN_JSON_SERIALIZE_ENUM(SimulationConfig::Report::Scaling,
                              {{SimulationConfig::Report::Scaling::invalid, nullptr},
                               {SimulationConfig::Report::Scaling::none, "none"},
@@ -1146,14 +1147,40 @@ class SimulationConfig::Parser
             const std::string debugStr = "report " + it.key();
 
             parseOptional(valueIt, "cells", report.cells, parseNodeSet());
-            parseOptional(valueIt, "sections", report.sections, {Report::Sections::soma});
             parseMandatory(valueIt, "type", debugStr, report.type);
-            parseOptional(valueIt, "scaling", report.scaling, {Report::Scaling::area});
+            if (report.type == Report::Type::compartment_set) {
+                parseMandatory(valueIt, "compartment_set", debugStr, report.compartment_set);
+                if (valueIt.find("sections") != valueIt.end()) {
+                    throw SonataError(
+                        "Field 'sections' is not allowed for reports of type 'compartment_set'.");
+                }
+                if (report.type == Report::Type::compartment_set &&
+                    valueIt.find("compartments") != valueIt.end()) {
+                    throw SonataError(
+                        "Field 'compartments' is not allowed for reports of type "
+                        "'compartment_set'.");
+                }
+            } else {
+                if (valueIt.find("compartment_set") != valueIt.end()) {
+                    throw SonataError(
+                        "Field 'compartment_set' is not allowed for reports of type "
+                        "'compartment_set'.");
+                }
+            }
+            parseOptional(valueIt,
+                          "sections",
+                          report.sections,
+                          {report.type == Report::Type::compartment_set ? Report::Sections::invalid
+                                                                        : Report::Sections::soma});
             parseOptional(valueIt,
                           "compartments",
                           report.compartments,
-                          {report.sections == Report::Sections::soma ? Report::Compartments::center
-                                                                     : Report::Compartments::all});
+                          {report.type == Report::Type::compartment_set
+                               ? Report::Compartments::invalid
+                               : (report.sections == Report::Sections::soma
+                                      ? Report::Compartments::center
+                                      : Report::Compartments::all)});
+            parseOptional(valueIt, "scaling", report.scaling, {Report::Scaling::area});
             parseMandatory(valueIt, "variable_name", debugStr, report.variableName);
             parseOptional(valueIt, "unit", report.unit, {"mV"});
             parseMandatory(valueIt, "dt", debugStr, report.dt);

tests/data/compartment_sets.json

@@ -4,7 +4,7 @@
       "compartment_set": [
         [0, 10, 0.1],
         [0, 10, 0.2],
-        [0, 10, 0.1],
+        [0, 10, 0.3],
         [2, 3, 0.1],
         [3, 6, 0.3]
       ]