Fix example, rename tool to ☕ Arbata (#30)

☕ Arbata 🎉 

Fixes #27, #28.

.github/workflows/sonata.yml

@@ -45,4 +45,4 @@ jobs:
       - name: Run unit test
         run: ./build/bin/unit
       - name: Run the example
-        run: ./build/bin/example example/simulation_config.json
+        run: ./build/bin/arbata example/simulation_config.json

.gitignore

@@ -0,0 +1,4 @@
+build/
+install-arbor/build/
+install-arbor/install/
+*.h5

CMakeLists.txt

@@ -22,4 +22,4 @@ install(TARGETS sonata-public-deps EXPORT sonata-targets)
 
 add_subdirectory(sonata)
 add_subdirectory(test)
-add_subdirectory(example)
+add_subdirectory(arbata)

arbata/CMakeLists.txt

@@ -0,0 +1,2 @@
+add_executable(arbata arbata.cpp)
+target_link_libraries(arbata PRIVATE sonata arbor::arbor arbor::arborenv ${HDF5_C_LIBRARIES})

example/example.cpp → arbata/arbata.cpp

@@ -134,8 +134,13 @@ int main(int argc, char **argv)
         std::cout << summary << "\n";
 
     }
+    catch (arb::arbor_exception& e) {
+        std::cerr << "Arbor exception caught in SONATA miniapp: " << e.what() << "\n";
+        std::cerr << "the error occured here: " << e.where << "\n";
+        return 1;
+    }
     catch (std::exception& e) {
-        std::cerr << "exception caught in SONATA miniapp: " << e.what() << "\n";
+        std::cerr << "runtime_error caught in SONATA miniapp: " << e.what() << "\n";
         return 1;
     }
     return 0;

install-arbor/scripts/build_arbor.sh

@@ -39,6 +39,8 @@ cmake_args="$cmake_args -DARB_ARCH=$arb_arch"
 cmake_args="$cmake_args -DARB_VECTORIZE=$arb_vectorize"
 cmake_args="$cmake_args -DARB_WITH_PROFILING=OFF"
 cmake_args="$cmake_args -DARB_USE_BUNDLED_LIBS=ON"
+# cmake_args="$cmake_args -DARB_BACKTRACE=ON"
+# cmake_args="$cmake_args -DCMAKE_BUILD_TYPE=debug"
 msg "ARBOR: cmake $arb_repo_path $cmake_args"
 cmake "$arb_repo_path" $cmake_args >> "$out" 2>&1
 [ $? != 0 ] && exit_on_error "see ${out}"

readme.md

@@ -1,6 +1,8 @@
-# SONATA for Arbor
+# SONATA for Arbor: ☕ Arbata
 
-### Build the example and unit tests
+This project builds ☕ Arbata: it executes Sonata models with the Arbor simulator. Point `arbata` at a ``simulation_config.json` and let us know what happened 😄.
+
+### Build the runner and unit tests
 
 #### Optional: Clone and build the arbor library in the `install-arbor` directory
 
@@ -46,8 +48,7 @@ $ python generate/gen_spikes.py
 Make sure you have an hdf5 development package installed, e.g. `libhdf5-dev`.
 ```
 $ mkdir build && cd build
-$ cmake .. -Darbor_DIR=/absolute/path/to/install-arbor/install/lib/cmake/arbor -DCMAKE_BUILD_TYPE=release
-$ make
+$ cmake .. -Darbor_DIR=/absolute/path/to/install-arbor/install/lib/cmake/arbor -DCMAKE_BUILD_TYPE=release && make -j4
 ```
 
 #### Run the unit tests
@@ -57,9 +58,36 @@ $ ./build/bin/unit
 
 #### Run the example
 ```
-$ ./build/bin/example example/simulation_config.json
+$ ./build/bin/arbata example/simulation_config.json
 ```
 ##### Expected outcome of running the example:
 * **2492** spikes generated
 * Output spikes report: `output_spikes.h5`
 * Voltage and current probe reports: `voltage_report.h5` and `current_report.h5`
+
+### Code org
+
+- `sonata/`
+  - contains the library, also libsonata.
+- `test/`
+  - Unit tests for libsonata
+- `arbata/`
+  - contains the ☕ Arbata application.
+- `example/`
+  - contains an example Sonata simulation.
+
+### Developer notes
+
+- class `sonata_recipe`: public arb::recipe, in `sonata_recipe.hpp`
+  - takes  as arg
+  - main entry point to setting up a Sonata simulation, see `example.cpp`.
+
+- class `sonata_cell` -> arb::cable_cell, in `sonata_cell.hpp`
+  - not meant for direct use.
+
+- class `model_desc`, in `data_management_lib.hpp`
+  - member of `sonata_recipe`
+  - place where most parsing takes place.
+
+- `arbata.cpp`
+  - An example user of libsonata and the end-user tool for Sonata-on-Arbor simulation.

sonata/csv_lib.cpp

@@ -150,7 +150,6 @@ std::unordered_map<std::string, variable_map> csv_node_record::dynamic_params(ty
 
 std::unordered_map<section_kind, std::vector<arb::mechanism_desc>> csv_node_record::density_mech_desc(
         type_pop_id id, std::unordered_map<std::string, variable_map> overrides) {
-    std::unordered_map<section_kind, std::vector<arb::mechanism_desc>> ret;
 
     std::unordered_map<std::string, mech_groups> density_mechs = density_params_[id];
 
@@ -166,14 +165,12 @@ std::unordered_map<section_kind, std::vector<arb::mechanism_desc>> csv_node_reco
         }
     }
 
-    // For every mech_id
-    for (auto mech: density_mechs) {
-        auto mech_id = mech.first;
-        auto mech_gp = mech.second;
-
-        mech_gp.apply_variables();
+    std::unordered_map<section_kind, std::vector<arb::mechanism_desc>> ret;
+    for (const auto& [mech_id,mech_gp]: density_mechs) {
+        auto desc = mech_gp;
+        desc.apply_variables();
 
-        for (auto mech_instance: mech_gp.mech_details) {
+        for (auto mech_instance: desc.mech_details) {
             ret[mech_instance.section].push_back(mech_instance.mech);
         }
     }

sonata/data_management_lib.cpp

@@ -262,8 +262,7 @@ void model_desc::get_connections(cell_gid_type gid, std::vector<arb::cell_connec
                     if (loc != source_maps_[source_gid].end()) {
                         if (*loc == src_rng[s]) {
                             unsigned index = loc - source_maps_[source_gid].begin();
-                            // TODO not sure about "synapse"
-                            sources.emplace_back(source_gid, std::string{"synapse@"} + std::to_string(index));
+                            sources.emplace_back(source_gid, std::string{"detector@"} + std::to_string(index));
                         } else {
                             throw sonata_exception("source maps initialized incorrectly");
                         }
@@ -284,8 +283,7 @@ void model_desc::get_connections(cell_gid_type gid, std::vector<arb::cell_connec
                     if (loc != target_maps_[gid].end()) {
                         if ((*loc).second == edges_.globalize({edge_pop_name, (cell_gid_type) t})) {
                             unsigned index = loc - target_maps_[gid].begin();
-                             // TODO not sure about "detector"
-                            targets.emplace_back(std::string{"detector@"} + std::to_string(index));
+                            targets.emplace_back(std::string{"synapse@"} + std::to_string(index));
                         } else {
                             throw sonata_exception("target maps initialized incorrectly");
                         }

sonata/include/sonata/density_mech_helper.hpp

@@ -17,14 +17,14 @@ enum section_kind {
 
 using variable_map = std::unordered_map<std::string, double>;
 
-// mech_param is a descrtiption of a mechanism `mech` and related attributes:
+// mech_param is a descrription of a mechanism `mech` and related attributes:
 // `section` is the name of the sections where `mech` is located (eg. soma, dend etc);
 // some of the parameters of `mech` aren't set yet because they aren't assigned a numerical value
 // instead, the parameters are assigned aliases, which they should get their numerical values from
-// these alieases are stored in `param_alias`
+// these aliases are stored in `param_alias`
 struct mech_params {
     section_kind section;
-    std::unordered_map<std::string, std::string> param_alias;     // Map from mechansim parameter to alias
+    std::unordered_map<std::string, std::string> param_alias;     // Map from mechanism parameter to alias
     arb::mechanism_desc mech;
 
     mech_params(std::string sec, std::unordered_map<std::string, std::string> map, arb::mechanism_desc full_mech) :

sonata/include/sonata/sonata_cell.hpp

@@ -44,7 +44,26 @@ using namespace arborio::literals;
 
 namespace sonata {
 // Generate a cell.
-arb::cable_cell dummy_cell(
+
+arb::density fix_my_mech (const arb::mechanism_catalogue& cat, const arb::mechanism_desc& desc) {
+    std::string name = desc.name();
+    std::unordered_map<std::string, double> params;
+    const auto& info = cat[name];
+    std::string sep = "/";
+    for (const auto& [param, value]: desc.values()) {
+        if (info.globals.count(param)) {
+            name += sep + param + "=" + std::to_string(value);
+            sep = ",";
+        }
+        else {
+            params[param] = value;
+        }
+    }
+    return {name, std::move(params)};
+}
+
+arb::cable_cell sonata_cell(
+        const arb::mechanism_catalogue& cat,
         arb::decor dec,
         arb::morphology morph,
         std::unordered_map<section_kind, std::vector<arb::mechanism_desc>> mechs,
@@ -59,13 +78,13 @@ arb::cable_cell dummy_cell(
     ld.set("dend", join(tagged(3), tagged(4)));
 
     for (auto mech: mechs[section_kind::soma]) {
-        dec.paint("soma"_lab, arb::density(mech));
+        dec.paint("soma"_lab, fix_my_mech(cat,mech));
     }
     for (auto mech: mechs[section_kind::dend]) {
-        dec.paint("dend"_lab, arb::density(mech));
+        dec.paint("dend"_lab, fix_my_mech(cat,mech));
     }
     for (auto mech: mechs[section_kind::axon]) {
-        dec.paint("axon"_lab, arb::density(mech));
+        dec.paint("axon"_lab, fix_my_mech(cat,mech));
     }
 
     // Add spike threshold detector at the soma.

sonata/include/sonata/sonata_recipe.hpp

@@ -40,6 +40,7 @@ using arb::time_type;
 namespace sonata {
 class sonata_recipe: public arb::recipe {
 public:
+    arb::cable_cell_global_properties gprop;
     sonata_recipe(sonata_params params):
             model_desc_(params.network.nodes,
                        params.network.edges,
@@ -52,7 +53,14 @@ class sonata_recipe: public arb::recipe {
             run_params_(params.run),
             sim_cond_(params.conditions),
             probe_info_(params.probes_info),
-            num_cells_(model_desc_.num_cells()) {}
+            num_cells_(model_desc_.num_cells()) {
+                gprop.default_parameters = arb::neuron_parameter_defaults;
+                gprop.default_parameters.axial_resistivity = 100;
+                gprop.default_parameters.temperature_K = sim_cond_.temp_c + 273.15;
+                gprop.default_parameters.init_membrane_potential = sim_cond_.v_init;
+                gprop.default_parameters.reversal_potential_method["k"] = "nernst/k";
+                gprop.default_parameters.reversal_potential_method["na"] = "nernst/na";
+            }
 
     cell_size_type num_cells() const override {
         return num_cells_;
@@ -89,12 +97,12 @@ class sonata_recipe: public arb::recipe {
                 decor.place(s.stim_loc, stim, std::string{"i_clamp"} + std::to_string(i));
             }
 
-            return dummy_cell(decor, morph, mechs, src_types, tgt_types);
+            return sonata_cell(gprop.catalogue, decor, morph, mechs, src_types, tgt_types);
         }
         else if (get_cell_kind(gid) == cell_kind::spike_source) {
             std::lock_guard<std::mutex> l(mtx_);
             std::vector<double> time_sequence = io_desc_.get_spikes(gid);
-            return arb::util::unique_any(arb::spike_source_cell{"det@0",arb::explicit_schedule(time_sequence)});
+            return arb::util::unique_any(arb::spike_source_cell{"detector@0",arb::explicit_schedule(time_sequence)});
         }
         return {};
     }
@@ -104,16 +112,6 @@ class sonata_recipe: public arb::recipe {
         return model_desc_.get_cell_kind(gid);
     }
 
-    // cell_size_type num_sources(cell_gid_type gid) const override {
-    //     std::lock_guard<std::mutex> l(mtx_);
-    //     return model_desc_.num_sources(gid);
-    // }
-
-    // cell_size_type num_targets(cell_gid_type gid) const override {
-    //     std::lock_guard<std::mutex> l(mtx_);
-    //     return model_desc_.num_targets(gid);
-    // }
-
     std::vector<arb::cell_connection> connections_on(cell_gid_type gid) const override {
         std::vector<arb::cell_connection> conns;
 
@@ -151,14 +149,6 @@ class sonata_recipe: public arb::recipe {
     }
 
     std::any get_global_properties(cell_kind k) const override {
-        arb::cable_cell_global_properties gprop;
-        gprop.default_parameters = arb::neuron_parameter_defaults;
-        gprop.default_parameters.axial_resistivity = 100;
-        gprop.default_parameters.temperature_K = sim_cond_.temp_c + 273.15;
-        gprop.default_parameters.init_membrane_potential = sim_cond_.v_init;
-        gprop.default_parameters.reversal_potential_method["k"] = "nernst/k";
-        gprop.default_parameters.reversal_potential_method["na"] = "nernst/na";
-
         return gprop;
     }
 

test/unit/test_csv.cpp

@@ -54,8 +54,7 @@ TEST(csv_node_record, constructor) {
 
     EXPECT_EQ(m0.num_branches(), m1.num_branches());
 
-    // EXPECT_TRUE(false) << m0.num_branches();
-    EXPECT_EQ(3, m0.num_branches()); //was 1. should be 2 looking at the swc/segment_tree. ???
+    EXPECT_EQ(3, m0.num_branches());
 
     EXPECT_THROW(r.morph(t2), sonata_exception);
 

test/unit/test_model_desc.cpp

@@ -172,56 +172,56 @@ TEST(model_desc, connections) {
     md.get_connections(0, conns);
     EXPECT_EQ(1, conns.size());
     EXPECT_EQ(5,conns[0].source.gid);
-    EXPECT_EQ("synapse@"+arb::cell_tag_type{"0"},conns[0].source.label.tag);
+    EXPECT_EQ("detector@"+arb::cell_tag_type{"0"},conns[0].source.label.tag);
     // EXPECT_EQ(0,conns[0].dest.gid);
-    EXPECT_EQ("detector@"+arb::cell_tag_type{"0"},conns[0].target.tag);
+    EXPECT_EQ("synapse@"+arb::cell_tag_type{"0"},conns[0].target.tag);
     EXPECT_NEAR(0.01,conns[0].weight,1e-5);
     EXPECT_NEAR(0.1,conns[0].delay,1e-5);
 
     conns.clear();
     md.get_connections(1, conns);
     EXPECT_EQ(1, conns.size());
     EXPECT_EQ(4,conns[0].source.gid);
-    EXPECT_EQ("synapse@"+arb::cell_tag_type{"0"},conns[0].source.label.tag);
+    EXPECT_EQ("detector@"+arb::cell_tag_type{"0"},conns[0].source.label.tag);
     // EXPECT_EQ(1,conns[0].dest.gid);
-    EXPECT_EQ("detector@"+arb::cell_tag_type{"0"},conns[0].target.tag);
+    EXPECT_EQ("synapse@"+arb::cell_tag_type{"0"},conns[0].target.tag);
     EXPECT_NEAR(-0.02,conns[0].weight,1e-5);
     EXPECT_NEAR(0.1,conns[0].delay,1e-5);
 
     conns.clear();
     md.get_connections(2, conns);
     EXPECT_EQ(1, conns.size());
     EXPECT_EQ(5,conns[0].source.gid);
-    EXPECT_EQ("synapse@"+arb::cell_tag_type{"0"},conns[0].source.label.tag);
+    EXPECT_EQ("detector@"+arb::cell_tag_type{"0"},conns[0].source.label.tag);
     // EXPECT_EQ(2,conns[0].dest.gid);
-    EXPECT_EQ("detector@"+arb::cell_tag_type{"0"},conns[0].target.tag);
+    EXPECT_EQ("synapse@"+arb::cell_tag_type{"0"},conns[0].target.tag);
     EXPECT_NEAR(0.01,conns[0].weight,1e-5);
     EXPECT_NEAR(0.1,conns[0].delay,1e-5);
 
     conns.clear();
     md.get_connections(3, conns);
     EXPECT_EQ(1, conns.size());
     EXPECT_EQ(1,conns[0].source.gid);
-    EXPECT_EQ("synapse@"+arb::cell_tag_type{"0"},conns[0].source.label.tag);
+    EXPECT_EQ("detector@"+arb::cell_tag_type{"0"},conns[0].source.label.tag);
     // EXPECT_EQ(3,conns[0].dest.gid);
-    EXPECT_EQ("detector@"+arb::cell_tag_type{"0"},conns[0].target.tag);
+    EXPECT_EQ("synapse@"+arb::cell_tag_type{"0"},conns[0].target.tag);
     EXPECT_NEAR(0.05,conns[0].weight,1e-5);
     EXPECT_NEAR(0.2,conns[0].delay,1e-5);
 
     conns.clear();
     md.get_connections(4, conns);
     EXPECT_EQ(2, conns.size());
     EXPECT_EQ(0,conns[0].source.gid);
-    EXPECT_EQ("synapse@"+arb::cell_tag_type{"0"},conns[0].source.label.tag);
+    EXPECT_EQ("detector@"+arb::cell_tag_type{"0"},conns[0].source.label.tag);
     // EXPECT_EQ(4,conns[0].dest.gid);
-    EXPECT_EQ("detector@"+arb::cell_tag_type{"0"},conns[0].target.tag);
+    EXPECT_EQ("synapse@"+arb::cell_tag_type{"0"},conns[0].target.tag);
     EXPECT_NEAR(0.0235,conns[0].weight,1e-5);
     EXPECT_NEAR(0.3,conns[0].delay,1e-5);
 
     EXPECT_EQ(2,conns[1].source.gid);
-    EXPECT_EQ("synapse@"+arb::cell_tag_type{"0"},conns[1].source.label.tag);
+    EXPECT_EQ("detector@"+arb::cell_tag_type{"0"},conns[1].source.label.tag);
     // EXPECT_EQ(4,conns[1].dest.gid);
-    EXPECT_EQ("detector@"+arb::cell_tag_type{"1"},conns[1].target.tag);
+    EXPECT_EQ("synapse@"+arb::cell_tag_type{"1"},conns[1].target.tag);
     EXPECT_NEAR(0.04,conns[1].weight,1e-5);
     EXPECT_NEAR(0.3,conns[1].delay,1e-5);