Update to arbor/master (0.10.0)

CMakeLists.txt

@@ -3,9 +3,11 @@ cmake_minimum_required(VERSION 3.19)
 find_package(Git)
 
 project(the-arbor-gui
-  VERSION 0.8.1
+  VERSION 0.10.1
   LANGUAGES C CXX)
 
+set(EXPORT_COMPILE_COMMANDS ON)
+
 set(CMAKE_CXX_STANDARD 20)
 include(GNUInstallDirs)
 
@@ -141,6 +143,7 @@ target_include_directories(imgui PRIVATE 3rd-party/imgui)
 target_include_directories(imgui PUBLIC 3rd-party/implot)
 target_include_directories(imgui PUBLIC 3rd-party/ImGuizmo/)
 target_link_libraries(imgui PRIVATE glfw)
+target_compile_definitions(arbor-gui-deps INTERFACE SPDLOG_FMT_EXTERNAL)
 
 # Get commit hashes from git and inject into config
 execute_process(COMMAND ${GIT_EXECUTABLE} -C ${PROJECT_SOURCE_DIR} describe --always --dirty --abbrev=0

src/gui_state.cpp

@@ -23,6 +23,7 @@
 #include <arbor/simulation.hpp>
 #include <arbor/context.hpp>
 #include <arbor/load_balance.hpp>
+#include <arbor/units.hpp>
 
 #include "gui.hpp"
 #include "utils.hpp"
@@ -36,6 +37,7 @@ extern float delta_zoom;
 extern glm::vec2 mouse;
 
 using namespace std::literals;
+namespace U = arb::units;
 
 namespace {
   inline void gui_read_morphology(gui_state& state, bool& open);
@@ -852,20 +854,22 @@ namespace {
         i_clamp.frequency = item.frequency;
         i_clamp.phase     = item.phase;
         std::sort(item.envelope.begin(), item.envelope.end());
-        for (const auto& [t, i]: item.envelope) i_clamp.envelope.emplace_back(arb::i_clamp::envelope_point{t, i});
+        for (const auto& [t, i]: item.envelope) {
+          i_clamp.envelope.emplace_back(arb::i_clamp::envelope_point{t * U::ms , i * U::nA });
+        }
         decor.place(locset, i_clamp, item.tag);
       }
       for (const auto child: state.detectors.get_children(id)) {
         auto item = state.detectors[child];
-        decor.place(locset, arb::threshold_detector{item.threshold}, item.tag);
+        decor.place(locset, arb::threshold_detector{item.threshold * U::mV}, item.tag);
       }
     }
 
     auto param = state.parameter_defaults;
-    if (param.RL) decor.set_default(arb::axial_resistivity{param.RL.value()});
-    if (param.Cm) decor.set_default(arb::membrane_capacitance{param.Cm.value()});
-    if (param.TK) decor.set_default(arb::temperature_K{param.TK.value()});
-    if (param.Vm) decor.set_default(arb::init_membrane_potential{param.Vm.value()});
+    if (param.RL) decor.set_default(arb::axial_resistivity{param.RL.value() * U::Ohm * U::cm});
+    if (param.Cm) decor.set_default(arb::membrane_capacitance{param.Cm.value() * U::F / U::m2});
+    if (param.TK) decor.set_default(arb::temperature{param.TK.value() * U::Kelvin });
+    if (param.Vm) decor.set_default(arb::init_membrane_potential{param.Vm.value() * U::mV});
 
     for (const auto& ion: state.ions) {
       const auto& data = state.ion_defaults[ion];
@@ -878,31 +882,31 @@ namespace {
 
       if (state.presets.ion_data.contains(name)) {
         auto p = state.presets.ion_data.at(name);
-        decor.set_default(arb::init_int_concentration{name,  data.Xi.value_or(p.init_int_concentration.value())});
-        decor.set_default(arb::init_ext_concentration{name,  data.Xo.value_or(p.init_ext_concentration.value())});
-        decor.set_default(arb::init_reversal_potential{name, data.Er.value_or(p.init_reversal_potential.value())});
+        decor.set_default(arb::init_int_concentration{name,  data.Xi.value_or(p.init_int_concentration.value()) * U::mM});
+        decor.set_default(arb::init_ext_concentration{name,  data.Xo.value_or(p.init_ext_concentration.value()) * U::mM});
+        decor.set_default(arb::init_reversal_potential{name, data.Er.value_or(p.init_reversal_potential.value())* U::mV});
       } else {
-        decor.set_default(arb::init_int_concentration{name,  data.Xi.value()});
-        decor.set_default(arb::init_ext_concentration{name,  data.Xo.value()});
-        decor.set_default(arb::init_reversal_potential{name, data.Er.value()});
+        decor.set_default(arb::init_int_concentration{name,  data.Xi.value() * U::mM});
+        decor.set_default(arb::init_ext_concentration{name,  data.Xo.value() * U::mM});
+        decor.set_default(arb::init_reversal_potential{name, data.Er.value() * U::mV});
       }
     }
 
     for (const auto& id: state.regions) {
       auto rg = state.region_defs[id];
       if (!rg.data) continue;
       auto param  = state.parameter_defs[id];
-      if (param.RL) decor.paint(rg.data.value(), arb::axial_resistivity{param.RL.value()});
-      if (param.Cm) decor.paint(rg.data.value(), arb::membrane_capacitance{param.Cm.value()});
-      if (param.TK) decor.paint(rg.data.value(), arb::temperature_K{param.TK.value()});
-      if (param.Vm) decor.paint(rg.data.value(), arb::init_membrane_potential{param.Vm.value()});
+      if (param.RL) decor.paint(rg.data.value(), arb::axial_resistivity{param.RL.value() * U::Ohm * U::cm});
+      if (param.Cm) decor.paint(rg.data.value(), arb::membrane_capacitance{param.Cm.value() * U::F / U::m2});
+      if (param.TK) decor.paint(rg.data.value(), arb::temperature{param.TK.value() * U::Kelvin});
+      if (param.Vm) decor.paint(rg.data.value(), arb::init_membrane_potential{param.Vm.value() * U::mV});
 
       for (const auto& ion: state.ions) {
         const auto& data = state.ion_par_defs[{id, ion}];
         const auto& name = state.ion_defs[ion].name;
-        if (data.Xi) decor.paint(rg.data.value(), arb::init_int_concentration{name,  data.Xi.value()});
-        if (data.Xo) decor.paint(rg.data.value(), arb::init_ext_concentration{name,  data.Xo.value()});
-        if (data.Er) decor.paint(rg.data.value(), arb::init_reversal_potential{name, data.Er.value()});
+        if (data.Xi) decor.paint(rg.data.value(), arb::init_int_concentration{name,  data.Xi.value() * U::mM});
+        if (data.Xo) decor.paint(rg.data.value(), arb::init_ext_concentration{name,  data.Xo.value() * U::mM});
+        if (data.Er) decor.paint(rg.data.value(), arb::init_reversal_potential{name, data.Er.value() * U::mV});
       }
 
       for (const auto child: state.mechanisms.get_children(id)) {
@@ -940,7 +944,7 @@ namespace {
       if (ImGui::BeginChild("TracePlot", {-180.0f, 0.0f})) {
         if (to_plot) {
           auto probe = to_plot.value();
-          auto trace = state.sim.traces.at(probe);
+          auto trace = state.sim.traces.at(probe.value);
           const auto& [lo, hi] = std::accumulate(trace.values.begin(),
                                                  trace.values.end(),
                                                  std::make_pair(std::numeric_limits<float>::max(), std::numeric_limits<float>::min()),
@@ -1014,9 +1018,9 @@ void gui_state::deserialize(const std::filesystem::path& fn) {
   struct ls_visitor {
     gui_state* state;
     id_type locset;
-    std::string tag;
+    arb::hash_type tag;
 
-    ls_visitor(gui_state* s, const arb::locset& l, const std::string& t): state{s}, tag{t} {
+    ls_visitor(gui_state* s, const arb::locset& l, const arb::hash_type& t): state{s}, tag{t} {
       auto ls = to_string(l);
       auto res = std::find_if(state->locsets.begin(), state->locsets.end(),
                               [&](const auto& id) {
@@ -1082,33 +1086,33 @@ void gui_state::deserialize(const std::filesystem::path& fn) {
       return id;
     }
 
-    void operator()(const arb::init_membrane_potential& t) { state->parameter_defs[region].Vm = t.value.get_scalar(); }
-    void operator()(const arb::axial_resistivity& t)       { state->parameter_defs[region].RL = t.value.get_scalar(); }
-    void operator()(const arb::temperature_K& t)           { state->parameter_defs[region].TK = t.value.get_scalar(); }
-    void operator()(const arb::membrane_capacitance& t)    { state->parameter_defs[region].Cm = t.value.get_scalar(); }
+    void operator()(const arb::init_membrane_potential& t) { state->parameter_defs[region].Vm = t.value; }
+    void operator()(const arb::axial_resistivity& t)       { state->parameter_defs[region].RL = t.value; }
+    void operator()(const arb::temperature& t)             { state->parameter_defs[region].TK = t.value; }
+    void operator()(const arb::membrane_capacitance& t)    { state->parameter_defs[region].Cm = t.value; }
     void operator()(const arb::init_int_concentration& t) {
       auto ion = std::find_if(state->ions.begin(), state->ions.end(),
                               [&](const auto& id) { return state->ion_defs[id].name == t.ion; });
       if (ion == state->ions.end()) log_error("Unknown ion");
-      state->ion_par_defs[{region, *ion}].Xi = t.value.get_scalar();
+      state->ion_par_defs[{region, *ion}].Xi = t.value;
     }
     void operator()(const arb::init_ext_concentration& t) {
       auto ion = std::find_if(state->ions.begin(), state->ions.end(),
                               [&](const auto& id) { return state->ion_defs[id].name == t.ion; });
       if (ion == state->ions.end()) log_error("Unknown ion");
-      state->ion_par_defs[{region, *ion}].Xo = t.value.get_scalar();
+      state->ion_par_defs[{region, *ion}].Xo = t.value;
     }
     void operator()(const arb::init_reversal_potential& t) {
       auto ion = std::find_if(state->ions.begin(), state->ions.end(),
                               [&](const auto& id) { return state->ion_defs[id].name == t.ion; });
       if (ion == state->ions.end()) log_error("Unknown ion");
-      state->ion_par_defs[{region, *ion}].Er = t.value.get_scalar();
+      state->ion_par_defs[{region, *ion}].Er = t.value;
     }
     void operator()(const arb::ion_diffusivity& t) {
       auto ion = std::find_if(state->ions.begin(), state->ions.end(),
                               [&](const auto& id) { return state->ion_defs[id].name == t.ion; });
       if (ion == state->ions.end()) log_error("Unknown ion");
-      state->ion_par_defs[{region, *ion}].D = t.value.get_scalar();
+      state->ion_par_defs[{region, *ion}].D = t.value;
     }
     void operator()(const arb::scaled_mechanism<arb::density>& s) {
       auto id = make_density(s.t_mech);
@@ -1426,15 +1430,15 @@ void gui_state::run_simulation() {
       auto p = presets.ion_data.at(name);
       prop.add_ion(name,
                    data.charge,
-                   def.Xi.value_or(p.init_int_concentration.value()),
-                   def.Xo.value_or(p.init_ext_concentration.value()),
-                   def.Er.value_or(p.init_reversal_potential.value()));
+                   def.Xi.value_or(p.init_int_concentration.value()) * U::mM,
+                   def.Xo.value_or(p.init_ext_concentration.value()) * U::mM,
+                   def.Er.value_or(p.init_reversal_potential.value())* U::mV);
     } else {
       prop.add_ion(name,
                    data.charge,
-                   def.Xi.value(),
-                   def.Xo.value(),
-                   def.Er.value());
+                   def.Xi.value() * U::mM,
+                   def.Xo.value() * U::mM,
+                   def.Er.value() * U::mV);
     }
   }
   auto rec = make_recipe(prop, cell);
@@ -1444,33 +1448,42 @@ void gui_state::run_simulation() {
       const auto& where = locset_defs[ls];
       if (!where.data) continue;
       auto loc = where.data.value();
+      // TODO this is quite crude...
+      auto tag = std::to_string(pb.value);
       if (data.kind == "Voltage") {
-        rec.probes.emplace_back(arb::cable_probe_membrane_voltage{loc}, pb.value);
+        rec.probes.emplace_back(arb::cable_probe_membrane_voltage{loc}, tag);
       } else if (data.kind == "Axial Current") {
-        rec.probes.emplace_back(arb::cable_probe_axial_current{loc}, pb.value);
+        rec.probes.emplace_back(arb::cable_probe_axial_current{loc}, tag);
       } else if (data.kind == "Membrane Current") {
-        rec.probes.emplace_back(arb::cable_probe_total_ion_current_density{loc}, pb.value);
+        rec.probes.emplace_back(arb::cable_probe_total_ion_current_density{loc}, tag);
       }
       // TODO Finish
     }
   }
   // Make simulation
-  auto sm  = arb::simulation(rec);
+  auto sm = arb::simulation(rec);
   sim.traces.clear();
+  sim.tag_to_id.clear();
   sm.add_sampler(arb::all_probes,
-                 arb::regular_schedule(this->sim.dt),
-                 [&](arb::probe_metadata pm, std::size_t n, const arb::sample_record* samples) {
-                    auto loc = arb::util::any_cast<const arb::mlocation*>(pm.meta); 
-                    trace t{(size_t)pm.tag, pm.index, loc->pos, loc->branch, {}, {}};
+                 arb::regular_schedule(this->sim.dt * U::ms),
+                 [&](const arb::probe_metadata pm, std::size_t n, const arb::sample_record* samples) {
+                    auto loc = arb::util::any_cast<const arb::mlocation*>(pm.meta);
+                    auto tag = pm.id.tag;
+                    if (sim.tag_to_id.count(tag) == 0) {
+                      id_type id = {sim.traces.size()};
+                      sim.tag_to_id[tag] = id;
+                      sim.traces.emplace_back(tag, id, pm.index, loc->pos, loc->branch);
+                    }
+                    auto id = sim.tag_to_id[tag];
+                    auto& t = sim.traces.at(id.value);
                     for (std::size_t i = 0; i<n; ++i) {
                       const double* value = arb::util::any_cast<const double*>(samples[i].data);
                       t.times.push_back(samples[i].time);
                       t.values.push_back(*value);
                     }
-                    sim.traces[t.id] = t;                    
                   });
   try {
-    sm.run(sim.until, sim.dt);
+    sm.run(sim.until * U::ms, sim.dt * U::ms);
   } catch (...) {
       ImGui::OpenPopup("Error");
   }

src/loader.cpp

@@ -14,6 +14,7 @@ namespace io {
 
 loaded_morphology load_swc(const std::filesystem::path &fn,
                            std::function<arb::morphology(const std::vector<arborio::swc_record> &)> swc_to_morph) {
+
     return { swc_to_morph(arborio::parse_swc(slurp(fn)).records()),
              {{"soma",   "(tag 1)"},
               {"axon",   "(tag 2)"},
@@ -22,8 +23,21 @@ loaded_morphology load_swc(const std::filesystem::path &fn,
              {}};
 }
 
-loaded_morphology load_neuron_swc(const std::filesystem::path &fn) { return load_swc(fn, arborio::load_swc_neuron); }
-loaded_morphology load_arbor_swc(const std::filesystem::path &fn)  { return load_swc(fn, arborio::load_swc_arbor); }
+loaded_morphology load_neuron_swc(const std::filesystem::path &fn) {
+    auto loaded = arborio::load_swc_neuron(fn);
+    loaded_morphology res{.morph=loaded.morphology};
+    for (const auto& [k, v]: loaded.labels.regions()) res.regions.emplace_back(k, to_string(v));
+    for (const auto& [k, v]: loaded.labels.locsets()) res.locsets.emplace_back(k, to_string(v));
+    return res;
+}
+
+loaded_morphology load_arbor_swc(const std::filesystem::path &fn) {
+    auto loaded = arborio::load_swc_arbor(fn);
+    loaded_morphology res{.morph=loaded.morphology};
+    for (const auto& [k, v]: loaded.labels.regions()) res.regions.emplace_back(k, to_string(v));
+    for (const auto& [k, v]: loaded.labels.locsets()) res.locsets.emplace_back(k, to_string(v));
+    return res;
+}
 
 loaded_morphology load_neuroml_morph(const std::filesystem::path &fn) {
     arborio::neuroml nml(slurp(fn));
@@ -33,8 +47,8 @@ loaded_morphology load_neuroml_morph(const std::filesystem::path &fn) {
     if (!morph_data) log_error("Invalid morphology id {} in NML file.");
     auto morph = morph_data.value();
     loaded_morphology result{.morph=morph.morphology};
-    for (const auto& [k, v]: morph.groups.regions()) result.regions.emplace_back(k, to_string(v));
-    for (const auto& [k, v]: morph.groups.locsets()) result.locsets.emplace_back(k, to_string(v));
+    for (const auto& [k, v]: morph.labels.regions()) result.regions.emplace_back(k, to_string(v));
+    for (const auto& [k, v]: morph.labels.locsets()) result.locsets.emplace_back(k, to_string(v));
     return result;
 }
 
@@ -46,8 +60,8 @@ loaded_morphology load_neuroml_cell(const std::filesystem::path &fn) {
     if (!morph_data) log_error("Invalid cell id {} in NML file.");
     auto morph = morph_data.value();
     loaded_morphology result{.morph=morph.morphology};
-    for (const auto& [k, v]: morph.groups.regions()) result.regions.emplace_back(k, to_string(v));
-    for (const auto& [k, v]: morph.groups.locsets()) result.locsets.emplace_back(k, to_string(v));
+    for (const auto& [k, v]: morph.labels.regions()) result.regions.emplace_back(k, to_string(v));
+    for (const auto& [k, v]: morph.labels.locsets()) result.locsets.emplace_back(k, to_string(v));
     return result;
 }
 

src/loader.hpp

@@ -2,7 +2,6 @@
 
 #include <vector>
 #include <string>
-#include <tuple>
 #include <filesystem>
 
 #include <arbor/morph/morphology.hpp>

src/recipe.hpp

@@ -22,7 +22,7 @@ struct recipe: arb::recipe {
     std::any get_global_properties(arb::cell_kind) const override { return properties; }
 };
 
-
+inline
 recipe make_recipe(const arb::cable_cell_global_properties& properties,
                    const arb::cable_cell& cell) {
     recipe result;

src/simulation.hpp

@@ -5,16 +5,20 @@
 #include <string>
 
 #include "id.hpp"
-#include "utils.hpp"
 
 struct trace {
+    std::string tag;
     id_type id;
     size_t index;
     double location;
     size_t branch;
     bool show = true;
     std::vector<float> times;
     std::vector<float> values;
+
+    trace(const std::string t, const id_type i, size_t x, const double l, const size_t b):
+        tag{std::move(t)}, id{i}, index{x}, location{l}, branch{b}
+    {}
 };
 
 
@@ -25,7 +29,8 @@ struct simulation {
     bool should_run = false;
     bool show_trace = false;
 
-    std::unordered_map<id_type, trace> traces;
+    std::unordered_map<std::string, id_type> tag_to_id;
+    std::vector<trace> traces;
 };
 
 void gui_sim(simulation&);