diff --git a/arbor/include/arbor/mechcat.hpp b/arbor/include/arbor/mechcat.hpp index bbe80e6e7ad6c89a1fe722c15311bda731b3478b..17cd7ef962917e55e913af6d66c154a4a7c8ef56 100644 --- a/arbor/include/arbor/mechcat.hpp +++ b/arbor/include/arbor/mechcat.hpp @@ -75,6 +75,8 @@ public: const std::vector<std::pair<std::string, double>>& global_params, const std::vector<std::pair<std::string, std::string>>& ion_remap = {}); + void derive(const std::string& name, const std::string& parent); + // Remove mechanism from catalogue, together with any derivations of it. void remove(const std::string& name); diff --git a/arbor/mechcat.cpp b/arbor/mechcat.cpp index 8d2a9c19ab2cd5a2220d02ea289b55f3ac9d4296..59b45df7943c0f666c925ce7d0bd63d3fa4ff468 100644 --- a/arbor/mechcat.cpp +++ b/arbor/mechcat.cpp @@ -536,6 +536,10 @@ void mechanism_catalogue::derive(const std::string& name, const std::string& par state_->bind(name, value(state_->derive(name, parent, global_params, ion_remap_vec))); } +void mechanism_catalogue::derive(const std::string& name, const std::string& parent) { + state_->bind(name, value(state_->derive(parent))); +} + void mechanism_catalogue::remove(const std::string& name) { if (!has(name)) { throw no_such_mechanism(name); diff --git a/arbor/morph/region.cpp b/arbor/morph/region.cpp index 3893ab11871faa946d5f45a77be9e8ec92600fc1..90c30690c75d2f410f8a8961939bbcf71a2cf013 100644 --- a/arbor/morph/region.cpp +++ b/arbor/morph/region.cpp @@ -416,7 +416,7 @@ struct radius_ge_ { }; region radius_ge(region reg, double val) { - return region(radius_gt_{reg, val}); + return region(radius_ge_{reg, val}); } mextent thingify_(const radius_ge_& r, const mprovider& p) { diff --git a/python/CMakeLists.txt b/python/CMakeLists.txt index 0b90b6d0e47e5667457adf2a38ca1dcc5b78e060..21fe5be9123db9ec8f9c784a99a4007401770933 100644 --- a/python/CMakeLists.txt +++ b/python/CMakeLists.txt @@ -23,6 +23,7 @@ set(pyarb_source event_generator.cpp flat_cell_builder.cpp identifiers.cpp + mechanism.cpp morphology.cpp morph_parse.cpp mpi.cpp diff --git a/python/cells.cpp b/python/cells.cpp index e0ae99ac8b9c7af8b4d69e0345acd6e875218a5e..89fe1fc696267858268f66ab3b10a30f576327b7 100644 --- a/python/cells.cpp +++ b/python/cells.cpp @@ -12,6 +12,7 @@ #include <arbor/util/any.hpp> #include <arbor/util/unique_any.hpp> +#include "cells.hpp" #include "conversion.hpp" #include "error.hpp" #include "morph_parse.hpp" @@ -142,16 +143,59 @@ struct label_dict_proxy { std::string s; s += "(label_dict"; for (auto& x: dict.regions()) { - s += util::pprintf(" (region '{}' {})", x.first, x.second); + s += util::pprintf(" (region \"{}\" {})", x.first, x.second); } for (auto& x: dict.locsets()) { - s += util::pprintf(" (locset '{}' {})", x.first, x.second); + s += util::pprintf(" (locset \"{}\" {})", x.first, x.second); } s += ")"; return s; } }; +global_props_shim::global_props_shim(): + cat(arb::global_default_catalogue()) +{ + props.catalogue = &cat; +} + +// This isn't pretty. Partly because the information in the global parameters +// is all over the place. +std::string to_string(const global_props_shim& G) { + std::string s = "{arbor.cable_global_properties"; + + auto nan_is_none = [](double x) {return x==x? std::to_string(x): "None";}; + const auto& P = G.props; + const auto& D = P.default_parameters; + const auto& I = D.ion_data; + // name, valence, int_con, ext_con, rev_pot, rev_pot_method + s += "\n ions: {"; + for (auto& ion: P.ion_species) { + if (!I.count(ion.first)) { + s += util::pprintf("\n {name: '{}', valence: {}, int_con: None, ext_con: None, rev_pot: None, rev_pot_method: None}", + ion.first, ion.second); + } + else { + auto& props = I.at(ion.first); + std::string method = D.reversal_potential_method.count(ion.first)? + "'"+D.reversal_potential_method.at(ion.first).name()+"'": "None"; + s += util::pprintf("\n {name: '{}', valence: {}, int_con: {}, ext_con: {}, rev_pot: {}, rev_pot_method: {}}", + ion.first, ion.second, + nan_is_none(props.init_int_concentration), + nan_is_none(props.init_ext_concentration), + nan_is_none(props.init_reversal_potential), + method); + } + } + s += "}\n"; + s += util::pprintf(" parameters: {Vm: {}, cm: {}, rL: {}, tempK: {}}\n", + D.init_membrane_potential, D.membrane_capacitance, + D.axial_resistivity, D.temperature_K); + s += "}"; + return s; +} + + // // string printers // @@ -164,7 +208,7 @@ std::string lif_str(const arb::lif_cell& c){ std::string mechanism_desc_str(const arb::mechanism_desc& md) { - return util::pprintf("<arbor.mechanism: name '{}', parameters {}", + return util::pprintf("mechanism('{}', {})", md.name(), util::dictionary_csv(md.values())); } @@ -277,77 +321,6 @@ void register_cells(pybind11::module& m) { .def("__repr__", [](const label_dict_proxy& d){return d.to_string();}) .def("__str__", [](const label_dict_proxy& d){return d.to_string();}); - // Data structures used to describe mechanisms, electrical properties, - // gap junction properties, etc. - - // arb::cable_cell_ion_data - - pybind11::class_<arb::initial_ion_data> ion_data(m, "ion", - "For setting ion properties (internal and external concentration and reversal potential) on cells and regions."); - ion_data - .def(pybind11::init( - [](const char* name, - optional<double> int_con, optional<double> ext_con, - optional<double> rev_pot) - { - arb::initial_ion_data x; - x.ion = name; - if (int_con) x.initial.init_int_concentration = *int_con; - if (ext_con) x.initial.init_int_concentration = *ext_con; - if (rev_pot) x.initial.init_reversal_potential = *rev_pot; - return x; - } - ), - "ion_name"_a, - pybind11::arg_v("int_con", pybind11::none(), "Intial internal concentration [mM]"), - pybind11::arg_v("ext_con", pybind11::none(), "Intial external concentration [mM]"), - pybind11::arg_v("rev_pot", pybind11::none(), "Intial reversal potential [mV]"), - "If concentrations or reversal potential are specified as 'None'," - " cell default or global default value will be used, in that order if set."); - - // arb::mechanism_desc - - pybind11::class_<arb::mechanism_desc> mechanism_desc(m, "mechanism"); - mechanism_desc - .def(pybind11::init([](const char* n) {return arb::mechanism_desc{n};})) - // allow construction of a description with parameters provided in a dictionary: - // mech = arbor.mechanism('mech_name', {'param1': 1.2, 'param2': 3.14}) - .def(pybind11::init( - [](const char* name, std::unordered_map<std::string, double> params) { - arb::mechanism_desc md(name); - for (const auto& p: params) md.set(p.first, p.second); - return md; - }), - "name"_a, "The name of the mechanism", - "params"_a, "A dictionary of parameter values, with parameter name as key.", - "Example usage setting pararmeters:\n" - " m = arbor.mechanism('expsyn', {'tau': 1.4})\n" - "will create parameters for the 'expsyn' mechanism, with the provided value\n" - "for 'tau' overrides the default. If a parameter is not set, the default\n" - "(as defined in NMODL) is used.\n\n" - "Example overriding a global parameter:\n" - " m = arbor.mechanism('nernst/R=8.3145,F=96485')") - .def("set", - [](arb::mechanism_desc& md, std::string name, double value) { - md.set(name, value); - }, - "name"_a, "value"_a, "Set parameter value.") - .def_property_readonly("name", - [](const arb::mechanism_desc& md) { - return md.name(); - }, - "The name of the mechanism.") - .def_property_readonly("values", - [](const arb::mechanism_desc& md) { - return md.values(); - }, "A dictionary of parameter values with parameter name as key.") - .def("__repr__", - [](const arb::mechanism_desc& md) { - return util::pprintf("<arbor.mechanism: name '{}', parameters {}", md.name(), util::dictionary_csv(md.values())); }) - .def("__str__", - [](const arb::mechanism_desc& md) { - return util::pprintf("('{}' {})", md.name(), util::dictionary_csv(md.values())); }); - // arb::gap_junction_site pybind11::class_<arb::gap_junction_site> gjsite(m, "gap_junction", @@ -390,6 +363,102 @@ void register_cells(pybind11::module& m) { .def("__str__", [](const arb::threshold_detector& d){ return util::pprintf("(threshold_detector {})", d.threshold);}); + // arb::cable_cell_ion_data + + pybind11::class_<arb::initial_ion_data> ion_data(m, "ion", + "For setting ion properties (internal and external concentration and reversal potential) on cells and regions."); + ion_data + .def(pybind11::init( + [](const char* name, + optional<double> int_con, + optional<double> ext_con, + optional<double> rev_pot) + { + arb::initial_ion_data x; + x.ion = name; + if (int_con) x.initial.init_int_concentration = *int_con; + if (ext_con) x.initial.init_int_concentration = *ext_con; + if (rev_pot) x.initial.init_reversal_potential = *rev_pot; + return x; + } + ), + "ion_name"_a, + pybind11::arg_v("int_con", pybind11::none(), "Intial internal concentration [mM]"), + pybind11::arg_v("ext_con", pybind11::none(), "Intial external concentration [mM]"), + pybind11::arg_v("rev_pot", pybind11::none(), "Intial reversal potential [mV]"), + "If concentrations or reversal potential are specified as 'None'," + " cell default or global default value will be used, in that order if set."); + + // arb::cable_cell_global_properties + + pybind11::class_<global_props_shim> gprop(m, "cable_global_properties"); + gprop + .def(pybind11::init<>()) + .def(pybind11::init<const global_props_shim&>()) + .def("check", [](const global_props_shim& shim) { + arb::check_global_properties(shim.props);}, + "Test whether all default parameters and ion specids properties have been set.") + // set cable properties + .def("set_properties", + [](global_props_shim& G, + optional<double> Vm, optional<double> cm, + optional<double> rL, optional<double> tempK) + { + if (Vm) G.props.default_parameters.init_membrane_potential = Vm; + if (cm) G.props.default_parameters.membrane_capacitance=cm; + if (rL) G.props.default_parameters.axial_resistivity=rL; + if (tempK) G.props.default_parameters.temperature_K=tempK; + }, + "Vm"_a=pybind11::none(), "cm"_a=pybind11::none(), "rL"_a=pybind11::none(), "tempK"_a=pybind11::none(), + "Set global default values for cable and cell properties.\n" + " Vm: initial membrane voltage [mV].\n" + " cm: membrane capacitance [F/m²].\n" + " rL: axial resistivity [Ω·cm].\n" + " tempK: temperature [Kelvin].") + .def("foo", [](global_props_shim& G, double x, pybind11::object method) { + std::cout << "foo(" << x << ", " << method << ")\n";}, + "x"_a, "method"_a=pybind11::none()) + // add/modify ion species + .def("set_ion", + [](global_props_shim& G, const char* ion, + optional<double> int_con, optional<double> ext_con, + optional<double> rev_pot, pybind11::object method) + { + auto& data = G.props.default_parameters.ion_data[ion]; + if (int_con) data.init_int_concentration = *int_con; + if (ext_con) data.init_ext_concentration = *ext_con; + if (rev_pot) data.init_reversal_potential = *rev_pot; + // allow rev_pot_method to be specified with string or mechanism_desc + if (!method.is_none()) { + if (auto m=try_cast<std::string>(method)) { + G.props.default_parameters.reversal_potential_method[ion] = *m; + } + else if (auto m=try_cast<arb::mechanism_desc>(method)) { + G.props.default_parameters.reversal_potential_method[ion] = *m; + } + else { + throw std::runtime_error(util::pprintf("invalid rev_pot_method: {}", method)); + } + } + }, + "ion"_a, + "int_con"_a=pybind11::none(), + "ext_con"_a=pybind11::none(), + "rev_pot"_a=pybind11::none(), + "rev_pot_method"_a=pybind11::none(), + "Set the global default propoerties of ion species named 'ion'.\n" + "Species concentrations and reversal potential can be overridden on\n" + "specific regions using the paint interface, while the method for calculating\n" + "reversal potential is global for all compartments in the cell, and can't be\n" + "overriden locally.\n" + " ion: name of ion species.\n" + " int_con: initial internal concentration [mM].\n" + " ext_con: initial external concentration [mM].\n" + " rev_pot: initial reversal potential [mV].\n" + " rev_pot_method: method for calculating reversal potential.") + .def_readwrite("catalogue", &global_props_shim::cat, "The mechanism catalogue.") + .def("__str__", [](const global_props_shim& p){return to_string(p);}); + // arb::cable_cell pybind11::class_<arb::cable_cell> cable_cell(m, "cable_cell", @@ -480,10 +549,10 @@ void register_cells(pybind11::module& m) { }, "region"_a, "Vm"_a=pybind11::none(), "cm"_a=pybind11::none(), "rL"_a=pybind11::none(), "tempK"_a=pybind11::none(), "Set cable properties on a region.\n" - " region: initial membrane voltage [mV].\n" - " cm: membrane capacitance [F/m²].\n" - " rL: axial resistivity [Ω·cm].\n" - " tempK: temperature [Kelvin].") + " Vm: initial membrane voltage [mV].\n" + " cm: membrane capacitance [F/m²].\n" + " rL: axial resistivity [Ω·cm].\n" + " tempK: temperature [Kelvin].") // Paint ion species initial conditions on a region. .def("paint", @@ -529,6 +598,7 @@ void register_cells(pybind11::module& m) { .def("locations", [](arb::cable_cell& c, const char* label) {return c.concrete_locset(label);}, "label"_a, "The locations of the cell morphology for a locset label.") + // Get cables associated with a region label. .def("cables", [](arb::cable_cell& c, const char* label) {return c.concrete_region(label).cables();}, "label"_a, "The cable segments of the cell morphology for a region label.") diff --git a/python/cells.hpp b/python/cells.hpp index 3dbe72b43c5cd8bae949214313c81df4ca143e50..f4fd31759eedb981e4fee25f7f23928125b6dd12 100644 --- a/python/cells.hpp +++ b/python/cells.hpp @@ -2,8 +2,17 @@ #include <pybind11/pybind11.h> +#include <arbor/cable_cell_param.hpp> +#include <arbor/mechcat.hpp> #include <arbor/util/unique_any.hpp> namespace pyarb { arb::util::unique_any convert_cell(pybind11::object o); + +struct global_props_shim { + arb::mechanism_catalogue cat; + arb::cable_cell_global_properties props; + global_props_shim(); +}; + } diff --git a/python/conversion.hpp b/python/conversion.hpp index 5d58748adcc479fc4e5529390a3593e474e7448d..eda94ba044a36997dbe6b58dc1a303fdfe4c985d 100644 --- a/python/conversion.hpp +++ b/python/conversion.hpp @@ -78,6 +78,8 @@ arb::util::optional<T> py2optional(pybind11::object o, const char* msg) { // to cast is not an error. template <typename T> arb::util::optional<T> try_cast(pybind11::object o) { + if (o.is_none()) return arb::util::nullopt; + try { return o.cast<T>(); } diff --git a/python/mechanism.cpp b/python/mechanism.cpp new file mode 100644 index 0000000000000000000000000000000000000000..0df5d95da1394ca7712f3a1e91075ee6d26e24d2 --- /dev/null +++ b/python/mechanism.cpp @@ -0,0 +1,180 @@ +#include <pybind11/pybind11.h> +#include "pybind11/pytypes.h" +#include <pybind11/stl.h> + +#include <arbor/cable_cell_param.hpp> +#include <arbor/mechanism.hpp> +#include <arbor/mechcat.hpp> +#include <arbor/util/optional.hpp> +#include <stdexcept> + +#include "arbor/mechinfo.hpp" + +#include "conversion.hpp" +#include "strprintf.hpp" + +namespace pyarb { + +void apply_derive(arb::mechanism_catalogue& m, + const std::string& name, + const std::string& parent, + const std::unordered_map<std::string, double>& globals, + const std::unordered_map<std::string, std::string>& ions) +{ + if (globals.empty() && ions.empty()) { + m.derive(name, parent); + return; + } + std::vector<std::pair<std::string, double>> G; + for (auto& g: globals) { + G.push_back({g.first, g.second}); + } + std::vector<std::pair<std::string, std::string>> I; + for (auto& i: ions) { + I.push_back({i.first, i.second}); + } + m.derive(name, parent, G, I); +} + +void register_mechanisms(pybind11::module& m) { + using arb::util::optional; + using namespace pybind11::literals; + + pybind11::class_<arb::mechanism_field_spec> field_spec(m, "mechanism_field", + "Basic information about a mechanism field."); + field_spec + .def(pybind11::init<const arb::mechanism_field_spec&>()) + .def_readonly("units", &arb::mechanism_field_spec::units) + .def_readonly("default", &arb::mechanism_field_spec::default_value) + .def_readonly("min", &arb::mechanism_field_spec::lower_bound) + .def_readonly("max", &arb::mechanism_field_spec::upper_bound) + .def("__repr__", + [](const arb::mechanism_field_spec& spec) { + return util::pprintf("{units: '{}', default: {}, min: {}, max: {}}", + (spec.units.size()? spec.units.c_str(): "1"), spec.default_value, + spec.lower_bound, spec.upper_bound); }) + .def("__str__", + [](const arb::mechanism_field_spec& spec) { + return util::pprintf("{units: '{}', default: {}, min: {}, max: {}}", + (spec.units.size()? spec.units.c_str(): "1"), spec.default_value, + spec.lower_bound, spec.upper_bound); }); + + pybind11::class_<arb::ion_dependency> ion_dep(m, "ion_dependency", + "Information about a mechanism's dependence on an ion species."); + ion_dep + .def(pybind11::init<const arb::ion_dependency&>()) + .def_readonly("write_int_con", &arb::ion_dependency::write_concentration_int) + .def_readonly("write_ext_con", &arb::ion_dependency::write_concentration_ext) + .def_readonly("write_rev_pot", &arb::ion_dependency::write_reversal_potential) + .def_readonly("read_rev_pot", &arb::ion_dependency::read_reversal_potential) + .def("__repr__", + [](const arb::ion_dependency& dep) { + auto tf = [](bool x) {return x? "True": "False";}; + return util::pprintf("{write_int_con: {}, write_ext_con: {}, write_rev_pot: {}, read_rev_pot: {}}", + tf(dep.write_concentration_int), tf(dep.write_concentration_ext), + tf(dep.write_reversal_potential), tf(dep.read_reversal_potential)); }) + .def("__str__", + [](const arb::ion_dependency& dep) { + auto tf = [](bool x) {return x? "True": "False";}; + return util::pprintf("{write_int_con: {}, write_ext_con: {}, write_rev_pot: {}, read_rev_pot: {}}", + tf(dep.write_concentration_int), tf(dep.write_concentration_ext), + tf(dep.write_reversal_potential), tf(dep.read_reversal_potential)); }) + ; + + pybind11::class_<arb::mechanism_info> mech_inf(m, "mechanism_info", + "Meta data about a mechanism's fields and ion dependendencies."); + mech_inf + .def(pybind11::init<const arb::mechanism_info&>()) + .def_readonly("globals", &arb::mechanism_info::globals, + "Global fields have one value common to an instance of a mechanism, are constant in time and set at instantiation.") + .def_readonly("parameters", &arb::mechanism_info::parameters, + "Parameter fields may vary across the extent of a mechanism, but are constant in time and set at instantiation.") + .def_readonly("state", &arb::mechanism_info::state, + "State fields vary in time and across the extent of a mechanism, and potentially can be sampled at run-time.") + .def_readonly("ions", &arb::mechanism_info::ions, + "Ion dependencies.") + .def_readonly("linear", &arb::mechanism_info::linear, + "True if a synapse mechanism has linear current contributions so that multiple instances on the same compartment can be coalesed.") + .def("__repr__", + [](const arb::mechanism_info& inf) { + return util::pprintf("(arbor.mechanism_info)"); }) + .def("__str__", + [](const arb::mechanism_info& inf) { + return util::pprintf("(arbor.mechanism_info)"); }); + + pybind11::class_<arb::mechanism_catalogue> cat(m, "mechanism_catalogue"); + cat + .def(pybind11::init<const arb::mechanism_catalogue&>()) + .def("has", &arb::mechanism_catalogue::has, + "name"_a, "Is 'name' in the catalogue?") + .def("is_derived", &arb::mechanism_catalogue::is_derived, + "name"_a, "Is 'name' a derived mechanism or can it be implicitly derived?") + .def("__getitem__", + [](arb::mechanism_catalogue& c, const char* name) { + try { + return c[name]; + } + catch (...) { + throw std::runtime_error(util::pprintf("\nKeyError: '{}'", name)); + } + }) + .def("derive", &apply_derive, + "name"_a, "parent"_a, + "globals"_a=std::unordered_map<std::string, double>{}, + "ions"_a=std::unordered_map<std::string, std::string>{}) + .def("__repr__", + [](const arb::mechanism_catalogue& cat) { + return util::pprintf("<arbor.mechanism_catalogue>"); }) + .def("__str__", + [](const arb::mechanism_catalogue& cat) { + return util::pprintf("<arbor.mechanism_catalogue>"); }); + + m.def("default_catalogue", [](){return arb::global_default_catalogue();}); + + // arb::mechanism_desc + // For specifying a mechanism in the cable_cell interface. + + pybind11::class_<arb::mechanism_desc> mechanism_desc(m, "mechanism"); + mechanism_desc + .def(pybind11::init([](const char* n) {return arb::mechanism_desc{n};})) + // allow construction of a description with parameters provided in a dictionary: + // mech = arbor.mechanism('mech_name', {'param1': 1.2, 'param2': 3.14}) + .def(pybind11::init( + [](const char* name, std::unordered_map<std::string, double> params) { + arb::mechanism_desc md(name); + for (const auto& p: params) md.set(p.first, p.second); + return md; + }), + "name"_a, "The name of the mechanism", + "params"_a, "A dictionary of parameter values, with parameter name as key.", + "Example usage setting pararmeters:\n" + " m = arbor.mechanism('expsyn', {'tau': 1.4})\n" + "will create parameters for the 'expsyn' mechanism, with the provided value\n" + "for 'tau' overrides the default. If a parameter is not set, the default\n" + "(as defined in NMODL) is used.\n\n" + "Example overriding a global parameter:\n" + " m = arbor.mechanism('nernst/R=8.3145,F=96485')") + .def("set", + [](arb::mechanism_desc& md, std::string name, double value) { + md.set(name, value); + }, + "name"_a, "value"_a, "Set parameter value.") + .def_property_readonly("name", + [](const arb::mechanism_desc& md) { + return md.name(); + }, + "The name of the mechanism.") + .def_property_readonly("values", + [](const arb::mechanism_desc& md) { + return md.values(); + }, "A dictionary of parameter values with parameter name as key.") + .def("__repr__", + [](const arb::mechanism_desc& md) { + return util::pprintf("<arbor.mechanism: name '{}', parameters {}", md.name(), util::dictionary_csv(md.values())); }) + .def("__str__", + [](const arb::mechanism_desc& md) { + return util::pprintf("('{}' {})", md.name(), util::dictionary_csv(md.values())); }); + +} + +} // namespace pyarb diff --git a/python/morphology.cpp b/python/morphology.cpp index 443a1642de2dd0371f8918d5da234db46c4e8deb..ea2ea8101b9ea64c632782f16444dab0fb5cd1fd 100644 --- a/python/morphology.cpp +++ b/python/morphology.cpp @@ -32,10 +32,10 @@ void register_morphology(pybind11::module& m) { pyarb::assert_throw(arb::test_invariants(mloc), "invalid location"); return mloc; }), - "branch"_a, "position"_a, + "branch"_a, "pos"_a, "Construct a location specification holding:\n" " branch: The id of the branch.\n" - " position: The relative position (from 0., proximal, to 1., distal) on the branch.\n") + " pos: The relative position (from 0., proximal, to 1., distal) on the branch.\n") .def_readonly("branch", &arb::mlocation::branch, "The id of the branch.") .def_readonly("pos", &arb::mlocation::pos, @@ -45,28 +45,8 @@ void register_morphology(pybind11::module& m) { .def("__repr__", [](arb::mlocation l) { return util::pprintf("(location {} {})", l.branch, l.pos); }); - // arb::mpoint - pybind11::class_<arb::mpoint> mpoint(m, "mpoint"); - mpoint - .def(pybind11::init( - [](double x, double y, double z, double r) { - return arb::mpoint{x,y,z,r}; - }), - "x"_a, "y"_a, "z"_a, "radius"_a, "All values in μm.") - .def_readonly("x", &arb::mpoint::x, "X coordinate [μm].") - .def_readonly("y", &arb::mpoint::y, "Y coordinate [μm].") - .def_readonly("z", &arb::mpoint::z, "Z coordinate [μm].") - .def_readonly("radius", &arb::mpoint::radius, - "Radius of cable at sample location centered at coordinates [μm].") - .def("__str__", - [](const arb::mpoint& p) { - return util::pprintf("<arbor.mpoint: x {}, y {}, z {}, radius {}>", p.x, p.y, p.z, p.radius); - }) - .def("__repr__", - [](const arb::mpoint& p) {return util::pprintf("{}>", p);}); - // arb::msample - pybind11::class_<arb::msample> msample(m, "msample"); + pybind11::class_<arb::msample> msample(m, "sample"); msample .def(pybind11::init( [](arb::mpoint loc, int tag){ @@ -77,8 +57,16 @@ void register_morphology(pybind11::module& m) { return arb::msample{{x,y,z,r}, tag};}), "x"_a, "y"_a, "z"_a, "radius"_a, "tag"_a, "spatial values {x, y, z, radius} are in μm.") - .def_readonly("loc", &arb::msample::loc, - "Location of sample.") + .def_readonly("x", &arb::msample::loc) + .def_property_readonly("x", + [](const arb::msample& s) {return s.loc.x;}, + "X coordinate [μm].") + .def_property_readonly("y", + [](const arb::msample& s) {return s.loc.y;}, + "Y coordinate [μm].") + .def_property_readonly("z", + [](const arb::msample& s) {return s.loc.z;}, + "Z coordinate [μm].") .def_readonly("tag", &arb::msample::tag, "Property tag of sample. " "If loaded from standard SWC file the following tags are used: soma=1, axon=2, dendrite=3, apical dendrite=4, however arbitrary tags can be used.") @@ -93,22 +81,21 @@ void register_morphology(pybind11::module& m) { pybind11::class_<arb::mcable> cable(m, "cable"); cable .def(pybind11::init( - [](arb::msize_t bid, double prox, double dist) { - arb::mcable c{bid, prox, dist}; - if (!test_invariants(c)) { - throw pyarb_error("Invalid cable description. Cable segments must have proximal and distal end points in the range [0,1]."); - } - return c; - }), - "branch_id"_a, "prox"_a, "dist"_a) + [](arb::msize_t bid, double prox, double dist) { + arb::mcable c{bid, prox, dist}; + if (!test_invariants(c)) { + throw pyarb_error("Invalid cable description. Cable segments must have proximal and distal end points in the range [0,1]."); + } + return c; + }), + "branch"_a, "prox"_a, "dist"_a) + .def_readonly("branch", &arb::mcable::branch, + "The id of the branch on which the cable lies.") .def_readonly("prox", &arb::mcable::prox_pos, "The relative position of the proximal end of the cable on its branch ∈ [0,1].") .def_readonly("dist", &arb::mcable::dist_pos, "The relative position of the distal end of the cable on its branch ∈ [0,1].") - .def_readonly("branch", &arb::mcable::branch, - "The id of the branch on which the cable lies.") - .def("__str__", [](const arb::mcable& c) { - return util::pprintf("<arbor.cable: branch {}, prox {}, dist {}", c.branch, c.prox_pos, c.dist_pos); }) + .def("__str__", [](const arb::mcable& c) { return util::pprintf("{}", c); }) .def("__repr__", [](const arb::mcable& c) { return util::pprintf("{}", c); }); // @@ -190,7 +177,7 @@ void register_morphology(pybind11::module& m) { // be implemented as read-only properties. .def_property_readonly("empty", [](const arb::morphology& m){return m.empty();}, - "A list with the parent index of each sample.") + "Whether the morphology is empty.") .def_property_readonly("spherical_root", [](const arb::morphology& m){return m.spherical_root();}, "Whether the root of the morphology is spherical.") diff --git a/python/pyarb.cpp b/python/pyarb.cpp index ae2dc9a0843350567741874f42ebe10820f55249..a3d97d5136b68dc4446ac695166906e7c96eeb56 100644 --- a/python/pyarb.cpp +++ b/python/pyarb.cpp @@ -13,6 +13,7 @@ void register_domain_decomposition(pybind11::module& m); void register_event_generators(pybind11::module& m); void register_flat_builder(pybind11::module& m); void register_identifiers(pybind11::module& m); +void register_mechanisms(pybind11::module& m); void register_morphology(pybind11::module& m); void register_profiler(pybind11::module& m); void register_recipe(pybind11::module& m); @@ -38,6 +39,7 @@ PYBIND11_MODULE(_arbor, m) { pyarb::register_event_generators(m); pyarb::register_flat_builder(m); pyarb::register_identifiers(m); + pyarb::register_mechanisms(m); pyarb::register_morphology(m); pyarb::register_profiler(m); pyarb::register_recipe(m); diff --git a/python/s_expr.cpp b/python/s_expr.cpp index e879c17d1a02c7a81f8749e54950fd5d7fb85fcd..129d10801266277d017ee0c05f1471aa848d543a 100644 --- a/python/s_expr.cpp +++ b/python/s_expr.cpp @@ -38,6 +38,9 @@ std::ostream& operator<<(std::ostream& o, const tok& t) { } std::ostream& operator<<(std::ostream& o, const token& t) { + if (t.kind==tok::string) { + return o << util::pprintf("\"{}\"", t.spelling); + } return o << util::pprintf("{}", t.spelling); } @@ -111,6 +114,9 @@ private: character(); continue; // skip to next character + case ';': + eat_comment(); + continue; case '(': token_ = {loc_, tok::lparen, {character()}}; return; @@ -154,6 +160,14 @@ private: return; } + // Consumes characters in the stream until end of stream or a new line. + // Assumes that the current_ location is the `;` that starts the comment. + void eat_comment() { + while (*current_ && *current_!='\n') { + character(); + } + } + // Parse alphanumeric sequence that starts with an alphabet character, // and my contain alphabet, numeric or underscore '_' characters. // diff --git a/python/s_expr.hpp b/python/s_expr.hpp index f24e5f7cdd1f8b3fadb93bcc008356d69e3220fa..37975ee77275dd32c83904fc46d1af0f25d1cd59 100644 --- a/python/s_expr.hpp +++ b/python/s_expr.hpp @@ -90,7 +90,7 @@ struct s_expr { // with a std::shared_ptr. using pair_type = s_pair<value_wrapper<s_expr>>; - arb::util::either<token, pair_type> state; + arb::util::either<token, pair_type> state = token{-1, tok::nil, "nil"}; s_expr(const s_expr& s): state(s.state) {} s_expr() = default; diff --git a/python/sampling.cpp b/python/sampling.cpp index 50d68a42c1585ba20d437a3f6dbdd49831979bd3..6220f625f1cbb9bfc33d8b0de8d1004777b19da1 100644 --- a/python/sampling.cpp +++ b/python/sampling.cpp @@ -14,7 +14,7 @@ namespace pyarb { // TODO: trace entry of different types/container (e.g. vector of doubles to get all samples of a cell) -struct trace_entry { +struct trace_sample { arb::time_type t; double v; }; @@ -22,9 +22,9 @@ struct trace_entry { // A helper struct (state) ensuring that only one thread can write to the probe_buffers holding the trace entries (mapped by probe id) struct sampler_state { std::mutex mutex; - std::unordered_map<arb::cell_member_type, std::vector<trace_entry>> probe_buffers; + std::unordered_map<arb::cell_member_type, std::vector<trace_sample>> probe_buffers; - std::vector<trace_entry>& probe_buffer(arb::cell_member_type pid) { + std::vector<trace_sample>& probe_buffer(arb::cell_member_type pid) { // lock the mutex, s.t. other threads cannot write std::lock_guard<std::mutex> lock(mutex); // return or create entry @@ -37,20 +37,20 @@ struct sampler_state { } // helper function to push back to locked vector - void push_back(arb::cell_member_type pid, trace_entry value) { + void push_back(arb::cell_member_type pid, trace_sample value) { auto& v = probe_buffer(pid); v.push_back(std::move(value)); } // Access the probe buffers - const std::unordered_map<arb::cell_member_type, std::vector<trace_entry>>& samples() const { + const std::unordered_map<arb::cell_member_type, std::vector<trace_sample>>& samples() const { return probe_buffers; } }; // A functor that models arb::sampler_function. -// Holds a shared pointer to the trace_entry used to store the samples, so that if -// the trace_entry in sampler is garbage collected in Python, stores will +// Holds a shared pointer to the trace_sample used to store the samples, so that if +// the trace_sample in sampler is garbage collected in Python, stores will // not seg fault. struct sample_callback { @@ -89,7 +89,7 @@ struct sampler { return sample_callback(sample_store); } - const std::vector<trace_entry>& samples(arb::cell_member_type pid) const { + const std::vector<trace_sample>& samples(arb::cell_member_type pid) const { if (!sample_store->has_pid(pid)) { throw std::runtime_error(util::pprintf("probe id {} does not exist", pid)); } @@ -117,7 +117,7 @@ std::shared_ptr<sampler> attach_sampler(arb::simulation& sim, arb::time_type int return r; } -std::string sample_str(const trace_entry& s) { +std::string sample_str(const trace_sample& s) { return util::pprintf("<arbor.sample: time {} ms, \tvalue {}>", s.t, s.v); } @@ -125,10 +125,10 @@ void register_sampling(pybind11::module& m) { using namespace pybind11::literals; // Sample - pybind11::class_<trace_entry> sample(m, "sample"); - sample - .def_readonly("time", &trace_entry::t, "The sample time [ms] at a specific probe.") - .def_readonly("value", &trace_entry::v, "The sample record at a specific probe.") + pybind11::class_<trace_sample> trace_sample(m, "trace_sample"); + trace_sample + .def_readonly("time", &trace_sample::t, "The sample time [ms] at a specific probe.") + .def_readonly("value", &trace_sample::v, "The sample record at a specific probe.") .def("__str__", &sample_str) .def("__repr__", &sample_str); diff --git a/python/single_cell_model.cpp b/python/single_cell_model.cpp index dfba5158ff45faa0aa829cb66582f27cc0a976cc..4e406d5afe6044cdec08da8a4e7ad08982a901fb 100644 --- a/python/single_cell_model.cpp +++ b/python/single_cell_model.cpp @@ -10,6 +10,7 @@ #include <arbor/simulation.hpp> #include "error.hpp" +#include "cells.hpp" namespace pyarb { @@ -133,7 +134,6 @@ class single_cell_model { arb::cable_cell cell_; arb::context ctx_; bool run_ = false; - arb::cable_cell_global_properties gprop_; std::vector<probe_site> probes_; std::unique_ptr<arb::simulation> sim_; @@ -142,10 +142,15 @@ class single_cell_model { std::vector<trace> traces_; public: + // Make gprop public to make it possible to expose it as a field in Python: + // model.properties.catalogue = my_custom_cat + //arb::cable_cell_global_properties gprop; + global_props_shim gprop; + single_cell_model(arb::cable_cell c): cell_(std::move(c)), ctx_(arb::make_context()) { - gprop_.default_parameters = arb::neuron_parameter_defaults; + gprop.props.default_parameters = arb::neuron_parameter_defaults; } // example use: @@ -167,12 +172,8 @@ public: } } - void add_ion(const std::string& ion, double valence, double int_con, double ext_con, double rev_pot) { - gprop_.add_ion(ion, valence, int_con, ext_con, rev_pot); - } - void run(double tfinal, double dt) { - single_cell_recipe rec(cell_, probes_, gprop_); + single_cell_recipe rec(cell_, probes_, gprop.props); auto domdec = arb::partition_load_balance(rec, ctx_); @@ -254,20 +255,14 @@ void register_single_cell(pybind11::module& m) { " what: Name of the variable to record (currently only 'voltage').\n" " where: Location on cell morphology at which to sample the variable.\n" " frequency: The target frequency at which to sample [Hz].") - .def("add_ion", &single_cell_model::add_ion, - "ion"_a, "valence"_a, "int_con"_a, "ext_con"_a, "rev_pot"_a, - "Add a new ion species to the model.\n" - " ion: name of the ion species.\n" - " valence: valence of the ion species.\n" - " int_con: initial internal concentration [mM].\n" - " ext_con: initial external concentration [mM].\n" - " rev_pot: reversal potential [mV].") .def_property_readonly("spikes", [](const single_cell_model& m) { return m.spike_times();}, "Holds spike times [ms] after a call to run().") .def_property_readonly("traces", [](const single_cell_model& m) { - return m.traces();}, "Holds sample traces after a call to run().") + return m.traces();}, + "Holds sample traces after a call to run().") + .def_readwrite("properties", &single_cell_model::gprop, "Global properties.") .def("__repr__", [](const single_cell_model&){return "<arbor.single_cell_model>";}) .def("__str__", [](const single_cell_model&){return "<arbor.single_cell_model>";}); } diff --git a/python/strprintf.hpp b/python/strprintf.hpp index e8f96100fbfff0f74320eec7dc5eb751bdb645c1..9e4687129fe1bbee314b3b02489117c122181dc0 100644 --- a/python/strprintf.hpp +++ b/python/strprintf.hpp @@ -13,6 +13,8 @@ #include <arbor/util/optional.hpp> +#include "s_expr.hpp" + namespace pyarb { namespace util { @@ -202,13 +204,15 @@ impl::sepval_lim<Seq> csv(const Seq& seq, unsigned n) { // Print dictionary: this could be done easily with range adaptors in C++17 template <typename Key, typename T> std::string dictionary_csv(const std::unordered_map<Key, T>& dict) { - constexpr bool string_key = std::is_same<std::string, std::decay_t<Key>>::value; + constexpr bool string_key = std::is_same<std::string, std::decay_t<Key>>::value; + constexpr bool string_value = std::is_same<std::string, std::decay_t<T>>::value; + std::string format = pprintf("{}: {}", string_key? "\"{}\"": "{}", string_value? "\"{}\"": "{}"); std::string s = "{"; bool first = true; for (auto& p: dict) { if (!first) s += ", "; first = false; - s += pprintf(string_key? "'{}': {}": "{}: {}", p.first, p.second); + s += pprintf(format.c_str(), p.first, p.second); } s += "}"; return s; @@ -216,10 +220,13 @@ std::string dictionary_csv(const std::unordered_map<Key, T>& dict) { } // namespace util +} + +namespace arb { +namespace util { template <typename T> std::ostream& operator<<(std::ostream& o, const arb::util::optional<T>& x) { - return o << (x? util::to_string(*x): "None"); + return o << (x? pyarb::util::to_string(*x): "None"); +} +} } - -} // namespace pyarb - diff --git a/python/test/cpp/s_expr.cpp b/python/test/cpp/s_expr.cpp index 772103e507a15b036b8ca515619d041ebae18d65..62300624a8f7c0f2ac8576fc1a59fc62afb29603 100644 --- a/python/test/cpp/s_expr.cpp +++ b/python/test/cpp/s_expr.cpp @@ -5,6 +5,7 @@ #include <morph_parse.hpp> #include <s_expr.hpp> +#include <strprintf.hpp> using namespace pyarb; using namespace std::string_literals; @@ -90,3 +91,22 @@ TEST(s_expr, parse) { EXPECT_EQ(util::pprintf("{}", join(lhs,rhs)), "(join (region \"dend\") (all))"); } + +TEST(s_expr, comments) { + auto round_trip_reg = [](const char* in) { + auto x = eval(parse(in)); + return util::pprintf("{}", arb::util::any_cast<arb::region>(*x)); + }; + + EXPECT_EQ("(all)", round_trip_reg("(all) ; a comment")); + const char *multi_line = + "; comment at start\n" + "(radius_lt\n" + " (join\n" + " (tag 3) ; end of line\n" + " ;comment on whole line\n" + " (tag 4))\n" + " 0.5) ; end of string"; + EXPECT_EQ("(radius_lt (join (tag 3) (tag 4)) 0.5)", + round_trip_reg(multi_line)); +} diff --git a/python/tutorial/example1.py b/python/tutorial/example1.py index 43bfd5c7ce4f2a275b1d26f82195ae0d344c6d0f..5e9b49cdb2767e42993326b3bf49ea4d17719baf 100644 --- a/python/tutorial/example1.py +++ b/python/tutorial/example1.py @@ -2,7 +2,7 @@ import arbor # Create a sample tree with a single sample of radius 3 μm tree = arbor.sample_tree() -tree.append(arbor.msample(x=0, y=0, z=0, radius=3, tag=2)) +tree.append(arbor.sample(x=0, y=0, z=0, radius=3, tag=2)) labels = arbor.label_dict({'soma': '(tag 2)', 'center': '(location 0 0.5)'}) @@ -12,20 +12,20 @@ cell = arbor.cable_cell(tree, labels) cell.set_properties(Vm=-40) # Put hh dynamics on soma, and passive properties on the dendrites. cell.paint('soma', 'hh') -# Attach stimuli with duration of 2 ms and current of 0.8 nA. -cell.place('center', arbor.iclamp( 10, 2, 0.8)) +# Attach stimuli with duration of 1 ms and current of 0.8 nA. +cell.place('center', arbor.iclamp( 10, duration=1, current=0.8)) # Add a spike detector with threshold of -10 mV. cell.place('center', arbor.spike_detector(-10)) # Make single cell model. m = arbor.single_cell_model(cell) -# Attach voltage probes, sampling at 10 kHz. -m.probe('voltage', 'center', 10000) +# Attach voltage probes, sampling at 1 MHz. +m.probe('voltage', 'center', 1000000) -# Run simulation for 100 ms of simulated activity. +# Run simulation for tfinal ms with time steps of 1 μs. tfinal=30 -m.run(tfinal) +m.run(tfinal, dt=0.001) # Print spike times. if len(m.spikes)>0: @@ -45,7 +45,7 @@ legend_labels = ['{}: {}'.format(s.variable, s.location) for s in m.traces] ax.legend(legend_labels) ax.set(xlabel='time (ms)', ylabel='voltage (mV)', title='cell builder demo') plt.xlim(0,tfinal) -plt.ylim(-80,50) +plt.ylim(-80,80) ax.grid() # Set to True to save the image to file instead of opening a plot window. diff --git a/test/unit/test_morph_expr.cpp b/test/unit/test_morph_expr.cpp index af9ef3c697902f57956d2bfac045c9cbd2caf954..08e52d1e919ebc46b3f66ab4c966c8837ce63b76 100644 --- a/test/unit/test_morph_expr.cpp +++ b/test/unit/test_morph_expr.cpp @@ -567,7 +567,7 @@ TEST(region, thingify_moderate_morphologies) { EXPECT_TRUE(region_eq(mp, radius_le(all(), 2), cl{{0,0,0.55}, {1,0,0.325}, {2,1,1}, {3,0.375,0.75}})); EXPECT_TRUE(region_eq(mp, radius_le(all(), 3), cl{{0,0,1}, {1,0,0.55}, {2,6.0/9.0,1}, {3,0.25,1}})); EXPECT_TRUE(region_eq(mp, radius_ge(all(), 2), cl{{0,0.55,1}, {1,0.325,1}, {2,0,1}, {3,0,0.375}, {3,0.75,1}})); - EXPECT_TRUE(region_eq(mp, radius_ge(all(), 3), cl{{1,0.55,1}, {2,0,6.0/9.0}, {3,0,0.25}})); + EXPECT_TRUE(region_eq(mp, radius_ge(all(), 3), cl{{0,1,1}, {1,0.55,1}, {2,0,6.0/9.0}, {3,0,0.25}, {3,1,1}})); EXPECT_TRUE(region_eq(mp, radius_lt(reg_a_, 2), cl{{0,0.1,0.4},{3,0.375,0.4}})); EXPECT_TRUE(region_eq(mp, radius_gt(reg_a_, 2), cl{{2,0,1},{3,0.1,0.375}}));