diff --git a/arbor/cable_cell.cpp b/arbor/cable_cell.cpp index 617a7c6e455e4f5eec56c0f81337ea11c640a5de..35268b7e234a188aef7784d5f47215154eefa9d7 100644 --- a/arbor/cable_cell.cpp +++ b/arbor/cable_cell.cpp @@ -44,6 +44,9 @@ struct cable_cell_impl { // Track number of point assignments by type for lid/target numbers. dynamic_typed_map<constant_type<cell_lid_type>::type> placed_count; + // The label dictionary. + const label_dict dictionary; + // The decorations on the cell. decor decorations; @@ -52,6 +55,7 @@ struct cable_cell_impl { cable_cell_impl(const arb::morphology& m, const label_dict& labels, const decor& decorations): provider(m, labels), + dictionary(labels), decorations(decorations) { init(decorations); @@ -169,6 +173,10 @@ cable_cell::cable_cell(const cable_cell& other): impl_(make_impl(new cable_cell_impl(*other.impl_))) {} +const label_dict& cable_cell::labels() const { + return impl_->dictionary; +} + const concrete_embedding& cable_cell::embedding() const { return impl_->provider.embedding(); } diff --git a/arbor/cable_cell_param.cpp b/arbor/cable_cell_param.cpp index 0868f410e321e90e93665ffbab7732e5f45cf923..c8666e451b863b3e56150b00f0f4c5d264670264 100644 --- a/arbor/cable_cell_param.cpp +++ b/arbor/cable_cell_param.cpp @@ -6,9 +6,9 @@ #include <arbor/cable_cell.hpp> #include <arbor/cable_cell_param.hpp> +#include <arbor/s_expr.hpp> #include "util/maputil.hpp" -#include "s_expr.hpp" namespace arb { diff --git a/arbor/include/arbor/cable_cell.hpp b/arbor/include/arbor/cable_cell.hpp index 3185cf0f240f86cba756432da3af1b6e68555273..cf267d79738c0c379161694c69a1e6adb9208acc 100644 --- a/arbor/include/arbor/cable_cell.hpp +++ b/arbor/include/arbor/cable_cell.hpp @@ -257,6 +257,9 @@ public: cable_cell(m, {}, {}) {} + /// Access to labels + const label_dict& labels() const; + /// Access to morphology and embedding const concrete_embedding& embedding() const; const arb::morphology& morphology() const; diff --git a/arbor/include/arbor/cable_cell_param.hpp b/arbor/include/arbor/cable_cell_param.hpp index e7f5de4cf1cd3a69e83f430d5e9b8b67232eda0b..c21649a6a5e772dfdbc4e7b76f1cc76649931bf6 100644 --- a/arbor/include/arbor/cable_cell_param.hpp +++ b/arbor/include/arbor/cable_cell_param.hpp @@ -41,8 +41,8 @@ struct cable_cell_ion_data { // * The time points must be monotonically increasing. // * Onset and initial amplitude is given by the first point. // * The amplitude for time after the last time point is that of the last -// amplitgude point; an explicit zero amplitude point must be provided if the -// envelope is intended to have finite support. +// amplitude point; an explicit zero amplitude point must be provided if the +// envelope is intended to have finite support. // // Periodic envelopes are not supported, but may well be a feature worth // considering in the future. @@ -195,14 +195,14 @@ struct ion_reversal_potential_method { }; using paintable = - std::variant<mechanism_desc, - init_membrane_potential, + std::variant<init_membrane_potential, axial_resistivity, temperature_K, membrane_capacitance, init_int_concentration, init_ext_concentration, - init_reversal_potential>; + init_reversal_potential, + mechanism_desc>; using placeable = std::variant<mechanism_desc, diff --git a/arbor/include/arbor/morph/label_parse.hpp b/arbor/include/arbor/morph/label_parse.hpp index ecc82d533d0b33333cd2d1e01a723529b3ae4e6a..db3b21c18e0930645b7e0a1147b8259390b27f74 100644 --- a/arbor/include/arbor/morph/label_parse.hpp +++ b/arbor/include/arbor/morph/label_parse.hpp @@ -2,8 +2,9 @@ #include <any> -#include <arbor/morph/region.hpp> #include <arbor/arbexcept.hpp> +#include <arbor/morph/region.hpp> +#include <arbor/s_expr.hpp> #include <arbor/util/expected.hpp> namespace arb { @@ -16,6 +17,7 @@ template <typename T> using parse_hopefully = arb::util::expected<T, label_parse_error>; parse_hopefully<std::any> parse_label_expression(const std::string&); +parse_hopefully<std::any> parse_label_expression(const s_expr&); parse_hopefully<arb::region> parse_region_expression(const std::string& s); parse_hopefully<arb::locset> parse_locset_expression(const std::string& s); diff --git a/arbor/s_expr.hpp b/arbor/include/arbor/s_expr.hpp similarity index 87% rename from arbor/s_expr.hpp rename to arbor/include/arbor/s_expr.hpp index c3d1309076b6204a92ac6e25482d909815d510e7..e6ac0d5156f279605dc833e4265ca602b47afc37 100644 --- a/arbor/s_expr.hpp +++ b/arbor/include/arbor/s_expr.hpp @@ -48,6 +48,15 @@ struct token { std::ostream& operator<<(std::ostream&, const token&); +struct symbol { + std::string str; + operator std::string() const { return str; } +}; + +inline symbol operator"" _symbol(const char* chars, size_t size) { + return {chars}; +} + struct s_expr { template <typename U> struct s_pair { @@ -112,7 +121,9 @@ struct s_expr { s_expr_iterator_impl(reference e): inner_(&e) { - if (inner_->is_atom()) { + // We can't iterate over an atom, unless the atom is + // nil, which is both an atom and an empty list. + if (inner_->is_atom() && inner_->atom().kind!=tok::nil) { throw std::runtime_error("Attempt to create s_expr_iterator on an atom."); } if (finished()) inner_ = nullptr; @@ -193,7 +204,7 @@ struct s_expr { // with a std::unique_ptr via value_wrapper. using pair_type = s_pair<value_wrapper<s_expr>>; - std::variant<token, pair_type> state = token{{0,0}, tok::nil, "nil"}; + std::variant<token, pair_type> state = token{{0,0}, tok::nil, "()"}; s_expr(const s_expr& s): state(s.state) {} s_expr() = default; @@ -202,6 +213,17 @@ struct s_expr { state(pair_type(std::move(l), std::move(r))) {} + explicit s_expr(std::string s): + s_expr(token{{0,0}, tok::string, std::move(s)}) {} + explicit s_expr(const char* s): + s_expr(token{{0,0}, tok::string, s}) {} + s_expr(double x): + s_expr(token{{0,0}, tok::real, std::to_string(x)}) {} + s_expr(int x): + s_expr(token{{0,0}, tok::integer, std::to_string(x)}) {} + s_expr(symbol s): + s_expr(token{{0,0}, tok::symbol, s}) {} + bool is_atom() const; const token& atom() const; @@ -223,10 +245,13 @@ struct s_expr { friend std::ostream& operator<<(std::ostream& o, const s_expr& x); }; -std::size_t length(const s_expr& l); -src_location location(const s_expr& l); - +// Build s-expr from string s_expr parse_s_expr(const std::string& line); +// Length of the s-expr +std::size_t length(const s_expr& l); + +// Location of the head of the s-expr +src_location location(const s_expr& l); } // namespace arb diff --git a/arbor/morph/label_parse.cpp b/arbor/morph/label_parse.cpp index a40bb160a89fe698260cb44d3d96f9d0a175a91c..855e6e7a14ac0ee0986c2c985d5f301774241481 100644 --- a/arbor/morph/label_parse.cpp +++ b/arbor/morph/label_parse.cpp @@ -2,14 +2,12 @@ #include <limits> #include <arbor/arbexcept.hpp> -#include <arbor/morph/region.hpp> -#include <arbor/morph/locset.hpp> #include <arbor/morph/label_parse.hpp> -#include <arbor/morph/region.hpp> #include <arbor/morph/locset.hpp> +#include <arbor/morph/region.hpp> +#include <arbor/s_expr.hpp> +#include <arbor/util/expected.hpp> -#include "arbor/util/expected.hpp" -#include "s_expr.hpp" #include "util/strprintf.hpp" namespace arb { @@ -351,7 +349,7 @@ parse_hopefully<std::any> eval(const s_expr& e) { // This must be a function evaluation, where head is the function name, and // tail is a list of arguments. - // Evaluate the arguments, and return error state if an error ocurred. + // Evaluate the arguments, and return error state if an error occurred. auto args = eval_args(e.tail()); if (!args) { return args.error(); @@ -387,6 +385,9 @@ parse_hopefully<std::any> eval(const s_expr& e) { parse_hopefully<std::any> parse_label_expression(const std::string& e) { return eval(parse_s_expr(e)); } +parse_hopefully<std::any> parse_label_expression(const s_expr& s) { + return eval(s); +} parse_hopefully<arb::region> parse_region_expression(const std::string& s) { if (auto e = eval(parse_s_expr(s))) { diff --git a/arbor/morph/region.cpp b/arbor/morph/region.cpp index 487360fdc6bec2aa70063cd27b6b107b4a5e61d5..64234f554fb3cd00cd36fe7394b4a9bb247c6457 100644 --- a/arbor/morph/region.cpp +++ b/arbor/morph/region.cpp @@ -10,8 +10,8 @@ #include <arbor/morph/morphexcept.hpp> #include <arbor/morph/mprovider.hpp> #include <arbor/morph/region.hpp> +#include <arbor/s_expr.hpp> -#include "s_expr.hpp" #include "util/mergeview.hpp" #include "util/span.hpp" #include "util/strprintf.hpp" diff --git a/arbor/s_expr.cpp b/arbor/s_expr.cpp index f916700b74649225d6ef4b30d8a3028d82afc3d5..d6e0ef58b3681946dfd046d2470686c881cdfc32 100644 --- a/arbor/s_expr.cpp +++ b/arbor/s_expr.cpp @@ -8,9 +8,10 @@ #include <vector> #include <arbor/arbexcept.hpp> +#include <arbor/s_expr.hpp> #include "util/strprintf.hpp" -#include "s_expr.hpp" + namespace arb { inline bool is_alphanumeric(char c) { @@ -390,19 +391,38 @@ s_expr::operator bool() const { return !(is_atom() && atom().kind==tok::nil); } -std::ostream& operator<<(std::ostream& o, const s_expr& x) { - if (x.is_atom()) return o << x.atom(); -#if 1 +// Assume that stream indented and ready to go at location to start printing. +std::ostream& print(std::ostream& o, const s_expr& x, int indent) { + std::string in(std::string::size_type(2*indent), ' '); + if (x.is_atom()) { + return o << x.atom(); + } + auto it = std::begin(x); + auto end = std::end(x); + bool first=true; o << "("; - bool first = true; - for (auto& e: x) { - o << (first? "": " ") << e; + while (it!=end) { + if (!first && !it->is_atom() && length(*it)>=0) { + o << "\n" << in; + print(o, *it, indent+1); + ++it; + if (it!=end && it->is_atom()) { + o << "\n" << in; + } + } + else { + print(o, *it, indent+1); + if (++it!=end) { + o << " "; + } + } first = false; } return o << ")"; -#else - return o << "(" << x.head() << " . " << x.tail() << ")"; -#endif +} + +std::ostream& operator<<(std::ostream& o, const s_expr& x) { + return print(o, x, 1); } std::size_t length(const s_expr& l) { @@ -496,5 +516,4 @@ s_expr parse_s_expr(const std::string& line) { } return result; } - } // namespace arb diff --git a/arborio/CMakeLists.txt b/arborio/CMakeLists.txt index 14522ad6253fad35e13f61b5715999e2644f2a3a..63855f42a7029aa50681f925d4d6ba2a775564b0 100644 --- a/arborio/CMakeLists.txt +++ b/arborio/CMakeLists.txt @@ -2,6 +2,7 @@ set(arborio-sources asc_lexer.cpp neurolucida.cpp swcio.cpp + cableio.cpp ) if(ARB_WITH_NEUROML) list(APPEND arborio-sources @@ -17,10 +18,14 @@ endif() add_library(arborio ${arborio-sources}) add_library(arborio-public-headers INTERFACE) +add_library(arborio-private-headers INTERFACE) + target_include_directories(arborio-public-headers INTERFACE $<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}/include> $<INSTALL_INTERFACE:include> ) +target_include_directories(arborio-private-headers INTERFACE + "$<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}>") if(ARB_WITH_NEUROML) target_link_libraries(arborio PUBLIC arbor arborio-public-headers LibXml2::LibXml2) diff --git a/arborio/cableio.cpp b/arborio/cableio.cpp new file mode 100644 index 0000000000000000000000000000000000000000..9a7ff5af09f4abb9cad353791fc32a3b66ef6c13 --- /dev/null +++ b/arborio/cableio.cpp @@ -0,0 +1,860 @@ +#include <algorithm> +#include <sstream> +#include <unordered_map> + +#include <arbor/morph/label_parse.hpp> +#include <arbor/s_expr.hpp> +#include <arbor/util/pp_util.hpp> +#include <arbor/util/any_visitor.hpp> + +#include <arborio/cableio.hpp> + +#include "parse_s_expr.hpp" + +namespace arborio { + +using namespace arb; + +std::string acc_version() {return "0.1-dev";} + +cableio_parse_error::cableio_parse_error(const std::string& msg, const arb::src_location& loc): + arb::arbor_exception(msg+" at :"+ + std::to_string(loc.line)+ ":"+ + std::to_string(loc.column)) +{} + +cableio_morphology_error::cableio_morphology_error(const unsigned bid): + arb::arbor_exception("Invalid morphology: branch `" +std::to_string(bid) + + "` only has one child branch, making it an invalid branch specification") +{} + +cableio_version_error::cableio_version_error(const std::string& version): + arb::arbor_exception("Unsupported cable-cell format version `" + version + "`") +{} + +struct nil_tag {}; + +// Define s-expr makers for various types +s_expr mksexp(const init_membrane_potential& p) { + return slist("membrane-potential"_symbol, p.value); +} +s_expr mksexp(const axial_resistivity& r) { + return slist("axial-resistivity"_symbol, r.value); +} +s_expr mksexp(const temperature_K& t) { + return slist("temperature-kelvin"_symbol, t.value); +} +s_expr mksexp(const membrane_capacitance& c) { + return slist("membrane-capacitance"_symbol, c.value); +} +s_expr mksexp(const init_int_concentration& c) { + return slist("ion-internal-concentration"_symbol, s_expr(c.ion), c.value); +} +s_expr mksexp(const init_ext_concentration& c) { + return slist("ion-external-concentration"_symbol, s_expr(c.ion), c.value); +} +s_expr mksexp(const init_reversal_potential& c) { + return slist("ion-reversal-potential"_symbol, s_expr(c.ion), c.value); +} +s_expr mksexp(const i_clamp& c) { + std::vector<s_expr> evlps; + std::transform(c.envelope.begin(), c.envelope.end(), std::back_inserter(evlps), + [](const auto& x){return slist(x.t, x.amplitude);}); + auto envelope = slist("envelope"_symbol, slist_range(evlps)); + return slist("current-clamp"_symbol, envelope, c.frequency); +} +s_expr mksexp(const threshold_detector& d) { + return slist("threshold-detector"_symbol, d.threshold); +} +s_expr mksexp(const gap_junction_site& s) { + return slist("gap-junction-site"_symbol); +} +s_expr mksexp(const mechanism_desc& d) { + std::vector<s_expr> mech; + mech.push_back(s_expr(d.name())); + std::transform(d.values().begin(), d.values().end(), std::back_inserter(mech), + [](const auto& x){return slist(s_expr(x.first), x.second);}); + return s_expr{"mechanism"_symbol, slist_range(mech)}; +} +s_expr mksexp(const ion_reversal_potential_method& e) { + return slist("ion-reversal-potential-method"_symbol, s_expr(e.ion), mksexp(e.method)); +} +s_expr mksexp(const mpoint& p) { + return slist("point"_symbol, p.x, p.y, p.z, p.radius); +} +s_expr mksexp(const msegment& seg) { + return slist("segment"_symbol, (int)seg.id, mksexp(seg.prox), mksexp(seg.dist), seg.tag); +} +// This can be removed once cv_policy is removed from the decor. +s_expr mksexp(const cv_policy& c) { + return s_expr(); +} +s_expr mksexp(const decor& d) { + auto round_trip = [](auto& x) { + std::stringstream s; + s << x; + return parse_s_expr(s.str()); + }; + std::vector<s_expr> decorations; + for (const auto& p: d.defaults().serialize()) { + decorations.push_back(std::visit([&](auto& x) + { return slist("default"_symbol, mksexp(x)); }, p)); + } + for (const auto& p: d.paintings()) { + decorations.push_back(std::visit([&](auto& x) + { return slist("paint"_symbol, round_trip(p.first), mksexp(x)); }, p.second)); + } + for (const auto& p: d.placements()) { + decorations.push_back(std::visit([&](auto& x) + { return slist("place"_symbol, round_trip(p.first), mksexp(x)); }, p.second)); + } + return {"decor"_symbol, slist_range(decorations)}; +} +s_expr mksexp(const label_dict& dict) { + auto round_trip = [](auto& x) { + std::stringstream s; + s << x; + return parse_s_expr(s.str()); + }; + auto defs = slist(); + for (auto& r: dict.locsets()) { + defs = s_expr(slist("locset-def"_symbol, s_expr(r.first), round_trip(r.second)), std::move(defs)); + } + for (auto& r: dict.regions()) { + defs = s_expr(slist("region-def"_symbol, s_expr(r.first), round_trip(r.second)), std::move(defs)); + } + return {"label-dict"_symbol, std::move(defs)}; +} +s_expr mksexp(const morphology& morph) { + // s-expression representation of branch i in the morphology + auto make_branch = [&morph](int i) { + std::vector<s_expr> segments; + const auto& segs = morph.branch_segments(i); + std::transform(segs.begin(), segs.end(), std::back_inserter(segments), + [](const auto& x){return mksexp(x);}); + return s_expr{"branch"_symbol, {i, {(int)morph.branch_parent(i), slist_range(segments)}}}; + }; + std::vector<s_expr> branches; + for (msize_t i = 0; i < morph.num_branches(); ++i) { + branches.push_back(make_branch(i)); + } + return s_expr{"morphology"_symbol, slist_range(branches)}; +} +s_expr mksexp(const meta_data& meta) { + return slist("meta-data"_symbol, slist("version"_symbol, s_expr(meta.version))); +} + +// Implement public facing s-expr writers +std::ostream& write_component(std::ostream& o, const decor& x, const meta_data& m) { + if (m.version != acc_version()) { + throw cableio_version_error(m.version); + } + return o << s_expr{"arbor-component"_symbol, slist(mksexp(m), mksexp(x))}; +} +std::ostream& write_component(std::ostream& o, const label_dict& x, const meta_data& m) { + if (m.version != acc_version()) { + throw cableio_version_error(m.version); + } + return o << s_expr{"arbor-component"_symbol, slist(mksexp(m), mksexp(x))}; +} +std::ostream& write_component(std::ostream& o, const morphology& x, const meta_data& m) { + if (m.version != acc_version()) { + throw cableio_version_error(m.version); + } + return o << s_expr{"arbor-component"_symbol, slist(mksexp(m), mksexp(x))}; +} +std::ostream& write_component(std::ostream& o, const cable_cell& x, const meta_data& m) { + if (m.version != acc_version()) { + throw cableio_version_error(m.version); + } + auto cell = s_expr{"cable-cell"_symbol, slist(mksexp(x.morphology()), mksexp(x.labels()), mksexp(x.decorations()))}; + return o << s_expr{"arbor-component"_symbol, slist(mksexp(m), cell)}; +} +std::ostream& write_component(std::ostream& o, const cable_cell_component& x) { + if (x.meta.version != acc_version()) { + throw cableio_version_error(x.meta.version); + } + std::visit([&](auto&& c){write_component(o, c, x.meta);}, x.component); + return o; +} + +// Anonymous namespace containing helper functions and types for parsing s-expr +namespace { +// Test whether a value wrapped in std::any can be converted to a target type +template <typename T> +bool match(const std::type_info& info) { + return info == typeid(T); +} +template <> +bool match<double>(const std::type_info& info) { + return info == typeid(double) || info == typeid(int); +} + +// Convert a value wrapped in a std::any to target type. +template <typename T> +T eval_cast(std::any arg) { + return std::move(std::any_cast<T&>(arg)); +} +template <> +double eval_cast<double>(std::any arg) { + if (arg.type()==typeid(int)) return std::any_cast<int>(arg); + return std::any_cast<double>(arg); +} + +// Convert a value wrapped in a std::any to an optional std::variant type +template <typename T, std::size_t I=0> +std::optional<T> eval_cast_variant(const std::any& a) { + if constexpr (I<std::variant_size_v<T>) { + using var_type = std::variant_alternative_t<I, T>; + return match<var_type>(a.type())? eval_cast<var_type>(a): eval_cast_variant<T, I+1>(a); + } + return std::nullopt; +} + +// Useful tuple aliases +using envelope_tuple = std::tuple<double,double>; +using pulse_tuple = std::tuple<double,double,double>; +using param_tuple = std::tuple<std::string,double>; +using branch_tuple = std::tuple<int,int,std::vector<arb::msegment>>; +using version_tuple = std::tuple<std::string>; + +// Define makers for defaultables, paintables, placeables +#define ARBIO_DEFINE_SINGLE_ARG(name) arb::name make_##name(double val) { return arb::name{val}; } +#define ARBIO_DEFINE_DOUBLE_ARG(name) arb::name make_##name(const std::string& ion, double val) { return arb::name{ion, val};} + +ARB_PP_FOREACH(ARBIO_DEFINE_SINGLE_ARG, init_membrane_potential, temperature_K, axial_resistivity, membrane_capacitance, threshold_detector) +ARB_PP_FOREACH(ARBIO_DEFINE_DOUBLE_ARG, init_int_concentration, init_ext_concentration, init_reversal_potential) + +std::vector<arb::i_clamp::envelope_point> make_envelope(const std::vector<std::variant<envelope_tuple>>& vec) { + std::vector<arb::i_clamp::envelope_point> envlp; + std::transform(vec.begin(), vec.end(), std::back_inserter(envlp), + [](const auto& x){ + auto t = std::get<envelope_tuple>(x); + return arb::i_clamp::envelope_point{std::get<0>(t), std::get<1>(t)}; + }); + return envlp; +} +arb::i_clamp make_i_clamp(const std::vector<arb::i_clamp::envelope_point>& envlp, double freq) { + return arb::i_clamp(envlp, freq); +} +pulse_tuple make_envelope_pulse(double delay, double duration, double amplitude) { + return pulse_tuple{delay, duration, amplitude}; +} +arb::i_clamp make_i_clamp_pulse(pulse_tuple p, double freq) { + return arb::i_clamp(std::get<0>(p), std::get<1>(p), std::get<2>(p), freq); +} +arb::gap_junction_site make_gap_junction_site() { + return arb::gap_junction_site{}; +} +arb::ion_reversal_potential_method make_ion_reversal_potential_method(const std::string& ion, const arb::mechanism_desc& mech) { + return ion_reversal_potential_method{ion, mech}; +} +#undef ARBIO_DEFINE_SINGLE_ARG +#undef ARBIO_DEFINE_DOUBLE_ARG + +// Define makers for placeable pairs, paintable pairs, defaultables and decors +using place_pair = std::pair<arb::locset, arb::placeable>; +using paint_pair = std::pair<arb::region, arb::paintable>; +place_pair make_place(locset where, placeable what) { + return place_pair{where, what}; +} +paint_pair make_paint(region where, paintable what) { + return paint_pair{where, what}; +} +defaultable make_default(defaultable what) { + return what; +} +decor make_decor(const std::vector<std::variant<place_pair, paint_pair, defaultable>>& args) { + decor d; + for(const auto& a: args) { + auto decor_visitor = arb::util::overload( + [&](const place_pair & p) { d.place(p.first, p.second); }, + [&](const paint_pair & p) { d.paint(p.first, p.second); }, + [&](const defaultable & p){ d.set_default(p); }); + std::visit(decor_visitor, a); + } + return d; +} + +// Define maker for locset pairs, region pairs and label_dicts +using locset_pair = std::pair<std::string, locset>; +using region_pair = std::pair<std::string, region>; +locset_pair make_locset_pair(std::string name, locset desc) { + return locset_pair{name, desc}; +} +region_pair make_region_pair(std::string name, region desc) { + return region_pair{name, desc}; +} +label_dict make_label_dict(const std::vector<std::variant<locset_pair, region_pair>>& args) { + label_dict d; + for(const auto& a: args) { + std::visit([&](auto&& x){d.set(x.first, x.second);}, a); + } + return d; +} +// Define makers for mpoints and msegments and morphologies +arb::mpoint make_point(double x, double y, double z, double r) { + return arb::mpoint{x, y, z, r}; +} +arb::msegment make_segment(unsigned id, arb::mpoint prox, arb::mpoint dist, int tag) { + return arb::msegment{id, prox, dist, tag}; +} +morphology make_morphology(const std::vector<std::variant<branch_tuple>>& args) { + segment_tree tree; + std::vector<unsigned> branch_final_seg(args.size()); + std::vector<unsigned> branch_children(args.size(), 0); + std::vector<std::pair<msegment, int>> segs; + for (const auto& br: args) { + auto b = std::get<branch_tuple>(br); + auto b_id = std::get<0>(b); + auto b_pid = std::get<1>(b); + auto b_segments = std::get<2>(b); + + if (b_pid!=-1) branch_children[b_pid]++; + + auto s_pid = b_pid==-1? arb::mnpos: branch_final_seg[b_pid]; + for (const auto& s: b_segments) { + segs.emplace_back(s, s_pid); + s_pid = s.id; + } + branch_final_seg[b_id] = s_pid; + } + // A branch should have either 0 or >1 children. + auto it = std::find_if(branch_children.begin(), branch_children.end(), [](unsigned i){return i==1;}); + if (it != branch_children.end()) { + throw cableio_morphology_error(it - branch_children.begin()); + } + + // Append segments according to id order. + std::sort(segs.begin(), segs.end(), [](const auto& lhs, const auto& rhs){return lhs.first.id < rhs.first.id;}); + for (const auto& [seg, s_pid]: segs) { + tree.append(s_pid, seg.prox, seg.dist, seg.tag); + } + return morphology(tree); +} + +// Define cable-cell maker +// Accepts the morphology, decor and label_dict arguments in any order as a vector +cable_cell make_cable_cell(const std::vector<std::variant<morphology, label_dict, decor>>& args) { + decor dec; + label_dict dict; + morphology morpho; + for(const auto& a: args) { + auto cable_cell_visitor = arb::util::overload( + [&](const morphology & p) { morpho = p; }, + [&](const label_dict & p) { dict = p; }, + [&](const decor & p){ dec = p; }); + std::visit(cable_cell_visitor, a); + } + return cable_cell(morpho, dict, dec); +} +version_tuple make_version(std::string v) { + return version_tuple{v}; +} +meta_data make_meta_data(version_tuple v) { + return meta_data{std::get<0>(v)}; +} +template <typename T> +cable_cell_component make_component(const meta_data& m, const T& d) { + return cable_cell_component{m, d}; +} + +// Evaluator: member of make_call, make_arg_vec_call, make_mech_call, make_branch_call, make_unordered_call +struct evaluator { + using any_vec = std::vector<std::any>; + using eval_fn = std::function<std::any(any_vec)>; + using args_fn = std::function<bool(const any_vec&)>; + + eval_fn eval; + args_fn match_args; + const char* message; + + evaluator(eval_fn f, args_fn a, const char* m): + eval(std::move(f)), + match_args(std::move(a)), + message(m) + {} + + std::any operator()(any_vec args) { + return eval(std::move(args)); + } +}; + +// Test whether a list of arguments passed as a std::vector<std::any> can be converted +// to the types in Args. +// +// For example, the following would return true: +// +// call_match<int, int, string>(vector<any(4), any(12), any(string("hello"))>) +template <typename... Args> +struct call_match { + template <std::size_t I, typename T, typename Q, typename... Rest> + bool match_args_impl(const std::vector<std::any>& args) const { + return match<T>(args[I].type()) && match_args_impl<I+1, Q, Rest...>(args); + } + + template <std::size_t I, typename T> + bool match_args_impl(const std::vector<std::any>& args) const { + return match<T>(args[I].type()); + } + + template <std::size_t I> + bool match_args_impl(const std::vector<std::any>& args) const { + return true; + } + + bool operator()(const std::vector<std::any>& args) const { + const auto nargs_in = args.size(); + const auto nargs_ex = sizeof...(Args); + return nargs_in==nargs_ex? match_args_impl<0, Args...>(args): false; + } +}; +// Evaluate a call to a function where the arguments are provided as a std::vector<std::any>. +// The arguments are expanded and converted to the correct types, as specified by Args. +template <typename... Args> +struct call_eval { + using ftype = std::function<std::any(Args...)>; + ftype f; + call_eval(ftype f): f(std::move(f)) {} + + template<std::size_t... I> + std::any expand_args_then_eval(const std::vector<std::any>& args, std::index_sequence<I...>) { + return f(eval_cast<Args>(std::move(args[I]))...); + } + + std::any operator()(const std::vector<std::any>& args) { + return expand_args_then_eval(std::move(args), std::make_index_sequence<sizeof...(Args)>()); + } +}; +// Wrap call_match and call_eval in an evaluator +template <typename... Args> +struct make_call { + evaluator state; + + template <typename F> + make_call(F&& f, const char* msg="call"): + state(call_eval<Args...>(std::forward<F>(f)), call_match<Args...>(), msg) + {} + operator evaluator() const { + return state; + } +}; + +// Test whether a list of arguments passed as a std::vector<std::any> can be converted +// to a std::vector<std::variant<Args...>>. +// +// For example, the following would return true: +// +// call_match<int, string>(vector<any(4), any(12), any(string("hello"))>) +template <typename... Args> +struct arg_vec_match { + template <typename T, typename Q, typename... Rest> + bool match_args_impl(const std::any& arg) const { + return match<T>(arg.type()) || match_args_impl<Q, Rest...>(arg); + } + + template <typename T> + bool match_args_impl(const std::any& arg) const { + return match<T>(arg.type()); + } + + bool operator()(const std::vector<std::any>& args) const { + for (const auto& a: args) { + if (!match_args_impl<Args...>(a)) return false; + } + return true; + } +}; +// Evaluate a call to a function where the arguments are provided as a std::vector<std::any>. +// The arguments are converted to std::variant<Args...> and passed to the function as a std::vector. +template <typename... Args> +struct arg_vec_eval { + using ftype = std::function<std::any(std::vector<std::variant<Args...>>)>; + ftype f; + arg_vec_eval(ftype f): f(std::move(f)) {} + + std::any operator()(const std::vector<std::any>& args) { + std::vector<std::variant<Args...>> vars; + for (const auto& a: args) { + vars.push_back(eval_cast_variant<std::variant<Args...>>(a).value()); + } + return f(vars); + } +}; +// Wrap arg_vec_match and arg_vec_eval in an evaluator +template <typename... Args> +struct make_arg_vec_call { + evaluator state; + + template <typename F> + make_arg_vec_call(F&& f, const char* msg="argument vector"): + state(arg_vec_eval<Args...>(std::forward<F>(f)), arg_vec_match<Args...>(), msg) + {} + operator evaluator() const { + return state; + } +}; + +// Test whether a list of arguments passed as a std::vector<std::any> can be converted +// to a string followed by any number of std::pair<std::string, double> +struct mech_match { + bool operator()(const std::vector<std::any>& args) const { + // First argument is the mech name + if (args.size() < 1) return false; + if (!match<std::string>(args.front().type())) return false; + + // The rest of the arguments should be param_tuples + for (auto it = args.begin()+1; it != args.end(); ++it) { + if (!match<param_tuple>(it->type())) return false; + } + return true; + } +}; +// Create a mechanism_desc from a std::vector<std::any>. +struct mech_eval { + arb::mechanism_desc operator()(const std::vector<std::any>& args) { + auto name = eval_cast<std::string>(args.front()); + arb::mechanism_desc mech(name); + for (auto it = args.begin()+1; it != args.end(); ++it) { + auto p = eval_cast<param_tuple>(*it); + mech.set(std::get<0>(p), std::get<1>(p)); + } + return mech; + } +}; +// Wrap mech_match and mech_eval in an evaluator +struct make_mech_call { + evaluator state; + make_mech_call(const char* msg="mechanism"): + state(mech_eval(), mech_match(), msg) + {} + operator evaluator() const { + return state; + } +}; + +// Test whether a list of arguments passed as a std::vector<std::any> can be converted +// to 2 integers followed by at least 1 msegment +struct branch_match { + bool operator()(const std::vector<std::any>& args) const { + if (args.size() < 2) return false; + auto it = args.begin(); + if (!match<int>(it++->type())) return false; + if (!match<int>(it++->type())) return false; + if (it == args.end()) return false; + for (; it != args.end(); ++it) { + if(!match<arb::msegment>(it->type())) return false; + } + return true; + } +}; +// Create a `branch` from a std::vector<std::any>. +struct branch_eval { + branch_tuple operator()(const std::vector<std::any>& args) { + std::vector<msegment> segs; + auto it = args.begin(); + auto id = eval_cast<int>(*it++); + auto parent = eval_cast<int>(*it++); + for (; it != args.end(); ++it) { + segs.push_back(eval_cast<msegment>(*it)); + } + return branch_tuple{id, parent, segs}; + } +}; +// Wrap branch_match and branch_eval in an evaluator +struct make_branch_call { + evaluator state; + make_branch_call(const char* msg="branch"): + state(branch_eval(), branch_match(), msg) + {} + operator evaluator() const { + return state; + } +}; + +// Test whether a list of arguments passed as a std::vector<std::any> `args` can be +// converted to a std::vector<std::variant<Args...>>. +// - `args` must have the same size as Args... +// - no more than one element in `args` can match a given template parameter of Args... +// +// For example, the following would return true: +// +// call_match<int, string>(vector<any(4), any(string("hello"))>) +// call_match<int, string>(vector<any(string("hello")), any(4)>) +// +// And the following would return false: +// +// call_match<int, string>(vector<any(4), any(string("hello")), any(string("bye"))>) +// call_match<int, string>(vector<any(4), any(2)>) +// +// Not an efficient implementation, but should be okay for a few arguments. +template <typename... Args> +struct unordered_match { + template <typename T, typename Q, typename... Rest> + bool match_args_impl(const std::vector<std::any>& args) const { + bool found_match = false; + for (const auto& a: args) { + auto new_match = match<T>(a.type()); + if (new_match && found_match) return false; + found_match = new_match; + } + return found_match || match_args_impl<Q, Rest...>(args); + } + + template <typename T> + bool match_args_impl(const std::vector<std::any>& args) const { + bool found_match = false; + for (const auto& a: args) { + auto new_match = match<T>(a.type()); + if (new_match && found_match) return false; + found_match = new_match; + } + return found_match; + } + + bool operator()(const std::vector<std::any>& args) const { + const auto nargs_in = args.size(); + const auto nargs_ex = sizeof...(Args); + return (nargs_in == nargs_ex) && match_args_impl<Args...>(args); + } +}; +// Wrap unordered_match and arg_vec_eval in an evaluator +template <typename... Args> +struct make_unordered_call { + evaluator state; + + template <typename F> + make_unordered_call(F&& f, const char* msg="call"): + state(arg_vec_eval<Args...>(std::forward<F>(f)), unordered_match<Args...>(), msg) + {} + operator evaluator() const { + return state; + } +}; +} // anonymous namespace + +using eval_map = std::unordered_multimap<std::string, evaluator>; +using eval_vec = std::vector<evaluator>; + +// Parse s-expression into std::any given a function evaluation map and a tuple evaluation vector. +parse_hopefully<std::any> eval(const arb::s_expr&, const eval_map&, const eval_vec&); + +parse_hopefully<std::vector<std::any>> eval_args(const s_expr& e, const eval_map& map, const eval_vec& vec) { + if (!e) return {std::vector<std::any>{}}; + std::vector<std::any> args; + for (auto& h: e) { + if (auto arg=eval(h, map, vec)) { + args.push_back(std::move(*arg)); + } + else { + return util::unexpected(std::move(arg.error())); + } + } + return args; +} + +parse_hopefully<std::any> eval(const s_expr& e, const eval_map& map, const eval_vec& vec) { + if (e.is_atom()) { + auto& t = e.atom(); + switch (t.kind) { + case tok::integer: + return {std::stoi(t.spelling)}; + case tok::real: + return {std::stod(t.spelling)}; + case tok::nil: + return {nil_tag()}; + case tok::string: + return std::any{std::string(t.spelling)}; + // An arbitrary symbol in a region/locset expression is an error, and is + // often a result of not quoting a label correctly. + case tok::symbol: + return util::unexpected(cableio_parse_error("Unexpected symbol "+e.atom().spelling, location(e))); + case tok::error: + return util::unexpected(cableio_parse_error("Unexpected term "+e.atom().spelling, location(e))); + default: + return util::unexpected(cableio_parse_error("Unexpected term "+e.atom().spelling, location(e))); + } + } + if (e.head().is_atom()) { + // If this is not a symbol, parse as a tuple + if (e.head().atom().kind != tok::symbol) { + auto args = eval_args(e, map, vec); + if (!args) { + return util::unexpected(args.error()); + } + for (auto& e: vec) { + if (e.match_args(*args)) { // found a match: evaluate and return. + return e.eval(*args); + } + } + + // Unable to find a match: try to return a helpful error message. + const auto nc = vec.size(); + std::string msg = "No matches for found for unnamed tuple with "+std::to_string(args->size())+" arguments.\n" + "There are "+std::to_string(nc)+" potential candiates"+(nc?":":"."); + int count = 0; + for (auto& e: vec) { + msg += "\n Candidate "+std::to_string(++count)+": "+e.message; + } + return util::unexpected(cableio_parse_error(msg, location(e))); + }; + + // Otherwise this must be a function evaluation, where head is the function name, + // and tail is a list of arguments. + // Evaluate the arguments, and return error state if an error occurred. + auto args = eval_args(e.tail(), map, vec); + if (!args) { + return util::unexpected(args.error()); + } + + // Find all candidate functions that match the name of the function. + auto& name = e.head().atom().spelling; + auto matches = map.equal_range(name); + + // Search for a candidate that matches the argument list. + for (auto i=matches.first; i!=matches.second; ++i) { + if (i->second.match_args(*args)) { // found a match: evaluate and return. + return i->second.eval(*args); + } + } + + // If it's not in the provided map, maybe it's a label expression + // the corresponding parser is provided by the arbor lib + if (auto l = parse_label_expression(e)) { + if (match<region>(l->type())) return eval_cast<region>(l.value()); + if (match<locset>(l->type())) return eval_cast<locset>(l.value()); + } + + // Unable to find a match: try to return a helpful error message. + const auto nc = std::distance(matches.first, matches.second); + std::string msg = "No matches for found for "+name+" with "+std::to_string(args->size())+" arguments.\n" + "There are "+std::to_string(nc)+" potential candiates"+(nc?":":"."); + int count = 0; + for (auto i=matches.first; i!=matches.second; ++i) { + msg += "\n Candidate "+std::to_string(++count)+": "+i->second.message; + } + return util::unexpected(cableio_parse_error(msg, location(e))); + } + return util::unexpected(cableio_parse_error("Expression is not integer, real expression of the form (op <args>) nor tuple of the form (e0 e1 ... en)", location(e))); +} + +eval_map named_evals{ + {"membrane-potential", make_call<double>(make_init_membrane_potential, + "'membrane-potential' with 1 argument (val:real)")}, + {"temperature-kelvin", make_call<double>(make_temperature_K, + "'temperature-kelvin' with 1 argument (val:real)")}, + {"axial-resistivity", make_call<double>(make_axial_resistivity, + "'axial-resistivity' with 1 argument (val:real)")}, + {"membrane-capacitance", make_call<double>(make_membrane_capacitance, + "'membrane-capacitance' with 1 argument (val:real)")}, + {"ion-internal-concentration", make_call<std::string, double>(make_init_int_concentration, + "'ion_internal_concentration' with 2 arguments (ion:string val:real)")}, + {"ion-external-concentration", make_call<std::string, double>(make_init_ext_concentration, + "'ion_external_concentration' with 2 arguments (ion:string val:real)")}, + {"ion-reversal-potential", make_call<std::string, double>(make_init_reversal_potential, + "'ion_reversal_potential' with 2 arguments (ion:string val:real)")}, + {"envelope", make_arg_vec_call<envelope_tuple>(make_envelope, + "`envelope` with one or more pairs of start time and amplitude (start:real amplitude:real)")}, + {"envelope-pulse", make_call<double, double, double>(make_envelope_pulse, + "'envelope-pulse' with 3 arguments (delay:real duration:real amplitude:real)")}, + {"current-clamp", make_call<std::vector<arb::i_clamp::envelope_point>, double>(make_i_clamp, + "`current-clamp` with 2 arguments (env:envelope freq:real)")}, + {"current-clamp", make_call<pulse_tuple, double>(make_i_clamp_pulse, + "`current-clamp` with 2 arguments (env:envelope_pulse freq:real)")}, + {"threshold-detector", make_call<double>(make_threshold_detector, + "'threshold-detector' with 1 argument (threshold:real)")}, + {"gap-junction-site", make_call<>(make_gap_junction_site, + "'gap-junction-site' with 0 arguments")}, + {"ion-reversal-potential-method", make_call<std::string, arb::mechanism_desc>(make_ion_reversal_potential_method, + "'ion-reversal-potential-method' with 2 ""arguments (ion:string mech:mechanism)")}, + {"mechanism", make_mech_call("'mechanism' with a name argument, and 0 or more parameter settings" + "(name:string (param:string val:real))")}, + + {"place", make_call<locset, gap_junction_site>(make_place, "'place' with 2 arguments (locset gap-junction-site)")}, + {"place", make_call<locset, i_clamp>(make_place, "'place' with 2 arguments (locset current-clamp)")}, + {"place", make_call<locset, threshold_detector>(make_place, "'place' with 2 arguments (locset threshold-detector)")}, + {"place", make_call<locset, mechanism_desc>(make_place, "'place' with 2 arguments (locset mechanism)")}, + + {"paint", make_call<region, init_membrane_potential>(make_paint, "'paint' with 2 arguments (region membrane-potential)")}, + {"paint", make_call<region, temperature_K>(make_paint, "'paint' with 2 arguments (region temperature-kelvin)")}, + {"paint", make_call<region, membrane_capacitance>(make_paint, "'paint' with 2 arguments (region membrane-capacitance)")}, + {"paint", make_call<region, axial_resistivity>(make_paint, "'paint' with 2 arguments (region axial-resistivity)")}, + {"paint", make_call<region, init_int_concentration>(make_paint, "'paint' with 2 arguments (region ion-internal-concentration)")}, + {"paint", make_call<region, init_ext_concentration>(make_paint, "'paint' with 2 arguments (region ion-external-concentration)")}, + {"paint", make_call<region, init_reversal_potential>(make_paint, "'paint' with 2 arguments (region ion-reversal-potential)")}, + {"paint", make_call<region, mechanism_desc>(make_paint, "'paint' with 2 arguments (region mechanism)")}, + + {"default", make_call<init_membrane_potential>(make_default, "'default' with 1 argument (membrane-potential)")}, + {"default", make_call<temperature_K>(make_default, "'default' with 1 argument (temperature-kelvin)")}, + {"default", make_call<membrane_capacitance>(make_default, "'default' with 1 argument (membrane-capacitance)")}, + {"default", make_call<axial_resistivity>(make_default, "'default' with 1 argument (axial-resistivity)")}, + {"default", make_call<init_int_concentration>(make_default, "'default' with 1 argument (ion-internal-concentration)")}, + {"default", make_call<init_ext_concentration>(make_default, "'default' with 1 argument (ion-external-concentration)")}, + {"default", make_call<init_reversal_potential>(make_default, "'default' with 1 argument (ion-reversal-potential)")}, + {"default", make_call<ion_reversal_potential_method>(make_default, "'default' with 1 argument (ion-reversal-potential-method)")}, + + {"locset-def", make_call<std::string, locset>(make_locset_pair, + "'locset-def' with 2 arguments (name:string ls:locset)")}, + {"region-def", make_call<std::string, region>(make_region_pair, + "'region-def' with 2 arguments (name:string reg:region)")}, + + {"point", make_call<double, double, double, double>(make_point, + "'point' with 4 arguments (x:real y:real z:real radius:real)")}, + {"segment", make_call<int, mpoint, mpoint, int>(make_segment, + "'segment' with 4 arguments (parent:int prox:point dist:point tag:int)")}, + {"branch", make_branch_call( + "'branch' with 2 integers and 1 or more segment arguments (id:int parent:int s0:segment s1:segment ..)")}, + + {"decor", make_arg_vec_call<place_pair, paint_pair, defaultable>(make_decor, + "'decor' with 1 or more `paint`, `place` or `default` arguments")}, + {"label-dict", make_arg_vec_call<locset_pair, region_pair>(make_label_dict, + "'label-dict' with 1 or more `locset-def` or `region-def` arguments")}, + {"morphology", make_arg_vec_call<branch_tuple>(make_morphology, + "'morphology' 1 or more `branch` arguments")}, + + {"cable-cell", make_unordered_call<morphology, label_dict, decor>(make_cable_cell, + "'cable-cell' with 3 arguments: `morphology`, `label-dict`, and `decor` in any order")}, + + {"version", make_call<std::string>(make_version, "'version' with one argment (val:std::string)")}, + {"meta-data", make_call<version_tuple>(make_meta_data, "'meta-data' with one argument (v:version)")}, + + { "arbor-component", make_call<meta_data, decor>(make_component<decor>, "'arbor-component' with 2 arguments (m:meta_data p:decor)")}, + { "arbor-component", make_call<meta_data, label_dict>(make_component<label_dict>, "'arbor-component' with 2 arguments (m:meta_data p:label_dict)")}, + { "arbor-component", make_call<meta_data, morphology>(make_component<morphology>, "'arbor-component' with 2 arguments (m:meta_data p:morphology)")}, + { "arbor-component", make_call<meta_data, cable_cell>(make_component<cable_cell>, "'arbor-component' with 2 arguments (m:meta_data p:cable_cell)")} +}; + +eval_vec unnamed_evals{ + make_call<std::string, double>(std::make_tuple<std::string, double>, "tuple<std::string, double>"), + make_call<double, double>(std::make_tuple<double, double>, "tuple<double, double>") +}; + +inline parse_hopefully<std::any> parse(const arb::s_expr& s) { + return eval(std::move(s), named_evals, unnamed_evals); +} + +parse_hopefully<std::any> parse_expression(const std::string& s) { + return parse(parse_s_expr(s)); +} + +// Read s-expr +parse_hopefully<cable_cell_component> parse_component(const std::string& s) { + auto sexp = parse_s_expr(s); + auto try_parse = parse(sexp); + if (!try_parse) { + return util::unexpected(cableio_parse_error(try_parse.error())); + } + if (!match<cable_cell_component>(try_parse.value().type())) { + return util::unexpected(cableio_parse_error("Expected arbor-component", location(sexp))); + } + auto comp = eval_cast<cable_cell_component>(try_parse.value()); + if (comp.meta.version != acc_version()) { + return util::unexpected(cableio_parse_error("Unsupported cable-cell format version "+ comp.meta.version, location(sexp))); + } + return comp; +}; + +parse_hopefully<cable_cell_component> parse_component(std::istream& s) { + return parse_component(std::string(std::istreambuf_iterator<char>(s), {})); +} +} // namespace arborio \ No newline at end of file diff --git a/arborio/include/arborio/cableio.hpp b/arborio/include/arborio/cableio.hpp new file mode 100644 index 0000000000000000000000000000000000000000..7502e7510bf9b114a35ed92c3fd41af235a1adb3 --- /dev/null +++ b/arborio/include/arborio/cableio.hpp @@ -0,0 +1,42 @@ +#pragma once + +#include <arbor/cable_cell.hpp> +#include <arbor/s_expr.hpp> + +namespace arborio { +std::string acc_version(); + +struct cableio_parse_error: arb::arbor_exception { + explicit cableio_parse_error(const std::string& msg, const arb::src_location& loc); +}; + +struct cableio_morphology_error: arb::arbor_exception { + explicit cableio_morphology_error(const unsigned bid); +}; + +struct cableio_version_error: arb::arbor_exception { + explicit cableio_version_error(const std::string& version); +}; + +template <typename T> +using parse_hopefully = arb::util::expected<T, cableio_parse_error>; +using cable_cell_variant = std::variant<arb::morphology, arb::label_dict, arb::decor, arb::cable_cell>; + +struct meta_data { + std::string version = acc_version(); +}; +struct cable_cell_component { + meta_data meta; + cable_cell_variant component; +}; + +std::ostream& write_component(std::ostream&, const cable_cell_component&); +std::ostream& write_component(std::ostream&, const arb::decor& x, const meta_data& m = {}); +std::ostream& write_component(std::ostream&, const arb::label_dict& x, const meta_data& m = {}); +std::ostream& write_component(std::ostream&, const arb::morphology& x, const meta_data& m = {}); +std::ostream& write_component(std::ostream&, const arb::cable_cell& x, const meta_data& m = {}); + +parse_hopefully<cable_cell_component> parse_component(const std::string&); +parse_hopefully<cable_cell_component> parse_component(std::istream&); + +} // namespace arborio \ No newline at end of file diff --git a/arborio/parse_s_expr.hpp b/arborio/parse_s_expr.hpp new file mode 100644 index 0000000000000000000000000000000000000000..b240897dd65490bf5b9c852f5dbb9d01eb289bbd --- /dev/null +++ b/arborio/parse_s_expr.hpp @@ -0,0 +1,45 @@ +#pragma once + +#include <arbor/s_expr.hpp> +#include <arborio/cableio.hpp> + +namespace arborio { +using namespace arb; + +// Helper function for programmatically building lists +// slist(1, 2, "hello world", "banjax@cat/3"_symbol); +// Produces the following s-expression: +// (1 2 "hello world" banjax@cat/3) +// Can be nested: +// slist(1, slist(2, 3), 4, 5 ); +// Produces: +// (1 (2 3) 4 5) + +template <typename T> +s_expr slist(T v) { + return {v, {}}; +} + +template <typename T, typename... Args> +s_expr slist(T v, Args... args) { + return {v, slist(args...)}; +} + +inline s_expr slist() { + return {}; +} + +template <typename I, typename S> +s_expr slist_range(I b, S e) { + return b==e ? s_expr{} + : s_expr{*b, slist_range(++b,e)}; +} + +template <typename Range> +s_expr slist_range(const Range& range) { + return slist_range(std::begin(range), std::end(range)); +} + +parse_hopefully<std::any> parse_expression(const std::string&); + +} // namespace arborio diff --git a/doc/concepts/decor.rst b/doc/concepts/decor.rst index 2851e931ec54466e6df6476e1e6b55da7072a62e..93f6c0b53ae5da08a0151911ecfdcb4d566a08ad 100644 --- a/doc/concepts/decor.rst +++ b/doc/concepts/decor.rst @@ -287,7 +287,7 @@ See :ref:`modelgapjunctions`. 4. Stimuli ~~~~~~~~~~ -A current stimulus is a DC or sinusoidal current of fixed frequemcy with a time-varying amplitude +A current stimulus is a DC or sinusoidal current of fixed frequency with a time-varying amplitude governed by a piecewise-linear envelope. The stimulus is described by two parameters: a frequency in Hertz, where a value of zero denotes DC; diff --git a/doc/concepts/labels.rst b/doc/concepts/labels.rst index eae960b78b23bf7cd1f16afdae826c980076ea40..d99a499c05fb38c32e88c9d3eb6333c64eaec518 100644 --- a/doc/concepts/labels.rst +++ b/doc/concepts/labels.rst @@ -197,6 +197,7 @@ dendritic tree where the radius first is less than or equal to 0.2 μm. thing could be achieved using hoc in NEURON, and whether it would be free of bugs and applicable to arbitrary morphologies. +.. _labels-locset-expr: Locset expressions ~~~~~~~~~~~~~~~~~~ @@ -345,6 +346,8 @@ Locset expressions (join (location 1 0.5) (location 2 0.1) (location 1 0.2) (location 1 0.5) (location 4 0)) +.. _labels-region-expr: + Region expressions ~~~~~~~~~~~~~~~~~~ diff --git a/doc/cpp/cable_cell.rst b/doc/cpp/cable_cell.rst index eab88668f01dcf78401c8b06608d2d3635c478db..ae46477f4b8fc1e9d91e1932397f63a4b2f4e7fc 100644 --- a/doc/cpp/cable_cell.rst +++ b/doc/cpp/cable_cell.rst @@ -8,6 +8,7 @@ Cable cells morphology probe_sample + cable_cell_format .. Warning:: The interface for building and modifying cable cell objects diff --git a/doc/cpp/cable_cell_format.rst b/doc/cpp/cable_cell_format.rst new file mode 100644 index 0000000000000000000000000000000000000000..d0dc5b82ce457c48da2055c16f9b808d0f1e7418 --- /dev/null +++ b/doc/cpp/cable_cell_format.rst @@ -0,0 +1,88 @@ +.. _cppcablecellformat: + +Description Format +================== + +Arbor provides readers and writers for describing :ref:`label dictionaries <labels>`, +:ref:`decoration objects <cablecell-decoration>`, :ref:`morphologies <morph>` and +:ref:`cable cells <cablecell>`, referred to here as *arbor-components*. + +A detailed description of the s-expression format used to describe each of these components +can be found :ref:`here <formatcablecell>`. + +Reading and writing of the arbor-component description format is delegated to the ``arborio`` +library and the responsible classes and functions are present in the ``arborio`` namespace. + +The arbor-components and meta-data +---------------------------------- + +.. cpp:type:: cable_cell_variant = std::variant<arb::morphology, arb::label_dict, arb::decor, arb::cable_cell> + +.. cpp:type:: template <typename T> parse_hopefully = arb::util::expected<T, cableio_parse_error> + + ``arborio::cableio_parse_error`` is derived from ``arb::arbor_exception`` which in turn is derived + from ``std::runtime_error``. It contains information about the ``line`` and ``position`` of an encountered + error in a document. + + ``arb::util::expected`` contains either an object of type ``T`` or an error object. + +.. cpp:class:: meta_data + + .. cpp:member:: std::string version + + Stores the version of the format being used. + +.. cpp:class:: cable_cell_component + + .. cpp:member:: meta_data meta + + Stores meta-data pertaining to the description of a cable cell component. + + .. cpp:member:: cable_cell_variant component + + Stores one of :cpp:class:`decor`, :cpp:class:`label_dict`, :cpp:class:`morphology` or :cpp:class:`cable_cell`. + +Reading arbor-components +------------------------ + +.. cpp:function:: parse_hopefully<cable_cell_component> parse_component(const std::string&) + + This function will attempt to construct a :cpp:class:`cable_cell_component` object by parsing the + contents of a string. It will return a :cpp:type:`parse_hopefully` containing the constructed object, + or, if parsing fails, a helpful ``cableio_parse_error``. + +.. cpp:function:: parse_hopefully<cable_cell_component> parse_component(std::istream&) + + Performs the same functionality as ``parse_component`` above, but starting from + ``std::istream``. + +Writing arbor-components +------------------------ + +.. cpp:function:: std::ostream& write_component(std::ostream&, const cable_cell_component&) + + Writes the :cpp:class:`cable_cell_component` object to the given ``std::ostream``. + +.. cpp:function:: std::ostream& write_component(std::ostream& o, const arb::decor& x, const meta_data& m = {}) + + Constructs a :cpp:class:`cable_cell_component` from a :cpp:class:`decor` object, and optional + :cpp:class:`meta_data`. If no meta_data is provided, the most recent version of + the format is used to create it. The resulting object is written to the given ``std::ostream``. + +.. cpp:function:: std::ostream& write_component(std::ostream& o, const arb::label_dict& x, const meta_data& m = {}) + + Constructs a :cpp:class:`cable_cell_component` from a :cpp:class:`label_dict` object, and optional + :cpp:class:`meta_data`. If no meta_data is provided, the most recent version of + the format is used to create it. The resulting object is written to the given ``std::ostream``. + +.. cpp:function:: std::ostream& write_component(std::ostream& o, const arb::morphology& x, const meta_data& m = {}) + + Constructs a :cpp:class:`cable_cell_component` from a :cpp:class:`morphology` object, and optional + :cpp:class:`meta_data`. If no meta_data is provided, the most recent version of + the format is used to create it. The resulting object is written to the given ``std::ostream``. + +.. cpp:function:: std::ostream& write_component(std::ostream& o, const arb::cable_cell& x, const meta_data& m = {}) + + Constructs a :cpp:class:`cable_cell_component` from a :cpp:class:`cable_cell` object, and optional + :cpp:class:`meta_data`. If no meta_data is provided, the most recent version of + the format is used to create it. The resulting object is written to the given ``std::ostream``. \ No newline at end of file diff --git a/doc/fileformat/cable_cell.rst b/doc/fileformat/cable_cell.rst new file mode 100644 index 0000000000000000000000000000000000000000..5419912e01737136708210fb2c98ec5a81ac17b6 --- /dev/null +++ b/doc/fileformat/cable_cell.rst @@ -0,0 +1,440 @@ +.. _formatcablecell: + +Arbor Cable Cell format (ACC) +============================= + +We define an s-expression format for describing :ref:`cable cells <cablecell>`. +Cable cells are constructed from three components: a :ref:`label dictionary <labels>`, +a :ref:`decoration object <cablecell-decoration>` and a :ref:`morphology <morph>`. +The cable cell *format* is constructed in the same way. + +.. Note:: + + Extra line breaks and indentations in the s-expressions presented below have been + added for clarity. They are not a requirement of the format and will be treated as + whitespace. + +Label dictionary +---------------- + +The label dictionary stores :term:`region` and :term:`locset` s-expressions with +associated labels which can then be used to refer to the region or locset when +decorating the cable cell. + +Arbor provides many useful :ref:`region expressions <labels-region-expr>` and +:ref:`locset expressions <labels-locset-expr>` which are explained in detail at the +provided links. + +The components of the label dictionary are the following: + +.. label:: (region-def label:string reg:region) + + This defines a ``label`` which can be used to refer to the region ``reg``. + For example: + + .. code:: lisp + + (region-def "my_region" (branch 1)) + + This expression identifies the branch with id 1 as "my_region". + +.. label:: (locset-def label:string ls:locset) + + This defines a ``label`` which can be used to refer to the locset ``ls``. + For example: + + .. code:: lisp + + (locset-def "my_locset" (location 3 0.5)) + + This expression identifies the midpoint of the branch with id 3 as "my_locset". + +Any number of locset and region definitions can be grouped in a label dictionary as follows: + +.. label:: (label-dict [...def:region-def/locset-def]) + + This describes a label dictionary of zero or more region and locset definitons. + For example: + + .. code:: lisp + + (label-dict + (region-def "my_soma" (tag 1)) + (locset-def "root" (root)) + (region-def "all" (all)) + (region-def "my_region" (radius-ge (region "my_soma") 1.5)) + (locset-def "terminal" (terminal))) + +Decor +----- + +The decor of a cable cell describes the dynamics and properties of the cell which can be assigned on +:term:`regions <region>` or :term:`locsets <locset>`, or set as defaults on the entire cell. + +This table lists all supported dynamics and properties and whether they are *placeable* (i.e. they can +be placed on one or more locations on the cell described by a locset); *paintable* (i.e. they can be set +on an entire area of the cell described by a region) or *defaultable* (i.e. they are the default settings +of the cell): + +.. csv-table:: Property applicability. + :widths: 20, 10, 10, 10 + + , **placeable**, **paintable**, **defaultable** + initial membrane potential, --, ✓, ✓ + axial resistivity, --, ✓, ✓ + temperature, --, ✓, ✓ + membrane capacitance, --, ✓, ✓ + ion initial internal concentration, --, ✓, ✓ + ion initial external concentration, --, ✓, ✓ + ion initial reversal potential, --, ✓, ✓ + ion reversal potential method, --, --, ✓ + density mechanism, --, ✓, -- + point mechanism, ✓, --, -- + current clamp, ✓, --, -- + threshold detector, ✓, --, -- + gap junction site, ✓, --, -- + +The various properties and dynamics of the decor are described as follows: + +.. label:: (membrane-potential val:real) + + This describes an *initial membrane potential* object with value ``val`` (unit mV). + +.. label:: (axial-resistivity val:real) + + This describes an *axial resistivity* object with value ``val`` (unit Ω·cm). + +.. label:: (temperature-kelvin val:real) + + This describes a *temperature* object with value ``val`` (unit K). + +.. label:: (membrane-capacitance val:real) + + This describes a *membrane capacitance* object with value ``val`` (unit F/m²). + +.. label:: (ion-internal-concentration ion:string val:real) + + This describes an *initial internal concentration* object for ion ``ion`` with value ``val`` (unit mM). + +.. label:: (ion-external-concentration ion:string val:real) + + This describes an *initial external concentration* object for ion ``ion`` with value ``val`` (unit mM). + +.. label:: (ion-reversal-potential ion:string val:real) + + This describes an *initial reversal potential* object for ion ``ion`` with value ``val`` (unit mV). + +.. label:: (mechanism name:string [...(param:string val:real)]) + + This describes a (point or density) mechanism object of the mechanism called ``name``. This expression + accepts zero or more ``(param:string val:real)`` expressions. Each of these expressions sets the value of + parameter ``param`` to ``val``. + For example: + + .. code:: lisp + + (mechanism "hh" ("gl" 0.5) ("el" 2)) + + This expression creates an "hh" mechanism and sets the "gl" and "el" parameters of the mechanism to 0.5 + and 2 respectively (units depend on the :ref:`nmodl <formatnmodl>` mechanism). + +.. label:: (ion-reversal-potential-method ion:string method:mechanism) + + This creates a *reversal potential method* (able to modify the reversal potential) of ion ``ion`` from + mechanism ``method``. + For example: + + .. code:: lisp + + (ion-reversal-potential-method "ca" (mechanism "nernst/ca")) + +.. label:: (current-clamp (envelope-pulse delay:real duration:real amplitude:real) freq:real) + + This creates a *current clamp*. If the frequency ``freq`` (unit Hz) is zero, the current is a square + pulse with amplitude ``amplitude`` (unit nA) starting at ``delay`` (unit ms) and lasting for ``duration`` + (unit ms). If ``freq`` is non-zero, the current is sinusoidal with amplitude ``amplitude`` and frequency + ``freq`` from time ``delay`` and lasting for ``duration``. + (More information about current clamps can be found :ref:`here <cablecell-stimuli>`). + +.. label:: (current-clamp [...(envelope time:real amplitude:real)] freq:real) + + This creates a *current clamp* with an amplitude governed by the given envelopes (``time`` unit ms and + ``amplitude`` unit nA). A frequency ``freq`` (unit Hz) of zero implies that the generated current simply + follows the envelope. A non-zero ``freq`` implies the current is sinusoidal with that frequency and amplitude + that varies according to the envelope. (More information about current clamps can be found + :ref:`here <cablecell-stimuli>`). + For example: + + .. code:: + + (current-clamp (envelope (0 10) (50 10) (50 0)) 40) + + This expression describes a sinusoidal current with amplitude 10nA and frequency 40Hz and that lasts + from t = 0ms to t = 50ms, finally leaving the current at 0nA (final amplitude in the envelope). + +.. label:: (threshold-detector val:real). + + This describes a *threshold-detector* object with value ``val`` (unit mV). + +.. label:: (gap-junction-site) + + This describes a *gap-junction-site*. + +*Paintable* and *placeable* properties and dynamics are placed on regions (generated from :ref:`region expressions +<labels-region-expr>`) and locsets (generated from :ref:`locset expressions <labels-locset-expr>`) respectively. +*Defaultable* properties and dynamics apply to an entire cell. + +.. label:: (paint reg:region prop:paintable) + + This applies the painatble property ``prop`` to region ``reg``. + For example: + + .. code:: lisp + + (paint (tag 1) (membrane-capacitance 0.02)) + + This expression sets the membrane capacitance of the region tagged ``1`` to 0.02 F/m². + + +.. label:: (place ls:locset prop:placeable) + + This places the property ``prop`` on locset ``ls``. + For example: + + .. code:: lisp + + (place (locset "mylocset") (threshold-detector 10)) + + This expression places a 10 mV threshold detector on the locset labeled ``mylocset``. + (The definition of ``mylocset`` should be provided in a label dictionary associated + with the decor). + +.. label:: (default prop:defaultable) + + This sets the property ``prop`` as default for the entire cell. (This default property can be overridden on region + using a ``paint`` expression). + For example: + + .. code:: lisp + + (default (membrane-potential -65)) + + This expression sets the default membrane potential of the cell to -65 mV. + +Any number of paint, place and default expressions can be used to create a decor as follows: + +.. label:: (decor [...def:paint/place/default]) + + This describes a decor object with zero or more paint, place or default expressions in any order. + For example: + + .. code:: lisp + + (decor + (default (membrane-potential -55.000000)) + (paint (region "custom") (temperature-kelvin 270)) + (paint (region "soma") (membrane-potential -50.000000)) + (paint (all) (mechanism "pas")) + (paint (tag 4) (mechanism "Ih" ("gbar" 0.001))) + (place (locset "root") (mechanism "expsyn")) + (place (terminal) (gap-junction-site))) + +Morphology +---------- + +The morphology of a cable cell can be described in terms of points, tagged segments and branches. + +.. label:: (point x:real y:real z:real radius:real) + + This describes a 3D *point* in space with ``x``, ``y``, and ``z`` coordinates and a radius ``r`` (unit µm). + +.. label:: (segment id:int prox:point dist:point tag:int) + + This describes a tapered segment from point ``prox`` to point ``dist`` with a tag ``tag`` and id ``id``. + For example: + + .. code:: lisp + + (segment 3 (point 0 0 0 5) (point 0 0 10 2) 1) + + This expression creates a segment with id 3, with a radius that tapers linearly from 5 to 2 µm, which has a + a tag of 1. + +.. label:: (branch id:int parent_id:int seg:segment [...seg:segment]) + + This describes a branch with a given ``id`` which has as a parent the branch with id ``parent_id`` (a + ``parent_id`` equal to -1 means the branch is at the root of the morphology). The branch is composed of 1 or + more contiguous segments ``seg``. + + +.. label:: (morphology [...b:branch]) + + This creates the morphology from a set of branches. There exists more than one valid s-expression to + describe the same morphology. + + For example, the shown morphology can be represented using the following s-expression. If we change + any of the branch or segment ids, we would obtain an identical morphology. + + .. figure:: ../gen-images/label_morph.svg + :width: 600 + :align: center + + On the left the morphology visualized using its segments, on the right using its branches. + Python code to generate this cable cell is in the :class:`segment_tree<arbor.segment_tree>` + documentation :ref:`here <morph-label-seg-code>`. + + .. code:: lisp + + (morphology + (branch 0 -1 + (segment 0 (point 0 0 0 2) (point 4 0 0 2) 1) + (segment 1 (point 4 0 0 0.8) (point 8 0 0 0.8) 3) + (segment 2 (point 8 0 0 0.8) (point 12 -0.5 0 0.8) 3)) + (branch 1 0 + (segment 3 (point 12 -0.5 0 0.8) (point 20 4 0 0.4) 3) + (segment 4 (point 20 4 0 0.4) (point 26 6 0 0.2) 3)) + (branch 2 0 + (segment 5 (point 12 -0.5 0 0.5) (point 19 -3 0 0.5) 3)) + (branch 3 2 + (segment 6 (point 19 -3 0 0.5) (point 24 -7 0 0.2) 3)) + (branch 4 2 + (segment 7 (point 19 -3 0 0.5) (point 23 -1 0 0.2) 3) + (segment 8 (point 23 -1 0 0.3) (point 26 -2 0 0.2) 3)) + (branch 5 -1 + (segment 9 (point 0 0 0 2) (point -7 0 0 0.4) 2) + (segment 10 (point -7 0 0 0.4) (point -10 0 0 0.4) 2))) + +Cable cell +---------- + +The entire cable-cell can then be constructed given the 3 previously described component +expressions. + +.. label:: (cable-cell morph:morphology dec:decor dict:label-dict) + + The arguments of the cable-cell can be in any order, as long as all 3 components are listed. + For example: + + .. code:: lisp + + (cable-cell + (label-dict + (region-def "my_soma" (tag 1)) + (locset-def "root" (root)) + (region-def "all" (all)) + (region-def "my_region" (radius-ge (region "my_soma") 1.5)) + (locset-def "terminal" (terminal))) + (decor + (default (membrane-potential -55.000000)) + (paint (region "my_soma") (temperature-kelvin 270)) + (paint (region "my_region") (membrane-potential -50.000000)) + (paint (tag 4) (mechanism "Ih" ("gbar" 0.001))) + (place (locset "root") (mechanism "expsyn")) + (place (location 1 0.2) (gap-junction-site))) + (morphology + (branch 0 -1 + (segment 0 (point 0 0 0 2) (point 4 0 0 2) 1) + (segment 1 (point 4 0 0 0.8) (point 8 0 0 0.8) 3) + (segment 2 (point 8 0 0 0.8) (point 12 -0.5 0 0.8) 3)) + (branch 1 0 + (segment 3 (point 12 -0.5 0 0.8) (point 20 4 0 0.4) 3) + (segment 4 (point 20 4 0 0.4) (point 26 6 0 0.2) 3)) + (branch 2 0 + (segment 5 (point 12 -0.5 0 0.5) (point 19 -3 0 0.5) 3)) + (branch 3 2 + (segment 6 (point 19 -3 0 0.5) (point 24 -7 0 0.2) 3)) + (branch 4 2 + (segment 7 (point 19 -3 0 0.5) (point 23 -1 0 0.2) 3) + (segment 8 (point 23 -1 0 0.3) (point 26 -2 0 0.2) 3)) + (branch 5 -1 + (segment 9 (point 0 0 0 2) (point -7 0 0 0.4) 2) + (segment 10 (point -7 0 0 0.4) (point -10 0 0 0.4) 2)))) + + This expression uses the *label-dictionary* in the *decoration* specification + to get the descriptions of regions and locsets specified using labels. + The *decor* is then applied on the provided *morphology*, creating a cable cell. + +Parsable arbor-components and meta-data +--------------------------------------- + +The formats described above can be used to generate a :ref:`label dictionary <labels>`, +:ref:`decoration <cablecell-decoration>`, :ref:`morphology <morph>`, or :ref:`cable cell <cablecell>` +object. These are denoted as arbor-components. Arbor-components need to be accompanied by *meta-data* +specifying the version of the format being used. The only version currently supported is ``0.1-dev``. + +.. label:: (version val:string) + + Specifies that the version of the component description format is ``val``. + +.. label:: (meta-data v:version) + + Add the version information ``v`` to the meta-data of the described component. + +.. label:: (arbor-component data:meta-data comp:decor/label-dict/morphology/cable-cell) + + Associates the component ``comp`` with meta-data ``data``. + +The final form of each arbor-component looks as follows: + +Label-dict +^^^^^^^^^^ + +.. code:: lisp + + (arbor-component + (meta-data (version "0.1-dev")) + (label-dict + (region-def "my_soma" (tag 1)) + (locset-def "root" (root)))) + +Decoration +^^^^^^^^^^ + +.. code:: lisp + + (arbor-component + (meta-data (version "0.1-dev")) + (decor + (default (membrane-potential -55.000000)) + (place (locset "root") (mechanism "expsyn")) + (paint (region "my_soma") (temperature-kelvin 270)))) + +Morphology +^^^^^^^^^^ + +.. code:: lisp + + (arbor-component + (meta-data (version "0.1-dev")) + (morphology + (branch 0 -1 + (segment 0 (point 0 0 0 2) (point 4 0 0 2) 1) + (segment 1 (point 4 0 0 0.8) (point 8 0 0 0.8) 3) + (segment 2 (point 8 0 0 0.8) (point 12 -0.5 0 0.8) 3)))) + +Cable-cell +^^^^^^^^^^ + +.. code:: lisp + + (arbor-component + (meta-data (version "0.1-dev")) + (cable-cell + (label-dict + (region-def "my_soma" (tag 1)) + (locset-def "root" (root))) + (decor + (default (membrane-potential -55.000000)) + (place (locset "root") (mechanism "expsyn")) + (paint (region "my_soma") (temperature-kelvin 270))) + (morphology + (branch 0 -1 + (segment 0 (point 0 0 0 2) (point 4 0 0 2) 1) + (segment 1 (point 4 0 0 0.8) (point 8 0 0 0.8) 3) + (segment 2 (point 8 0 0 0.8) (point 12 -0.5 0 0.8) 3))))) + +API +--- + +* :ref:`Python <pycablecellformat>` +* :ref:`C++ <cppcablecellformat>` \ No newline at end of file diff --git a/doc/index.rst b/doc/index.rst index 944841201e1aace4b947799a80de9b72af827235..a181dd529c977a14b49f9e12d9c14937c2980836 100644 --- a/doc/index.rst +++ b/doc/index.rst @@ -109,6 +109,7 @@ Arbor is an `eBrains project <https://ebrains.eu/service/arbor/>`_. fileformat/neuroml fileformat/asc fileformat/nmodl + fileformat/cable_cell .. toctree:: :caption: API reference: diff --git a/doc/python/cable_cell.rst b/doc/python/cable_cell.rst index 75a6144b5a56b9e5c6ec65959c264f93b2409bc7..86fe945d289fe5df54bbd40ebcfe9c2ac0e0156a 100644 --- a/doc/python/cable_cell.rst +++ b/doc/python/cable_cell.rst @@ -3,6 +3,16 @@ Cable cells =========== +.. toctree:: + :maxdepth: 1 + + morphology + labels + mechanisms + decor + probe_sample + cable_cell_format + .. currentmodule:: arbor .. py:class:: cable_cell @@ -68,14 +78,3 @@ Cable cells .. py:class:: ion properties of an ionic species. - -.. _pycablecell-probes: - -.. toctree:: - :maxdepth: 1 - - morphology - labels - mechanisms - decor - probe_sample diff --git a/doc/python/cable_cell_format.rst b/doc/python/cable_cell_format.rst new file mode 100644 index 0000000000000000000000000000000000000000..a8a10822ca27185140f2f6b2f80665c121793689 --- /dev/null +++ b/doc/python/cable_cell_format.rst @@ -0,0 +1,80 @@ +.. _pycablecellformat: + +Description Format +================== + +Arbor provides readers and writers for describing :ref:`label dictionaries <labels>`, +:ref:`decoration objects <cablecell-decoration>`, :ref:`morphologies <morph>` and +:ref:`cable cells <cablecell>`, referred to here as *arbor-components*. + +A detailed description of the s-expression format used to describe each of these components +can be found :ref:`here <formatcablecell>`. + +The arbor-components and meta-data +---------------------------------- +.. currentmodule:: arbor + +.. py:class:: meta_data + + .. py:attribute:: string version + + Stores the version of the format being used. + +.. py:class:: cable_cell_component + + .. py:attribute:: meta_data meta + + Stores meta-data pertaining to the description of a cable cell component. + + .. py:attribute:: component + + Stores one of :class:`decor`, :class:`label_dict`, :class:`morphology` or :class:`cable_cell`. + +Reading and writing arbor-components +------------------------------------ + +.. py:function:: load_component(filename) + + Load :class:`cable_cell_component` (decor, morphology, label_dict, cable_cell) from file. + + :param str filename: the name of the file containing the component description. + :rtype: :class:`cable_cell_component` + +.. py:function:: write_component(comp, filename) + + Write the :class:`cable_cell_component` to file. + + :param cable_cell_component comp: the component to be written to file. + :param str filename: the name of the file. + +.. py:function:: write_component(dec, filename) + :noindex: + + Write the :class:`decor` to file. Use the most recent version of the cable cell format to construct the meta-data. + + :param decor dec: the decor to be written to file. + :param str filename: the name of the file. + +.. py:function:: write_component(dict, filename) + :noindex: + + Write the :class:`label_dict` to file. Use the most recent version of the cable cell format to construct the meta-data. + + :param label_dict dict: the label dictionary to be written to file. + :param str filename: the name of the file. + +.. py:function:: write_component(morpho, filename) + :noindex: + + Write the :class:`morphology` to file. Use the most recent version of the cable cell format to construct the meta-data. + + :param morphology morpho: the morphology to be written to file. + :param str filename: the name of the file. + +.. py:function:: write_component(cell, filename) + :noindex: + + Write the :class:`cable_cell` to file. Use the most recent version of the cable cell format to construct the meta-data. + + :param cable_cell cell: the cable_cell to be written to file. + :param str filename: the name of the file. \ No newline at end of file diff --git a/python/CMakeLists.txt b/python/CMakeLists.txt index 8ecc727aa81cf1242d4b7bfb3c1d2e21c00b4c71..1c44f8d03ac12bfb3c1e8676a77212eefd6c44f0 100644 --- a/python/CMakeLists.txt +++ b/python/CMakeLists.txt @@ -20,6 +20,7 @@ else() endif() set(pyarb_source + cable_cell_io.cpp cable_probes.cpp cells.cpp config.cpp diff --git a/python/cable_cell_io.cpp b/python/cable_cell_io.cpp new file mode 100644 index 0000000000000000000000000000000000000000..b5c1272cf5cc171460a6d6f3d7d45054e3d4fe13 --- /dev/null +++ b/python/cable_cell_io.cpp @@ -0,0 +1,106 @@ +#include <pybind11/pybind11.h> +#include <pybind11/stl.h> + +#include <fstream> +#include <iomanip> + +#include <arbor/cable_cell.hpp> + +#include <arborio/cableio.hpp> + +#include "error.hpp" +#include "strprintf.hpp" + +namespace pyarb { + +arborio::cable_cell_component load_component(const std::string& fname) { + std::ifstream fid{fname}; + if (!fid.good()) { + throw pyarb_error("Can't open file '{}'" + fname); + } + auto component = arborio::parse_component(fid); + if (!component) { + throw pyarb_error("Error while trying to load component from \"" + fname + "\": " + component.error().what()); + } + return component.value(); +}; + +template<typename T> +void write_component(const T& component, const std::string& fname) { + std::ofstream fid(fname); + arborio::write_component(fid, component, arborio::meta_data{}); +} + +void write_component(const arborio::cable_cell_component& component, const std::string& fname) { + std::ofstream fid(fname); + arborio::write_component(fid, component); +} + +void register_cable_loader(pybind11::module& m) { + m.def("load_component", + &load_component, + pybind11::arg_v("filename", "the name of the file."), + "Load arbor-component (decor, morphology, label_dict, cable_cell) from file."); + + m.def("write_component", + [](const arborio::cable_cell_component& d, const std::string& fname) { + return write_component(d, fname); + }, + pybind11::arg_v("object", "the cable_component object."), + pybind11::arg_v("filename", "the name of the file."), + "Write cable_component to file."); + + m.def("write_component", + [](const arb::decor& d, const std::string& fname) { + return write_component<arb::decor>(d, fname); + }, + pybind11::arg_v("object", "the decor object."), + pybind11::arg_v("filename", "the name of the file."), + "Write decor to file."); + + m.def("write_component", + [](const arb::label_dict& d, const std::string& fname) { + return write_component<arb::label_dict>(d, fname); + }, + pybind11::arg_v("object", "the label_dict object."), + pybind11::arg_v("filename", "the name of the file."), + "Write label_dict to file."); + + m.def("write_component", + [](const arb::morphology& d, const std::string& fname) { + return write_component<arb::morphology>(d, fname); + }, + pybind11::arg_v("object", "the morphology object."), + pybind11::arg_v("filename", "the name of the file."), + "Write morphology to file."); + + m.def("write_component", + [](const arb::cable_cell& d, const std::string& fname) { + return write_component<arb::cable_cell>(d, fname); + }, + pybind11::arg_v("object", "the cable_cell object."), + pybind11::arg_v("filename", "the name of the file."), + "Write cable_cell to file."); + + // arborio::meta_data + pybind11::class_<arborio::meta_data> component_meta_data(m, "component_meta_data"); + component_meta_data + .def_readwrite("version", &arborio::meta_data::version, "cable-cell component version."); + + // arborio::cable_cell_component + pybind11::class_<arborio::cable_cell_component> cable_component(m, "cable_component"); + cable_component + .def_readwrite("meta_data", &arborio::cable_cell_component::meta, "cable-cell component meta-data.") + .def_readwrite("component", &arborio::cable_cell_component::component, "cable-cell component.") + .def("__repr__", [](const arborio::cable_cell_component& comp) { + std::stringstream stream; + arborio::write_component(stream, comp); + return "<arbor.cable_component>\n"+stream.str(); + }) + .def("__str__", [](const arborio::cable_cell_component& comp) { + std::stringstream stream; + arborio::write_component(stream, comp); + return "<arbor.cable_component>\n"+stream.str(); + }); +} +} diff --git a/python/pyarb.cpp b/python/pyarb.cpp index 6d4146c01f844892833f89884d68a990f54851da..4c4280885ad8c8029ec2686c8f1d62fb7e33c90e 100644 --- a/python/pyarb.cpp +++ b/python/pyarb.cpp @@ -11,6 +11,7 @@ // types and functions to be exposed to Python. namespace pyarb { +void register_cable_loader(pybind11::module& m); void register_cable_probes(pybind11::module& m, pyarb_global_ptr); void register_cells(pybind11::module& m); void register_config(pybind11::module& m); @@ -43,6 +44,7 @@ PYBIND11_MODULE(_arbor, m) { m.doc() = "arbor: multicompartment neural network models."; m.attr("__version__") = ARB_VERSION; + pyarb::register_cable_loader(m); pyarb::register_cable_probes(m, global_ptr); pyarb::register_cells(m); pyarb::register_config(m); diff --git a/test/unit/CMakeLists.txt b/test/unit/CMakeLists.txt index 7333b9ae262ceaa3aa714ff54982086767c331ae..8f9d5296e88ddee1ef6e1f64d5310a5f04f7f8e0 100644 --- a/test/unit/CMakeLists.txt +++ b/test/unit/CMakeLists.txt @@ -239,4 +239,4 @@ target_compile_options(unit PRIVATE ${ARB_CXXOPT_ARCH}) target_compile_definitions(unit PRIVATE "-DDATADIR=\"${CMAKE_CURRENT_SOURCE_DIR}/swc\"") target_compile_definitions(unit PRIVATE "-DLIBDIR=\"${PROJECT_BINARY_DIR}/lib\"") target_include_directories(unit PRIVATE "${CMAKE_CURRENT_BINARY_DIR}") -target_link_libraries(unit PRIVATE gtest arbor arborenv arborio arbor-private-headers arbor-sup) +target_link_libraries(unit PRIVATE gtest arbor arborenv arborio arborio-private-headers arbor-private-headers arbor-sup) diff --git a/test/unit/test_cable_cell.cpp b/test/unit/test_cable_cell.cpp index df2fc66ad91d06eb5ebb441a6e0db626d77c4091..d57a0caa59135b938f152a1f1f0d3ea29c0d62fb 100644 --- a/test/unit/test_cable_cell.cpp +++ b/test/unit/test_cable_cell.cpp @@ -1,8 +1,6 @@ #include "../gtest.h" #include "common_cells.hpp" -#include "s_expr.hpp" - #include <arbor/cable_cell.hpp> #include <arbor/cable_cell_param.hpp> #include <arbor/string_literals.hpp> diff --git a/test/unit/test_s_expr.cpp b/test/unit/test_s_expr.cpp index 13bd7b7b303cf385d8cce21f8e2f8fe8cc988b56..543f77796ad482ece56f4232ab46692b9eae4854 100644 --- a/test/unit/test_s_expr.cpp +++ b/test/unit/test_s_expr.cpp @@ -1,13 +1,16 @@ #include <any> +#include <typeinfo> #include "../test/gtest.h" #include <arbor/morph/region.hpp> #include <arbor/morph/locset.hpp> #include <arbor/morph/label_parse.hpp> -#include <typeinfo> +#include <arbor/s_expr.hpp> + +#include <arborio/cableio.hpp> -#include "s_expr.hpp" +#include "parse_s_expr.hpp" #include "util/strprintf.hpp" using namespace arb; @@ -75,6 +78,92 @@ TEST(s_expr, atoms_in_parens) { } } +TEST(s_expr, list) { + using namespace arborio; + auto to_string = [](const s_expr& obj) { + std::stringstream s; + s << obj; + return s.str(); + }; + { + auto l = slist(); + EXPECT_EQ("()", to_string(l)); + EXPECT_EQ(0u, std::distance(l.begin(), l.end())); + } + { + auto l = slist(1); + EXPECT_EQ("(1)", to_string(l)); + EXPECT_EQ(1u, std::distance(l.begin(), l.end())); + } + { + auto l = slist(1, 2.3); + EXPECT_EQ("(1 2.300000)", to_string(l)); + EXPECT_EQ(2u, std::distance(l.begin(), l.end())); + } + { + auto l = slist(1, 2.3, s_expr("hello")); + EXPECT_EQ("(1 2.300000 \"hello\")", to_string(l)); + EXPECT_EQ(3u, std::distance(l.begin(), l.end())); + } + { + auto l = slist(1, 2.3, "hello"_symbol); + EXPECT_EQ("(1 2.300000 hello)", to_string(l)); + EXPECT_EQ(3u, std::distance(l.begin(), l.end())); + EXPECT_EQ(tok::symbol, (l.begin()+2)->atom().kind); + } + { + auto l = slist(1, slist(1, 2), 3); + EXPECT_EQ("(1 \n (1 2)\n 3)", to_string(l)); + EXPECT_EQ(3u, std::distance(l.begin(), l.end())); + EXPECT_EQ(tok::integer, (l.begin()+2)->atom().kind); + auto l1 = *(l.begin()+1); + EXPECT_EQ(2u, std::distance(l1.begin(), l1.end())); + } +} + +TEST(s_expr, list_range) { + using namespace arborio; + auto to_string = [](const s_expr& obj) { + std::stringstream s; + s << obj; + return s.str(); + }; + auto s0 = slist_range(std::vector{1,2,3}); + EXPECT_EQ("(1 2 3)", to_string(s0)); + EXPECT_EQ(3u, length(s0)); + + auto s1 = slist_range(std::vector<int>{}); + EXPECT_EQ("()", to_string(s1)); + EXPECT_EQ(0u, length(s1)); + + auto s2 = slist_range(std::vector{12.1, 0.1}); + EXPECT_EQ("(12.100000 0.100000)", to_string(s2)); + EXPECT_EQ(2u, length(s2)); + + auto s3 = slist_range(slist(1, 2, s_expr("hello world"))); + EXPECT_EQ("(1 2 \"hello world\")", to_string(s3)); + EXPECT_EQ(3u, length(s3)); +} + +TEST(s_expr, iterate) { + using namespace arborio; + { + auto l = slist(); + EXPECT_EQ(0u, std::distance(l.begin(), l.end())); + } + { + auto l = slist(1, 2, 3., s_expr("hello")); + auto b = l.begin(); + auto e = l.end(); + EXPECT_EQ(4u, std::distance(b, e)); + EXPECT_EQ(tok::integer, b++->atom().kind); + EXPECT_EQ(tok::integer, b++->atom().kind); + EXPECT_EQ(tok::real, b++->atom().kind); + EXPECT_EQ(tok::string, b->atom().kind); + EXPECT_EQ("hello", b->atom().spelling); + } +} + template <typename L> std::string round_trip_label(const char* in) { if (auto x = parse_label_expression(in)) { @@ -205,3 +294,701 @@ TEST(regloc, errors) { EXPECT_FALSE(parse_label_expression(expr)); } } + +namespace arb { +namespace cable_s_expr { +template <typename T, std::size_t I = 0> +std::optional<T> eval_cast_variant(const std::any& a) { + if constexpr (I<std::variant_size_v<T>) { + using var_type = std::variant_alternative_t<I, T>; + return (typeid(var_type)==a.type())? std::any_cast<var_type>(a): eval_cast_variant<T, I+1>(a); + } + return std::nullopt; +} + +using branch = std::tuple<int, int, std::vector<arb::msegment>>; +using place_pair = std::pair<arb::locset, arb::placeable>; +using paint_pair = std::pair<arb::region, arb::paintable>; +using locset_pair = std::pair<std::string, locset>; +using region_pair = std::pair<std::string, region>; + +std::ostream& operator<<(std::ostream& o, const i_clamp& c) { + o << "(current-clamp (envelope"; + for (const auto& p: c.envelope) { + o << " (" << p.t << " " << p.amplitude << ')'; + } + return o << ") " << c.frequency << ')'; +} +std::ostream& operator<<(std::ostream& o, const threshold_detector& p) { + return o << "(threshold-detector " << p.threshold << ')'; +} +std::ostream& operator<<(std::ostream& o, const gap_junction_site& p) { + return o << "(gap-junction-site)"; +} +std::ostream& operator<<(std::ostream& o, const init_membrane_potential& p) { + return o << "(membrane-potential " << p.value << ')'; +} +std::ostream& operator<<(std::ostream& o, const temperature_K& p) { + return o << "(temperature-kelvin " << p.value << ')'; +} +std::ostream& operator<<(std::ostream& o, const axial_resistivity& p) { + return o << "(axial-resistivity " << p.value << ')'; +} +std::ostream& operator<<(std::ostream& o, const membrane_capacitance& p) { + return o << "(membrane-capacitance " << p.value << ')'; +} +std::ostream& operator<<(std::ostream& o, const init_int_concentration& p) { + return o << "(ion-internal-concentration \"" << p.ion << "\" " << p.value << ')'; +} +std::ostream& operator<<(std::ostream& o, const init_ext_concentration& p) { + return o << "(ion-external-concentration \"" << p.ion << "\" " << p.value << ')'; +} +std::ostream& operator<<(std::ostream& o, const init_reversal_potential& p) { + return o << "(ion-reversal-potential \"" << p.ion << "\" " << p.value << ')'; +} +std::ostream& operator<<(std::ostream& o, const mechanism_desc& m) { + o << "(mechanism \"" << m.name() << "\""; + for (const auto& p: m.values()) { + o << " (\"" << p.first << "\" " << p.second << ')'; + } + return o << ')'; +} +std::ostream& operator<<(std::ostream& o, const ion_reversal_potential_method& p) { + return o << "(ion-reversal-potential-method \"" << p.ion << "\" " << p.method << ')'; +} +std::ostream& operator<<(std::ostream& o, const cv_policy&) { + return o; +} +std::ostream& operator<<(std::ostream& o, const branch& b) { + o << "(branch " << std::to_string(std::get<0>(b)) << " " << std::to_string(std::get<1>(b)); + for (auto s: std::get<2>(b)) { + o << " " << s; + } + return o << ")"; +} +std::ostream& operator<<(std::ostream& o, const paint_pair& p) { + o << "(paint " << p.first << " "; + std::visit([&](auto&& x) {o << x;}, p.second); + return o << ")"; +} +std::ostream& operator<<(std::ostream& o, const place_pair& p) { + o << "(place " << p.first << " "; + std::visit([&](auto&& x) {o << x;}, p.second); + return o << ")"; +} +std::ostream& operator<<(std::ostream& o, const defaultable& p) { + o << "(default "; + std::visit([&](auto&& x) {o << x;}, p); + return o << ")"; +} +std::ostream& operator<<(std::ostream& o, const locset_pair& p) { + return o << "(locset-def \"" << p.first << "\" " << p.second << ")"; +} +std::ostream& operator<<(std::ostream& o, const region_pair& p) { + return o << "(region-def \"" << p.first << "\" " << p.second << ")"; +} + +template <typename T> +std::string to_string(const T& obj) { + std::stringstream s; + s << obj; + return s.str(); +} +std::string to_string(const arborio::cable_cell_component& c) { + std::stringstream s; + arborio::write_component(s, c); + return s.str(); +} +} // namespace +} // namespace arb + +template <typename T> +std::string round_trip(const char* in) { + using namespace cable_s_expr; + if (auto x = arborio::parse_expression(in)) { + return to_string(std::any_cast<T>(*x)); + } + else { + return x.error().what(); + } +} + +template <typename T> +std::string round_trip_variant(const char* in) { + using namespace cable_s_expr; + if (auto x = arborio::parse_expression(in)) { + std::string str; + std::visit([&](auto&& p){str = to_string(p);}, *(eval_cast_variant<T>(*x))); + return str; + } + else { + return x.error().what(); + } +} + +std::string round_trip_component(const char* in) { + using namespace cable_s_expr; + if (auto x = arborio::parse_component(in)) { + return to_string(x.value()); + } + else { + return x.error().what(); + } +} + +std::string round_trip_component(std::istream& stream) { + using namespace cable_s_expr; + if (auto x = arborio::parse_component(stream)) { + return to_string(x.value()); + } + else { + return x.error().what(); + } +} + + +TEST(decor_literals, round_tripping) { + auto paint_default_literals = { + "(membrane-potential -65.1)", + "(temperature-kelvin 301)", + "(axial-resistivity 102)", + "(membrane-capacitance 0.01)", + "(ion-internal-concentration \"ca\" 75.1)", + "(ion-external-concentration \"h\" -50.1)", + "(ion-reversal-potential \"na\" 30)"}; + auto paint_literals = { + "(mechanism \"hh\")", + "(mechanism \"pas\" (\"g\" 0.02))", + }; + auto default_literals = { + "(ion-reversal-potential-method \"ca\" (mechanism \"nernst/ca\"))"}; + auto place_literals = { + "(current-clamp (envelope (10 0.5) (110 0.5) (110 0)) 0)", + "(threshold-detector -10)", + "(gap-junction-site)", + "(mechanism \"expsyn\")"}; + for (auto l: paint_default_literals) { + EXPECT_EQ(l, round_trip_variant<defaultable>(l)); + EXPECT_EQ(l, round_trip_variant<paintable>(l)); + } + for (auto l: paint_literals) { + EXPECT_EQ(l, round_trip_variant<paintable>(l)); + } + for (auto l: default_literals) { + EXPECT_EQ(l, round_trip_variant<defaultable>(l)); + } + for (auto l: place_literals) { + EXPECT_EQ(l, round_trip_variant<placeable>(l)); + } + auto clamp_literal = "(current-clamp (envelope-pulse 10 5 0.1) 50)"; + EXPECT_EQ("(current-clamp (envelope (10 0.1) (15 0.1) (15 0)) 50)", round_trip_variant<placeable>(clamp_literal)); + + std::string mech_str = "(mechanism \"kamt\" (\"gbar\" 50) (\"zetam\" 0.1) (\"q10\" 5))"; + auto maybe_mech = arborio::parse_expression(mech_str); + EXPECT_TRUE(maybe_mech); + + auto any_mech = maybe_mech.value(); + EXPECT_TRUE(typeid(mechanism_desc) == any_mech.type()); + + auto mech = std::any_cast<mechanism_desc>(any_mech); + EXPECT_EQ("kamt", mech.name()); + EXPECT_EQ(3u, mech.values().size()); + + EXPECT_EQ(50, mech.values().at("gbar")); + EXPECT_EQ(0.1, mech.values().at("zetam")); + EXPECT_EQ(5, mech.values().at("q10")); + +} + +TEST(decor_expressions, round_tripping) { + using namespace cable_s_expr; + auto decorate_paint_literals = { + "(paint (region \"all\") (membrane-potential -65.1))", + "(paint (tag 1) (temperature-kelvin 301))", + "(paint (distal-interval (location 3 0)) (axial-resistivity 102))", + "(paint (join (region \"dend\") (all)) (membrane-capacitance 0.01))", + "(paint (radius-gt (tag 3) 1) (ion-internal-concentration \"ca\" 75.1))", + "(paint (intersect (cable 2 0 0.5) (region \"axon\")) (ion-external-concentration \"h\" -50.1))", + "(paint (region \"my_region\") (ion-reversal-potential \"na\" 30))", + "(paint (cable 2 0.1 0.4) (mechanism \"hh\"))", + "(paint (all) (mechanism \"pas\" (\"g\" 0.02)))" + }; + auto decorate_default_literals = { + "(default (membrane-potential -65.1))", + "(default (temperature-kelvin 301))", + "(default (axial-resistivity 102))", + "(default (membrane-capacitance 0.01))", + "(default (ion-internal-concentration \"ca\" 75.1))", + "(default (ion-external-concentration \"h\" -50.1))", + "(default (ion-reversal-potential \"na\" 30))", + "(default (ion-reversal-potential-method \"ca\" (mechanism \"nernst/ca\")))" + }; + auto decorate_place_literals = { + "(place (location 3 0.2) (current-clamp (envelope (10 0.5) (110 0.5) (110 0)) 0))", + "(place (terminal) (threshold-detector -10))", + "(place (root) (gap-junction-site))", + "(place (locset \"my!ls\") (mechanism \"expsyn\"))"}; + + for (auto l: decorate_paint_literals) { + EXPECT_EQ(l, round_trip<paint_pair>(l)); + } + for (auto l: decorate_place_literals) { + EXPECT_EQ(l, round_trip<place_pair>(l)); + } + for (auto l: decorate_default_literals) { + EXPECT_EQ(l, round_trip<defaultable>(l)); + } +} + +TEST(label_dict_expressions, round_tripping) { + using namespace cable_s_expr; + auto locset_def_literals = { + "(locset-def \"my! locset~\" (root))", + "(locset-def \"ls0\" (location 0 1))" + }; + auto region_def_literals = { + "(region-def \"my region\" (all))", + "(region-def \"reg42\" (cable 0 0.1 0.9))" + }; + for (auto l: locset_def_literals) { + EXPECT_EQ(l, round_trip<locset_pair>(l)); + } + for (auto l: region_def_literals) { + EXPECT_EQ(l, round_trip<region_pair>(l)); + } +} + +TEST(morphology_literals, round_tripping) { + using namespace cable_s_expr; + auto point = "(point 701.6 -3.1 -10 0.6)"; + auto segment = "(segment 5 (point 5 2 3 1) (point 5 2 5 6) 42)"; + auto branches = { + "(branch -1 5 (segment 5 (point 5 2 3 1) (point 5 2 3.1 0.5) 2))", + "(branch -1 5" + " (segment 2 (point 5 2 3 1) (point 5 2 5 6) 42)" + " (segment 3 (point 5 2 3 1) (point 5 2 3.1 0.5) 2)" + ")" + }; + + EXPECT_EQ(point, round_trip<mpoint>(point)); + EXPECT_EQ(segment, round_trip<msegment>(segment)); + for (auto l: branches) { + EXPECT_EQ(l, round_trip<branch>(l)); + } +} + +TEST(decor, round_tripping) { + std::string component_str = "(arbor-component \n" + " (meta-data \n" + " (version \"" + arborio::acc_version() +"\"))\n" + " (decor \n" + " (default \n" + " (axial-resistivity 100.000000))\n" + " (default \n" + " (ion-reversal-potential-method \"na\" \n" + " (mechanism \"nernst\")))\n" + " (paint \n" + " (region \"dend\")\n" + " (mechanism \"pas\"))\n" + " (paint \n" + " (region \"soma\")\n" + " (mechanism \"hh\"))\n" + " (paint \n" + " (join \n" + " (tag 1)\n" + " (tag 2))\n" + " (ion-internal-concentration \"ca\" 0.500000))\n" + " (place \n" + " (location 0 0)\n" + " (threshold-detector 10.000000))\n" + " (place \n" + " (location 0 0.5)\n" + " (mechanism \"expsyn\" \n" + " (\"tau\" 1.500000)))))"; + + EXPECT_EQ(component_str, round_trip_component(component_str.c_str())); +} + +TEST(label_dict, round_tripping) { + std::string component_str = "(arbor-component \n" + " (meta-data \n" + " (version \"" + arborio::acc_version() + "\"))\n" + " (label-dict \n" + " (region-def \"soma\" \n" + " (tag 1))\n" + " (region-def \"dend\" \n" + " (tag 3))\n" + " (locset-def \"root\" \n" + " (root))))"; + + EXPECT_EQ(component_str, round_trip_component(component_str.c_str())); +} + +TEST(morphology, round_tripping) { + std::string component_str = "(arbor-component \n" + " (meta-data \n" + " (version \"" + arborio::acc_version() +"\"))\n" + " (morphology \n" + " (branch 0 -1 \n" + " (segment 0 \n" + " (point 0.000000 0.000000 -6.307850 6.307850)\n" + " (point 0.000000 0.000000 6.307850 6.307850)\n" + " 1))\n" + " (branch 1 0 \n" + " (segment 1 \n" + " (point 0.000000 0.000000 6.307850 6.307850)\n" + " (point 0.000000 0.000000 72.974517 0.500000)\n" + " 3)\n" + " (segment 2 \n" + " (point 0.000000 0.000000 72.974517 0.500000)\n" + " (point 0.000000 0.000000 139.641183 0.500000)\n" + " 3)\n" + " (segment 3 \n" + " (point 0.000000 0.000000 139.641183 0.500000)\n" + " (point 0.000000 0.000000 206.307850 0.500000)\n" + " 3))\n" + " (branch 2 0 \n" + " (segment 4 \n" + " (point 0.000000 0.000000 6.307850 6.307850)\n" + " (point 0.000000 0.000000 72.974517 0.500000)\n" + " 3)\n" + " (segment 5 \n" + " (point 0.000000 0.000000 72.974517 0.500000)\n" + " (point 0.000000 0.000000 139.641183 0.500000)\n" + " 3)\n" + " (segment 6 \n" + " (point 0.000000 0.000000 139.641183 0.500000)\n" + " (point 0.000000 0.000000 206.307850 0.500000)\n" + " 3))\n" + " (branch 3 2 \n" + " (segment 7 \n" + " (point 0.000000 0.000000 206.307850 0.500000)\n" + " (point 0.000000 0.000000 257.974517 0.500000)\n" + " 3)\n" + " (segment 8 \n" + " (point 0.000000 0.000000 257.974517 0.500000)\n" + " (point 0.000000 0.000000 309.641183 0.500000)\n" + " 3)\n" + " (segment 9 \n" + " (point 0.000000 0.000000 309.641183 0.500000)\n" + " (point 0.000000 0.000000 361.307850 0.500000)\n" + " 3))\n" + " (branch 4 2 \n" + " (segment 10 \n" + " (point 0.000000 0.000000 206.307850 0.500000)\n" + " (point 0.000000 0.000000 257.974517 0.500000)\n" + " 3)\n" + " (segment 11 \n" + " (point 0.000000 0.000000 257.974517 0.500000)\n" + " (point 0.000000 0.000000 309.641183 0.500000)\n" + " 3)\n" + " (segment 12 \n" + " (point 0.000000 0.000000 309.641183 0.500000)\n" + " (point 0.000000 0.000000 361.307850 0.500000)\n" + " 3))\n" + " (branch 5 3 \n" + " (segment 13 \n" + " (point 0.000000 0.000000 361.307850 0.500000)\n" + " (point 0.000000 0.000000 416.307850 0.500000)\n" + " 3)\n" + " (segment 14 \n" + " (point 0.000000 0.000000 416.307850 0.500000)\n" + " (point 0.000000 0.000000 471.307850 0.500000)\n" + " 3)\n" + " (segment 21 \n" + " (point 0.000000 0.000000 471.307850 0.500000)\n" + " (point 0.000000 0.000000 503.807850 0.500000)\n" + " 3)\n" + " (segment 22 \n" + " (point 0.000000 0.000000 503.807850 0.500000)\n" + " (point 0.000000 0.000000 536.307850 0.500000)\n" + " 3)\n" + " (segment 29 \n" + " (point 0.000000 0.000000 536.307850 0.500000)\n" + " (point 0.000000 0.000000 556.307850 0.500000)\n" + " 3))\n" + " (branch 6 3 \n" + " (segment 15 \n" + " (point 0.000000 0.000000 361.307850 0.500000)\n" + " (point 0.000000 0.000000 416.307850 0.500000)\n" + " 3)\n" + " (segment 16 \n" + " (point 0.000000 0.000000 416.307850 0.500000)\n" + " (point 0.000000 0.000000 471.307850 0.500000)\n" + " 3)\n" + " (segment 23 \n" + " (point 0.000000 0.000000 471.307850 0.500000)\n" + " (point 0.000000 0.000000 503.807850 0.500000)\n" + " 3)\n" + " (segment 24 \n" + " (point 0.000000 0.000000 503.807850 0.500000)\n" + " (point 0.000000 0.000000 536.307850 0.500000)\n" + " 3))\n" + " (branch 7 4 \n" + " (segment 17 \n" + " (point 0.000000 0.000000 361.307850 0.500000)\n" + " (point 0.000000 0.000000 416.307850 0.500000)\n" + " 3)\n" + " (segment 18 \n" + " (point 0.000000 0.000000 416.307850 0.500000)\n" + " (point 0.000000 0.000000 471.307850 0.500000)\n" + " 3)\n" + " (segment 25 \n" + " (point 0.000000 0.000000 471.307850 0.500000)\n" + " (point 0.000000 0.000000 503.807850 0.500000)\n" + " 3)\n" + " (segment 26 \n" + " (point 0.000000 0.000000 503.807850 0.500000)\n" + " (point 0.000000 0.000000 536.307850 0.500000)\n" + " 3))\n" + " (branch 8 4 \n" + " (segment 19 \n" + " (point 0.000000 0.000000 361.307850 0.500000)\n" + " (point 0.000000 0.000000 416.307850 0.500000)\n" + " 3)\n" + " (segment 20 \n" + " (point 0.000000 0.000000 416.307850 0.500000)\n" + " (point 0.000000 0.000000 471.307850 0.500000)\n" + " 3)\n" + " (segment 27 \n" + " (point 0.000000 0.000000 471.307850 0.500000)\n" + " (point 0.000000 0.000000 503.807850 0.500000)\n" + " 3)\n" + " (segment 28 \n" + " (point 0.000000 0.000000 503.807850 0.500000)\n" + " (point 0.000000 0.000000 536.307850 0.500000)\n" + " 3))))"; + + EXPECT_EQ(component_str, round_trip_component(component_str.c_str())); +} + +TEST(morphology, invalid) { + auto component_str = "(morphology\n" + " (branch 0 -1\n" + " (segment 0 (point 0 0 0 2) (point 4 0 0 2) 1))\n" + " (branch 1 0\n" + " (segment 1 (point 4 0 0 0.8) (point 8 0 0 0.8) 3)\n" + " (segment 2 (point 8 0 0 0.8) (point 12 -0.5 0 0.8) 3))\n" + ")"; + EXPECT_THROW(round_trip_component(component_str), arborio::cableio_morphology_error); +} + +TEST(cable_cell, round_tripping) { + std::string component_str = "(arbor-component \n" + " (meta-data \n" + " (version \"" + arborio::acc_version() +"\"))\n" + " (cable-cell \n" + " (morphology \n" + " (branch 0 -1 \n" + " (segment 0 \n" + " (point -6.300000 0.000000 0.000000 6.300000)\n" + " (point 6.300000 0.000000 0.000000 6.300000)\n" + " 1)\n" + " (segment 1 \n" + " (point 6.300000 0.000000 0.000000 0.500000)\n" + " (point 206.300000 0.000000 0.000000 0.200000)\n" + " 3)))\n" + " (label-dict \n" + " (region-def \"soma\" \n" + " (tag 1))\n" + " (region-def \"dend\" \n" + " (join \n" + " (join \n" + " (tag 3)\n" + " (tag 4))\n" + " (tag 42))))\n" + " (decor \n" + " (paint \n" + " (region \"dend\")\n" + " (mechanism \"pas\"))\n" + " (paint \n" + " (region \"soma\")\n" + " (mechanism \"hh\"))\n" + " (place \n" + " (location 0 1)\n" + " (mechanism \"exp2syn\")))))"; + + EXPECT_EQ(component_str, round_trip_component(component_str.c_str())); + + std::stringstream stream(component_str); + EXPECT_EQ(component_str, round_trip_component(stream)); +} + +TEST(cable_cell_literals, errors) { + for (auto expr: {"(membrane-potential \"56\")", // invalid argument + "(axial-resistivity 1 2)", // too many arguments to otherwise valid decor literal + "(membrane-capacitance ", // syntax error + "(mem-capacitance 3.5)", // invalid function + "(ion-initial-concentration ca 0.1)", // unquoted ion + "(mechanism \"hh\" (gl 3.5))", // unqouted parameter + "(mechanism \"pas\" ((\"g\" 0.5) (\"e\" 0.2)))", // paranthesis around params + "(mechanism \"pas\" (\"g\" 0.5 0.1) (\"e\" 0.2))", // too many values + "(gap-junction-site 0)", // too many arguments + "(current-clamp (envelope (10 0.5) (110 0.5) (110 0)))", // too few arguments + "(paint (region) (mechanism \"hh\"))", // invalid region + "(paint (tag 1) (mechanims hh))", // invalid painting + "(paint (terminal) (membrance-capacitance 0.2))", // can't paint a locset + "(paint (tag 3))", // too few arguments + "(place (locset) (gap-junction-site))", // invalid locset + "(place (gap-junction-site) (location 0 1))", // swapped argument order + "(region-def my_region (tag 3))", // unquoted region name + "(locset-def \"my_ls\" (tag 3))", // invalid locset + "(locset-def \"my_ls\")", // too few arguments + "(branch 0 1)", // branch with zero segments + "(segment -1 (point 1 2 3 4) 3)", // segment with 1 point + "(point 1 2 3)", // too few arguments + "(morphology (segment -1 (point 1 2 3 4) (point 2 3 4 5) 3))", // morphology with segment instead of branch + "(decor (region-def \"reg\" (tag 3)))", // decor with label definiton + "(cable-cell (paint (tag 3) (axial-resistivity 3.1)))" // cable-cell with paint + }) + { + // If an expression was passed and handled correctly the parse call will + // return without throwing, with the error stored in the return type. + // So it is sufficient to assert that it evaluates to false. + EXPECT_FALSE(arborio::parse_expression(expr)); + } + for (const std::string& expr: std::vector<std::string>{ + "(arbor-component (meta-data (version \"0.1\")) (decor))", // invalid component + "(arbor-component (morphology))", // arbor-component missing meta-data + "(arbor-component (meta-data (version \"" + arborio::acc_version() +"\")))", // arbor-component missing component + "(arbor-component (meta-data (version \"" + arborio::acc_version() +"\")) (membrane-potential 56))", // invalid component + "(arbor-component (meta-data (version \"" + arborio::acc_version() +"\")) (morphology (segment 1 (point 1 2 3 4) (point 2 3 4 5) 3)))", // morphology with segment instead of branch + "(arbor-component (meta-data (version \"" + arborio::acc_version() +"\")) (decor (region-def \"reg\" (tag 3))))", // decor with label definition + "(arbor-component (meta-data (version \"" + arborio::acc_version() +"\")) (cable-cell (paint (tag 3) (axial-resistivity 3.1))))", // cable-cell with paint + "(morphology (branch 0 -1 (segment 0 (point 0 1 2 3 ) (point 1 2 3 4) 3)))" // morphology without arbor-component + }) + { + // If an expression was passed and handled correctly the parse call will + // return without throwing, with the error stored in the return type. + // So it is sufficient to assert that it evaluates to false. + EXPECT_FALSE(arborio::parse_component(expr)); + } +} + +// Check that the examples used in the docs are valid (formats/cable_cell.rst) +TEST(doc_expressions, parse) { + // literals + for (auto expr: {"(region-def \"my_region\" (branch 1))", + "(locset-def \"my_locset\" (location 3 0.5))", + "(mechanism \"hh\" (\"gl\" 0.5) (\"el\" 2))", + "(ion-reversal-potential-method \"ca\" (mechanism \"nernst/ca\"))", + "(current-clamp (envelope (0 10) (50 10) (50 0)) 40)", + "(paint (tag 1) (membrane-capacitance 0.02))", + "(place (locset \"mylocset\") (threshold-detector 10))", + "(default (membrane-potential -65))", + "(segment 3 (point 0 0 0 5) (point 0 0 10 2) 1)"}) + { + EXPECT_TRUE(arborio::parse_expression(expr)); + } + + // objects + for (auto expr: {"(label-dict" + " (region-def \"my_soma\" (tag 1))\n" + " (locset-def \"root\" (root))\n" + " (region-def \"all\" (all))\n" + " (region-def \"my_region\" (radius-ge (region \"my_soma\") 1.5))\n" + " (locset-def \"terminal\" (terminal)))", + "(decor\n" + " (default (membrane-potential -55.000000))\n" + " (paint (region \"custom\") (temperature-kelvin 270))\n" + " (paint (region \"soma\") (membrane-potential -50.000000))\n" + " (paint (all) (mechanism \"pas\"))\n" + " (paint (tag 4) (mechanism \"Ih\" (\"gbar\" 0.001)))\n" + " (place (locset \"root\") (mechanism \"expsyn\"))\n" + " (place (terminal) (gap-junction-site)))", + "(morphology\n" + " (branch 0 -1\n" + " (segment 0 (point 0 0 0 2) (point 4 0 0 2) 1)\n" + " (segment 1 (point 4 0 0 0.8) (point 8 0 0 0.8) 3)\n" + " (segment 2 (point 8 0 0 0.8) (point 12 -0.5 0 0.8) 3))\n" + " (branch 1 0\n" + " (segment 3 (point 12 -0.5 0 0.8) (point 20 4 0 0.4) 3)\n" + " (segment 4 (point 20 4 0 0.4) (point 26 6 0 0.2) 3))\n" + " (branch 2 0\n" + " (segment 5 (point 12 -0.5 0 0.5) (point 19 -3 0 0.5) 3))\n" + " (branch 3 2\n" + " (segment 6 (point 19 -3 0 0.5) (point 24 -7 0 0.2) 3))\n" + " (branch 4 2\n" + " (segment 7 (point 19 -3 0 0.5) (point 23 -1 0 0.2) 3)\n" + " (segment 8 (point 23 -1 0 0.3) (point 26 -2 0 0.2) 3))\n" + " (branch 5 -1\n" + " (segment 9 (point 0 0 0 2) (point -7 0 0 0.4) 2)\n" + " (segment 10 (point -7 0 0 0.4) (point -10 0 0 0.4) 2)))", + "(cable-cell\n" + " (label-dict\n" + " (region-def \"my_soma\" (tag 1))\n" + " (locset-def \"root\" (root))\n" + " (region-def \"all\" (all))\n" + " (region-def \"my_region\" (radius-ge (region \"my_soma\") 1.5))\n" + " (locset-def \"terminal\" (terminal)))\n" + " (decor\n" + " (default (membrane-potential -55.000000))\n" + " (paint (region \"my_soma\") (temperature-kelvin 270))\n" + " (paint (region \"my_region\") (membrane-potential -50.000000))\n" + " (paint (tag 4) (mechanism \"Ih\" (\"gbar\" 0.001)))\n" + " (place (locset \"root\") (mechanism \"expsyn\"))\n" + " (place (location 1 0.2) (gap-junction-site)))\n" + " (morphology\n" + " (branch 0 -1\n" + " (segment 0 (point 0 0 0 2) (point 4 0 0 2) 1)\n" + " (segment 1 (point 4 0 0 0.8) (point 8 0 0 0.8) 3)\n" + " (segment 2 (point 8 0 0 0.8) (point 12 -0.5 0 0.8) 3))\n" + " (branch 1 0\n" + " (segment 3 (point 12 -0.5 0 0.8) (point 20 4 0 0.4) 3)\n" + " (segment 4 (point 20 4 0 0.4) (point 26 6 0 0.2) 3))\n" + " (branch 2 0\n" + " (segment 5 (point 12 -0.5 0 0.5) (point 19 -3 0 0.5) 3))\n" + " (branch 3 2\n" + " (segment 6 (point 19 -3 0 0.5) (point 24 -7 0 0.2) 3))\n" + " (branch 4 2\n" + " (segment 7 (point 19 -3 0 0.5) (point 23 -1 0 0.2) 3)\n" + " (segment 8 (point 23 -1 0 0.3) (point 26 -2 0 0.2) 3))\n" + " (branch 5 -1\n" + " (segment 9 (point 0 0 0 2) (point -7 0 0 0.4) 2)\n" + " (segment 10 (point -7 0 0 0.4) (point -10 0 0 0.4) 2))))"}) + { + auto t = arborio::parse_expression(expr); + if (!t) { + std::cout << t.error().what() << std::endl; + } + EXPECT_TRUE(arborio::parse_expression(expr)); + } + + // components + for (std::string expr: {"(arbor-component\n" + " (meta-data (version \"" + arborio::acc_version() +"\"))\n" + " (label-dict\n" + " (region-def \"my_soma\" (tag 1))\n" + " (locset-def \"root\" (root))))", + "(arbor-component\n" + " (meta-data (version \"" + arborio::acc_version() +"\"))\n" + " (decor\n" + " (default (membrane-potential -55.000000))\n" + " (place (locset \"root\") (mechanism \"expsyn\"))\n" + " (paint (region \"my_soma\") (temperature-kelvin 270))))", + "(arbor-component\n" + " (meta-data (version \"" + arborio::acc_version() +"\"))\n" + " (morphology\n" + " (branch 0 -1\n" + " (segment 0 (point 0 0 0 2) (point 4 0 0 2) 1)\n" + " (segment 1 (point 4 0 0 0.8) (point 8 0 0 0.8) 3)\n" + " (segment 2 (point 8 0 0 0.8) (point 12 -0.5 0 0.8) 3))))\n", + "(arbor-component\n" + " (meta-data (version \"" + arborio::acc_version() +"\"))\n" + " (cable-cell\n" + " (label-dict\n" + " (region-def \"my_soma\" (tag 1))\n" + " (locset-def \"root\" (root)))\n" + " (decor\n" + " (default (membrane-potential -55.000000))\n" + " (place (locset \"root\") (mechanism \"expsyn\"))\n" + " (paint (region \"my_soma\") (temperature-kelvin 270)))\n" + " (morphology\n" + " (branch 0 -1\n" + " (segment 0 (point 0 0 0 2) (point 4 0 0 2) 1)\n" + " (segment 1 (point 4 0 0 0.8) (point 8 0 0 0.8) 3)\n" + " (segment 2 (point 8 0 0 0.8) (point 12 -0.5 0 0.8) 3)))))"}) + { + EXPECT_TRUE(arborio::parse_expression(expr)); + } +} \ No newline at end of file